displacementLayeredMotionMotionSolver.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) 2011-2022 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 
29 #include "pointFields.H"
30 #include "PointEdgeWave.H"
31 #include "syncTools.H"
32 #include "Table.H"
33 #include "pointConstraints.H"
34 #include "polyTopoChangeMap.H"
35 
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 
38 namespace Foam
39 {
40  defineTypeNameAndDebug(displacementLayeredMotionMotionSolver, 0);
41 
43  (
44  motionSolver,
45  displacementLayeredMotionMotionSolver,
46  dictionary
47  );
48 }
49 
50 
51 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
52 
53 void Foam::displacementLayeredMotionMotionSolver::calcZoneMask
54 (
55  const label cellZoneI,
56  PackedBoolList& isZonePoint,
57  PackedBoolList& isZoneEdge
58 ) const
59 {
60  if (cellZoneI == -1)
61  {
62  isZonePoint.setSize(mesh().nPoints());
63  isZonePoint = 1;
64 
65  isZoneEdge.setSize(mesh().nEdges());
66  isZoneEdge = 1;
67  }
68  else
69  {
70  const cellZone& cz = mesh().cellZones()[cellZoneI];
71 
72  label nPoints = 0;
73  forAll(cz, i)
74  {
75  const labelList& cPoints = mesh().cellPoints(cz[i]);
76  forAll(cPoints, cPointi)
77  {
78  if (!isZonePoint[cPoints[cPointi]])
79  {
80  isZonePoint[cPoints[cPointi]] = 1;
81  nPoints++;
82  }
83  }
84  }
86  (
87  mesh(),
88  isZonePoint,
89  orEqOp<unsigned int>(),
90  0
91  );
92 
93 
94  // Mark edge inside cellZone
95  label nEdges = 0;
96  forAll(cz, i)
97  {
98  const labelList& cEdges = mesh().cellEdges(cz[i]);
99  forAll(cEdges, cEdgeI)
100  {
101  if (!isZoneEdge[cEdges[cEdgeI]])
102  {
103  isZoneEdge[cEdges[cEdgeI]] = 1;
104  nEdges++;
105  }
106  }
107  }
109  (
110  mesh(),
111  isZoneEdge,
112  orEqOp<unsigned int>(),
113  0
114  );
115 
116  if (debug)
117  {
118  Info<< "On cellZone " << cz.name()
119  << " marked " << returnReduce(nPoints, sumOp<label>())
120  << " points and " << returnReduce(nEdges, sumOp<label>())
121  << " edges." << endl;
122  }
123  }
124 }
125 
126 
127 // Find distance to starting point
128 void Foam::displacementLayeredMotionMotionSolver::walkStructured
129 (
130  const label cellZoneI,
131  const PackedBoolList& isZonePoint,
132  const PackedBoolList& isZoneEdge,
133  const labelList& seedPoints,
134  const vectorField& seedData,
135  scalarField& distance,
136  vectorField& data
137 ) const
138 {
139  List<pointEdgeStructuredWalk> seedInfo(seedPoints.size());
140 
141  forAll(seedPoints, i)
142  {
143  seedInfo[i] = pointEdgeStructuredWalk
144  (
145  points0()[seedPoints[i]], // location of data
146  points0()[seedPoints[i]], // previous location
147  0.0,
148  seedData[i]
149  );
150  }
151 
152  // Current info on points
153  List<pointEdgeStructuredWalk> allPointInfo(mesh().nPoints());
154 
155  // Mark points inside cellZone.
156  // Note that we use points0, not mesh.points()
157  // so as not to accumulate errors.
158  forAll(isZonePoint, pointi)
159  {
160  if (isZonePoint[pointi])
161  {
162  allPointInfo[pointi] = pointEdgeStructuredWalk
163  (
164  points0()[pointi], // location of data
165  vector::max, // not valid
166  0.0,
167  Zero // passive data
168  );
169  }
170  }
171 
172  // Current info on edges
173  List<pointEdgeStructuredWalk> allEdgeInfo(mesh().nEdges());
174 
175  // Mark edges inside cellZone
176  forAll(isZoneEdge, edgeI)
177  {
178  if (isZoneEdge[edgeI])
179  {
180  allEdgeInfo[edgeI] = pointEdgeStructuredWalk
181  (
182  mesh().edges()[edgeI].centre(points0()), // location of data
183  vector::max, // not valid
184  0.0,
185  Zero
186  );
187  }
188  }
189 
190  // Walk
191  PointEdgeWave<pointEdgeStructuredWalk> wallCalc
192  (
193  mesh(),
194  seedPoints,
195  seedInfo,
196 
197  allPointInfo,
198  allEdgeInfo,
199  mesh().globalData().nTotalPoints() // max iterations
200  );
201 
202  // Extract distance and passive data
203  forAll(allPointInfo, pointi)
204  {
205  if (isZonePoint[pointi])
206  {
207  distance[pointi] = allPointInfo[pointi].dist();
208  data[pointi] = allPointInfo[pointi].data();
209  }
210  }
211 }
212 
213 
214 // Evaluate faceZone patch
216 Foam::displacementLayeredMotionMotionSolver::faceZoneEvaluate
217 (
218  const faceZone& fz,
219  const labelList& meshPoints,
220  const dictionary& dict,
221  const PtrList<pointVectorField>& patchDisp,
222  const label patchi
223 ) const
224 {
225  tmp<vectorField> tfld(new vectorField(meshPoints.size()));
226  vectorField& fld = tfld.ref();
227 
228  const word type(dict.lookup("type"));
229 
230  if (type == "fixedValue")
231  {
232  fld = vectorField("value", dict, meshPoints.size());
233  }
234  else if (type == "timeVaryingUniformFixedValue")
235  {
236  Function1s::Table<vector> timeSeries(word::null, dict);
237 
238  fld = timeSeries.value(mesh().time().userTimeValue());
239  }
240  else if (type == "slip")
241  {
242  if ((patchi % 2) != 1)
243  {
245  << "FaceZone:" << fz.name()
246  << exit(FatalError);
247  }
248  // Use field set by previous bc
249  fld = vectorField(patchDisp[patchi - 1], meshPoints);
250  }
251  else if (type == "follow")
252  {
253  // Only on boundary faces - follow boundary conditions
254  fld = vectorField(pointDisplacement_, meshPoints);
255  }
256  else if (type == "uniformFollow")
257  {
258  // Reads name of name of patch. Then get average point displacement on
259  // patch. That becomes the value of fld.
260  const word patchName(dict.lookup("patch"));
261  label patchID = mesh().boundaryMesh().findPatchID(patchName);
262  pointField pdf
263  (
264  pointDisplacement_.boundaryField()[patchID].patchInternalField()
265  );
266  fld = gAverage(pdf);
267  }
268  else
269  {
271  << "Unknown faceZonePatch type " << type << " for faceZone "
272  << fz.name() << exit(FatalError);
273  }
274  return tfld;
275 }
276 
277 
278 void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
279 (
280  const label cellZoneI,
281  const dictionary& zoneDict
282 )
283 {
284  PackedBoolList isZonePoint(mesh().nPoints());
285  PackedBoolList isZoneEdge(mesh().nEdges());
286  calcZoneMask(cellZoneI, isZonePoint, isZoneEdge);
287 
288  const dictionary& patchesDict = zoneDict.subDict("boundaryField");
289 
290  if (patchesDict.size() != 2)
291  {
293  << "Two faceZones (patches) must be specified per cellZone. "
294  << " cellZone:" << cellZoneI
295  << " patches:" << patchesDict.toc()
296  << exit(FatalError);
297  }
298 
299  PtrList<scalarField> patchDist(patchesDict.size());
300  PtrList<pointVectorField> patchDisp(patchesDict.size());
301 
302  // Allocate the fields
303  label patchi = 0;
304  forAllConstIter(dictionary, patchesDict, patchiter)
305  {
306  const word& faceZoneName = patchiter().keyword();
307  label zoneI = mesh().faceZones().findZoneID(faceZoneName);
308  if (zoneI == -1)
309  {
311  << "Cannot find faceZone " << faceZoneName
312  << endl << "Valid zones are " << mesh().faceZones().names()
313  << exit(FatalError);
314  }
315 
316  // Determine the points of the faceZone within the cellZone
317  const faceZone& fz = mesh().faceZones()[zoneI];
318 
319  patchDist.set(patchi, new scalarField(mesh().nPoints()));
320  patchDisp.set
321  (
322  patchi,
323  new pointVectorField
324  (
325  IOobject
326  (
327  mesh().cellZones()[cellZoneI].name() + "_" + fz.name(),
328  mesh().time().timeName(),
329  mesh(),
332  false
333  ),
334  pointDisplacement_ // to inherit the boundary conditions
335  )
336  );
337 
338  patchi++;
339  }
340 
341 
342 
343  // 'correctBoundaryConditions'
344  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
345  // Loops over all the faceZones and walks their boundary values
346 
347  // Make sure we can pick up bc values from field
348  pointDisplacement_.correctBoundaryConditions();
349 
350  patchi = 0;
351  forAllConstIter(dictionary, patchesDict, patchiter)
352  {
353  const word& faceZoneName = patchiter().keyword();
354  const dictionary& faceZoneDict = patchiter().dict();
355 
356  // Determine the points of the faceZone within the cellZone
357  const faceZone& fz = mesh().faceZones()[faceZoneName];
358  const labelList& fzMeshPoints = fz().meshPoints();
359  DynamicList<label> meshPoints(fzMeshPoints.size());
360  forAll(fzMeshPoints, i)
361  {
362  if (isZonePoint[fzMeshPoints[i]])
363  {
364  meshPoints.append(fzMeshPoints[i]);
365  }
366  }
367 
368  // Get initial value for all the faceZone points
369  tmp<vectorField> tseed = faceZoneEvaluate
370  (
371  fz,
372  meshPoints,
373  faceZoneDict,
374  patchDisp,
375  patchi
376  );
377 
378  if (debug)
379  {
380  Info<< "For cellZone:" << cellZoneI
381  << " for faceZone:" << fz.name()
382  << " nPoints:" << tseed().size()
383  << " have patchField:"
384  << " max:" << gMax(tseed())
385  << " min:" << gMin(tseed())
386  << " avg:" << gAverage(tseed())
387  << endl;
388  }
389 
390  // Set distance and transported value
391  walkStructured
392  (
393  cellZoneI,
394  isZonePoint,
395  isZoneEdge,
396 
397  meshPoints,
398  tseed,
399  patchDist[patchi],
400  patchDisp[patchi]
401  );
402 
403  // Implement real bc.
404  patchDisp[patchi].correctBoundaryConditions();
405 
406  patchi++;
407  }
408 
409 
410  // Solve
411  // ~~~~~
412 
413  if (debug)
414  {
415  // Normalised distance
417  (
418  IOobject
419  (
420  mesh().cellZones()[cellZoneI].name() + ":distance",
421  mesh().time().timeName(),
422  mesh(),
425  false
426  ),
427  pointMesh::New(mesh()),
429  );
430 
431  forAll(distance, pointi)
432  {
433  scalar d1 = patchDist[0][pointi];
434  scalar d2 = patchDist[1][pointi];
435  if (d1 + d2 > small)
436  {
437  scalar s = d1/(d1 + d2);
438  distance[pointi] = s;
439  }
440  }
441 
442  Info<< "Writing " << pointScalarField::typeName
443  << distance.name() << " to "
444  << mesh().time().timeName() << endl;
445  distance.write();
446  }
447 
448 
449  const word interpolationScheme = zoneDict.lookup("interpolationScheme");
450 
451  if (interpolationScheme == "oneSided")
452  {
453  forAll(pointDisplacement_, pointi)
454  {
455  if (isZonePoint[pointi])
456  {
457  pointDisplacement_[pointi] = patchDisp[0][pointi];
458  }
459  }
460  }
461  else if (interpolationScheme == "linear")
462  {
463  forAll(pointDisplacement_, pointi)
464  {
465  if (isZonePoint[pointi])
466  {
467  scalar d1 = patchDist[0][pointi];
468  scalar d2 = patchDist[1][pointi];
469  scalar s = d1/(d1 + d2 + vSmall);
470 
471  const vector& pd1 = patchDisp[0][pointi];
472  const vector& pd2 = patchDisp[1][pointi];
473 
474  pointDisplacement_[pointi] = (1 - s)*pd1 + s*pd2;
475  }
476  }
477  }
478  else
479  {
481  << "Invalid interpolationScheme: " << interpolationScheme
482  << ". Valid schemes are 'oneSided' and 'linear'"
483  << exit(FatalError);
484  }
485 }
486 
487 
488 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
489 
492 (
493  const word& name,
494  const polyMesh& mesh,
495  const dictionary& dict
496 )
497 :
498  displacementMotionSolver(name, mesh, dict, typeName)
499 {}
500 
501 
502 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
503 
506 {}
507 
508 
509 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
510 
513 {
514  tmp<pointField> tcurPoints
515  (
516  points0() + pointDisplacement_.primitiveField()
517  );
518 
519  return tcurPoints;
520 }
521 
522 
524 {
525  // The points have moved so before interpolation update the motionSolver
526  movePoints(mesh().points());
527 
528  // Apply boundary conditions
529  pointDisplacement_.boundaryFieldRef().updateCoeffs();
530 
531  // Solve motion on all regions (=cellZones)
532  const dictionary& regionDicts = coeffDict().subDict("regions");
533  forAllConstIter(dictionary, regionDicts, regionIter)
534  {
535  const word& cellZoneName = regionIter().keyword();
536  const dictionary& regionDict = regionIter().dict();
537 
538  label zoneI = mesh().cellZones().findZoneID(cellZoneName);
539 
540  Info<< "solving for zone: " << cellZoneName << endl;
541 
542  if (zoneI == -1)
543  {
545  << "Cannot find cellZone " << cellZoneName
546  << endl << "Valid zones are " << mesh().cellZones().names()
547  << exit(FatalError);
548  }
549 
550  cellZoneSolve(zoneI, regionDict);
551  }
552 
553  // Update pointDisplacement for solved values
554  const pointConstraints& pcs =
555  pointConstraints::New(pointDisplacement_.mesh());
556  pcs.constrainDisplacement(pointDisplacement_, false);
557 }
558 
559 
561 (
562  const polyTopoChangeMap& map
563 )
564 {
566  << "Probably inconsistent with points0MotionSolver" << nl
567  << " Needs to be updated and tested."
568  << exit(FatalError);
569 
571 
572  const vectorField displacement(this->newPoints() - points0_);
573 
574  forAll(points0_, pointi)
575  {
576  const label oldPointi = map.pointMap()[pointi];
577 
578  if (oldPointi >= 0)
579  {
580  label masterPointi = map.reversePointMap()[oldPointi];
581 
582  if ((masterPointi != pointi))
583  {
584  // newly inserted point in this cellZone
585 
586  // need to set point0 so that it represents the position that
587  // it would have had if it had existed for all time
588  points0_[pointi] -= displacement[pointi];
589  }
590  }
591  }
592 }
593 
594 
595 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:453
const labelListList & cellEdges() const
label findZoneID(const word &zoneName) const
Find zone index given a name.
Definition: MeshZones.C:341
Virtual base class for displacement motion solver.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
virtual Ostream & write(const char)=0
Write character.
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
bool set(const label) const
Is element set.
Definition: PtrListI.H:65
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
const meshCellZones & cellZones() const
Return cell zones.
Definition: polyMesh.H:501
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
Type gMin(const FieldField< Field, Type > &f)
virtual void topoChange(const polyTopoChangeMap &)
Update topology.
static const char *const typeName
Definition: Field.H:105
void constrainDisplacement(pointVectorField &displacement, const bool overrideValue=false) const
Apply boundary conditions (single-patch constraints),.
const labelList & pointMap() const
Old point map.
virtual tmp< pointField > curPoints() const
Return point location obtained from the current motion field.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
label findPatchID(const word &patchName) const
Find patch index given a name.
GeometricField< vector, pointPatchField, pointMesh > pointVectorField
scalar distance(const vector &p1, const vector &p2)
Definition: curveTools.C:12
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:372
fvMesh & mesh
Macros for easy insertion into run-time selection tables.
const dimensionSet dimLength
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:1002
wordList names() const
Return a list of zone names.
Definition: MeshZones.C:256
Application of (multi-)patch point constraints.
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
const pointField & points
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronise values on all mesh edges.
static word timeName(const scalar, const int precision=curPrecision_)
Return time name of given scalar time.
Definition: Time.C:666
A class for handling words, derived from string.
Definition: word.H:59
label nPoints
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static const word null
An empty word.
Definition: word.H:77
const labelList & reversePointMap() const
Reverse point map.
word timeName
Definition: getTimeIndex.H:3
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
virtual void topoChange(const polyTopoChangeMap &)
Update local data for topology changes.
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
static const char nl
Definition: Ostream.H:260
Type gMax(const FieldField< Field, Type > &f)
defineTypeNameAndDebug(combustionModel, 0)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
label patchi
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Type gAverage(const FieldField< Field, Type > &f)
const meshFaceZones & faceZones() const
Return face zones.
Definition: polyMesh.H:495
displacementLayeredMotionMotionSolver(const word &name, const polyMesh &, const dictionary &)
Construct from polyMesh and dictionary.
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronise values on all mesh points.
messageStream Info
Field< vector > vectorField
Specialisation of Field<T> for vector.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
A class for managing temporary objects.
Definition: PtrList.H:53
const labelListList & cellPoints() const
Namespace for OpenFOAM.