motionSmootherAlgo.H
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) 2011-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 Class
25  Foam::motionSmootherAlgo
26 
27 Description
28  Given a displacement moves the mesh by scaling the displacement back
29  until there are no more mesh errors.
30 
31  Holds displacement field (read upon construction since need boundary
32  conditions) and scaling factor and optional patch number on which to
33  scale back displacement.
34 
35  E.g.
36  \verbatim
37  // Construct iterative mesh mover.
38  motionSmoother meshMover(mesh, labelList(1, patchi));
39 
40  // Set desired displacement:
41  meshMover.displacement() = ..
42 
43  for (label iter = 0; iter < maxIter; iter++)
44  {
45  if (meshMover.scaleMesh(true))
46  {
47  Info<< "Successfully moved mesh" << endl;
48  return true;
49  }
50  }
51  \endverbatim
52 
53 Note
54  - Shared points (parallel): a processor can have points which are part of
55  pp on another processor but have no pp itself (i.e. it has points
56  and/or edges but no faces of pp). Hence we have to be careful when e.g.
57  synchronising displacements that the value from the processor which has
58  faces of pp get priority. This is currently handled in setDisplacement
59  by resetting the internal displacement to zero before doing anything
60  else. The combine operator used will give preference to non-zero
61  values.
62 
63  - Various routines take baffles. These are sets of boundary faces that
64  are treated as a single internal face. This is a hack used to apply
65  movement to internal faces.
66 
67  - Mesh constraints are looked up from the supplied dictionary. (uses
68  recursive lookup)
69 
70 SourceFiles
71  motionSmootherAlgo.C
72  motionSmootherAlgoTemplates.C
73 
74 \*---------------------------------------------------------------------------*/
75 
76 #ifndef motionSmootherAlgo_H
77 #define motionSmootherAlgo_H
78 
79 #include "pointFields.H"
80 #include "HashSet.H"
81 #include "PackedBoolList.H"
82 #include "indirectPrimitivePatch.H"
83 #include "className.H"
84 
85 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
86 
87 namespace Foam
88 {
89 
90 class polyMeshGeometry;
91 class faceSet;
92 
93 /*---------------------------------------------------------------------------*\
94  Class motionSmootherAlgo Declaration
95 \*---------------------------------------------------------------------------*/
96 
98 {
99  // Private class
100 
101  //- To synchronise displacements. We want max displacement since
102  // this is what is specified on pp and internal mesh will have
103  // zero displacement.
104  class maxMagEqOp
105  {
106 
107  public:
108 
109  void operator()(vector& x, const vector& y) const
110  {
111  for (direction i = 0; i < vector::nComponents; i++)
112  {
113  scalar magX = mag(x[i]);
114  scalar magY = mag(y[i]);
115 
116  if (magX < magY)
117  {
118  x[i] = y[i];
119  }
120  else if (magX == magY)
121  {
122  if (y[i] > x[i])
123  {
124  x[i] = y[i];
125  }
126  }
127  }
128  }
129  };
130 
131 
132  // Private Data
133 
134  //- Reference to polyMesh. Non-const since we move mesh.
135  polyMesh& mesh_;
136 
137  //- Reference to pointMesh
138  pointMesh& pMesh_;
139 
140  //- Reference to face subset of all adaptPatchIDs
142 
143  //- Displacement field
144  pointVectorField& displacement_;
145 
146  //- Scale factor for displacement
147  pointScalarField& scale_;
148 
149  //- Starting mesh position
150  pointField& oldPoints_;
151 
152 
153  // Internal data
154 
155  //- Indices of fixedValue patches that we're allowed to modify the
156  // displacement on.
157  const labelList adaptPatchIDs_;
158 
159  // Smoothing and checking parameters
160  dictionary paramDict_;
161 
162  //- Is mesh point on boundary or not
163  PackedBoolList isInternalPoint_;
164 
165  //- Is edge master (always except if on coupled boundary and on
166  // lower processor)
167  PackedBoolList isMasterEdge_;
168 
169 
170  // Private Member Functions
171 
172  //- Average of connected points.
173  template<class Type>
175  (
177  const scalarField& edgeWeight
178  ) const;
179 
180  //- Average position of connected points.
181  template<class Type>
183  (
185  const scalarField& edgeWeight
186  ) const;
187 
188  //- Check constraints
189  template<class Type>
190  static void checkConstraints
191  (
193  );
194 
195  //- Test synchronisation of generic field (not positions!) on points
196  template<class Type, class CombineOp>
197  void testSyncField
198  (
199  const Field<Type>&,
200  const CombineOp& cop,
201  const Type& zero,
202  const scalar maxMag
203  ) const;
204 
205  //- Test synchronisation of points
206  void testSyncPositions(const pointField&, const scalar maxMag) const;
207 
208  static void checkFld(const pointScalarField&);
209 
210  //- Get points used by given faces
211  labelHashSet getPoints(const labelHashSet&) const;
212 
213  //- Calculate per-edge weight
214  tmp<scalarField> calcEdgeWeights(const pointField&) const;
215 
216  //- Explicit smoothing and min on all affected internal points
217  void minSmooth
218  (
219  const scalarField& edgeWeights,
220  const PackedBoolList& isAffectedPoint,
221  const pointScalarField& fld,
222  pointScalarField& newFld
223  ) const;
224 
225  //- Same but only on selected points (usually patch points)
226  void minSmooth
227  (
228  const scalarField& edgeWeights,
229  const PackedBoolList& isAffectedPoint,
230  const labelList& meshPoints,
231  const pointScalarField& fld,
232  pointScalarField& newFld
233  ) const;
234 
235  //- Scale certain (internal) points of a field
236  void scaleField
237  (
238  const labelHashSet& pointLabels,
239  const scalar scale,
241  ) const;
242 
243  //- As above but points have to be in meshPoints as well
244  // (usually to scale patch points)
245  void scaleField
246  (
247  const labelList& meshPoints,
248  const labelHashSet& pointLabels,
249  const scalar scale,
251  ) const;
252 
253  //- Lower on internal points
254  void subtractField
255  (
256  const labelHashSet& pointLabels,
257  const scalar f,
259  ) const;
260 
261  //- As above but points have to be in meshPoints as well
262  // (usually to scale patch points)
263  void subtractField
264  (
265  const labelList& meshPoints,
266  const labelHashSet& pointLabels,
267  const scalar scale,
269  ) const;
270 
271  //- Helper function. Is point internal?
272  bool isInternalPoint(const label pointi) const;
273 
274  //- Given a set of faces that cause smoothing and a number of
275  // iterations determine the maximum set of points who are affected
276  // and the accordingly affected faces.
277  void getAffectedFacesAndPoints
278  (
279  const label nPointIter,
280  const faceSet& wrongFaces,
281 
282  labelList& affectedFaces,
283  PackedBoolList& isAffectedPoint
284  ) const;
285 
286 
287 public:
288 
289  ClassName("motionSmootherAlgo");
290 
291  // Constructors
292 
293  //- Construct from mesh, patches to work on and smoothing parameters.
295  (
296  polyMesh&,
297  pointMesh&,
298  indirectPrimitivePatch& pp, // 'outside' points
299  pointVectorField& displacement,
300  pointScalarField& scale,
301  pointField& oldPoints,
302  const labelList& adaptPatchIDs, // patches forming 'outside'
303  const dictionary& paramDict
304  );
305 
306  //- Disallow default bitwise copy construction
307  motionSmootherAlgo(const motionSmootherAlgo&) = delete;
308 
309 
310  //- Destructor
312 
313 
314  // Member Functions
315 
316  // Access
317 
318  //- Reference to mesh
319  const polyMesh& mesh() const;
320 
321  //- Reference to pointMesh
322  const pointMesh& pMesh() const;
323 
324  //- Reference to patch
325  const indirectPrimitivePatch& patch() const;
326 
327  //- Patch labels that are being adapted
328  const labelList& adaptPatchIDs() const;
329 
330  const dictionary& paramDict() const;
331 
332 
333 
334  // Edit
335 
336  //- Take over existing mesh position.
337  void correct();
338 
339 
340  //- Set patch fields on patchIDs to be consistent with
341  // all other boundary conditions
342  static void setDisplacementPatchFields
343  (
344  const labelList& patchIDs,
345  pointVectorField& pointDisplacement
346  );
347 
348  //- Set patch fields on displacement to be consistent with
349  // internal values.
351 
352  //- Set displacement field from displacement on patch points.
353  // Modify provided displacement to be consistent with actual
354  // boundary conditions on displacement. Note: resets the
355  // displacement to be 0 on coupled patches beforehand
356  // to make sure shared points
357  // partially on pp (on some processors) and partially not
358  // (on other processors) get the value from pp.
359  static void setDisplacement
360  (
361  const labelList& patchIDs,
362  const indirectPrimitivePatch& pp,
363  pointField& patchDisp,
364  pointVectorField& displacement
365  );
366 
367  void setDisplacement(pointField& patchDisp);
368 
369  //- Special correctBoundaryConditions which evaluates fixedValue
370  // patches first so they get overwritten with any constraint
371  // bc's.
373 
374  //- Apply optional point constraint (2d correction)
375  void modifyMotionPoints(pointField& newPoints) const;
376 
377  //- Get the current points (oldPoints+scale*displacement)
378  tmp<pointField> curPoints() const;
379 
380  //- Set the errorReduction (by how much to scale the displacement
381  // at error locations) parameter. Returns the old value.
382  // Set to 0 (so revert to old mesh) grows out one cell layer
383  // from error faces.
384  scalar setErrorReduction(const scalar);
385 
386  //- Move mesh with given scale. Return true if mesh ok or has
387  // less than nAllow errors, false
388  // otherwise and locally update scale. Smoothmesh=false means only
389  // patch points get moved.
390  // Parallel ok (as long as displacement field is consistent
391  // across patches)
392  bool scaleMesh
393  (
394  labelList& checkFaces,
395  const bool smoothMesh = true,
396  const label nAllow = 0
397  );
398 
399  //- Move mesh (with baffles) with given scale.
400  bool scaleMesh
401  (
402  labelList& checkFaces,
403  const List<labelPair>& baffles,
404  const bool smoothMesh = true,
405  const label nAllow = 0
406  );
407 
408  //- Move mesh with externally provided mesh constraints
409  bool scaleMesh
410  (
411  labelList& checkFaces,
412  const List<labelPair>& baffles,
413  const dictionary& paramDict,
414  const dictionary& meshQualityDict,
415  const bool smoothMesh = true,
416  const label nAllow = 0
417  );
418 
419 
420  //- Update for new mesh geometry
421  void movePoints();
422 
423  //- Update for new mesh topology
424  void updateMesh();
425 
426 
427  //- Check mesh with mesh settings in dict. Collects incorrect faces
428  // in set. Returns true if one or more faces in error.
429  // Parallel ok.
430  static bool checkMesh
431  (
432  const bool report,
433  const polyMesh& mesh,
434  const dictionary& dict,
435  labelHashSet& wrongFaces
436  );
437 
438  //- Check (subset of mesh) with mesh settings in dict.
439  // Collects incorrect faces in set. Returns true if one
440  // or more faces in error. Parallel ok.
441  static bool checkMesh
442  (
443  const bool report,
444  const polyMesh& mesh,
445  const dictionary& dict,
446  const labelList& checkFaces,
447  labelHashSet& wrongFaces
448  );
449 
450  //- Check (subset of mesh including baffles) with mesh settings
451  // in dict. Collects incorrect faces in set. Returns true if one
452  // or more faces in error. Parallel ok.
453  static bool checkMesh
454  (
455  const bool report,
456  const polyMesh& mesh,
457  const dictionary& dict,
458  const labelList& checkFaces,
459  const List<labelPair>& baffles,
460  labelHashSet& wrongFaces
461  );
462 
463  //- Check part of mesh with mesh settings in dict.
464  // Collects incorrect faces in set. Returns true if one or
465  // more faces in error. Parallel ok.
466  static bool checkMesh
467  (
468  const bool report,
469  const dictionary& dict,
470  const polyMeshGeometry&,
471  const labelList& checkFaces,
472  labelHashSet& wrongFaces
473  );
474 
475  //- Check part of mesh including baffles with mesh settings in dict.
476  // Collects incorrect faces in set. Returns true if one or
477  // more faces in error. Parallel ok.
478  static bool checkMesh
479  (
480  const bool report,
481  const dictionary& dict,
482  const polyMeshGeometry&,
483  const labelList& checkFaces,
484  const List<labelPair>& baffles,
485  labelHashSet& wrongFaces
486  );
487 
488  // Helper functions to manipulate displacement vector.
489 
490  //- Fully explicit smoothing of fields (not positions)
491  // of internal points with varying diffusivity.
492  template<class Type>
493  void smooth
494  (
496  const scalarField& edgeWeight,
498  ) const;
499 
500 
501  // Member Operators
502 
503  //- Disallow default bitwise assignment
504  void operator=(const motionSmootherAlgo&) = delete;
505 };
506 
507 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
508 
509 } // End namespace Foam
510 
511 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
512 
513 #ifdef NoRepository
515 #endif
516 
517 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
518 
519 #endif
520 
521 // ************************************************************************* //
static bool checkMesh(const bool report, const polyMesh &mesh, const dictionary &dict, labelHashSet &wrongFaces)
Check mesh with mesh settings in dict. Collects incorrect faces.
dictionary dict
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
const polyMesh & mesh() const
Reference to mesh.
A list of face labels.
Definition: faceSet.H:48
motionSmootherAlgo(polyMesh &, pointMesh &, indirectPrimitivePatch &pp, pointVectorField &displacement, pointScalarField &scale, pointField &oldPoints, const labelList &adaptPatchIDs, const dictionary &paramDict)
Construct from mesh, patches to work on and smoothing parameters.
void operator=(const motionSmootherAlgo &)=delete
Disallow default bitwise assignment.
labelList pointLabels(nPoints, -1)
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
static void setDisplacement(const labelList &patchIDs, const indirectPrimitivePatch &pp, pointField &patchDisp, pointVectorField &displacement)
Set displacement field from displacement on patch points.
void correctBoundaryConditions(pointVectorField &) const
Special correctBoundaryConditions which evaluates fixedValue.
bool scaleMesh(labelList &checkFaces, const bool smoothMesh=true, const label nAllow=0)
Move mesh with given scale. Return true if mesh ok or has.
uint8_t direction
Definition: direction.H:45
const pointMesh & pMesh() const
Reference to pointMesh.
const dictionary & paramDict() const
Generic GeometricField class.
void movePoints()
Update for new mesh geometry.
Mesh representing a set of points created from polyMesh.
Definition: pointMesh.H:48
void correct()
Take over existing mesh position.
scalar y
A list of faces which address into the list of points.
static const direction nComponents
Number of components in this vector space.
Definition: VectorSpace.H:99
scalar setErrorReduction(const scalar)
Set the errorReduction (by how much to scale the displacement.
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){ const word &name=lagrangianScalarNames[i];IOField< scalar > fld(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
void setDisplacementPatchFields()
Set patch fields on displacement to be consistent with.
Updateable mesh geometry and checking routines.
void modifyMotionPoints(pointField &newPoints) const
Apply optional point constraint (2d correction)
void smooth(const GeometricField< Type, pointPatchField, pointMesh > &fld, const scalarField &edgeWeight, GeometricField< Type, pointPatchField, pointMesh > &newFld) const
Fully explicit smoothing of fields (not positions)
const indirectPrimitivePatch & patch() const
Reference to patch.
labelList f(nPoints)
A bit-packed bool list.
tmp< pointField > curPoints() const
Get the current points (oldPoints+scale*displacement)
ClassName("motionSmootherAlgo")
Macro definitions for declaring ClassName(), NamespaceName(), etc.
A class representing the concept of 0 used to avoid unnecessary manipulations for objects that are kn...
Definition: zero.H:49
dimensioned< scalar > mag(const dimensioned< Type > &)
void updateMesh()
Update for new mesh topology.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
const labelList & adaptPatchIDs() const
Patch labels that are being adapted.
A class for managing temporary objects.
Definition: PtrList.H:53
Given a displacement moves the mesh by scaling the displacement back until there are no more mesh err...
Namespace for OpenFOAM.