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-2023 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 faceSet;
91 
92 /*---------------------------------------------------------------------------*\
93  Class motionSmootherAlgo Declaration
94 \*---------------------------------------------------------------------------*/
95 
97 {
98  // Private class
99 
100  //- To synchronise displacements. We want max displacement since
101  // this is what is specified on pp and internal mesh will have
102  // zero displacement.
103  class maxMagEqOp
104  {
105 
106  public:
107 
108  void operator()(vector& x, const vector& y) const
109  {
110  for (direction i = 0; i < vector::nComponents; i++)
111  {
112  scalar magX = mag(x[i]);
113  scalar magY = mag(y[i]);
114 
115  if (magX < magY)
116  {
117  x[i] = y[i];
118  }
119  else if (magX == magY)
120  {
121  if (y[i] > x[i])
122  {
123  x[i] = y[i];
124  }
125  }
126  }
127  }
128  };
129 
130 
131  // Private Data
132 
133  //- Reference to polyMesh. Non-const since we move mesh.
134  polyMesh& mesh_;
135 
136  //- Reference to pointMesh
137  pointMesh& pMesh_;
138 
139  //- Reference to face subset of all adaptPatchIDs
141 
142  //- Displacement field
143  pointVectorField& displacement_;
144 
145  //- Scale factor for displacement
146  pointScalarField& scale_;
147 
148  //- Starting mesh position
149  pointField& oldPoints_;
150 
151 
152  // Internal data
153 
154  //- Indices of fixedValue patches that we're allowed to modify the
155  // displacement on.
156  const labelList adaptPatchIndices_;
157 
158  // Smoothing and checking parameters
159  dictionary paramDict_;
160 
161  //- Is mesh point on boundary or not
162  PackedBoolList isInternalPoint_;
163 
164  //- Is edge master (always except if on coupled boundary and on
165  // lower processor)
166  PackedBoolList isMasterEdge_;
167 
168 
169  // Private Member Functions
170 
171  //- Average of connected points.
172  template<class Type>
174  (
175  const PointField<Type>& fld,
176  const scalarField& edgeWeight
177  ) const;
178 
179  //- Average position of connected points.
180  template<class Type>
181  tmp<PointField<Type>> avgPositions
182  (
183  const PointField<Type>& fld,
184  const scalarField& edgeWeight
185  ) const;
186 
187  //- Check constraints
188  template<class Type>
189  static void checkConstraints
190  (
192  );
193 
194  //- Test synchronisation of generic field (not positions!) on points
195  template<class Type, class CombineOp>
196  void testSyncField
197  (
198  const Field<Type>&,
199  const CombineOp& cop,
200  const Type& zero,
201  const scalar maxMag
202  ) const;
203 
204  //- Test synchronisation of points
205  void testSyncPositions(const pointField&, const scalar maxMag) const;
206 
207  static void checkFld(const pointScalarField&);
208 
209  //- Get points used by given faces
210  labelHashSet getPoints(const labelHashSet&) const;
211 
212  //- Calculate per-edge weight
213  tmp<scalarField> calcEdgeWeights(const pointField&) const;
214 
215  //- Explicit smoothing and min on all affected internal points
216  void minSmooth
217  (
218  const scalarField& edgeWeights,
219  const PackedBoolList& isAffectedPoint,
220  const pointScalarField& fld,
221  pointScalarField& newFld
222  ) const;
223 
224  //- Same but only on selected points (usually patch points)
225  void minSmooth
226  (
227  const scalarField& edgeWeights,
228  const PackedBoolList& isAffectedPoint,
229  const labelList& meshPoints,
230  const pointScalarField& fld,
231  pointScalarField& newFld
232  ) const;
233 
234  //- Scale certain (internal) points of a field
235  void scaleField
236  (
237  const labelHashSet& pointLabels,
238  const scalar scale,
240  ) const;
241 
242  //- As above but points have to be in meshPoints as well
243  // (usually to scale patch points)
244  void scaleField
245  (
246  const labelList& meshPoints,
247  const labelHashSet& pointLabels,
248  const scalar scale,
250  ) const;
251 
252  //- Lower on internal points
253  void subtractField
254  (
255  const labelHashSet& pointLabels,
256  const scalar f,
258  ) const;
259 
260  //- As above but points have to be in meshPoints as well
261  // (usually to scale patch points)
262  void subtractField
263  (
264  const labelList& meshPoints,
265  const labelHashSet& pointLabels,
266  const scalar scale,
268  ) const;
269 
270  //- Helper function. Is point internal?
271  bool isInternalPoint(const label pointi) const;
272 
273  //- Given a set of faces that cause smoothing and a number of
274  // iterations determine the maximum set of points who are affected
275  // and the accordingly affected faces.
276  void getAffectedFacesAndPoints
277  (
278  const label nPointIter,
279  const faceSet& wrongFaces,
280 
281  labelList& affectedFaces,
282  PackedBoolList& isAffectedPoint
283  ) const;
284 
285 
286 public:
287 
288  ClassName("motionSmootherAlgo");
289 
290  // Constructors
291 
292  //- Construct from mesh, patches to work on and smoothing parameters.
294  (
295  polyMesh&,
296  pointMesh&,
297  indirectPrimitivePatch& pp, // 'outside' points
298  pointVectorField& displacement,
299  pointScalarField& scale,
300  pointField& oldPoints,
301  const labelList& adaptPatchIDs, // patches forming 'outside'
302  const dictionary& paramDict
303  );
304 
305  //- Disallow default bitwise copy construction
306  motionSmootherAlgo(const motionSmootherAlgo&) = delete;
307 
308 
309  //- Destructor
311 
312 
313  // Member Functions
314 
315  // Access
316 
317  //- Reference to mesh
318  const polyMesh& mesh() const;
319 
320  //- Reference to pointMesh
321  const pointMesh& pMesh() const;
322 
323  //- Reference to patch
324  const indirectPrimitivePatch& patch() const;
325 
326  //- Patch labels that are being adapted
327  const labelList& adaptPatchIndices() const;
328 
329  const dictionary& paramDict() const;
330 
331 
332 
333  // Edit
334 
335  //- Take over existing mesh position.
336  void correct();
337 
338 
339  //- Set patch fields on patchIDs to be consistent with
340  // all other boundary conditions
341  static void setDisplacementPatchFields
342  (
343  const labelList& patchIDs,
344  pointVectorField& pointDisplacement
345  );
346 
347  //- Set patch fields on displacement to be consistent with
348  // internal values.
350 
351  //- Set displacement field from displacement on patch points.
352  // Modify provided displacement to be consistent with actual
353  // boundary conditions on displacement. Note: resets the
354  // displacement to be 0 on coupled patches beforehand
355  // to make sure shared points
356  // partially on pp (on some processors) and partially not
357  // (on other processors) get the value from pp.
358  static void setDisplacement
359  (
360  const labelList& patchIDs,
361  const indirectPrimitivePatch& pp,
362  pointField& patchDisp,
363  pointVectorField& displacement
364  );
365 
366  void setDisplacement(pointField& patchDisp);
367 
368  //- Special correctBoundaryConditions which evaluates fixedValue
369  // patches first so they get overwritten with any constraint
370  // bc's.
372 
373  //- Apply optional point constraint (2d correction)
374  void modifyMotionPoints(pointField& newPoints) const;
375 
376  //- Get the current points (oldPoints+scale*displacement)
377  tmp<pointField> curPoints() const;
378 
379  //- Set the errorReduction (by how much to scale the displacement
380  // at error locations) parameter. Returns the old value.
381  // Set to 0 (so revert to old mesh) grows out one cell layer
382  // from error faces.
383  scalar setErrorReduction(const scalar);
384 
385  //- Move mesh with given scale. Return true if mesh ok or has
386  // less than nAllow errors, false
387  // otherwise and locally update scale. Smoothmesh=false means only
388  // patch points get moved.
389  // Parallel ok (as long as displacement field is consistent
390  // across patches)
391  bool scaleMesh
392  (
393  const labelList& checkFaces,
394  const bool smoothMesh = true,
395  const label nAllow = 0
396  );
397 
398  //- Move mesh (with baffles) with given scale.
399  bool scaleMesh
400  (
401  const labelList& checkFaces,
402  const List<labelPair>& baffles,
403  const bool smoothMesh = true,
404  const label nAllow = 0
405  );
406 
407  //- Move mesh with externally provided mesh constraints
408  bool scaleMesh
409  (
410  const labelList& checkFaces,
411  const List<labelPair>& baffles,
412  const dictionary& paramDict,
413  const dictionary& meshQualityDict,
414  const bool smoothMesh = true,
415  const label nAllow = 0
416  );
417 
418 
419  //- Update for new mesh geometry
420  void movePoints();
421 
422  //- Update for new mesh topology
423  void topoChange();
424 
425 
426  // Helper functions to manipulate displacement vector.
427 
428  //- Fully explicit smoothing of fields (not positions)
429  // of internal points with varying diffusivity.
430  template<class Type>
431  void smooth
432  (
433  const PointField<Type>& fld,
434  const scalarField& edgeWeight,
435  PointField<Type>& newFld
436  ) const;
437 
438 
439  // Member Operators
440 
441  //- Disallow default bitwise assignment
442  void operator=(const motionSmootherAlgo&) = delete;
443 };
444 
445 
446 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
447 
448 } // End namespace Foam
449 
450 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
451 
452 #ifdef NoRepository
454 #endif
455 
456 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
457 
458 #endif
459 
460 // ************************************************************************* //
scalar y
A bit-packed bool list.
static const direction nComponents
Number of components in this vector space.
Definition: VectorSpace.H:105
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
A list of face labels.
Definition: faceSet.H:51
Given a displacement moves the mesh by scaling the displacement back until there are no more mesh err...
void setDisplacementPatchFields()
Set patch fields on displacement to be consistent with.
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.
tmp< pointField > curPoints() const
Get the current points (oldPoints+scale*displacement)
void topoChange()
Update for new mesh topology.
scalar setErrorReduction(const scalar)
Set the errorReduction (by how much to scale the displacement.
ClassName("motionSmootherAlgo")
static void setDisplacement(const labelList &patchIDs, const indirectPrimitivePatch &pp, pointField &patchDisp, pointVectorField &displacement)
Set displacement field from displacement on patch points.
void correct()
Take over existing mesh position.
const polyMesh & mesh() const
Reference to mesh.
void smooth(const PointField< Type > &fld, const scalarField &edgeWeight, PointField< Type > &newFld) const
Fully explicit smoothing of fields (not positions)
bool scaleMesh(const labelList &checkFaces, const bool smoothMesh=true, const label nAllow=0)
Move mesh with given scale. Return true if mesh ok or has.
void movePoints()
Update for new mesh geometry.
const indirectPrimitivePatch & patch() const
Reference to patch.
const pointMesh & pMesh() const
Reference to pointMesh.
const labelList & adaptPatchIndices() const
Patch labels that are being adapted.
const dictionary & paramDict() const
void correctBoundaryConditions(pointVectorField &) const
Special correctBoundaryConditions which evaluates fixedValue.
void modifyMotionPoints(pointField &newPoints) const
Apply optional point constraint (2d correction)
void operator=(const motionSmootherAlgo &)=delete
Disallow default bitwise assignment.
Mesh representing a set of points created from polyMesh.
Definition: pointMesh.H:53
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
A class for managing temporary objects.
Definition: tmp.H:55
A class representing the concept of 0 used to avoid unnecessary manipulations for objects that are kn...
Definition: zero.H:50
Macro definitions for declaring ClassName(), NamespaceName(), etc.
gmvFile<< "tracers "<< particles.size()<< nl;{ pointField positions(particles.size());label particlei=0;forAllConstIter(Cloud< passiveParticle >, particles, iter) { positions[particlei++]=iter().position(mesh);} for(i=0;i< pTraits< point >::nComponents;i++) { forAll(positions, particlei) { gmvFile<< component(positions[particlei], i)<< ' ';} gmvFile<< nl;}}forAll(lagrangianScalarNames, i){ const word &name=lagrangianScalarNames[i];IOField< scalar > fld(IOobject(name, runTime.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Namespace for OpenFOAM.
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
dimensioned< scalar > mag(const dimensioned< Type > &)
uint8_t direction
Definition: direction.H:45
labelList f(nPoints)
labelList pointLabels(nPoints, -1)