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  motionDict_.subDict("right").lookup<scalar>("minThickness"),
213  motionDict_.subDict("right").lookup<scalar>("maxThickness")
214  );
215  nMods++;
216 
217  tm[nMods] = new layerAdditionRemoval
218  (
219  "left",
220  nMods,
221  topoChanger_,
222  "leftExtrusionFaces",
223  motionDict_.subDict("left").lookup<scalar>("minThickness"),
224  motionDict_.subDict("left").lookup<scalar>("maxThickness")
225  );
226  nMods++;
227  tm.setSize(nMods);
228 
229  Info<< "Adding " << nMods << " mesh modifiers" << endl;
230  topoChanger_.addTopologyModifiers(tm);
231 
232  write();
233 }
234 
235 
236 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
237 
239 :
240  topoChangerFvMesh(io),
241  motionDict_(dynamicMeshDict().optionalSubDict(typeName + "Coeffs")),
242  motionVelAmplitude_(motionDict_.lookup("motionVelAmplitude")),
243  motionVelPeriod_(motionDict_.lookup<scalar>("motionVelPeriod")),
244  curMotionVel_
245  (
246  motionVelAmplitude_*sin(time().value()*pi/motionVelPeriod_)
247  ),
248  leftEdge_(motionDict_.lookup<scalar>("leftEdge")),
249  curLeft_(motionDict_.lookup<scalar>("leftObstacleEdge")),
250  curRight_(motionDict_.lookup<scalar>("rightObstacleEdge"))
251 {
252  Pout<< "Initial time:" << time().value()
253  << " Initial curMotionVel_:" << curMotionVel_
254  << endl;
255 
256  addZonesAndModifiers();
257 
258  curLeft_ = average
259  (
260  faceZones()
261  [
262  faceZones().findZoneID("leftExtrusionFaces")
263  ]().localPoints()
264  ).x() - small;
265 
266  curRight_ = average
267  (
268  faceZones()
269  [
270  faceZones().findZoneID("rightExtrusionFaces")
271  ]().localPoints()
272  ).x() + small;
273 
274  motionMask_ = vertexMarkup
275  (
276  points(),
277  curLeft_,
278  curRight_
279  );
280 }
281 
282 
283 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
284 
286 {}
287 
288 
289 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
290 
292 {
293  // Do mesh changes (use inflation - put new points in topoChangeMap)
294  autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh(true);
295 
296  // Calculate the new point positions depending on whether the
297  // topological change has happened or not
298  pointField newPoints;
299 
300  vector curMotionVel_ =
301  motionVelAmplitude_*sin(time().value()*pi/motionVelPeriod_);
302 
303  Pout<< "time:" << time().value() << " curMotionVel_:" << curMotionVel_
304  << " curLeft:" << curLeft_ << " curRight:" << curRight_
305  << endl;
306 
307  if (topoChangeMap.valid())
308  {
309  Info<< "Topology change. Calculating motion points" << endl;
310 
311  if (topoChangeMap().hasMotionPoints())
312  {
313  Info<< "Topology change. Has premotion points" << endl;
314 
315  motionMask_ =
316  vertexMarkup
317  (
318  topoChangeMap().preMotionPoints(),
319  curLeft_,
320  curRight_
321  );
322 
323  // Move points inside the motionMask
324  newPoints =
325  topoChangeMap().preMotionPoints()
326  + (
327  pos0(0.5 - mag(motionMask_)) // cells above the body
328  )*curMotionVel_*time().deltaT().value();
329  }
330  else
331  {
332  Info<< "Topology change. Already set mesh points" << endl;
333 
334  motionMask_ =
335  vertexMarkup
336  (
337  points(),
338  curLeft_,
339  curRight_
340  );
341 
342  // Move points inside the motionMask
343  newPoints =
344  points()
345  + (
346  pos0(0.5 - mag(motionMask_)) // cells above the body
347  )*curMotionVel_*time().deltaT().value();
348  }
349  }
350  else
351  {
352  Info<< "No topology change" << endl;
353  // Set the mesh motion
354  newPoints =
355  points()
356  + (
357  pos0(0.5 - mag(motionMask_)) // cells above the body
358  )*curMotionVel_*time().deltaT().value();
359  }
360 
361  // The mesh now contains the cells with zero volume
362  Info << "Executing mesh motion" << endl;
363  movePoints(newPoints);
364 
365  // The mesh now has got non-zero volume cells
366 
367  curLeft_ = average
368  (
369  faceZones()
370  [
371  faceZones().findZoneID("leftExtrusionFaces")
372  ]().localPoints()
373  ).x() - small;
374 
375  curRight_ = average
376  (
377  faceZones()
378  [
379  faceZones().findZoneID("rightExtrusionFaces")
380  ]().localPoints()
381  ).x() + small;
382 
383  return true;
384 }
385 
386 
387 // ************************************************************************* //
#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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
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
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.
void write(std::ostream &os, const bool binary, List< floatScalar > &fField)
Write floats ascii or binary.
List< label > labelList
A List of labels.
Definition: labelList.H:56
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
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
const meshFaceZones & faceZones() const
Return face zones.
Definition: polyMesh.H:476
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.