Memosa-FVM  0.2
MovingMeshModel.h
Go to the documentation of this file.
1 // This file os part of FVM
2 // Copyright (c) 2012 FVM Authors
3 // See LICENSE file for terms.
4 
5 #ifndef _MOVINGMESHMODEL_H_
6 #define _MOVINGMESHMODEL_H_
7 
8 #include "Model.h"
9 #include "GeomFields.h"
10 #include "FlowFields.h"
11 #include "NumType.h"
12 #include "Mesh.h"
13 #include "Array.h"
14 #include "Field.h"
15 #include "Vector.h"
16 #include "CRConnectivity.h"
17 #include "StorageSite.h"
18 #include "MovingMeshBC.h"
19 
20 
21 template<class T>
22 class MovingMeshModel : public Model
23 {
24 public:
25  typedef Array<T> TArray;
28 
29  void advance()
30  {
31  const int numMeshes = _meshes.size();
32  for(int n=0;n<numMeshes;n++)
33  {
34  const Mesh& mesh = *_meshes[n];
35  const StorageSite& nodes = mesh.getNodes();
36  const int nNodes = nodes.getCount();
37  shared_ptr<Array<int> > GlobalToLocalPtr = mesh.createAndGetBNglobalToLocal();
38  Array<int>& GlobalToLocal = *GlobalToLocalPtr;
39  const StorageSite& boundaryNodes = mesh.getBoundaryNodes();
40  VectorT3Array& nodeCoordinate =
41  dynamic_cast<VectorT3Array&> (_geomFields.coordinate[nodes]);
42  VectorT3Array& nodeDisplacement =
43  dynamic_cast<VectorT3Array&> (_geomFields.nodeDisplacement[nodes]);
44  VectorT3Array& dirichletNodeDisplacement =
45  dynamic_cast<VectorT3Array&> (_geomFields.dirichletNodeDisplacement[nodes]);
46  VectorT3Array& nodeNormal =
47  dynamic_cast<VectorT3Array&> (_geomFields.boundaryNodeNormal[boundaryNodes]);
48  const Array<int>& displacementOptions =
49  dynamic_cast<Array<int>& > (_geomFields.displacementOptions[nodes]);
50 
51  const CRConnectivity& cellNodes = mesh.getCellNodes();
52  shared_ptr<CRConnectivity> nodeCellsPtr = cellNodes.getTranspose();
53  shared_ptr<CRConnectivity> nodeNodesPtr = nodeCellsPtr->multiply(cellNodes,false);
54  CRConnectivity& nodeNodes = *nodeNodesPtr;
55  nodeDisplacement.zero();
56  const T one(1.0);
57  const T underrelaxation = _options["underrelaxation"];
58  const T small(1e-10);
59 
60  for(int i=0;i<_options.nNodeDisplacementSweeps;i++)
61  {
62 
63  int nDirichlet =0;
64  T averageDirichletDisplacement(0.);
65 
66  shared_ptr<VectorT3Array> _previousNodeDisplacementPtr =
67  dynamic_pointer_cast<VectorT3Array>(nodeDisplacement.newCopy());
68 
69  for(int j=0;j<nNodes;j++)
70  {
72  T weight(0.0);
73  for(int k=0;k<nodeNodes.getCount(j);k++)
74  {
75  const int num = nodeNodes(j,k);
76  if(num!=j)
77  {
78  const VectorT3 ds = nodeCoordinate[num]-nodeCoordinate[j];
79  if(mag(ds)!=T(0.0))
80  {
81  dr += nodeDisplacement[num]/mag(ds);
82  weight += one/mag(ds);
83  }
84  else
85  {
86  dr += nodeDisplacement[num]/small;
87  weight += one/small;
88  }
89  }
90  }
91  dr = dr/weight;
92  if(displacementOptions[j] == 0)
93  {
94  nodeDisplacement[j] = T(0.0);
95  nodeCoordinate[j] += nodeDisplacement[j] - (*(_previousNodeDisplacementPtr))[j];
96  }
97  else if(displacementOptions[j] == 1)
98  {
99  averageDirichletDisplacement += mag((dirichletNodeDisplacement)[j]);
100  nDirichlet ++;
101 
102  nodeDisplacement[j] = (dirichletNodeDisplacement)[j];
103  nodeCoordinate[j] += nodeDisplacement[j] - (*(_previousNodeDisplacementPtr))[j];
104  }
105  else if(displacementOptions[j] == 2)
106  {
107  T temp = dot(dr,nodeNormal[GlobalToLocal[j]]);
108  nodeDisplacement[j] = dr - temp*nodeNormal[GlobalToLocal[j]];
109  //if(dot(nodeDisplacement[j],nodeNormal[GlobalToLocal[j]])!=zero)
110  // nodeDisplacement[j] = dr + temp*nodeNormal[GlobalToLocal[j]];
111  nodeDisplacement[j] = (*(_previousNodeDisplacementPtr))[j] +
112  underrelaxation*(nodeDisplacement[j]-(*(_previousNodeDisplacementPtr))[j]);
113  nodeCoordinate[j] += nodeDisplacement[j] - (*(_previousNodeDisplacementPtr))[j];
114  }
115  else if(displacementOptions[j] == 3)
116  {
117  nodeDisplacement[j] = dr;
118  nodeDisplacement[j] = (*(_previousNodeDisplacementPtr))[j] +
119  underrelaxation*(nodeDisplacement[j]-(*(_previousNodeDisplacementPtr))[j]);
120  nodeCoordinate[j] += nodeDisplacement[j] - (*(_previousNodeDisplacementPtr))[j];
121  }
122  }
123 
124  if (nDirichlet > 0)
125  averageDirichletDisplacement /= nDirichlet;
126  else
127  averageDirichletDisplacement = T(1.);
128 
129  T maxChangeInDisplacement = T(0.0);
130  for(int j=0;j<nNodes;j++)
131  {
132  const T changeInDisplacement = (mag(nodeDisplacement[j]-
133  (*(_previousNodeDisplacementPtr))[j]));
134  if (maxChangeInDisplacement < changeInDisplacement)
135  maxChangeInDisplacement = changeInDisplacement;
136  }
137 
138  T maxChangeRelative = maxChangeInDisplacement / averageDirichletDisplacement;
139  //cout<<"\nsweep "<<i<<" max change is "<< maxChangeInDisplacement
140  // <<" and ratio is "<< maxChangeRelative<<"\n";
141  if((maxChangeInDisplacement<=_options.absTolerance)||
142  (maxChangeRelative<=_options.relativeTolerance))
143  {
144 
145  return;
146  }
147  }
148  }
149  }
150 
151  MovingMeshModel(const MeshList& meshes, GeomFields& geomFields,
152  FlowFields& flowFields) :
153  Model(meshes),
154  _geomFields(geomFields),
155  _flowFields(flowFields),
156  _meshes(meshes)
157  {
158  logCtor();
159  }
160 
161  virtual ~MovingMeshModel() {}
162 
164 
165  void volChange()
166  {
167  const int numMeshes = _meshes.size();
168  for (int m=0;m<numMeshes;m++)
169  {
170  const Mesh& mesh = *_meshes[m];
171  const StorageSite& nodes = mesh.getNodes();
172  const StorageSite& faces = mesh.getFaces();
173  const StorageSite& cells = mesh.getCells();
174  const int nNodes = nodes.getCount();
175  const int nFaces = faces.getCount();
176  const int nCells = cells.getCount();
177  const CRConnectivity& faceNodes = mesh.getAllFaceNodes();
178  const CRConnectivity& faceCells = mesh.getAllFaceCells();
179  const TArray& rho =
180  dynamic_cast<const TArray&>(_flowFields.density[cells]);
181  const VectorT3Array& faceArea =
182  dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
183  VectorT3Array& faceAreaN1 =
184  dynamic_cast<VectorT3Array&>(_geomFields.areaN1[faces]);
185  const VectorT3Array& nodeCoord =
186  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[nodes]);
187  VectorT3Array& nodeCoordN1 =
188  dynamic_cast<VectorT3Array&>(_geomFields.coordinateN1[nodes]);
189  TArray& cellVolume =
190  dynamic_cast<TArray&>(_geomFields.volume[cells]);
191  TArray& gridFlux =
192  dynamic_cast<TArray&>(_geomFields.gridFlux[faces]);
193  VectorT3Array& faceVel =
194  dynamic_cast<VectorT3Array&>(_geomFields.faceVel[faces]);
195  TArray& sweptVolDot =
196  dynamic_cast<TArray&>(_geomFields.sweptVolDot[faces]);
197  shared_ptr<VectorT3Array> gridVelPtr(new VectorT3Array(nNodes));
198  VectorT3Array& gridVel = *gridVelPtr;
199  TArray volChangeDot(cells.getCount());
200 
201  const T deltaT = _options["timeStep"];
202 
203  if (mesh.getDimension() == 2)
204  {
205 
206  faceVel.zero();
207  gridVel.zero();
208  volChangeDot.zero();
209  sweptVolDot.zero();
210  gridFlux.zero();
211  // const T deltaT = _options["timeStep"];
212  const T half(0.5);
213  const T onepointfive(1.5);
214  for(int i=0;i<nNodes;i++)
215  {
216  VectorT3 dr = (nodeCoord[i] - nodeCoordN1[i]);
217  gridVel[i] = dr/deltaT;
218  }
219  for (int f=0;f<nFaces;f++)
220  {
221  const int numNodes = faceNodes.getCount(f);
222  for (int n=0;n<numNodes;n++)
223  {
224  const int num = faceNodes(f,n);
225  faceVel[f] += gridVel[num];
226  }
227  faceVel[f] /= T(numNodes);
228  T temp = dot((half*(faceArea[f]+faceAreaN1[f])),faceVel[f]);
229  sweptVolDot[f] = temp;
230  const int c0 = faceCells(f,0);
231  volChangeDot[c0] += temp;
232  const int c1 = faceCells(f,1);
233  volChangeDot[c1] -= temp;
235  {
236  TArray& sweptVolDotN1 =
237  dynamic_cast<TArray&>(_geomFields.sweptVolDotN1[faces]);
238  gridFlux[f] = onepointfive*sweptVolDot[f] -
239  half*sweptVolDotN1[f];
240  }
241  else
242  gridFlux[f] = sweptVolDot[f];
243  gridFlux[f] *= half*(rho[c0]+rho[c1]);
244  }
245  T volumeChange(0.0);
246  for(int c=0;c<cells.getSelfCount();c++)
247  volumeChange += volChangeDot[c];
248  volumeChange *= deltaT;
249  cout << "volume change for Mesh " << mesh.getID() << " = " << volumeChange << endl;
250  }
251  else if (mesh.getDimension() == 3)
252  {
253 
254  faceVel.zero();
255  gridVel.zero();
256  volChangeDot.zero();
257  sweptVolDot.zero();
258  gridFlux.zero();
259  // const T deltaT = _options["timeStep"];
260  const T half(0.5);
261  const T onepointfive(1.5);
262  const T onesixth(1./6.);
263  const T onethird(1./3.);
264  T temp(0.);
265  for(int i=0;i<nNodes;i++)
266  {
267  VectorT3 dr = (nodeCoord[i] - nodeCoordN1[i]);
268  gridVel[i] = dr/deltaT;
269  }
270  for (int f=0;f<nFaces;f++)
271  {
272  const int numNodes = faceNodes.getCount(f);
273  for (int n=0;n<numNodes;n++)
274  {
275  const int num = faceNodes(f,n);
276  faceVel[f] += gridVel[num];
277  }
278  faceVel[f] /= T(numNodes);
279  if (numNodes == 3)
280  {
281  const int n0 = faceNodes(f,0);
282  const int n1 = faceNodes(f,1);
283  const int n2 = faceNodes(f,2);
284  VectorT3 dr10 = nodeCoord[n1]-nodeCoord[n0];
285  VectorT3 dr20 = nodeCoord[n2]-nodeCoord[n0];
286  VectorT3 dr10N1 = nodeCoordN1[n1]-nodeCoordN1[n0];
287  VectorT3 dr20N1 = nodeCoordN1[n2]-nodeCoordN1[n0];
288  VectorT3 eta = onesixth*(cross(dr10,dr20)+cross(dr10N1,dr20N1)
289  +half*(cross(dr10,dr20N1)+cross(dr10N1,dr20)));
290  VectorT3 del0 = nodeCoord[n0]-nodeCoordN1[n0];
291  VectorT3 del1 = nodeCoord[n1]-nodeCoordN1[n1];
292  VectorT3 del2 = nodeCoord[n2]-nodeCoordN1[n2];
293  VectorT3 avg = onethird*(del0 + del1 + del2);
294  temp = dot(avg,eta)/deltaT;
295  }
296  else if (numNodes == 4)
297  {
298  const int n0 = faceNodes(f,0);
299  const int n1 = faceNodes(f,1);
300  const int n2 = faceNodes(f,2);
301  const int n3 = faceNodes(f,3);
302 
303  VectorT3 del0 = nodeCoord[n0]-nodeCoordN1[n0];
304  VectorT3 del1 = nodeCoord[n1]-nodeCoordN1[n1];
305  VectorT3 del2 = nodeCoord[n2]-nodeCoordN1[n2];
306  VectorT3 del3 = nodeCoord[n3]-nodeCoordN1[n3];
307 
308  VectorT3 dr10 = nodeCoord[n1]-nodeCoord[n0];
309  VectorT3 dr20 = nodeCoord[n2]-nodeCoord[n0];
310  VectorT3 dr10N1 = nodeCoordN1[n1]-nodeCoordN1[n0];
311  VectorT3 dr20N1 = nodeCoordN1[n2]-nodeCoordN1[n0];
312  VectorT3 eta1 = onesixth*(cross(dr10,dr20)+cross(dr10N1,dr20N1)
313  +half*(cross(dr10,dr20N1)+cross(dr10N1,dr20)));
314  VectorT3 avg1 = onethird*(del0+del1+del2);
315  T temp1 = dot(avg1,eta1)/deltaT;
316 
317  // VectorT3 dr20 = nodeCoord[n2]-nodeCoord[n0];
318  VectorT3 dr30 = nodeCoord[n3]-nodeCoord[n0];
319  // VectorT3 dr20N1 = nodeCoordN1[n2]-nodeCoordN1[n0];
320  VectorT3 dr30N1 = nodeCoordN1[n3]-nodeCoordN1[n0];
321  VectorT3 eta2 = onesixth*(cross(dr20,dr30)+cross(dr20N1,dr30N1)
322  +half*(cross(dr20,dr30N1)+cross(dr20N1,dr30)));
323  VectorT3 avg2 = onethird*(del0+del2+del3);
324  T temp2 = dot(avg2,eta2)/deltaT;
325 
326  temp = temp1+temp2;
327  }
328  sweptVolDot[f] = temp;
329  const int c0 = faceCells(f,0);
330  volChangeDot[c0] += temp;
331  const int c1 = faceCells(f,1);
332  volChangeDot[c1] -= temp;
334  {
335  TArray& sweptVolDotN1 =
336  dynamic_cast<TArray&>(_geomFields.sweptVolDotN1[faces]);
337  gridFlux[f] = onepointfive*sweptVolDot[f] -
338  half*sweptVolDotN1[f];
339  }
340  else
341  gridFlux[f] = sweptVolDot[f];
342  gridFlux[f] *= half*(rho[c0]+rho[c1]);
343  }
344  T volumeChange(0.0);
345  for(int c=0;c<cells.getSelfCount();c++)
346  volumeChange += volChangeDot[c];
347  volumeChange *= deltaT;
348  cout << "volume change for Mesh " << mesh.getID() << " = " << volumeChange << endl;
349 
350  }
351 
352  // update boundary cells with adjacent interior cells values
353  foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
354  {
355  const FaceGroup& fg = *fgPtr;
356  if(fg.groupType != "interior")
357  {
358  const StorageSite& bfaces = fg.site;
359  const CRConnectivity& bfaceCells = mesh.getFaceCells(bfaces);
360  const int faceCount = bfaces.getCount();
361  for(int f=0; f<faceCount; f++)
362  {
363  const int c0 = bfaceCells(f,0);
364  const int c1 = bfaceCells(f,1);
365  volChangeDot[c1] = volChangeDot[c0];
366  }
367  }
368  }
369  for (int c=0;c<nCells;c++)
370  cellVolume[c] += volChangeDot[c]*deltaT;
371  }
372  }
373 
374  void updateTime()
375  {
376  const int numMeshes = _meshes.size();
377  for (int n=0;n<numMeshes;n++)
378  {
379  const Mesh& mesh = *_meshes[n];
380  const StorageSite& cells = mesh.getCells();
381  const StorageSite& faces = mesh.getFaces();
382  const StorageSite& nodes = mesh.getNodes();
383  const TArray& cellVol =
384  dynamic_cast<TArray&>(_geomFields.volume[cells]);
385  TArray& cellVolN1 =
386  dynamic_cast<TArray&>(_geomFields.volumeN1[cells]);
387  VectorT3Array& nodeCoord =
388  dynamic_cast<VectorT3Array&>(_geomFields.coordinate[nodes]);
389  VectorT3Array& nodeCoordN1 =
390  dynamic_cast<VectorT3Array&>(_geomFields.coordinateN1[nodes]);
391  const VectorT3Array& faceArea =
392  dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
393  VectorT3Array& faceAreaN1 =
394  dynamic_cast<VectorT3Array&>(_geomFields.areaN1[faces]);
395  if (_options.timeDiscretizationOrder > 1)
396  {
397  TArray& cellVolN2 =
398  dynamic_cast<TArray&> (_geomFields.volumeN2[cells]);
399  cellVolN2 = cellVolN1;
400  TArray& sweptVolDot =
401  dynamic_cast<TArray&>(_geomFields.sweptVolDot[faces]);
402  TArray& sweptVolDotN1 =
403  dynamic_cast<TArray&>(_geomFields.sweptVolDotN1[faces]);
404  sweptVolDotN1 = sweptVolDot;
405  }
406  cellVolN1 = cellVol;
407  faceAreaN1 = faceArea;
408  nodeCoordN1 = nodeCoord;
409  }
410  }
411 
412  void init()
413  {
414  const int numMeshes = _meshes.size();
415  for (int n=0;n<numMeshes;n++)
416  {
417  const Mesh& mesh = *_meshes[n];
418  const StorageSite& faces = mesh.getFaces();
419  const StorageSite& cells = mesh.getCells();
420  const StorageSite& nodes = mesh.getNodes();
421  const TArray& cellVolume =
422  dynamic_cast<TArray&>(_geomFields.volume[cells]);
424  dynamic_pointer_cast<ArrayBase>(cellVolume.newCopy()));
425  VectorT3Array& nodeCoord =
426  dynamic_cast<VectorT3Array&>(_geomFields.coordinate[nodes]);
428  dynamic_pointer_cast<ArrayBase>(nodeCoord.newCopy()));
429  shared_ptr<TArray> sweptVolume(new TArray(faces.getCount()));
430  sweptVolume->zero();
431  _geomFields.sweptVolDot.addArray(faces,sweptVolume);
432  shared_ptr<TArray> gridFlux(new TArray(faces.getCount()));
433  gridFlux->zero();
434  _geomFields.gridFlux.addArray(faces,gridFlux);
435  shared_ptr<VectorT3Array> faceVel(new VectorT3Array(faces.getCount()));
436  faceVel->zero();
437  _geomFields.faceVel.addArray(faces,faceVel);
438  const VectorT3Array& faceArea =
439  dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
441  dynamic_pointer_cast<ArrayBase>(faceArea.newCopy()));
442  shared_ptr<Array<int> > displacementOptions(new Array<int>(nodes.getCount()));
443  *displacementOptions = 3;
444  _geomFields.displacementOptions.addArray(nodes,displacementOptions);
445  shared_ptr<VectorT3Array> dirichletNodeDisplacement
446  (new VectorT3Array(nodes.getCount()));
447  dirichletNodeDisplacement->zero();
448  _geomFields.dirichletNodeDisplacement.addArray(nodes,dirichletNodeDisplacement);
449  shared_ptr<VectorT3Array> nodeDisplacement
450  (new VectorT3Array(nodes.getCount()));
451  nodeDisplacement->zero();
452  _geomFields.nodeDisplacement.addArray(nodes,nodeDisplacement);
453  if (_options.timeDiscretizationOrder > 1)
454  {
456  dynamic_pointer_cast<ArrayBase>(cellVolume.newCopy()));
458  dynamic_pointer_cast<ArrayBase>(sweptVolume->newCopy()));
459  }
460  }
461  }
462 
463 private:
468 };
469 
470 
471 #endif
472 
const CRConnectivity & getAllFaceNodes() const
Definition: Mesh.cpp:368
int getCount(const int i) const
MovingMeshModelOptions< T > & getOptions()
virtual void zero()
Definition: Array.h:281
Field displacementOptions
Definition: GeomFields.h:37
shared_ptr< FaceGroup > FaceGroupPtr
Definition: Mesh.h:46
int getSelfCount() const
Definition: StorageSite.h:40
Field coordinate
Definition: GeomFields.h:19
bool hasArray(const StorageSite &s) const
Definition: Field.cpp:37
Definition: Mesh.h:28
virtual ~MovingMeshModel()
const FaceGroupList & getAllFaceGroups() const
Definition: Mesh.h:193
const StorageSite & getNodes() const
Definition: Mesh.h:110
Field areaN1
Definition: GeomFields.h:24
Field faceVel
Definition: GeomFields.h:32
Field volumeN2
Definition: GeomFields.h:28
Field coordinateN1
Definition: GeomFields.h:20
Field sweptVolDot
Definition: GeomFields.h:29
Field nodeDisplacement
Definition: GeomFields.h:33
Definition: Mesh.h:49
T mag(const Vector< T, 3 > &a)
Definition: Vector.h:260
Array< VectorT3 > VectorT3Array
Field gridFlux
Definition: GeomFields.h:31
MovingMeshModel(const MeshList &meshes, GeomFields &geomFields, FlowFields &flowFields)
Array< T > TArray
#define logCtor()
Definition: RLogInterface.h:26
string groupType
Definition: Mesh.h:42
MovingMeshModelOptions< T > _options
const CRConnectivity & getAllFaceCells() const
Definition: Mesh.cpp:378
GeomFields & _geomFields
const MeshList _meshes
const StorageSite & getBoundaryNodes() const
Definition: Mesh.cpp:325
const StorageSite & getFaces() const
Definition: Mesh.h:108
Field volumeN1
Definition: GeomFields.h:27
Definition: Model.h:13
const StorageSite & getCells() const
Definition: Mesh.h:109
Vector< T, 3 > cross(const Vector< T, 3 > &a, const Vector< T, 3 > &b)
Definition: Vector.h:242
Vector< T, 3 > VectorT3
Field volume
Definition: GeomFields.h:26
void addArray(const StorageSite &, shared_ptr< ArrayBase > a)
Definition: Field.cpp:72
const CRConnectivity & getFaceCells(const StorageSite &site) const
Definition: Mesh.cpp:388
Definition: Array.h:14
FlowFields & _flowFields
Field density
Definition: FlowFields.h:22
Field boundaryNodeNormal
Definition: GeomFields.h:35
shared_ptr< Array< int > > createAndGetBNglobalToLocal() const
Definition: Mesh.cpp:288
int getCount() const
Definition: StorageSite.h:39
virtual shared_ptr< IContainer > newCopy() const
Definition: Array.h:483
Field area
Definition: GeomFields.h:23
int getDimension() const
Definition: Mesh.h:105
T dot(const Vector< T, 3 > &a, const Vector< T, 3 > &b)
Definition: Vector.h:253
int getID() const
Definition: Mesh.h:106
const CRConnectivity & getCellNodes() const
Definition: Mesh.cpp:426
vector< Mesh * > MeshList
Definition: Mesh.h:439
Field dirichletNodeDisplacement
Definition: GeomFields.h:36
StorageSite site
Definition: Mesh.h:40
Field sweptVolDotN1
Definition: GeomFields.h:30