projectFace.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 
26 #include "projectFace.H"
27 #include "blockDescriptor.H"
29 #include "OBJstream.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 namespace blockFaces
37 {
40 }
41 }
42 
43 
44 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
45 
46 const Foam::searchableSurface& Foam::blockFaces::projectFace::lookupSurface
47 (
48  const searchableSurfaces& geometry,
49  Istream& is
50 ) const
51 {
52  const word name(is);
53 
54  forAll(geometry, i)
55  {
56  if (geometry[i].name() == name)
57  {
58  return geometry[i];
59  }
60  }
61 
63  << "Cannot find surface " << name << " in geometry"
64  << exit(FatalIOError);
65 
66  return geometry[0];
67 }
68 
69 
70 Foam::label Foam::blockFaces::projectFace::index
71 (
72  const labelPair& n,
73  const labelPair& coord
74 ) const
75 {
76  return coord.first() + coord.second()*n.first();
77 }
78 
79 
80 void Foam::blockFaces::projectFace::calcLambdas
81 (
82  const labelPair& n,
83  const pointField& points,
84  scalarField& lambdaI,
85  scalarField& lambdaJ
86 ) const
87 {
88  lambdaI.setSize(points.size());
89  lambdaI = 0;
90  lambdaJ.setSize(points.size());
91  lambdaJ = 0;
92 
93  for (label i = 1; i < n.first(); i++)
94  {
95  for (label j = 1; j < n.second(); j++)
96  {
97  const label ij = index(n, labelPair(i, j));
98  const label iMin1j = index(n, labelPair(i-1, j));
99  lambdaI[ij] = lambdaI[iMin1j] + mag(points[ij] - points[iMin1j]);
100 
101  const label ijMin1 = index(n, labelPair(i, j-1));
102  lambdaJ[ij] = lambdaJ[ijMin1] + mag(points[ij] - points[ijMin1]);
103  }
104  }
105 
106  for (label i = 1; i < n.first(); i++)
107  {
108  const label ijLast = index(n, labelPair(i, n.second()-1));
109  for (label j = 1; j < n.second(); j++)
110  {
111  label ij = index(n, labelPair(i, j));
112  lambdaJ[ij] /= lambdaJ[ijLast];
113  }
114  }
115  for (label j = 1; j < n.second(); j++)
116  {
117  const label iLastj = index(n, labelPair(n.first()-1, j));
118  for (label i = 1; i < n.first(); i++)
119  {
120  label ij = index(n, labelPair(i, j));
121  lambdaI[ij] /= lambdaI[iLastj];
122  }
123  }
124 }
125 
126 
127 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
128 
130 (
131  const dictionary& dict,
132  const label index,
133  const searchableSurfaces& geometry,
134  Istream& is
135 )
136 :
137  blockFace(dict, index, is),
138  surface_(lookupSurface(geometry, is))
139 {}
140 
141 
142 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
143 
145 (
146  const blockDescriptor& desc,
147  const label blockFacei,
148  pointField& points
149 ) const
150 {
151  // For debugging to tag the output
152  static label fIter = 0;
153 
154  autoPtr<OBJstream> debugStr;
155  if (debug)
156  {
157  debugStr.reset
158  (
159  new OBJstream("projectFace_" + Foam::name(fIter++) + ".obj")
160  );
161  Info<< "Face:" << blockFacei << " on block:" << desc.blockShape()
162  << " with verts:" << desc.vertices()
163  << " writing lines from start points"
164  << " to projected points to " << debugStr().name() << endl;
165  }
166 
167 
168  // Find out range of vertices in face
169  labelPair n(-1, -1);
170  switch (blockFacei)
171  {
172  case 0:
173  case 1:
174  {
175  n.first() = desc.density()[1] + 1;
176  n.second() = desc.density()[2] + 1;
177  }
178  break;
179 
180  case 2:
181  case 3:
182  {
183  n.first() = desc.density()[0] + 1;
184  n.second() = desc.density()[2] + 1;
185  }
186  break;
187 
188  case 4:
189  case 5:
190  {
191  n.first() = desc.density()[0] + 1;
192  n.second() = desc.density()[1] + 1;
193  }
194  break;
195  }
196 
197 
198  // Calculate initial normalised edge lengths (= u,v coordinates)
199  scalarField lambdaI(points.size(), 0);
200  scalarField lambdaJ(points.size(), 0);
201  calcLambdas(n, points, lambdaI, lambdaJ);
202 
203 
204  // Upper limit for number of iterations
205  const label maxIter = 10;
206 
207  // Residual tolerance
208  const scalar relTol = 0.1;
209  const scalar absTol = 1e-4;
210 
211  scalar initialResidual = 0;
212 
213  for (label iter = 0; iter < maxIter; iter++)
214  {
215  // Do projection
216  {
217  List<pointIndexHit> hits;
218  const scalarField nearestDistSqr
219  (
220  points.size(),
221  magSqr(points[0] - points[points.size()-1])
222  );
223  surface_.findNearest(points, nearestDistSqr, hits);
224 
225  forAll(hits, i)
226  {
227  if (hits[i].hit())
228  {
229  const point& hitPt = hits[i].hitPoint();
230  if (debugStr.valid())
231  {
232  debugStr().write(linePointRef(points[i], hitPt));
233  }
234  points[i] = hitPt;
235  }
236  }
237  }
238 
239  // Predict along i
240  vectorField residual(points.size(), vector::zero);
241 
242  // Work arrays for interpolation
243  labelList indices;
244  scalarField weights;
245  for (label j = 1; j < n.second()-1; j++)
246  {
247  // Calculate actual lamdba along constant j line
248  scalarField projLambdas(n.first());
249  projLambdas[0] = 0;
250  for (label i = 1; i < n.first(); i++)
251  {
252  const label ij = index(n, labelPair(i, j));
253  const label iMin1j = index(n, labelPair(i-1, j));
254  projLambdas[i] =
255  projLambdas[i-1]
256  + mag(points[ij] - points[iMin1j]);
257  }
258  projLambdas /= projLambdas.last();
259 
260  linearInterpolationWeights interpolator(projLambdas);
261 
262  for (label i = 1; i < n.first()-1; i++)
263  {
264  const label ij = index(n, labelPair(i, j));
265 
266  interpolator.valueWeights(lambdaI[ij], indices, weights);
267 
268  point predicted = vector::zero;
269  forAll(indices, indexi)
270  {
271  const label ptIndex =
272  index(n, labelPair(indices[indexi], j));
273  predicted += weights[indexi]*points[ptIndex];
274  }
275  residual[ij] = predicted - points[ij];
276  }
277  }
278 
279  if (debugStr.valid())
280  {
281  forAll(points, i)
282  {
283  const point predicted(points[i] + residual[i]);
284  debugStr().write(linePointRef(points[i], predicted));
285  }
286  }
287 
288  scalar scalarResidual = sum(mag(residual));
289 
290  // Update points before doing j. Note: is this needed? Complicates
291  // residual checking.
292  points += residual;
293 
294 
295  // Predict along j
296  residual = vector::zero;
297  for (label i = 1; i < n.first()-1; i++)
298  {
299  // Calculate actual lamdba along constant i line
300  scalarField projLambdas(n.second());
301  projLambdas[0] = 0;
302  for (label j = 1; j < n.second(); j++)
303  {
304  const label ij = index(n, labelPair(i, j));
305  const label ijMin1 = index(n, labelPair(i, j-1));
306  projLambdas[j] =
307  projLambdas[j-1]
308  + mag(points[ij] - points[ijMin1]);
309  }
310 
311  projLambdas /= projLambdas.last();
312 
313  linearInterpolationWeights interpolator(projLambdas);
314 
315  for (label j = 1; j < n.second()-1; j++)
316  {
317  const label ij = index(n, labelPair(i, j));
318 
319  interpolator.valueWeights(lambdaJ[ij], indices, weights);
320 
321  point predicted = vector::zero;
322  forAll(indices, indexi)
323  {
324  const label ptIndex =
325  index(n, labelPair(i, indices[indexi]));
326  predicted += weights[indexi]*points[ptIndex];
327  }
328  residual[ij] = predicted - points[ij];
329  }
330  }
331 
332  if (debugStr.valid())
333  {
334  forAll(points, i)
335  {
336  const point predicted(points[i] + residual[i]);
337  debugStr().write(linePointRef(points[i], predicted));
338  }
339  }
340 
341  scalarResidual += sum(mag(residual));
342 
343  if (debug)
344  {
345  Pout<< "Iter:" << iter << " initialResidual:" << initialResidual
346  << " scalarResidual:" << scalarResidual << endl;
347  }
348 
349  if (scalarResidual < absTol*lambdaI.size())
350  {
351  break;
352  }
353  else if (iter == 0)
354  {
355  initialResidual = scalarResidual;
356  }
357  else if (scalarResidual/initialResidual < relTol)
358  {
359  break;
360  }
361 
362  points += residual;
363  }
364 }
365 
366 
367 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
const Vector< label > & density() const
Return the mesh density (number of cells) in the i,j,k directions.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
defineTypeNameAndDebug(projectFace, 0)
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
Definition: Istream.H:57
Takes the description of the block and the list of curved edges and creates a list of points on edges...
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
Base class of (analytical or triangulated) surface. Encapsulates all the search routines. WIP.
Macros for easy insertion into run-time selection tables.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
const cellShape & blockShape() const
Return the block shape.
virtual bool valueWeights(const scalar t, labelList &indices, scalarField &weights) const
Calculate weights and indices to calculate t from samples.
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
A class for handling words, derived from string.
Definition: word.H:59
virtual void project(const blockDescriptor &, const label blockFacei, pointField &points) const
Project the given points onto the surface.
Definition: projectFace.C:145
Container for searchableSurfaces.
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
bool valid() const
Return true if the autoPtr valid (ie, the pointer is set)
Definition: autoPtrI.H:83
dimensioned< scalar > magSqr(const dimensioned< Type > &)
OFstream which keeps track of vertices.
Definition: OBJstream.H:53
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
const pointField & vertices() const
Reference to point field defining the block mesh.
Define a curved face.
Definition: blockFace.H:55
const Type & second() const
Return second.
Definition: Pair.H:110
void setSize(const label)
Reset size of List.
Definition: List.C:281
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:335
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Projects the given set of face points onto the selected surface of the geometry provided as a searcha...
Definition: projectFace.H:51
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
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
addToRunTimeSelectionTable(blockFace, projectFace, Istream)
const Type & first() const
Return first.
Definition: Pair.H:98
projectFace(const dictionary &dict, const label index, const searchableSurfaces &geometry, Istream &)
Construct from Istream setting pointsList.
Definition: projectFace.C:130
Namespace for OpenFOAM.
IOerror FatalIOError