movingConeTopoFvMesh.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-2019 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 #include "mathematicalConstants.H"
34 
35 using namespace Foam::constant::mathematical;
36 
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41  defineTypeNameAndDebug(movingConeTopoFvMesh, 0);
42 
44  (
45  topoChangerFvMesh,
46  movingConeTopoFvMesh,
47  IOobject
48  );
49 }
50 
51 
52 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
53 
54 Foam::tmp<Foam::scalarField> Foam::movingConeTopoFvMesh::vertexMarkup
55 (
56  const pointField& p,
57  const scalar curLeft,
58  const scalar curRight
59 ) const
60 {
61  Info<< "Updating vertex markup. curLeft: "
62  << curLeft << " curRight: " << curRight << endl;
63 
64  tmp<scalarField> tvertexMarkup(new scalarField(p.size()));
65  scalarField& vertexMarkup = tvertexMarkup.ref();
66 
67  forAll(p, pI)
68  {
69  if (p[pI].x() < curLeft - small)
70  {
71  vertexMarkup[pI] = -1;
72  }
73  else if (p[pI].x() > curRight + small)
74  {
75  vertexMarkup[pI] = 1;
76  }
77  else
78  {
79  vertexMarkup[pI] = 0;
80  }
81  }
82 
83  return tvertexMarkup;
84 }
85 
86 
87 void Foam::movingConeTopoFvMesh::addZonesAndModifiers()
88 {
89  // Add zones and modifiers for motion action
90 
91  if
92  (
93  pointZones().size()
94  || faceZones().size()
95  || cellZones().size()
96  || topoChanger_.size()
97  )
98  {
100  << "Zones and modifiers already present. Skipping."
101  << endl;
102 
103  return;
104  }
105 
106  Info<< "Time = " << time().timeName() << endl
107  << "Adding zones and modifiers to the mesh" << endl;
108 
109  const vectorField& fc = faceCentres();
110  const vectorField& fa = faceAreas();
111 
112  labelList zone1(fc.size());
113  boolList flipZone1(fc.size(), false);
114  label nZoneFaces1 = 0;
115 
116  labelList zone2(fc.size());
117  boolList flipZone2(fc.size(), false);
118  label nZoneFaces2 = 0;
119 
120  forAll(fc, facei)
121  {
122  if
123  (
124  fc[facei].x() > -0.003501
125  && fc[facei].x() < -0.003499
126  )
127  {
128  if ((fa[facei] & vector(1, 0, 0)) < 0)
129  {
130  flipZone1[nZoneFaces1] = true;
131  }
132 
133  zone1[nZoneFaces1] = facei;
134  Info<< "face " << facei << " for zone 1. Flip: "
135  << flipZone1[nZoneFaces1] << endl;
136  nZoneFaces1++;
137  }
138  else if
139  (
140  fc[facei].x() > -0.00701
141  && fc[facei].x() < -0.00699
142  )
143  {
144  zone2[nZoneFaces2] = facei;
145 
146  if ((fa[facei] & vector(1, 0, 0)) > 0)
147  {
148  flipZone2[nZoneFaces2] = true;
149  }
150 
151  Info<< "face " << facei << " for zone 2. Flip: "
152  << flipZone2[nZoneFaces2] << endl;
153  nZoneFaces2++;
154  }
155  }
156 
157  zone1.setSize(nZoneFaces1);
158  flipZone1.setSize(nZoneFaces1);
159 
160  zone2.setSize(nZoneFaces2);
161  flipZone2.setSize(nZoneFaces2);
162 
163  Info<< "zone: " << zone1 << endl;
164  Info<< "zone: " << zone2 << endl;
165 
166  List<pointZone*> pz(0);
167  List<faceZone*> fz(2);
168  List<cellZone*> cz(0);
169 
170  label nFz = 0;
171 
172  fz[nFz] =
173  new faceZone
174  (
175  "rightExtrusionFaces",
176  zone1,
177  flipZone1,
178  nFz,
179  faceZones()
180  );
181  nFz++;
182 
183  fz[nFz] =
184  new faceZone
185  (
186  "leftExtrusionFaces",
187  zone2,
188  flipZone2,
189  nFz,
190  faceZones()
191  );
192  nFz++;
193 
194  fz.setSize(nFz);
195 
196  Info<< "Adding mesh zones." << endl;
197  addZones(pz, fz, cz);
198 
199 
200  // Add layer addition/removal interfaces
201 
202  List<polyMeshModifier*> tm(2);
203  label nMods = 0;
204 
205  tm[nMods] =
206  new layerAdditionRemoval
207  (
208  "right",
209  nMods,
210  topoChanger_,
211  "rightExtrusionFaces",
212  readScalar
213  (
214  motionDict_.subDict("right").lookup("minThickness")
215  ),
216  readScalar
217  (
218  motionDict_.subDict("right").lookup("maxThickness")
219  )
220  );
221  nMods++;
222 
223  tm[nMods] = new layerAdditionRemoval
224  (
225  "left",
226  nMods,
227  topoChanger_,
228  "leftExtrusionFaces",
229  readScalar
230  (
231  motionDict_.subDict("left").lookup("minThickness")
232  ),
233  readScalar
234  (
235  motionDict_.subDict("left").lookup("maxThickness")
236  )
237  );
238  nMods++;
239  tm.setSize(nMods);
240 
241  Info<< "Adding " << nMods << " mesh modifiers" << endl;
242  topoChanger_.addTopologyModifiers(tm);
243 
244  write();
245 }
246 
247 
248 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
249 
251 :
252  topoChangerFvMesh(io),
253  motionDict_(dynamicMeshDict().optionalSubDict(typeName + "Coeffs")),
254  motionVelAmplitude_(motionDict_.lookup("motionVelAmplitude")),
255  motionVelPeriod_(readScalar(motionDict_.lookup("motionVelPeriod"))),
256  curMotionVel_
257  (
258  motionVelAmplitude_*sin(time().value()*pi/motionVelPeriod_)
259  ),
260  leftEdge_(readScalar(motionDict_.lookup("leftEdge"))),
261  curLeft_(readScalar(motionDict_.lookup("leftObstacleEdge"))),
262  curRight_(readScalar(motionDict_.lookup("rightObstacleEdge")))
263 {
264  Pout<< "Initial time:" << time().value()
265  << " Initial curMotionVel_:" << curMotionVel_
266  << endl;
267 
268  addZonesAndModifiers();
269 
270  curLeft_ = average
271  (
272  faceZones()
273  [
274  faceZones().findZoneID("leftExtrusionFaces")
275  ]().localPoints()
276  ).x() - small;
277 
278  curRight_ = average
279  (
280  faceZones()
281  [
282  faceZones().findZoneID("rightExtrusionFaces")
283  ]().localPoints()
284  ).x() + small;
285 
286  motionMask_ = vertexMarkup
287  (
288  points(),
289  curLeft_,
290  curRight_
291  );
292 }
293 
294 
295 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
296 
298 {}
299 
300 
301 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
302 
304 {
305  // Do mesh changes (use inflation - put new points in topoChangeMap)
306  autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh(true);
307 
308  // Calculate the new point positions depending on whether the
309  // topological change has happened or not
310  pointField newPoints;
311 
312  vector curMotionVel_ =
313  motionVelAmplitude_*sin(time().value()*pi/motionVelPeriod_);
314 
315  Pout<< "time:" << time().value() << " curMotionVel_:" << curMotionVel_
316  << " curLeft:" << curLeft_ << " curRight:" << curRight_
317  << endl;
318 
319  if (topoChangeMap.valid())
320  {
321  Info<< "Topology change. Calculating motion points" << endl;
322 
323  if (topoChangeMap().hasMotionPoints())
324  {
325  Info<< "Topology change. Has premotion points" << endl;
326 
327  motionMask_ =
328  vertexMarkup
329  (
330  topoChangeMap().preMotionPoints(),
331  curLeft_,
332  curRight_
333  );
334 
335  // Move points inside the motionMask
336  newPoints =
337  topoChangeMap().preMotionPoints()
338  + (
339  pos0(0.5 - mag(motionMask_)) // cells above the body
340  )*curMotionVel_*time().deltaT().value();
341  }
342  else
343  {
344  Info<< "Topology change. Already set mesh points" << endl;
345 
346  motionMask_ =
347  vertexMarkup
348  (
349  points(),
350  curLeft_,
351  curRight_
352  );
353 
354  // Move points inside the motionMask
355  newPoints =
356  points()
357  + (
358  pos0(0.5 - mag(motionMask_)) // cells above the body
359  )*curMotionVel_*time().deltaT().value();
360  }
361  }
362  else
363  {
364  Info<< "No topology change" << endl;
365  // Set the mesh motion
366  newPoints =
367  points()
368  + (
369  pos0(0.5 - mag(motionMask_)) // cells above the body
370  )*curMotionVel_*time().deltaT().value();
371  }
372 
373  // The mesh now contains the cells with zero volume
374  Info << "Executing mesh motion" << endl;
375  movePoints(newPoints);
376 
377  // The mesh now has got non-zero volume cells
378 
379  curLeft_ = average
380  (
381  faceZones()
382  [
383  faceZones().findZoneID("leftExtrusionFaces")
384  ]().localPoints()
385  ).x() - small;
386 
387  curRight_ = average
388  (
389  faceZones()
390  [
391  faceZones().findZoneID("rightExtrusionFaces")
392  ]().localPoints()
393  ).x() + small;
394 
395  return true;
396 }
397 
398 
399 // ************************************************************************* //
#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
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:476
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
bool movePoints()
Do what is necessary if the mesh has moved.
movingConeTopoFvMesh(const IOobject &io)
Construct from database.
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.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:239
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:1131
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
void write(Ostream &, const label, const dictionary &)
Write with dictionary lookup.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
stressControl lookup("compactNormalStress") >> compactNormalStress
mathematical constants.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Abstract base class for a topology changing fvMesh.
const Type & value() const
Return const reference to value.
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 successful.
Definition: doubleScalar.H:68
bool valid() const
Return true if the autoPtr valid (ie, the pointer is set)
Definition: autoPtrI.H:83
dimensionedScalar pos0(const dimensionedScalar &ds)
dimensionedScalar sin(const dimensionedScalar &ds)
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
A class for managing temporary objects.
Definition: PtrList.H:53
polyTopoChanger topoChanger_
virtual ~movingConeTopoFvMesh()
Destructor.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
dimensionedScalar deltaT() const
Return time step.
Definition: TimeStateI.H:53
Namespace for OpenFOAM.
#define InfoInFunction
Report an information message using Foam::Info.