movingConeTopoFvMesh.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) 2011 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 
26 #include "movingConeTopoFvMesh.H"
27 #include "Time.H"
28 #include "mapPolyMesh.H"
29 #include "layerAdditionRemoval.H"
31 #include "meshTools.H"
32 #include "OFstream.H"
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38  defineTypeNameAndDebug(movingConeTopoFvMesh, 0);
39 
41  (
42  topoChangerFvMesh,
43  movingConeTopoFvMesh,
44  IOobject
45  );
46 }
47 
48 
49 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
50 
51 Foam::tmp<Foam::scalarField> Foam::movingConeTopoFvMesh::vertexMarkup
52 (
53  const pointField& p,
54  const scalar curLeft,
55  const scalar curRight
56 ) const
57 {
58  Info<< "Updating vertex markup. curLeft: "
59  << curLeft << " curRight: " << curRight << endl;
60 
61  tmp<scalarField> tvertexMarkup(new scalarField(p.size()));
62  scalarField& vertexMarkup = tvertexMarkup();
63 
64  forAll(p, pI)
65  {
66  if (p[pI].x() < curLeft - SMALL)
67  {
68  vertexMarkup[pI] = -1;
69  }
70  else if (p[pI].x() > curRight + SMALL)
71  {
72  vertexMarkup[pI] = 1;
73  }
74  else
75  {
76  vertexMarkup[pI] = 0;
77  }
78  }
79 
80  return tvertexMarkup;
81 }
82 
83 
84 void Foam::movingConeTopoFvMesh::addZonesAndModifiers()
85 {
86  // Add zones and modifiers for motion action
87 
88  if
89  (
90  pointZones().size()
91  || faceZones().size()
92  || cellZones().size()
93  || topoChanger_.size()
94  )
95  {
96  Info<< "void movingConeTopoFvMesh::addZonesAndModifiers() : "
97  << "Zones and modifiers already present. Skipping."
98  << endl;
99 
100  return;
101  }
102 
103  Info<< "Time = " << time().timeName() << endl
104  << "Adding zones and modifiers to the mesh" << endl;
105 
106  const vectorField& fc = faceCentres();
107  const vectorField& fa = faceAreas();
108 
109  labelList zone1(fc.size());
110  boolList flipZone1(fc.size(), false);
111  label nZoneFaces1 = 0;
112 
113  labelList zone2(fc.size());
114  boolList flipZone2(fc.size(), false);
115  label nZoneFaces2 = 0;
116 
117  forAll(fc, faceI)
118  {
119  if
120  (
121  fc[faceI].x() > -0.003501
122  && fc[faceI].x() < -0.003499
123  )
124  {
125  if ((fa[faceI] & vector(1, 0, 0)) < 0)
126  {
127  flipZone1[nZoneFaces1] = true;
128  }
129 
130  zone1[nZoneFaces1] = faceI;
131  Info<< "face " << faceI << " for zone 1. Flip: "
132  << flipZone1[nZoneFaces1] << endl;
133  nZoneFaces1++;
134  }
135  else if
136  (
137  fc[faceI].x() > -0.00701
138  && fc[faceI].x() < -0.00699
139  )
140  {
141  zone2[nZoneFaces2] = faceI;
142 
143  if ((fa[faceI] & vector(1, 0, 0)) > 0)
144  {
145  flipZone2[nZoneFaces2] = true;
146  }
147 
148  Info<< "face " << faceI << " for zone 2. Flip: "
149  << flipZone2[nZoneFaces2] << endl;
150  nZoneFaces2++;
151  }
152  }
153 
154  zone1.setSize(nZoneFaces1);
155  flipZone1.setSize(nZoneFaces1);
156 
157  zone2.setSize(nZoneFaces2);
158  flipZone2.setSize(nZoneFaces2);
159 
160  Info<< "zone: " << zone1 << endl;
161  Info<< "zone: " << zone2 << endl;
162 
163  List<pointZone*> pz(0);
164  List<faceZone*> fz(2);
165  List<cellZone*> cz(0);
166 
167  label nFz = 0;
168 
169  fz[nFz] =
170  new faceZone
171  (
172  "rightExtrusionFaces",
173  zone1,
174  flipZone1,
175  nFz,
176  faceZones()
177  );
178  nFz++;
179 
180  fz[nFz] =
181  new faceZone
182  (
183  "leftExtrusionFaces",
184  zone2,
185  flipZone2,
186  nFz,
187  faceZones()
188  );
189  nFz++;
190 
191  fz.setSize(nFz);
192 
193  Info<< "Adding mesh zones." << endl;
194  addZones(pz, fz, cz);
195 
196 
197  // Add layer addition/removal interfaces
198 
199  List<polyMeshModifier*> tm(2);
200  label nMods = 0;
201 
202  tm[nMods] =
203  new layerAdditionRemoval
204  (
205  "right",
206  nMods,
207  topoChanger_,
208  "rightExtrusionFaces",
209  readScalar
210  (
211  motionDict_.subDict("right").lookup("minThickness")
212  ),
213  readScalar
214  (
215  motionDict_.subDict("right").lookup("maxThickness")
216  )
217  );
218  nMods++;
219 
220  tm[nMods] = new layerAdditionRemoval
221  (
222  "left",
223  nMods,
224  topoChanger_,
225  "leftExtrusionFaces",
226  readScalar
227  (
228  motionDict_.subDict("left").lookup("minThickness")
229  ),
230  readScalar
231  (
232  motionDict_.subDict("left").lookup("maxThickness")
233  )
234  );
235  nMods++;
236  tm.setSize(nMods);
237 
238  Info<< "Adding " << nMods << " mesh modifiers" << endl;
239  topoChanger_.addTopologyModifiers(tm);
240 
241  write();
242 }
243 
244 
245 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
246 
247 // Construct from components
248 Foam::movingConeTopoFvMesh::movingConeTopoFvMesh(const IOobject& io)
249 :
250  topoChangerFvMesh(io),
251  motionDict_
252  (
254  (
255  IOobject
256  (
257  "dynamicMeshDict",
258  time().constant(),
259  *this,
260  IOobject::MUST_READ_IF_MODIFIED,
261  IOobject::NO_WRITE,
262  false
263  )
264  ).subDict(typeName + "Coeffs")
265  ),
266  motionVelAmplitude_(motionDict_.lookup("motionVelAmplitude")),
267  motionVelPeriod_(readScalar(motionDict_.lookup("motionVelPeriod"))),
268  curMotionVel_
269  (
270  motionVelAmplitude_*
271  Foam::sin(time().value()*M_PI/motionVelPeriod_)
272  ),
273  leftEdge_(readScalar(motionDict_.lookup("leftEdge"))),
274  curLeft_(readScalar(motionDict_.lookup("leftObstacleEdge"))),
275  curRight_(readScalar(motionDict_.lookup("rightObstacleEdge")))
276 {
277  Pout<< "Initial time:" << time().value()
278  << " Initial curMotionVel_:" << curMotionVel_
279  << endl;
280 
281  addZonesAndModifiers();
282 
283  curLeft_ = average
284  (
285  faceZones()
286  [
287  faceZones().findZoneID("leftExtrusionFaces")
288  ]().localPoints()
289  ).x() - SMALL;
290 
291  curRight_ = average
292  (
293  faceZones()
294  [
295  faceZones().findZoneID("rightExtrusionFaces")
296  ]().localPoints()
297  ).x() + SMALL;
298 
299  motionMask_ = vertexMarkup
300  (
301  points(),
302  curLeft_,
303  curRight_
304  );
305 }
306 
307 
308 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
309 
311 {}
312 
313 
314 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
315 
317 {
318  // Do mesh changes (use inflation - put new points in topoChangeMap)
319  autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh(true);
320 
321  // Calculate the new point positions depending on whether the
322  // topological change has happened or not
323  pointField newPoints;
324 
325  vector curMotionVel_ =
326  motionVelAmplitude_*
327  Foam::sin(time().value()*M_PI/motionVelPeriod_);
328 
329  Pout<< "time:" << time().value() << " curMotionVel_:" << curMotionVel_
330  << " curLeft:" << curLeft_ << " curRight:" << curRight_
331  << endl;
332 
333  if (topoChangeMap.valid())
334  {
335  Info<< "Topology change. Calculating motion points" << endl;
336 
337  if (topoChangeMap().hasMotionPoints())
338  {
339  Info<< "Topology change. Has premotion points" << endl;
340  //Info<< "preMotionPoints:" << topoChangeMap().preMotionPoints()
341  // << endl;
342 
343  //mkDir(time().timePath());
344  //{
345  // OFstream str(time().timePath()/"meshPoints.obj");
346  // Pout<< "Writing mesh with meshPoints to " << str.name()
347  // << endl;
348  //
349  // const pointField& currentPoints = points();
350  // label vertI = 0;
351  // forAll(currentPoints, pointI)
352  // {
353  // meshTools::writeOBJ(str, currentPoints[pointI]);
354  // vertI++;
355  // }
356  // forAll(edges(), edgeI)
357  // {
358  // const edge& e = edges()[edgeI];
359  // str << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
360  // }
361  //}
362  //{
363  // OFstream str(time().timePath()/"preMotionPoints.obj");
364  // Pout<< "Writing mesh with preMotionPoints to " << str.name()
365  // << endl;
366  //
367  // const pointField& newPoints =
368  // topoChangeMap().preMotionPoints();
369  // label vertI = 0;
370  // forAll(newPoints, pointI)
371  // {
372  // meshTools::writeOBJ(str, newPoints[pointI]);
373  // vertI++;
374  // }
375  // forAll(edges(), edgeI)
376  // {
377  // const edge& e = edges()[edgeI];
378  // str << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
379  // }
380  //}
381 
382 
383  motionMask_ =
384  vertexMarkup
385  (
386  topoChangeMap().preMotionPoints(),
387  curLeft_,
388  curRight_
389  );
390 
391  // Move points inside the motionMask
392  newPoints =
393  topoChangeMap().preMotionPoints()
394  + (
395  pos(0.5 - mag(motionMask_)) // cells above the body
396  )*curMotionVel_*time().deltaT().value();
397  }
398  else
399  {
400  Info<< "Topology change. Already set mesh points" << endl;
401 
402  motionMask_ =
403  vertexMarkup
404  (
405  points(),
406  curLeft_,
407  curRight_
408  );
409 
410  // Move points inside the motionMask
411  newPoints =
412  points()
413  + (
414  pos(0.5 - mag(motionMask_)) // cells above the body
415  )*curMotionVel_*time().deltaT().value();
416  }
417  }
418  else
419  {
420  Info<< "No topology change" << endl;
421  // Set the mesh motion
422  newPoints =
423  points()
424  + (
425  pos(0.5 - mag(motionMask_)) // cells above the body
426  )*curMotionVel_*time().deltaT().value();
427  }
428 
429  // The mesh now contains the cells with zero volume
430  Info << "Executing mesh motion" << endl;
431  movePoints(newPoints);
432  // The mesh now has got non-zero volume cells
433 
434  curLeft_ = average
435  (
436  faceZones()
437  [
438  faceZones().findZoneID("leftExtrusionFaces")
439  ]().localPoints()
440  ).x() - SMALL;
441 
442  curRight_ = average
443  (
444  faceZones()
445  [
446  faceZones().findZoneID("rightExtrusionFaces")
447  ]().localPoints()
448  ).x() + SMALL;
449 
450 
451  return true;
452 }
453 
454 
455 // ************************************************************************* //
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:463
polyTopoChanger topoChanger_
dimensioned< scalar > mag(const dimensioned< Type > &)
bool valid() const
Return true if the autoPtr valid (ie, the pointer is set).
Definition: autoPtrI.H:83
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
autoPtr< mapPolyMesh > changeMesh(const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
virtual bool update()
Update the mesh for both mesh motion and topology change.
messageStream Info
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Namespace for OpenFOAM.
runTime write()
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Abstract base class for a topology changing fvMesh.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
stressControl lookup("compactNormalStress") >> compactNormalStress
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1035
Field< vector > vectorField
Specialisation of Field<T> for vector.
#define forAll(list, i)
Definition: UList.H:421
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:53
Macros for easy insertion into run-time selection tables.
Constant dispersed-phase particle diameter model.
List< label > labelList
A List of labels.
Definition: labelList.H:56
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
dimensionedScalar pos(const dimensionedScalar &ds)
dimensionedScalar deltaT() const
Return time step.
Definition: TimeState.C:79
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:117
A class for managing temporary objects.
Definition: PtrList.H:118
bool movePoints()
Do what is neccessary if the mesh has moved.
const Type & value() const
Return const reference to value.
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
virtual ~movingConeTopoFvMesh()
Destructor.
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
dimensionedScalar sin(const dimensionedScalar &ds)