layerAdditionRemoval.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-2021 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 Description
25  Cell layer addition/removal mesh modifier
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "layerAdditionRemoval.H"
30 #include "polyTopoChanger.H"
31 #include "polyMesh.H"
32 #include "Time.H"
33 #include "primitiveMesh.H"
34 #include "polyTopoChange.H"
36 
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41  defineTypeNameAndDebug(layerAdditionRemoval, 0);
43  (
44  polyMeshModifier,
45  layerAdditionRemoval,
46  dictionary
47  );
48 }
49 
50 
51 const Foam::scalar Foam::layerAdditionRemoval::addDelta_ = 0.3;
52 const Foam::scalar Foam::layerAdditionRemoval::removeDelta_ = 0.1;
53 
54 
55 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
56 
57 void Foam::layerAdditionRemoval::checkDefinition()
58 {
59  if (!faceZoneID_.active())
60  {
62  << "Master face zone named " << faceZoneID_.name()
63  << " cannot be found."
64  << abort(FatalError);
65  }
66 
67  if
68  (
69  minLayerThickness_ < vSmall
70  || maxLayerThickness_ < minLayerThickness_
71  )
72  {
74  << "Incorrect layer thickness definition."
75  << abort(FatalError);
76  }
77 
78  label nFaces = topoChanger().mesh().faceZones()[faceZoneID_.index()].size();
79 
80  reduce(nFaces, sumOp<label>());
81 
82  if (nFaces == 0)
83  {
85  << "Face extrusion zone contains no faces. "
86  << "Please check your mesh definition."
87  << abort(FatalError);
88  }
89 
90  if (debug)
91  {
92  Pout<< "Cell layer addition/removal object " << name() << " :" << nl
93  << " faceZoneID: " << faceZoneID_ << endl;
94  }
95 }
96 
97 
98 Foam::scalar Foam::layerAdditionRemoval::readOldThickness
99 (
100  const dictionary& dict
101 )
102 {
103  return dict.lookupOrDefault("oldLayerThickness", -1.0);
104 }
105 
106 
107 void Foam::layerAdditionRemoval::clearAddressing() const
108 {
109  if (pointsPairingPtr_)
110  {
111  if (debug)
112  {
113  Pout<< "layerAdditionRemoval::clearAddressing()" << nl
114  << " clearing pointsPairingPtr_" << endl;
115  }
116 
117  deleteDemandDrivenData(pointsPairingPtr_);
118  }
119 
120  if (facesPairingPtr_)
121  {
122  if (debug)
123  {
124  Pout<< "layerAdditionRemoval::clearAddressing()" << nl
125  << " clearing facesPairingPtr_" << endl;
126  }
127 
128  deleteDemandDrivenData(facesPairingPtr_);
129  }
130 }
131 
132 
133 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
134 
136 (
137  const word& name,
138  const label index,
139  const polyTopoChanger& ptc,
140  const word& zoneName,
141  const scalar minThickness,
142  const scalar maxThickness,
143  const Switch thicknessFromVolume
144 )
145 :
146  polyMeshModifier(name, index, ptc, true),
147  faceZoneID_(zoneName, ptc.mesh().faceZones()),
148  minLayerThickness_(minThickness),
149  maxLayerThickness_(maxThickness),
150  thicknessFromVolume_(thicknessFromVolume),
151  oldLayerThickness_(-1.0),
152  pointsPairingPtr_(nullptr),
153  facesPairingPtr_(nullptr),
154  triggerRemoval_(-1),
155  triggerAddition_(-1)
156 {
157  checkDefinition();
158 }
159 
160 
162 (
163  const word& name,
164  const dictionary& dict,
165  const label index,
166  const polyTopoChanger& ptc
167 )
168 :
169  polyMeshModifier(name, index, ptc, Switch(dict.lookup("active"))),
170  faceZoneID_(dict.lookup("faceZoneName"), ptc.mesh().faceZones()),
171  minLayerThickness_(dict.lookup<scalar>("minLayerThickness")),
172  maxLayerThickness_(dict.lookup<scalar>("maxLayerThickness")),
173  thicknessFromVolume_
174  (
175  dict.lookupOrDefault<Switch>("thicknessFromVolume", true)
176  ),
177  oldLayerThickness_(readOldThickness(dict)),
178  pointsPairingPtr_(nullptr),
179  facesPairingPtr_(nullptr),
180  triggerRemoval_(-1),
181  triggerAddition_(-1)
182 {
183  checkDefinition();
184 }
185 
186 
187 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
188 
190 {
191  clearAddressing();
192 }
193 
194 
195 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
196 
198 {
199  // Protect from multiple calculation in the same time-step
200  if (triggerRemoval_ > -1 || triggerAddition_ > -1)
201  {
202  return true;
203  }
204 
205  // Go through all the cells in the master layer and calculate
206  // approximate layer thickness as the ratio of the cell volume and
207  // face area in the face zone.
208  // Layer addition:
209  // When the max thickness exceeds the threshold, trigger refinement.
210  // Layer removal:
211  // When the min thickness falls below the threshold, trigger removal.
212 
213  const polyMesh& mesh = topoChanger().mesh();
214 
215  const faceZone& fz = mesh.faceZones()[faceZoneID_.index()];
216  const labelList& mc = fz.masterCells();
217 
218  const scalarField& V = mesh.cellVolumes();
219  const vectorField& S = mesh.faceAreas();
220 
221  if (min(V) < -vSmall)
222  {
224  << "negative cell volume. Error in mesh motion before "
225  << "topological change.\n V: " << V
226  << abort(FatalError);
227  }
228 
229  scalar avgDelta = 0;
230  scalar minDelta = great;
231  scalar maxDelta = 0;
232  label nDelta = 0;
233 
234  if (thicknessFromVolume_)
235  {
236  // Thickness calculated from cell volume/face area
237  forAll(fz, facei)
238  {
239  scalar curDelta = V[mc[facei]]/mag(S[fz[facei]]);
240  avgDelta += curDelta;
241  minDelta = min(minDelta, curDelta);
242  maxDelta = max(maxDelta, curDelta);
243  }
244 
245  nDelta = fz.size();
246  }
247  else
248  {
249  // Thickness calculated from edges on layer
250  const Map<label>& meshZonePointMap = fz().meshPointMap();
251 
252  // Edges with only one point on zone
253  forAll(mc, facei)
254  {
255  const cell& cFaces = mesh.cells()[mc[facei]];
256  const edgeList cellEdges(cFaces.edges(mesh.faces()));
257 
258  forAll(cellEdges, i)
259  {
260  const edge& e = cellEdges[i];
261 
262  if (meshZonePointMap.found(e[0]))
263  {
264  if (!meshZonePointMap.found(e[1]))
265  {
266  scalar curDelta = e.mag(mesh.points());
267  avgDelta += curDelta;
268  nDelta++;
269  minDelta = min(minDelta, curDelta);
270  maxDelta = max(maxDelta, curDelta);
271  }
272  }
273  else
274  {
275  if (meshZonePointMap.found(e[1]))
276  {
277  scalar curDelta = e.mag(mesh.points());
278  avgDelta += curDelta;
279  nDelta++;
280  minDelta = min(minDelta, curDelta);
281  maxDelta = max(maxDelta, curDelta);
282  }
283  }
284  }
285  }
286  }
287 
288  reduce(minDelta, minOp<scalar>());
289  reduce(maxDelta, maxOp<scalar>());
290  reduce(avgDelta, sumOp<scalar>());
291  reduce(nDelta, sumOp<label>());
292 
293  avgDelta /= nDelta;
294 
295  if (debug)
296  {
297  Pout<< "bool layerAdditionRemoval::changeTopology() const "
298  << " for object " << name() << " : " << nl
299  << "Layer thickness: min: " << minDelta
300  << " max: " << maxDelta << " avg: " << avgDelta
301  << " old thickness: " << oldLayerThickness_ << nl
302  << "Removal threshold: " << minLayerThickness_
303  << " addition threshold: " << maxLayerThickness_ << endl;
304  }
305 
306  bool topologicalChange = false;
307 
308  // If the thickness is decreasing and crosses the min thickness,
309  // trigger removal
310  if (oldLayerThickness_ < 0)
311  {
312  if (debug)
313  {
314  Pout<< "First step. No addition/removal" << endl;
315  }
316 
317  // No topological changes allowed before first mesh motion
318  oldLayerThickness_ = avgDelta;
319 
320  topologicalChange = false;
321  }
322  else if (avgDelta < oldLayerThickness_)
323  {
324  // Layers moving towards removal
325  if (minDelta < minLayerThickness_)
326  {
327  // Check layer pairing
328  if (setLayerPairing())
329  {
330  // A mesh layer detected. Check that collapse is valid
331  if (validCollapse())
332  {
333  // At this point, info about moving the old mesh
334  // in a way to collapse the cells in the removed
335  // layer is available. Not sure what to do with
336  // it.
337 
338  if (debug)
339  {
340  Pout<< "bool layerAdditionRemoval::changeTopology() "
341  << " const for object " << name() << " : "
342  << "Triggering layer removal" << endl;
343  }
344 
345  triggerRemoval_ = mesh.time().timeIndex();
346 
347  // Old thickness looses meaning.
348  // Set it up to indicate layer removal
349  oldLayerThickness_ = great;
350 
351  topologicalChange = true;
352  }
353  else
354  {
355  // No removal, clear addressing
356  clearAddressing();
357  }
358  }
359  }
360  else
361  {
362  oldLayerThickness_ = avgDelta;
363  }
364  }
365  else
366  {
367  // Layers moving towards addition
368  if (maxDelta > maxLayerThickness_)
369  {
370  if (debug)
371  {
372  Pout<< "bool layerAdditionRemoval::changeTopology() const "
373  << " for object " << name() << " : "
374  << "Triggering layer addition" << endl;
375  }
376 
377  triggerAddition_ = mesh.time().timeIndex();
378 
379  // Old thickness looses meaning.
380  // Set it up to indicate layer removal
381  oldLayerThickness_ = 0;
382 
383  topologicalChange = true;
384  }
385  else
386  {
387  oldLayerThickness_ = avgDelta;
388  }
389  }
390 
391  return topologicalChange;
392 }
393 
394 
396 {
397  // Insert the layer addition/removal instructions
398  // into the topological change
399 
400  if (triggerRemoval_ == topoChanger().mesh().time().timeIndex())
401  {
402  removeCellLayer(ref);
403 
404  // Clear addressing. This also resets the addition/removal data
405  if (debug)
406  {
407  Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange&) "
408  << "for object " << name() << " : "
409  << "Clearing addressing after layer removal" << endl;
410  }
411 
412  triggerRemoval_ = -1;
413  clearAddressing();
414  }
415 
416  if (triggerAddition_ == topoChanger().mesh().time().timeIndex())
417  {
418  addCellLayer(ref);
419 
420  // Clear addressing. This also resets the addition/removal data
421  if (debug)
422  {
423  Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange&) "
424  << "for object " << name() << " : "
425  << "Clearing addressing after layer addition" << endl;
426  }
427 
428  triggerAddition_ = -1;
429  clearAddressing();
430  }
431 }
432 
433 
435 {
436  if (debug)
437  {
438  Pout<< "layerAdditionRemoval::updateMesh(const mapPolyMesh&) "
439  << "for object " << name() << " : "
440  << "Clearing addressing on external request";
441 
442  if (pointsPairingPtr_ || facesPairingPtr_)
443  {
444  Pout<< "Pointers set." << endl;
445  }
446  else
447  {
448  Pout<< "Pointers not set." << endl;
449  }
450  }
451 
452  // Mesh has changed topologically. Update local topological data
453  faceZoneID_.update(topoChanger().mesh().faceZones());
454 
455  clearAddressing();
456 }
457 
458 
460 {
461  if (t < vSmall || maxLayerThickness_ < t)
462  {
464  << "Incorrect layer thickness definition."
465  << abort(FatalError);
466  }
467 
468  minLayerThickness_ = t;
469 }
470 
471 
473 {
474  if (t < minLayerThickness_)
475  {
477  << "Incorrect layer thickness definition."
478  << abort(FatalError);
479  }
480 
481  maxLayerThickness_ = t;
482 }
483 
484 
486 {
487  os << nl << type() << nl
488  << name()<< nl
489  << faceZoneID_ << nl
490  << minLayerThickness_ << nl
491  << oldLayerThickness_ << nl
492  << maxLayerThickness_ << nl
493  << thicknessFromVolume_ << endl;
494 }
495 
496 
498 {
499  os << nl << name() << nl << token::BEGIN_BLOCK << nl
500  << " type " << type()
502  << " faceZoneName " << faceZoneID_.name()
504  << " minLayerThickness " << minLayerThickness_
506  << " maxLayerThickness " << maxLayerThickness_
508  << " thicknessFromVolume " << thicknessFromVolume_
510  << " oldLayerThickness " << oldLayerThickness_
512  << " active " << active()
514  << token::END_BLOCK << endl;
515 }
516 
517 
518 // ************************************************************************* //
#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
virtual void write(Ostream &) const
Write.
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
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
virtual bool changeTopology() const
Check for topology change.
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none/any.
Definition: Switch.H:60
const labelList & masterCells() const
Return labels of master cells (cells next to the master face.
Definition: faceZone.C:323
virtual void setRefinement(polyTopoChange &) const
Insert the layer addition/removal instructions.
virtual void updateMesh(const mapPolyMesh &)
Force recalculation of locally stored data on topological change.
layerAdditionRemoval(const word &name, const label index, const polyTopoChanger &ptc, const word &zoneName, const scalar minThickness, const scalar maxThickness, const Switch thicknessFromVolume=true)
Construct from components.
const cellList & cells() const
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:158
Macros for easy insertion into run-time selection tables.
void setMaxLayerThickness(const scalar t) const
Set max layer thickness which triggers removal.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1131
edgeList edges(const faceUList &) const
Return cell edges.
Definition: cell.C:120
virtual ~layerAdditionRemoval()
Destructor.
void setMinLayerThickness(const scalar t) const
Set min layer thickness which triggers removal.
dynamicFvMesh & mesh
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:113
List of mesh modifiers defining the mesh dynamics.
A class for handling words, derived from string.
Definition: word.H:59
tmp< fvMatrix< Type > > S(const Pair< tmp< volScalarField::Internal >> &, const GeometricField< Type, fvPatchField, volMesh > &)
scalar mag(const pointField &) const
Return scalar magnitude.
Definition: edgeI.H:181
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1156
Virtual base class for mesh modifiers.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:54
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
static const char nl
Definition: Ostream.H:260
const Time & time() const
Return time.
defineTypeNameAndDebug(combustionModel, 0)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
label timeIndex() const
Return current time index.
Definition: TimeStateI.H:35
virtual void writeDict(Ostream &) const
Write dictionary.
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
const meshFaceZones & faceZones() const
Return face zones.
Definition: polyMesh.H:476
const vectorField & faceAreas() const
Direct mesh changes based on v1.3 polyTopoChange syntax.
label timeIndex
Definition: getTimeIndex.H:4
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
dimensioned< scalar > mag(const dimensioned< Type > &)
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:64
rDeltaT ref()
const polyMesh & mesh() const
Return the mesh reference.
void deleteDemandDrivenData(DataPtr &dataPtr)
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:844
const scalarField & cellVolumes() const