projectEdge.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2016-2019 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
27 #include "projectEdge.H"
28 #include "unitConversion.H"
30 #include "pointConstraint.H"
31 #include "OBJstream.H"
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38  defineTypeNameAndDebug(projectEdge, 0);
39  addToRunTimeSelectionTable(blockEdge, projectEdge, Istream);
40 }
41 
42 
43 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
44 
45 void Foam::projectEdge::findNearest
46 (
47  const point& pt,
48  point& near,
49  pointConstraint& constraint
50 ) const
51 {
52  if (surfaces_.size())
53  {
54  const scalar distSqr = magSqr(points_[end_] - points_[start_]);
55 
56  pointField boundaryNear(1);
57  List<pointConstraint> boundaryConstraint(1);
58 
60  (
61  geometry_,
62  surfaces_,
63  pointField(1, pt),
64  scalarField(1, distSqr),
65  boundaryNear,
66  boundaryConstraint
67  );
68  near = boundaryNear[0];
69  constraint = boundaryConstraint[0];
70  }
71  else
72  {
73  near = pt;
74  constraint = pointConstraint();
75  }
76 }
77 
78 
79 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
80 
82 (
83  const dictionary& dict,
84  const label index,
85  const searchableSurfaces& geometry,
86  const pointField& points,
87  Istream& is
88 )
89 :
90  blockEdge(dict, index, points, is),
91  geometry_(geometry)
92 {
93  wordList names(is);
94  surfaces_.setSize(names.size());
95  forAll(names, i)
96  {
97  surfaces_[i] = geometry_.findSurfaceID(names[i]);
98 
99  if (surfaces_[i] == -1)
100  {
102  << "Cannot find surface " << names[i] << " in geometry"
103  << exit(FatalIOError);
104  }
105  }
106 }
107 
108 
109 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
110 
112 {
113  // Initial guess
114  const point start
115  (
116  points_[start_] + lambda*(points_[end_] - points_[start_])
117  );
118 
119  point near(start);
120 
121  if (lambda >= small && lambda < 1 - small)
122  {
123  pointConstraint constraint;
124  findNearest(start, near, constraint);
125  }
126 
127  return near;
128 }
129 
130 
133 {
134  // For debugging to tag the output
135  static label eIter = 0;
136 
137  autoPtr<OBJstream> debugStr;
138  if (debug)
139  {
140  debugStr.reset
141  (
142  new OBJstream("projectEdge_" + Foam::name(eIter++) + ".obj")
143  );
144  Info<< "Writing lines from straight-line start points"
145  << " to projected points to " << debugStr().name() << endl;
146  }
147 
148 
149  tmp<pointField> tpoints(new pointField(lambdas.size()));
150  pointField& points = tpoints.ref();
151 
152  const point& startPt = points_[start_];
153  const point& endPt = points_[end_];
154  const vector d = endPt-startPt;
155 
156  // Initial guess
157  forAll(lambdas, i)
158  {
159  points[i] = startPt+lambdas[i]*d;
160  }
161 
162 
163  // Upper limit for number of iterations
164  const label maxIter = 10;
165 
166  // Residual tolerance
167  const scalar relTol = 0.1;
168  const scalar absTol = 1e-4;
169 
170  scalar initialResidual = 0;
171 
172  for (label iter = 0; iter < maxIter; iter++)
173  {
174  // Do projection
175  {
176  List<pointConstraint> constraints(lambdas.size());
177  pointField start(points);
179  (
180  geometry_,
181  surfaces_,
182  start,
183  scalarField(start.size(), magSqr(d)),
184  points,
185  constraints
186  );
187 
188  // Reset start and end point
189  if (lambdas[0] < small)
190  {
191  points[0] = startPt;
192  }
193  if (lambdas.last() > 1 - small)
194  {
195  points.last() = endPt;
196  }
197 
198  if (debugStr.valid())
199  {
200  forAll(points, i)
201  {
202  debugStr().write(linePointRef(start[i], points[i]));
203  }
204  }
205  }
206 
207  // Calculate lambdas (normalised coordinate along edge)
208  scalarField projLambdas(points.size());
209  {
210  projLambdas[0] = 0;
211  for (label i = 1; i < points.size(); i++)
212  {
213  projLambdas[i] =
214  projLambdas[i-1] + mag(points[i] - points[i-1]);
215  }
216  projLambdas /= projLambdas.last();
217  }
218  linearInterpolationWeights interpolator(projLambdas);
219 
220  // Compare actual distances and move points (along straight line;
221  // not along surface)
222  vectorField residual(points.size(), vector::zero);
223  labelList indices;
224  scalarField weights;
225  for (label i = 1; i < points.size() - 1; i++)
226  {
227  interpolator.valueWeights(lambdas[i], indices, weights);
228 
229  point predicted = vector::zero;
230  forAll(indices, indexi)
231  {
232  predicted += weights[indexi]*points[indices[indexi]];
233  }
234  residual[i] = predicted - points[i];
235  }
236 
237  const scalar scalarResidual = sum(mag(residual));
238 
239  if (debug)
240  {
241  Pout<< "Iter:" << iter << " initialResidual:" << initialResidual
242  << " residual:" << scalarResidual << endl;
243  }
244 
245  if (scalarResidual < absTol*0.5*lambdas.size())
246  {
247  break;
248  }
249  else if (iter == 0)
250  {
251  initialResidual = scalarResidual;
252  }
253  else if (scalarResidual/initialResidual < relTol)
254  {
255  break;
256  }
257 
258  if (debugStr.valid())
259  {
260  forAll(points, i)
261  {
262  const point predicted(points[i] + residual[i]);
263  debugStr().write(linePointRef(points[i], predicted));
264  }
265  }
266 
267  points += residual;
268  }
269 
270  return tpoints;
271 }
272 
273 
274 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
void reset(T *=nullptr)
If object pointer already set, delete object and set to given.
Definition: autoPtrI.H:114
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
projectEdge(const dictionary &dict, const label index, const searchableSurfaces &geometry, const pointField &points, Istream &)
Construct from Istream setting pointsList.
Definition: projectEdge.C:82
static void findNearest(const PtrList< searchableSurface > &, const labelList &surfacesToTest, const pointField &, const scalarField &nearestDistSqr, labelList &surfaces, List< pointIndexHit > &)
Find nearest. Return -1 (and a miss()) or surface and nearest.
virtual point position(const scalar) const
Return the point positions corresponding to the curve parameters.
Definition: projectEdge.C:111
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
Unit conversion functions.
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
Definition: Istream.H:57
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
dimensionedScalar lambda(viscosity->lookup("lambda"))
Macros for easy insertion into run-time selection tables.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
virtual bool valueWeights(const scalar t, labelList &indices, scalarField &weights) const
Calculate weights and indices to calculate t from samples.
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Container for searchableSurfaces.
Accumulates point constraints through successive applications of the applyConstraint function...
bool valid() const
Return true if the autoPtr valid (ie, the pointer is set)
Definition: autoPtrI.H:83
dimensioned< scalar > magSqr(const dimensioned< Type > &)
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
OFstream which keeps track of vertices.
Definition: OBJstream.H:53
defineTypeNameAndDebug(combustionModel, 0)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
Define a curved edge that is parameterised for 0<lambda<1 between the start and end point...
Definition: blockEdge.H:56
vector point
Point is a vector.
Definition: point.H:41
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:318
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
A class for managing temporary objects.
Definition: PtrList.H:53
T & last()
Return the last element of the list.
Definition: UListI.H:128
Namespace for OpenFOAM.
IOerror FatalIOError