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-2016 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.ref();
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  {
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 
341  motionMask_ =
342  vertexMarkup
343  (
344  topoChangeMap().preMotionPoints(),
345  curLeft_,
346  curRight_
347  );
348 
349  // Move points inside the motionMask
350  newPoints =
351  topoChangeMap().preMotionPoints()
352  + (
353  pos(0.5 - mag(motionMask_)) // cells above the body
354  )*curMotionVel_*time().deltaT().value();
355  }
356  else
357  {
358  Info<< "Topology change. Already set mesh points" << endl;
359 
360  motionMask_ =
361  vertexMarkup
362  (
363  points(),
364  curLeft_,
365  curRight_
366  );
367 
368  // Move points inside the motionMask
369  newPoints =
370  points()
371  + (
372  pos(0.5 - mag(motionMask_)) // cells above the body
373  )*curMotionVel_*time().deltaT().value();
374  }
375  }
376  else
377  {
378  Info<< "No topology change" << endl;
379  // Set the mesh motion
380  newPoints =
381  points()
382  + (
383  pos(0.5 - mag(motionMask_)) // cells above the body
384  )*curMotionVel_*time().deltaT().value();
385  }
386 
387  // The mesh now contains the cells with zero volume
388  Info << "Executing mesh motion" << endl;
389  movePoints(newPoints);
390  // The mesh now has got non-zero volume cells
391 
392  curLeft_ = average
393  (
394  faceZones()
395  [
396  faceZones().findZoneID("leftExtrusionFaces")
397  ]().localPoints()
398  ).x() - SMALL;
399 
400  curRight_ = average
401  (
402  faceZones()
403  [
404  faceZones().findZoneID("rightExtrusionFaces")
405  ]().localPoints()
406  ).x() + SMALL;
407 
408 
409  return true;
410 }
411 
412 
413 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
bool movePoints()
Do what is neccessary if the mesh has moved.
const Type & value() const
Return const reference to value.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
virtual bool update()
Update the mesh for both mesh motion and topology change.
autoPtr< mapPolyMesh > changeMesh(const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Macros for easy insertion into run-time selection tables.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:979
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:53
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
dimensionedScalar pos(const dimensionedScalar &ds)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
stressControl lookup("compactNormalStress") >> compactNormalStress
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Abstract base class for a topology changing fvMesh.
List< label > labelList
A List of labels.
Definition: labelList.H:56
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
bool valid() const
Return true if the autoPtr valid (ie, the pointer is set).
Definition: autoPtrI.H:83
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
dimensionedScalar sin(const dimensionedScalar &ds)
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
dimensionedScalar deltaT() const
Return time step.
Definition: TimeStateI.H:53
Constant dispersed-phase particle diameter model.
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:463
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:53
A class for managing temporary objects.
Definition: PtrList.H:54
polyTopoChanger topoChanger_
virtual ~movingConeTopoFvMesh()
Destructor.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
runTime write()
Namespace for OpenFOAM.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
#define InfoInFunction
Report an information message using Foam::Info.