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-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 
26 #include "projectFace.H"
27 #include "unitConversion.H"
29 #include "blockDescriptor.H"
30 #include "OBJstream.H"
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37 namespace blockFaces
38 {
41 }
42 }
43 
44 
45 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
46 
47 const Foam::searchableSurface& Foam::blockFaces::projectFace::lookupSurface
48 (
49  const searchableSurfaces& geometry,
50  Istream& is
51 ) const
52 {
53  word name(is);
54 
55  forAll(geometry, i)
56  {
57  if (geometry[i].name() == name)
58  {
59  return geometry[i];
60  }
61  }
62 
64  << "Cannot find surface " << name << " in geometry"
65  << exit(FatalIOError);
66 
67  return geometry[0];
68 }
69 
70 
71 Foam::label Foam::blockFaces::projectFace::index
72 (
73  const labelPair& n,
74  const labelPair& coord
75 ) const
76 {
77  return coord.first() + coord.second()*n.first();
78 }
79 
80 
81 void Foam::blockFaces::projectFace::calcLambdas
82 (
83  const labelPair& n,
84  const pointField& points,
85  scalarField& lambdaI,
86  scalarField& lambdaJ
87 ) const
88 {
89  lambdaI.setSize(points.size());
90  lambdaI = 0.0;
91  lambdaJ.setSize(points.size());
92  lambdaJ = 0.0;
93 
94  for (label i = 1; i < n.first(); i++)
95  {
96  for (label j = 1; j < n.second(); j++)
97  {
98  label ij = index(n, labelPair(i, j));
99  label iMin1j = index(n, labelPair(i-1, j));
100  lambdaI[ij] = lambdaI[iMin1j] + mag(points[ij]-points[iMin1j]);
101 
102  label ijMin1 = index(n, labelPair(i, j-1));
103  lambdaJ[ij] = lambdaJ[ijMin1] + mag(points[ij]-points[ijMin1]);
104  }
105  }
106 
107  for (label i = 1; i < n.first(); i++)
108  {
109  label ijLast = index(n, labelPair(i, n.second()-1));
110  for (label j = 1; j < n.second(); j++)
111  {
112  label ij = index(n, labelPair(i, j));
113  lambdaJ[ij] /= lambdaJ[ijLast];
114  }
115  }
116  for (label j = 1; j < n.second(); j++)
117  {
118  label iLastj = index(n, labelPair(n.first()-1, j));
119  for (label i = 1; i < n.first(); i++)
120  {
121  label ij = index(n, labelPair(i, j));
122  lambdaI[ij] /= lambdaI[iLastj];
123  }
124  }
125 }
126 
127 
128 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
129 
130 Foam::blockFaces::projectFace::projectFace
131 (
132  const dictionary& dict,
133  const label index,
134  const searchableSurfaces& geometry,
135  Istream& is
136 )
137 :
138  blockFace(dict, index, is),
139  surface_(lookupSurface(geometry, is))
140 {}
141 
142 
143 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
144 
146 (
147  const blockDescriptor& desc,
148  const label blockFacei,
149  pointField& points
150 ) const
151 {
152  // For debugging to tag the output
153  static label fIter = 0;
154 
155  autoPtr<OBJstream> debugStr;
156  if (debug)
157  {
158  debugStr.reset
159  (
160  new OBJstream("projectFace_" + Foam::name(fIter++) + ".obj")
161  );
162  Info<< "Face:" << blockFacei << " on block:" << desc.blockShape()
163  << " with verts:" << desc.vertices()
164  << " writing lines from start points"
165  << " to projected points to " << debugStr().name() << endl;
166  }
167 
168 
169  // Find out range of vertices in face
170  labelPair n(-1, -1);
171  switch (blockFacei)
172  {
173  case 0:
174  case 1:
175  {
176  n.first() = desc.density()[1] + 1;
177  n.second() = desc.density()[2] + 1;
178  }
179  break;
180 
181  case 2:
182  case 3:
183  {
184  n.first() = desc.density()[0] + 1;
185  n.second() = desc.density()[2] + 1;
186  }
187  break;
188 
189  case 4:
190  case 5:
191  {
192  n.first() = desc.density()[0] + 1;
193  n.second() = desc.density()[1] + 1;
194  }
195  break;
196  }
197 
198 
199  // Calculate initial normalised edge lengths (= u,v coordinates)
200  scalarField lambdaI(points.size(), 0.0);
201  scalarField lambdaJ(points.size(), 0.0);
202  calcLambdas(n, points, lambdaI, lambdaJ);
203 
204 
205  // Upper limit for number of iterations
206  const label maxIter = 10;
207  // Residual tolerance
208  const scalar relTol = 0.1;
209 
210  scalar initialResidual = 0.0;
211  scalar iResidual = 0.0;
212  scalar jResidual = 0.0;
213 
214  for (label iter = 0; iter < maxIter; iter++)
215  {
216  // Do projection
217  {
218  List<pointIndexHit> hits;
219  scalarField nearestDistSqr
220  (
221  points.size(),
222  magSqr(points[0] - points[points.size()-1])
223  );
224  surface_.findNearest(points, nearestDistSqr, hits);
225 
226  forAll(hits, i)
227  {
228  if (hits[i].hit())
229  {
230  const point& hitPt = hits[i].hitPoint();
231  if (debugStr.valid())
232  {
233  debugStr().write(linePointRef(points[i], hitPt));
234  }
235  points[i] = hitPt;
236  }
237  }
238  }
239 
240  if (debug)
241  {
242  Pout<< "Iter:" << iter << " initialResidual:" << initialResidual
243  << " iResidual+jResidual:" << iResidual+jResidual << endl;
244  }
245 
246 
247  if (iter > 0 && (iResidual+jResidual)/initialResidual < relTol)
248  {
249  break;
250  }
251 
252 
253  // Predict along i
254  vectorField residual(points.size(), vector::zero);
255 
256  // Work arrays for interpolation
257  labelList indices;
258  scalarField weights;
259  for (label j = 1; j < n.second()-1; j++)
260  {
261  // Calculate actual lamdba along constant j line
262  scalarField projLambdas(n.first());
263  projLambdas[0] = 0.0;
264  for (label i = 1; i < n.first(); i++)
265  {
266  label ij = index(n, labelPair(i, j));
267  label iMin1j = index(n, labelPair(i-1, j));
268  projLambdas[i] =
269  projLambdas[i-1]
270  +mag(points[ij]-points[iMin1j]);
271  }
272  projLambdas /= projLambdas.last();
273 
274  linearInterpolationWeights interpolator(projLambdas);
275 
276  for (label i = 1; i < n.first()-1; i++)
277  {
278  label ij = index(n, labelPair(i, j));
279 
280  interpolator.valueWeights(lambdaI[ij], indices, weights);
281 
282  point predicted = vector::zero;
283  forAll(indices, indexi)
284  {
285  label ptIndex = index(n, labelPair(indices[indexi], j));
286  predicted += weights[indexi]*points[ptIndex];
287  }
288  residual[ij] = predicted-points[ij];
289  }
290  }
291 
292  if (debugStr.valid())
293  {
294  forAll(points, i)
295  {
296  const point predicted(points[i] + residual[i]);
297  debugStr().write(linePointRef(points[i], predicted));
298  }
299  }
300 
301  iResidual = sum(mag(residual));
302 
303  // Update points before doing j. Note: is this needed? Complicates
304  // residual checking.
305  points += residual;
306 
307 
308  // Predict along j
309  residual = vector::zero;
310  for (label i = 1; i < n.first()-1; i++)
311  {
312  // Calculate actual lamdba along constant i line
313  scalarField projLambdas(n.second());
314  projLambdas[0] = 0.0;
315  for (label j = 1; j < n.second(); j++)
316  {
317  label ij = index(n, labelPair(i, j));
318  label ijMin1 = index(n, labelPair(i, j-1));
319  projLambdas[j] =
320  projLambdas[j-1]
321  +mag(points[ij]-points[ijMin1]);
322  }
323 
324  projLambdas /= projLambdas.last();
325 
326  linearInterpolationWeights interpolator(projLambdas);
327 
328  for (label j = 1; j < n.second()-1; j++)
329  {
330  label ij = index(n, labelPair(i, j));
331 
332  interpolator.valueWeights(lambdaJ[ij], indices, weights);
333 
334  point predicted = vector::zero;
335  forAll(indices, indexi)
336  {
337  label ptIndex = index(n, labelPair(i, indices[indexi]));
338  predicted += weights[indexi]*points[ptIndex];
339  }
340  residual[ij] = predicted-points[ij];
341  }
342  }
343 
344  if (debugStr.valid())
345  {
346  forAll(points, i)
347  {
348  const point predicted(points[i] + residual[i]);
349  debugStr().write(linePointRef(points[i], predicted));
350  }
351  }
352 
353  jResidual = sum(mag(residual));
354 
355  if (iter == 0)
356  {
357  initialResidual = iResidual + jResidual;
358  }
359 
360  points += residual;
361  }
362 }
363 
364 
365 // ************************************************************************* //
#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
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:137
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:60
defineTypeNameAndDebug(projectFace, 0)
Unit conversion functions.
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:163
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
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:146
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:99
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:331
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
addToRunTimeSelectionTable(blockFace, projectFace, Istream)
const Type & first() const
Return first.
Definition: Pair.H:87
Namespace for OpenFOAM.
IOerror FatalIOError