projectCurveEdge.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 "projectCurveEdge.H"
28 #include "unitConversion.H"
30 #include "pointConstraint.H"
31 #include "OBJstream.H"
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39  defineTypeNameAndDebug(projectCurveEdge, 0);
40  addToRunTimeSelectionTable(blockEdge, projectCurveEdge, Istream);
41 }
42 
43 
44 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
45 
47 (
48  const dictionary& dict,
49  const label index,
50  const searchableSurfaces& geometry,
51  const pointField& points,
52  Istream& is
53 )
54 :
55  blockEdge(dict, index, points, is),
56  geometry_(geometry)
57 {
58  wordList names(is);
59  surfaces_.setSize(names.size());
60  forAll(names, i)
61  {
62  surfaces_[i] = geometry_.findSurfaceID(names[i]);
63 
64  if (surfaces_[i] == -1)
65  {
67  << "Cannot find surface " << names[i] << " in geometry"
68  << exit(FatalIOError);
69  }
70 
71  if (isA<searchableExtrudedCircle>(geometry_[surfaces_[i]]))
72  {
73  Info<< type() << " : Using curved surface "
74  << geometry_[surfaces_[i]].name()
75  << " to predict starting points." << endl;
76  }
77  }
78 }
79 
80 
81 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
82 
85 {
86  // For debugging to tag the output
87  static label eIter = 0;
88 
89  autoPtr<OBJstream> debugStr;
90  if (debug)
91  {
92  debugStr.reset
93  (
94  new OBJstream("projectCurveEdge_" + Foam::name(eIter++) + ".obj")
95  );
96  Info<< "Writing lines from straight-line start points"
97  << " to projected points to " << debugStr().name() << endl;
98  }
99 
100 
101  tmp<pointField> tpoints(new pointField(lambdas.size()));
102  pointField& points = tpoints.ref();
103 
104  const point& startPt = points_[start_];
105  const point& endPt = points_[end_];
106  const vector d = endPt-startPt;
107 
108  // Initial guess
109  forAll(lambdas, i)
110  {
111  points[i] = startPt+lambdas[i]*d;
112  }
113 
114  // Use special interpolation to keep initial guess on same position on
115  // surface
116  forAll(surfaces_, i)
117  {
118  if (isA<searchableExtrudedCircle>(geometry_[surfaces_[i]]))
119  {
120  const searchableExtrudedCircle& s =
121  refCast<const searchableExtrudedCircle>
122  (
123  geometry_[surfaces_[i]]
124  );
127  (
128  points[0],
129  points.last(),
130  scalarField(lambdas),
131  scalarField(points.size(), magSqr(d)),
132  nearInfo
133  );
134  forAll(nearInfo, i)
135  {
136  if (nearInfo[i].hit())
137  {
138  points[i] = nearInfo[i].hitPoint();
139  }
140  }
141 
142  break;
143  }
144  }
145 
146 
147 
148  // Upper limit for number of iterations
149  const label maxIter = 10;
150 
151  // Residual tolerance
152  const scalar relTol = 0.1;
153  const scalar absTol = 1e-4;
154 
155  scalar initialResidual = 0;
156 
157  for (label iter = 0; iter < maxIter; iter++)
158  {
159  // Do projection
160  {
161  List<pointConstraint> constraints(lambdas.size());
162  pointField start(points);
164  (
165  geometry_,
166  surfaces_,
167  start,
168  scalarField(start.size(), magSqr(d)),
169  points,
170  constraints
171  );
172 
173  // Reset start and end point
174  if (lambdas[0] < small)
175  {
176  points[0] = startPt;
177  }
178  if (lambdas.last() > 1 - small)
179  {
180  points.last() = endPt;
181  }
182 
183  if (debugStr.valid())
184  {
185  forAll(points, i)
186  {
187  debugStr().write(linePointRef(start[i], points[i]));
188  }
189  }
190  }
191 
192  // Calculate lambdas (normalised coordinate along edge)
193  scalarField projLambdas(points.size());
194  {
195  projLambdas[0] = 0;
196  for (label i = 1; i < points.size(); i++)
197  {
198  projLambdas[i] =
199  projLambdas[i-1] + mag(points[i] - points[i-1]);
200  }
201  projLambdas /= projLambdas.last();
202  }
203  linearInterpolationWeights interpolator(projLambdas);
204 
205  // Compare actual distances and move points (along straight line;
206  // not along surface)
207  vectorField residual(points.size(), vector::zero);
208  labelList indices;
209  scalarField weights;
210  for (label i = 1; i < points.size() - 1; i++)
211  {
212  interpolator.valueWeights(lambdas[i], indices, weights);
213 
214  point predicted = vector::zero;
215  forAll(indices, indexi)
216  {
217  predicted += weights[indexi]*points[indices[indexi]];
218  }
219  residual[i] = predicted - points[i];
220  }
221 
222  const scalar scalarResidual = sum(mag(residual));
223 
224  if (debug)
225  {
226  Pout<< "Iter:" << iter << " initialResidual:" << initialResidual
227  << " residual:" << scalarResidual << endl;
228  }
229 
230  if (scalarResidual < absTol*0.5*lambdas.size())
231  {
232  break;
233  }
234  else if (iter == 0)
235  {
236  initialResidual = scalarResidual;
237  }
238  else if (scalarResidual/initialResidual < relTol)
239  {
240  break;
241  }
242 
243 
244  if (debugStr.valid())
245  {
246  forAll(points, i)
247  {
248  const point predicted(points[i] + residual[i]);
249  debugStr().write(linePointRef(points[i], predicted));
250  }
251  }
252 
253  points += residual;
254  }
255 
256  return tpoints;
257 }
258 
259 
260 // ************************************************************************* //
#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
virtual point position(const scalar) const
Return the point positions corresponding to the curve parameters.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
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.
projectCurveEdge(const dictionary &dict, const label index, const searchableSurfaces &geometry, const pointField &points, Istream &)
Construct from Istream setting pointsList.
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
Tuple2< scalar, label > nearInfo
Private class for finding nearest.
Macros for easy insertion into run-time selection tables.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
virtual void findParametricNearest(const point &start, const point &end, const scalarField &lambdas, const scalarField &nearestDistSqr, List< pointIndexHit > &) const
Unique to parametric geometry: given points find.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
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.
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
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:318
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
Surface geometry with a tube shape, which can be used with snappyHexMesh. The geometry is formed from...
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