FreeFOAM The Cross-Platform CFD Toolkit
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) 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 \*---------------------------------------------------------------------------*/
25 
26 #include "movingConeTopoFvMesh.H"
27 #include <OpenFOAM/Time.H>
28 #include <OpenFOAM/mapPolyMesh.H>
31 #include <meshTools/meshTools.H>
32 #include <OpenFOAM/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,
261  IOobject::NO_WRITE
262  )
263  ).subDict(typeName + "Coeffs")
264  ),
265  motionVelAmplitude_(motionDict_.lookup("motionVelAmplitude")),
266  motionVelPeriod_(readScalar(motionDict_.lookup("motionVelPeriod"))),
267  curMotionVel_
268  (
269  motionVelAmplitude_*
270  Foam::sin(time().value()*M_PI/motionVelPeriod_)
271  ),
272  leftEdge_(readScalar(motionDict_.lookup("leftEdge"))),
273  curLeft_(readScalar(motionDict_.lookup("leftObstacleEdge"))),
274  curRight_(readScalar(motionDict_.lookup("rightObstacleEdge")))
275 {
276  Pout<< "Initial time:" << time().value()
277  << " Initial curMotionVel_:" << curMotionVel_
278  << endl;
279 
280  addZonesAndModifiers();
281 
282  curLeft_ = average
283  (
284  faceZones()
285  [
286  faceZones().findZoneID("leftExtrusionFaces")
287  ]().localPoints()
288  ).x() - SMALL;
289 
290  curRight_ = average
291  (
292  faceZones()
293  [
294  faceZones().findZoneID("rightExtrusionFaces")
295  ]().localPoints()
296  ).x() + SMALL;
297 
298  motionMask_ = vertexMarkup
299  (
300  points(),
301  curLeft_,
302  curRight_
303  );
304 }
305 
306 
307 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
308 
310 {}
311 
312 
313 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
314 
316 {
317  // Do mesh changes (use inflation - put new points in topoChangeMap)
318  autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh(true);
319 
320  // Calculate the new point positions depending on whether the
321  // topological change has happened or not
322  pointField newPoints;
323 
324  vector curMotionVel_ =
325  motionVelAmplitude_*
326  Foam::sin(time().value()*M_PI/motionVelPeriod_);
327 
328  Pout<< "time:" << time().value() << " curMotionVel_:" << curMotionVel_
329  << " curLeft:" << curLeft_ << " curRight:" << curRight_
330  << endl;
331 
332  if (topoChangeMap.valid())
333  {
334  Info << "Topology change. Calculating motion points" << endl;
335 
336  if (topoChangeMap().hasMotionPoints())
337  {
338  Info << "Topology change. Has premotion points" << endl;
339  //Info<< "preMotionPoints:" << topoChangeMap().preMotionPoints()
340  // << endl;
341 
342  //mkDir(time().timePath());
343  //{
344  // OFstream str(time().timePath()/"meshPoints.obj");
345  // Pout<< "Writing mesh with meshPoints to " << str.name()
346  // << endl;
347  //
348  // const pointField& currentPoints = points();
349  // label vertI = 0;
350  // forAll(currentPoints, pointI)
351  // {
352  // meshTools::writeOBJ(str, currentPoints[pointI]);
353  // vertI++;
354  // }
355  // forAll(edges(), edgeI)
356  // {
357  // const edge& e = edges()[edgeI];
358  // str << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
359  // }
360  //}
361  //{
362  // OFstream str(time().timePath()/"preMotionPoints.obj");
363  // Pout<< "Writing mesh with preMotionPoints to " << str.name()
364  // << endl;
365  //
366  // const pointField& newPoints =
367  // topoChangeMap().preMotionPoints();
368  // label vertI = 0;
369  // forAll(newPoints, pointI)
370  // {
371  // meshTools::writeOBJ(str, newPoints[pointI]);
372  // vertI++;
373  // }
374  // forAll(edges(), edgeI)
375  // {
376  // const edge& e = edges()[edgeI];
377  // str << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
378  // }
379  //}
380 
381 
382  motionMask_ =
383  vertexMarkup
384  (
385  topoChangeMap().preMotionPoints(),
386  curLeft_,
387  curRight_
388  );
389 
390  // Move points inside the motionMask
391  newPoints =
392  topoChangeMap().preMotionPoints()
393  + (
394  pos(0.5 - mag(motionMask_)) // cells above the body
395  )*curMotionVel_*time().deltaT().value();
396  }
397  else
398  {
399  Info << "Topology change. Already set mesh points" << endl;
400 
401  motionMask_ =
402  vertexMarkup
403  (
404  points(),
405  curLeft_,
406  curRight_
407  );
408 
409  // Move points inside the motionMask
410  newPoints =
411  points()
412  + (
413  pos(0.5 - mag(motionMask_)) // cells above the body
414  )*curMotionVel_*time().deltaT().value();
415  }
416  }
417  else
418  {
419  Info << "No topology change" << endl;
420  // Set the mesh motion
421  newPoints =
422  points()
423  + (
424  pos(0.5 - mag(motionMask_)) // cells above the body
425  )*curMotionVel_*time().deltaT().value();
426  }
427 
428  // The mesh now contains the cells with zero volume
429  Info << "Executing mesh motion" << endl;
430  movePoints(newPoints);
431  // The mesh now has got non-zero volume cells
432 
433  curLeft_ = average
434  (
435  faceZones()
436  [
437  faceZones().findZoneID("leftExtrusionFaces")
438  ]().localPoints()
439  ).x() - SMALL;
440 
441  curRight_ = average
442  (
443  faceZones()
444  [
445  faceZones().findZoneID("rightExtrusionFaces")
446  ]().localPoints()
447  ).x() + SMALL;
448 
449 
450  return true;
451 }
452 
453 
454 // ************************ vim: set sw=4 sts=4 et: ************************ //