5 #ifndef _MOVINGMESHMODEL_H_
6 #define _MOVINGMESHMODEL_H_
31 const int numMeshes =
_meshes.size();
32 for(
int n=0;n<numMeshes;n++)
52 shared_ptr<CRConnectivity> nodeCellsPtr = cellNodes.getTranspose();
53 shared_ptr<CRConnectivity> nodeNodesPtr = nodeCellsPtr->multiply(cellNodes,
false);
55 nodeDisplacement.
zero();
57 const T underrelaxation =
_options[
"underrelaxation"];
60 for(
int i=0;i<
_options.nNodeDisplacementSweeps;i++)
64 T averageDirichletDisplacement(0.);
66 shared_ptr<VectorT3Array> _previousNodeDisplacementPtr =
69 for(
int j=0;j<nNodes;j++)
73 for(
int k=0;k<nodeNodes.
getCount(j);k++)
75 const int num = nodeNodes(j,k);
78 const VectorT3 ds = nodeCoordinate[num]-nodeCoordinate[j];
81 dr += nodeDisplacement[num]/
mag(ds);
82 weight += one/
mag(ds);
86 dr += nodeDisplacement[num]/small;
92 if(displacementOptions[j] == 0)
94 nodeDisplacement[j] = T(0.0);
95 nodeCoordinate[j] += nodeDisplacement[j] - (*(_previousNodeDisplacementPtr))[j];
97 else if(displacementOptions[j] == 1)
99 averageDirichletDisplacement +=
mag((dirichletNodeDisplacement)[j]);
102 nodeDisplacement[j] = (dirichletNodeDisplacement)[j];
103 nodeCoordinate[j] += nodeDisplacement[j] - (*(_previousNodeDisplacementPtr))[j];
105 else if(displacementOptions[j] == 2)
107 T temp =
dot(dr,nodeNormal[GlobalToLocal[j]]);
108 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];
115 else if(displacementOptions[j] == 3)
117 nodeDisplacement[j] = dr;
118 nodeDisplacement[j] = (*(_previousNodeDisplacementPtr))[j] +
119 underrelaxation*(nodeDisplacement[j]-(*(_previousNodeDisplacementPtr))[j]);
120 nodeCoordinate[j] += nodeDisplacement[j] - (*(_previousNodeDisplacementPtr))[j];
125 averageDirichletDisplacement /= nDirichlet;
127 averageDirichletDisplacement = T(1.);
129 T maxChangeInDisplacement = T(0.0);
130 for(
int j=0;j<nNodes;j++)
132 const T changeInDisplacement = (
mag(nodeDisplacement[j]-
133 (*(_previousNodeDisplacementPtr))[j]));
134 if (maxChangeInDisplacement < changeInDisplacement)
135 maxChangeInDisplacement = changeInDisplacement;
138 T maxChangeRelative = maxChangeInDisplacement / averageDirichletDisplacement;
141 if((maxChangeInDisplacement<=
_options.absTolerance)||
142 (maxChangeRelative<=
_options.relativeTolerance))
167 const int numMeshes =
_meshes.size();
168 for (
int m=0;m<numMeshes;m++)
174 const int nNodes = nodes.
getCount();
175 const int nFaces = faces.
getCount();
176 const int nCells = cells.
getCount();
197 shared_ptr<VectorT3Array> gridVelPtr(
new VectorT3Array(nNodes));
201 const T deltaT =
_options[
"timeStep"];
213 const T onepointfive(1.5);
214 for(
int i=0;i<nNodes;i++)
216 VectorT3 dr = (nodeCoord[i] - nodeCoordN1[i]);
217 gridVel[i] = dr/deltaT;
219 for (
int f=0;f<nFaces;f++)
221 const int numNodes = faceNodes.
getCount(f);
222 for (
int n=0;n<numNodes;n++)
224 const int num = faceNodes(f,n);
225 faceVel[f] += gridVel[num];
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;
238 gridFlux[f] = onepointfive*sweptVolDot[f] -
239 half*sweptVolDotN1[f];
242 gridFlux[f] = sweptVolDot[f];
243 gridFlux[f] *= half*(rho[c0]+rho[c1]);
247 volumeChange += volChangeDot[c];
248 volumeChange *= deltaT;
249 cout <<
"volume change for Mesh " << mesh.
getID() <<
" = " << volumeChange << endl;
261 const T onepointfive(1.5);
262 const T onesixth(1./6.);
263 const T onethird(1./3.);
265 for(
int i=0;i<nNodes;i++)
267 VectorT3 dr = (nodeCoord[i] - nodeCoordN1[i]);
268 gridVel[i] = dr/deltaT;
270 for (
int f=0;f<nFaces;f++)
272 const int numNodes = faceNodes.
getCount(f);
273 for (
int n=0;n<numNodes;n++)
275 const int num = faceNodes(f,n);
276 faceVel[f] += gridVel[num];
278 faceVel[f] /= T(numNodes);
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];
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;
296 else if (numNodes == 4)
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);
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];
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];
313 +half*(
cross(dr10,dr20N1)+
cross(dr10N1,dr20)));
314 VectorT3 avg1 = onethird*(del0+del1+del2);
315 T temp1 =
dot(avg1,eta1)/deltaT;
318 VectorT3 dr30 = nodeCoord[n3]-nodeCoord[n0];
320 VectorT3 dr30N1 = nodeCoordN1[n3]-nodeCoordN1[n0];
322 +half*(
cross(dr20,dr30N1)+
cross(dr20N1,dr30)));
323 VectorT3 avg2 = onethird*(del0+del2+del3);
324 T temp2 =
dot(avg2,eta2)/deltaT;
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;
337 gridFlux[f] = onepointfive*sweptVolDot[f] -
338 half*sweptVolDotN1[f];
341 gridFlux[f] = sweptVolDot[f];
342 gridFlux[f] *= half*(rho[c0]+rho[c1]);
346 volumeChange += volChangeDot[c];
347 volumeChange *= deltaT;
348 cout <<
"volume change for Mesh " << mesh.
getID() <<
" = " << volumeChange << endl;
360 const int faceCount = bfaces.
getCount();
361 for(
int f=0; f<faceCount; f++)
363 const int c0 = bfaceCells(f,0);
364 const int c1 = bfaceCells(f,1);
365 volChangeDot[c1] = volChangeDot[c0];
369 for (
int c=0;c<nCells;c++)
370 cellVolume[c] += volChangeDot[c]*deltaT;
376 const int numMeshes =
_meshes.size();
377 for (
int n=0;n<numMeshes;n++)
395 if (
_options.timeDiscretizationOrder > 1)
399 cellVolN2 = cellVolN1;
404 sweptVolDotN1 = sweptVolDot;
407 faceAreaN1 = faceArea;
408 nodeCoordN1 = nodeCoord;
414 const int numMeshes =
_meshes.size();
415 for (
int n=0;n<numMeshes;n++)
421 const TArray& cellVolume =
424 dynamic_pointer_cast<ArrayBase>(cellVolume.
newCopy()));
428 dynamic_pointer_cast<ArrayBase>(nodeCoord.
newCopy()));
441 dynamic_pointer_cast<ArrayBase>(faceArea.
newCopy()));
443 *displacementOptions = 3;
445 shared_ptr<VectorT3Array> dirichletNodeDisplacement
447 dirichletNodeDisplacement->zero();
449 shared_ptr<VectorT3Array> nodeDisplacement
451 nodeDisplacement->zero();
453 if (
_options.timeDiscretizationOrder > 1)
456 dynamic_pointer_cast<ArrayBase>(cellVolume.
newCopy()));
458 dynamic_pointer_cast<ArrayBase>(sweptVolume->newCopy()));
const CRConnectivity & getAllFaceNodes() const
int getCount(const int i) const
MovingMeshModelOptions< T > & getOptions()
Field displacementOptions
shared_ptr< FaceGroup > FaceGroupPtr
bool hasArray(const StorageSite &s) const
virtual ~MovingMeshModel()
const FaceGroupList & getAllFaceGroups() const
const StorageSite & getNodes() const
T mag(const Vector< T, 3 > &a)
Array< VectorT3 > VectorT3Array
MovingMeshModel(const MeshList &meshes, GeomFields &geomFields, FlowFields &flowFields)
MovingMeshModelOptions< T > _options
const CRConnectivity & getAllFaceCells() const
const StorageSite & getBoundaryNodes() const
const StorageSite & getFaces() const
const StorageSite & getCells() const
Vector< T, 3 > cross(const Vector< T, 3 > &a, const Vector< T, 3 > &b)
void addArray(const StorageSite &, shared_ptr< ArrayBase > a)
const CRConnectivity & getFaceCells(const StorageSite &site) const
shared_ptr< Array< int > > createAndGetBNglobalToLocal() const
virtual shared_ptr< IContainer > newCopy() const
T dot(const Vector< T, 3 > &a, const Vector< T, 3 > &b)
const CRConnectivity & getCellNodes() const
vector< Mesh * > MeshList
Field dirichletNodeDisplacement