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-2018 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 
81 Foam::projectEdge::projectEdge
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(points_[start_] + lambda*(points_[end_]-points_[start_]));
115 
116  point near(start);
117 
118  if (lambda >= small && lambda < 1.0-small)
119  {
120  pointConstraint constraint;
121  findNearest(start, near, constraint);
122  }
123 
124  return near;
125 }
126 
127 
130 {
131  // For debugging to tag the output
132  static label eIter = 0;
133 
134  autoPtr<OBJstream> debugStr;
135  if (debug)
136  {
137  debugStr.reset
138  (
139  new OBJstream("projectEdge_" + Foam::name(eIter++) + ".obj")
140  );
141  Info<< "Writing lines from straight-line start points"
142  << " to projected points to " << debugStr().name() << endl;
143  }
144 
145 
146  tmp<pointField> tpoints(new pointField(lambdas.size()));
147  pointField& points = tpoints.ref();
148 
149  const point& startPt = points_[start_];
150  const point& endPt = points_[end_];
151  const vector d = endPt-startPt;
152 
153  // Initial guess
154  forAll(lambdas, i)
155  {
156  points[i] = startPt+lambdas[i]*d;
157  }
158 
159 
160  // Upper limit for number of iterations
161  const label maxIter = 10;
162  // Residual tolerance
163  const scalar relTol = 0.1;
164  const scalar absTol = 1e-4;
165 
166  scalar initialResidual = 0.0;
167 
168  for (label iter = 0; iter < maxIter; iter++)
169  {
170  // Do projection
171  {
172  List<pointConstraint> constraints(lambdas.size());
173  pointField start(points);
175  (
176  geometry_,
177  surfaces_,
178  start,
179  scalarField(start.size(), magSqr(d)),
180  points,
181  constraints
182  );
183 
184  // Reset start and end point
185  if (lambdas[0] < small)
186  {
187  points[0] = startPt;
188  }
189  if (lambdas.last() > 1.0-small)
190  {
191  points.last() = endPt;
192  }
193 
194  if (debugStr.valid())
195  {
196  forAll(points, i)
197  {
198  debugStr().write(linePointRef(start[i], points[i]));
199  }
200  }
201  }
202 
203  // Calculate lambdas (normalised coordinate along edge)
204  scalarField projLambdas(points.size());
205  {
206  projLambdas[0] = 0.0;
207  for (label i = 1; i < points.size(); i++)
208  {
209  projLambdas[i] = projLambdas[i-1] + mag(points[i]-points[i-1]);
210  }
211  projLambdas /= projLambdas.last();
212  }
213  linearInterpolationWeights interpolator(projLambdas);
214 
215  // Compare actual distances and move points (along straight line;
216  // not along surface)
217  vectorField residual(points.size(), vector::zero);
218  labelList indices;
219  scalarField weights;
220  for (label i = 1; i < points.size() - 1; i++)
221  {
222  interpolator.valueWeights(lambdas[i], indices, weights);
223 
224  point predicted = vector::zero;
225  forAll(indices, indexi)
226  {
227  predicted += weights[indexi]*points[indices[indexi]];
228  }
229  residual[i] = predicted-points[i];
230  }
231 
232  scalar scalarResidual = sum(mag(residual));
233 
234  if (debug)
235  {
236  Pout<< "Iter:" << iter << " initialResidual:" << initialResidual
237  << " residual:" << scalarResidual << endl;
238  }
239 
240  if (scalarResidual < absTol*0.5*lambdas.size())
241  {
242  break;
243  }
244  else if (iter == 0)
245  {
246  initialResidual = scalarResidual;
247  }
248  else if (scalarResidual/initialResidual < relTol)
249  {
250  break;
251  }
252 
253 
254  if (debugStr.valid())
255  {
256  forAll(points, i)
257  {
258  const point predicted(points[i] + residual[i]);
259  debugStr().write(linePointRef(points[i], predicted));
260  }
261  }
262 
263  points += residual;
264  }
265 
266  return tpoints;
267 }
268 
269 
270 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
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:137
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:174
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:163
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
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 parameterized 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:331
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
dimensionedScalar lambda(laminarTransport.lookup("lambda"))
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:98
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