FreeFOAM The Cross-Platform CFD Toolkit
layerAdditionRemoval.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
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"
31 #include <OpenFOAM/polyMesh.H>
32 #include <OpenFOAM/Time.H>
33 #include <OpenFOAM/primitiveMesh.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  (
63  "void Foam::layerAdditionRemoval::checkDefinition()"
64  ) << "Master face zone named " << faceZoneID_.name()
65  << " cannot be found."
66  << abort(FatalError);
67  }
68 
69  if
70  (
71  minLayerThickness_ < VSMALL
72  || maxLayerThickness_ < minLayerThickness_
73  )
74  {
76  (
77  "void Foam::layerAdditionRemoval::checkDefinition()"
78  ) << "Incorrect layer thickness definition."
79  << abort(FatalError);
80  }
81 
82  if (topoChanger().mesh().faceZones()[faceZoneID_.index()].empty())
83  {
85  (
86  "void Foam::layerAdditionRemoval::checkDefinition()"
87  ) << "Face extrusion zone contains no faces. "
88  << " Please check your mesh definition."
89  << abort(FatalError);
90  }
91 
92  if (debug)
93  {
94  Pout<< "Cell layer addition/removal object " << name() << " :" << nl
95  << " faceZoneID: " << faceZoneID_ << endl;
96  }
97 }
98 
99 Foam::scalar Foam::layerAdditionRemoval::readOldThickness
100 (
101  const dictionary& dict
102 )
103 {
104  return dict.lookupOrDefault("oldLayerThickness", -1.0);
105 }
106 
107 
108 void Foam::layerAdditionRemoval::clearAddressing() const
109 {
110  // Layer removal data
111  deleteDemandDrivenData(pointsPairingPtr_);
112  deleteDemandDrivenData(facesPairingPtr_);
113 }
114 
115 
116 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
117 
118 // Construct from components
119 Foam::layerAdditionRemoval::layerAdditionRemoval
120 (
121  const word& name,
122  const label index,
123  const polyTopoChanger& mme,
124  const word& zoneName,
125  const scalar minThickness,
126  const scalar maxThickness
127 )
128 :
129  polyMeshModifier(name, index, mme, true),
130  faceZoneID_(zoneName, mme.mesh().faceZones()),
131  minLayerThickness_(minThickness),
132  maxLayerThickness_(maxThickness),
133  oldLayerThickness_(-1.0),
134  pointsPairingPtr_(NULL),
135  facesPairingPtr_(NULL),
136  triggerRemoval_(-1),
137  triggerAddition_(-1)
138 {
139  checkDefinition();
140 }
141 
142 
143 // Construct from dictionary
144 Foam::layerAdditionRemoval::layerAdditionRemoval
145 (
146  const word& name,
147  const dictionary& dict,
148  const label index,
149  const polyTopoChanger& mme
150 )
151 :
152  polyMeshModifier(name, index, mme, Switch(dict.lookup("active"))),
153  faceZoneID_(dict.lookup("faceZoneName"), mme.mesh().faceZones()),
154  minLayerThickness_(readScalar(dict.lookup("minLayerThickness"))),
155  maxLayerThickness_(readScalar(dict.lookup("maxLayerThickness"))),
156  oldLayerThickness_(readOldThickness(dict)),
157  pointsPairingPtr_(NULL),
158  facesPairingPtr_(NULL),
159  triggerRemoval_(-1),
160  triggerAddition_(-1)
161 {
162  checkDefinition();
163 }
164 
165 
166 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
167 
169 {
170  clearAddressing();
171 }
172 
173 
174 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
175 
177 {
178  // Protect from multiple calculation in the same time-step
179  if (triggerRemoval_ > -1 || triggerAddition_ > -1)
180  {
181  return true;
182  }
183 
184  // Go through all the cells in the master layer and calculate
185  // approximate layer thickness as the ratio of the cell volume and
186  // face area in the face zone.
187  // Layer addition:
188  // When the max thickness exceeds the threshold, trigger refinement.
189  // Layer removal:
190  // When the min thickness falls below the threshold, trigger removal.
191 
192  const faceZone& fz = topoChanger().mesh().faceZones()[faceZoneID_.index()];
193  const labelList& mc = fz.masterCells();
194 
195  const scalarField& V = topoChanger().mesh().cellVolumes();
196  const vectorField& S = topoChanger().mesh().faceAreas();
197 
198  if (min(V) < -VSMALL)
199  {
200  FatalErrorIn("bool layerAdditionRemoval::changeTopology() const")
201  << "negative cell volume. Error in mesh motion before "
202  << "topological change.\n V: " << V
203  << abort(FatalError);
204  }
205 
206  scalar avgDelta = 0;
207  scalar minDelta = GREAT;
208  scalar maxDelta = 0;
209 
210  forAll (fz, faceI)
211  {
212  scalar curDelta = V[mc[faceI]]/mag(S[fz[faceI]]);
213  avgDelta += curDelta;
214  minDelta = min(minDelta, curDelta);
215  maxDelta = max(maxDelta, curDelta);
216  }
217 
218  avgDelta /= fz.size();
219 
221  //{
222  // // Edges on layer.
223  // const Map<label>& zoneMeshPointMap = fz().meshPointMap();
224  //
225  // label nDelta = 0;
226  //
227  // // Edges with only one point on zone
228  // const polyMesh& mesh = topoChanger().mesh();
229  //
230  // forAll(mc, faceI)
231  // {
232  // const cell& cFaces = mesh.cells()[mc[faceI]];
233  // const edgeList cellEdges(cFaces.edges(mesh.faces()));
234  //
235  // forAll(cellEdges, i)
236  // {
237  // const edge& e = cellEdges[i];
238  //
239  // if (zoneMeshPointMap.found(e[0]))
240  // {
241  // if (!zoneMeshPointMap.found(e[1]))
242  // {
243  // scalar curDelta = e.mag(mesh.points());
244  // avgDelta += curDelta;
245  // nDelta++;
246  // minDelta = min(minDelta, curDelta);
247  // maxDelta = max(maxDelta, curDelta);
248  // }
249  // }
250  // else
251  // {
252  // if (zoneMeshPointMap.found(e[1]))
253  // {
254  // scalar curDelta = e.mag(mesh.points());
255  // avgDelta += curDelta;
256  // nDelta++;
257  // minDelta = min(minDelta, curDelta);
258  // maxDelta = max(maxDelta, curDelta);
259  // }
260  // }
261  // }
262  // }
263  //
264  // avgDelta /= nDelta;
265  //}
266 
267  if (debug)
268  {
269  Pout<< "bool layerAdditionRemoval::changeTopology() const "
270  << " for object " << name() << " : " << nl
271  << "Layer thickness: min: " << minDelta
272  << " max: " << maxDelta << " avg: " << avgDelta
273  << " old thickness: " << oldLayerThickness_ << nl
274  << "Removal threshold: " << minLayerThickness_
275  << " addition threshold: " << maxLayerThickness_ << endl;
276  }
277 
278  bool topologicalChange = false;
279 
280  // If the thickness is decreasing and crosses the min thickness,
281  // trigger removal
282  if (oldLayerThickness_ < 0)
283  {
284  if (debug)
285  {
286  Pout << "First step. No addition/removal" << endl;
287  }
288 
289  // No topological changes allowed before first mesh motion
290  //
291  oldLayerThickness_ = avgDelta;
292 
293  topologicalChange = false;
294  }
295  else if (avgDelta < oldLayerThickness_)
296  {
297  // Layers moving towards removal
298  if (minDelta < minLayerThickness_)
299  {
300  // Check layer pairing
301  if (setLayerPairing())
302  {
303  // A mesh layer detected. Check that collapse is valid
304  if (validCollapse())
305  {
306  // At this point, info about moving the old mesh
307  // in a way to collapse the cells in the removed
308  // layer is available. Not sure what to do with
309  // it.
310 
311  if (debug)
312  {
313  Pout<< "bool layerAdditionRemoval::changeTopology() "
314  << " const for object " << name() << " : "
315  << "Triggering layer removal" << endl;
316  }
317 
318  triggerRemoval_ = topoChanger().mesh().time().timeIndex();
319 
320  // Old thickness looses meaning.
321  // Set it up to indicate layer removal
322  oldLayerThickness_ = GREAT;
323 
324  topologicalChange = true;
325  }
326  else
327  {
328  // No removal, clear addressing
329  clearAddressing();
330  }
331  }
332  }
333  else
334  {
335  oldLayerThickness_ = avgDelta;
336  }
337  }
338  else
339  {
340  // Layers moving towards addition
341  if (maxDelta > maxLayerThickness_)
342  {
343  if (debug)
344  {
345  Pout<< "bool layerAdditionRemoval::changeTopology() const "
346  << " for object " << name() << " : "
347  << "Triggering layer addition" << endl;
348  }
349 
350  triggerAddition_ = topoChanger().mesh().time().timeIndex();
351 
352  // Old thickness looses meaning.
353  // Set it up to indicate layer removal
354  oldLayerThickness_ = 0;
355 
356  topologicalChange = true;
357  }
358  else
359  {
360  oldLayerThickness_ = avgDelta;
361  }
362  }
363 
364  return topologicalChange;
365 }
366 
367 
369 {
370  // Insert the layer addition/removal instructions
371  // into the topological change
372 
373  if (triggerRemoval_ == topoChanger().mesh().time().timeIndex())
374  {
375  removeCellLayer(ref);
376 
377  // Clear addressing. This also resets the addition/removal data
378  if (debug)
379  {
380  Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange& ref) "
381  << " for object " << name() << " : "
382  << "Clearing addressing after layer removal. " << endl;
383  }
384 
385  triggerRemoval_ = -1;
386  clearAddressing();
387  }
388 
389  if (triggerAddition_ == topoChanger().mesh().time().timeIndex())
390  {
391  addCellLayer(ref);
392 
393  // Clear addressing. This also resets the addition/removal data
394  if (debug)
395  {
396  Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange& ref) "
397  << " for object " << name() << " : "
398  << "Clearing addressing after layer addition. " << endl;
399  }
400 
401  triggerAddition_ = -1;
402  clearAddressing();
403  }
404 }
405 
406 
408 {
409  if (debug)
410  {
411  Pout<< "layerAdditionRemoval::updateMesh(const mapPolyMesh&) "
412  << " for object " << name() << " : "
413  << "Clearing addressing on external request. ";
414 
415  if (pointsPairingPtr_ || facesPairingPtr_)
416  {
417  Pout << "Pointers set." << endl;
418  }
419  else
420  {
421  Pout << "Pointers not set." << endl;
422  }
423  }
424 
425  // Mesh has changed topologically. Update local topological data
426  faceZoneID_.update(topoChanger().mesh().faceZones());
427 
428  clearAddressing();
429 }
430 
431 
433 {
434  if
435  (
436  t < VSMALL
437  || maxLayerThickness_ < t
438  )
439  {
441  (
442  "void layerAdditionRemoval::setMinLayerThickness("
443  "const scalar t) const"
444  ) << "Incorrect layer thickness definition."
445  << abort(FatalError);
446  }
447 
448  minLayerThickness_ = t;
449 }
450 
451 
453 {
454  if
455  (
456  t < minLayerThickness_
457  )
458  {
460  (
461  "void layerAdditionRemoval::setMaxLayerThickness("
462  "const scalar t) const"
463  ) << "Incorrect layer thickness definition."
464  << abort(FatalError);
465  }
466 
467  maxLayerThickness_ = t;
468 }
469 
470 
472 {
473  os << nl << type() << nl
474  << name()<< nl
475  << faceZoneID_ << nl
476  << minLayerThickness_ << nl
477  << oldLayerThickness_ << nl
478  << maxLayerThickness_ << endl;
479 }
480 
481 
483 {
484  os << nl << name() << nl << token::BEGIN_BLOCK << nl
485  << " type " << type()
487  << " faceZoneName " << faceZoneID_.name()
489  << " minLayerThickness " << minLayerThickness_
491  << " maxLayerThickness " << maxLayerThickness_
493  << " oldLayerThickness " << oldLayerThickness_
495  << " active " << active()
497  << token::END_BLOCK << endl;
498 }
499 
500 
501 // ************************ vim: set sw=4 sts=4 et: ************************ //