Memosa-FVM  0.2
Public Types | Public Member Functions | Private Member Functions | Private Attributes | List of all members
MeshMetricsCalculator< T > Class Template Reference

#include <MeshMetricsCalculator.h>

Inheritance diagram for MeshMetricsCalculator< T >:
Inheritance graph
[legend]
Collaboration diagram for MeshMetricsCalculator< T >:
Collaboration graph
[legend]

Public Types

typedef Vector< T, 3 > VectorT3
 
typedef Array< T > TArray
 
typedef Array< int > IntArray
 
typedef Array< VectorT3 > VectorT3Array
 

Public Member Functions

 MeshMetricsCalculator (GeomFields &geomFields, const MeshList &meshes, bool transient=false)
 
virtual ~MeshMetricsCalculator ()
 
virtual void init ()
 
void createNodeDisplacement ()
 
void calculateBoundaryNodeNormal ()
 
void recalculate ()
 
void recalculate_deform ()
 
void computeIBInterpolationMatrices (const StorageSite &particles, const int option=0)
 
void computeIBInterpolationMatricesCells ()
 
void eraseIBInterpolationMatrices (const StorageSite &particles)
 
void computeSolidInterpolationMatrices (const StorageSite &particles)
 
void computeIBandSolidInterpolationMatrices (const StorageSite &particles)
 
void computeGridInterpolationMatrices (const StorageSite &grids, const StorageSite &faces)
 
void updateTime ()
 
- Public Member Functions inherited from Model
 Model (const MeshList &meshes)
 
virtual ~Model ()
 
 DEFINE_TYPENAME ("Model")
 
virtual map< string,
shared_ptr< ArrayBase > > & 		getPersistenceData ()
 
virtual void restart ()
 

Private Member Functions

virtual void calculateNodeCoordinates (const Mesh &mesh)
 
virtual void calculateFaceCentroids (const Mesh &mesh)
 
virtual void calculateCellCentroids (const Mesh &mesh)
 
virtual void calculateFaceAreas (const Mesh &mesh)
 
virtual void calculateFaceAreaMag (const Mesh &mesh)
 
virtual void calculateCellVolumes (const Mesh &mesh)
 
void computeIBInterpolationMatrices (const Mesh &mesh, const StorageSite &particles, const int option)
 
void computeIBInterpolationMatricesCells (const Mesh &mesh)
 
void computeSolidInterpolationMatrices (const Mesh &mesh, const StorageSite &particles)
 
void computeIBandSolidInterpolationMatrices (const Mesh &mesh, const StorageSite &particles)
 
void computeGridInterpolationMatrices (const Mesh &mesh, const StorageSite &grids, const StorageSite &faces)
 

Private Attributes

GeomFields & _geomFields
 
Field & _coordField
 
Field & _areaField
 
Field & _areaMagField
 
Field & _volumeField
 
Field & _nodeDisplacement
 
Field & _boundaryNodeNormal
 
bool _transient
 

Additional Inherited Members

- Protected Attributes inherited from Model
const MeshList _meshes
 
StorageSiteList _varSites
 
StorageSiteList _fluxSites
 
map< string, shared_ptr
< ArrayBase > > 		_persistenceData
 

Detailed Description

template<class T>
class MeshMetricsCalculator< T >

Definition at line 22 of file MeshMetricsCalculator.h.

Member Typedef Documentation

template<class T >
typedef Array<int> MeshMetricsCalculator< T >::IntArray

Definition at line 28 of file MeshMetricsCalculator.h.

template<class T >
typedef Array<T> MeshMetricsCalculator< T >::TArray

Definition at line 27 of file MeshMetricsCalculator.h.

template<class T >
typedef Vector<T,3> MeshMetricsCalculator< T >::VectorT3

Definition at line 26 of file MeshMetricsCalculator.h.

template<class T >
typedef Array<VectorT3> MeshMetricsCalculator< T >::VectorT3Array

Definition at line 29 of file MeshMetricsCalculator.h.

Constructor & Destructor Documentation

template<class T >
MeshMetricsCalculator< T >::MeshMetricsCalculator ( GeomFields &  geomFields,
const MeshList &  meshes,
bool  transient = false 
)

Definition at line 1928 of file MeshMetricsCalculator_impl.h.

References logCtor.

1929  :
1930  Model(meshes),
1931  _geomFields(geomFields),
1932  _coordField(geomFields.coordinate),
1933  _areaField(geomFields.area),
1934  _areaMagField(geomFields.areaMag),
1935  _volumeField(geomFields.volume),
1936  _nodeDisplacement(geomFields.nodeDisplacement),
1937  _boundaryNodeNormal(geomFields.boundaryNodeNormal),
1938  _transient(transient)
1939 {
1940  logCtor();
1941 }
GeomFields::coordinate
Field coordinate
Definition: GeomFields.h:19
Model::Model
Model(const MeshList &meshes)
Definition: Model.cpp:8
GeomFields::nodeDisplacement
Field nodeDisplacement
Definition: GeomFields.h:33
logCtor
#define logCtor()
Definition: RLogInterface.h:26
GeomFields::volume
Field volume
Definition: GeomFields.h:26
GeomFields::boundaryNodeNormal
Field boundaryNodeNormal
Definition: GeomFields.h:35
GeomFields::area
Field area
Definition: GeomFields.h:23
GeomFields::areaMag
Field areaMag
Definition: GeomFields.h:25
template<class T >
MeshMetricsCalculator< T >::~MeshMetricsCalculator ( )
virtual

Definition at line 2454 of file MeshMetricsCalculator_impl.h.

References logDtor.

2455 {
2456  logDtor();
2457 }
logDtor
#define logDtor()
Definition: RLogInterface.h:33

Member Function Documentation

template<class T >
void MeshMetricsCalculator< T >::calculateBoundaryNodeNormal ( )

Definition at line 308 of file MeshMetricsCalculator_impl.h.

References Mesh::createAndGetBNglobalToLocal(), Mesh::getAllFaceGroups(), Mesh::getBoundaryNodes(), StorageSite::getCount(), CRConnectivity::getCount(), Mesh::getFaceNodes(), FaceGroup::groupType, FaceGroup::site, and Array< T >::zero().

309 {
310  const int numMeshes = _meshes.size();
311  for (int n=0; n<numMeshes; n++)
312  {
313  const Mesh& mesh = *_meshes[n];
314  shared_ptr<Array<int> > GlobalToLocalPtr = mesh.createAndGetBNglobalToLocal();
315  Array<int>& GlobalToLocal = *GlobalToLocalPtr;
316  const StorageSite& boundaryNodes = mesh.getBoundaryNodes();
317  const int nBoundaryNodes = boundaryNodes.getCount();
318  Array<bool> nodeMark(nBoundaryNodes);
319  Array<bool> nodeMarked(nBoundaryNodes);
320  shared_ptr<VectorT3Array> nodeNormalPtr(new VectorT3Array(nBoundaryNodes));
321  VectorT3Array& nodeNormal = *nodeNormalPtr;
322  TArray number(nBoundaryNodes);
323  nodeNormal.zero();
324  number.zero();
325  const T one(1.0);
326  for(int i=0;i<nBoundaryNodes;i++)
327  nodeMarked[i] = false;
328  foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
329  {
330  const FaceGroup& fg = *fgPtr;
331  if (fg.groupType != "interior")
332  {
333  const StorageSite& faces = fg.site;
334  const int nFaces = faces.getCount();
335  const CRConnectivity& BoundaryFaceNodes = mesh.getFaceNodes(faces);
336  const VectorT3Array& faceArea = dynamic_cast<const VectorT3Array&>(_areaField[faces]);
337  const TArray& faceAreaMag = dynamic_cast<const TArray&>(_areaMagField[faces]);
338  for(int i=0;i<nBoundaryNodes;i++)
339  nodeMark[i] = false;
340  for(int i=0;i<nFaces;i++)
341  {
342  for(int j=0;j<BoundaryFaceNodes.getCount(i);j++)
343  {
344  const int num = BoundaryFaceNodes(i,j);
345  if(!nodeMark[GlobalToLocal[num]])
346  nodeMark[GlobalToLocal[num]] = true;
347  if((nodeMark[GlobalToLocal[num]])&&(!(nodeMarked[GlobalToLocal[num]])))
348  {
349  nodeNormal[GlobalToLocal[num]] += faceArea[i]/faceAreaMag[i];
350  number[GlobalToLocal[num]] += one;
351  }
352  }
353  }
354  for(int i=0;i<nBoundaryNodes;i++)
355  {
356  if((nodeMark[i])&&(!nodeMarked[i]))
357  {
358  nodeNormal[i][0] = nodeNormal[i][0]/number[i];
359  nodeNormal[i][1] = nodeNormal[i][1]/number[i];
360  nodeNormal[i][2] = nodeNormal[i][2]/number[i];
361  nodeMarked[i] = true;
362  }
363  }
364  }
365  }
366  _boundaryNodeNormal.addArray(boundaryNodes,nodeNormalPtr);
367  }
368 }
CRConnectivity::getCount
int getCount(const int i) const
Definition: CRConnectivity.h:165
Array::zero
virtual void zero()
Definition: Array.h:281
FaceGroupPtr
shared_ptr< FaceGroup > FaceGroupPtr
Definition: Mesh.h:46
MeshMetricsCalculator::VectorT3Array
Array< VectorT3 > VectorT3Array
Definition: MeshMetricsCalculator.h:29
FaceGroup
Definition: Mesh.h:28
Mesh::getAllFaceGroups
const FaceGroupList & getAllFaceGroups() const
Definition: Mesh.h:193
Mesh
Definition: Mesh.h:49
Mesh::getFaceNodes
const CRConnectivity & getFaceNodes(const StorageSite &site) const
Definition: Mesh.cpp:402
FaceGroup::groupType
string groupType
Definition: Mesh.h:42
Mesh::getBoundaryNodes
const StorageSite & getBoundaryNodes() const
Definition: Mesh.cpp:325
Model::_meshes
const MeshList _meshes
Definition: Model.h:29
Field::addArray
void addArray(const StorageSite &, shared_ptr< ArrayBase > a)
Definition: Field.cpp:72
Mesh::createAndGetBNglobalToLocal
shared_ptr< Array< int > > createAndGetBNglobalToLocal() const
Definition: Mesh.cpp:288
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
FaceGroup::site
StorageSite site
Definition: Mesh.h:40
template<class T >
void MeshMetricsCalculator< T >::calculateCellCentroids ( const Mesh &  mesh)
privatevirtual

calculate cell centroids. Needs to have the face centroid and areas computed already.

Definition at line 130 of file MeshMetricsCalculator_impl.h.

References dot(), Mesh::getAllFaceCells(), Mesh::getAllFaceGroups(), Mesh::getCells(), StorageSite::getCount(), CRConnectivity::getCount(), StorageSite::getCountLevel1(), Mesh::getFaceCells(), Mesh::getFaces(), Mesh::getParentFaceGroupSite(), StorageSite::getSelfCount(), FaceGroup::groupType, Mesh::isDoubleShell(), Mesh::isShell(), FaceGroup::site, and Array< T >::zero().

131 {
132  const StorageSite& cells = mesh.getCells();
133 
134 
135  const int cellCount = cells.getCountLevel1();
136  if (cellCount == 0)
137  return;
138 
139 
140  const int selfCellCount = cells.getSelfCount();
141 
142  shared_ptr<VectorT3Array> ccPtr(new VectorT3Array(cellCount));
143  VectorT3Array& cellCentroid = *ccPtr;
144  _coordField.addArray(cells,ccPtr);
145 
146 
147  // for shell mesh copy cell centroids from the faces it was created from
148  if (mesh.isShell())
149  {
150  const StorageSite& faces = mesh.getParentFaceGroupSite();
151  const int faceCount = faces.getCount();
152  const VectorT3Array& faceCentroid =
153  dynamic_cast<const VectorT3Array&>(_coordField[faces]);
154 
155  for(int f=0; f<faceCount; f++)
156  {
157  cellCentroid[f] = faceCentroid[f];
158 
159  //additional set of cells if double shell
160  if (mesh.isDoubleShell())
161  {
162  cellCentroid[f+faceCount] = faceCentroid[f];
163  }
164  }
165  return;
166  }
167 
168  const StorageSite& faces = mesh.getFaces();
169 
170  const VectorT3Array& faceCentroid =
171  dynamic_cast<const VectorT3Array&>(_coordField[faces]);
172  const TArray& faceAreaMag = dynamic_cast<const TArray&>(_areaMagField[faces]);
173  const int faceCount = faces.getCount();
174  const CRConnectivity& faceCells = mesh.getAllFaceCells();
175 
176  TArray weight(cellCount);
177 
178  weight.zero();
179  cellCentroid.zero();
180 
181  for(int f=0; f<faceCount; f++)
182  {
183  for(int nc=0; nc<faceCells.getCount(f); nc++)
184  {
185  const int c = faceCells(f,nc);
186  cellCentroid[c] += faceCentroid[f]*faceAreaMag[f];
187  weight[c] += faceAreaMag[f];
188  }
189  }
190 
191  for(int c=0; c<selfCellCount; c++)
192  cellCentroid[c] /= weight[c];
193 
194  // boundary cells have the corresponding face's centroid
195  foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
196  {
197  const FaceGroup& fg = *fgPtr;
198  const StorageSite& faces = fg.site;
199  const VectorT3Array& faceCentroid =
200  dynamic_cast<const VectorT3Array&>(_coordField[faces]);
201  const VectorT3Array& faceArea =
202  dynamic_cast<const VectorT3Array&>(_areaField[faces]);
203  const TArray& faceAreaMag =
204  dynamic_cast<const TArray&>(_areaMagField[faces]);
205  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
206  const int faceCount = faces.getCount();
207 
208  if ((fg.groupType!="interior") && (fg.groupType!="interface") &&
209  (fg.groupType!="dielectric interface"))
210  {
211  if (fg.groupType == "symmetry")
212  {
213  for(int f=0; f<faceCount; f++)
214  {
215  const int c0 = faceCells(f,0);
216  const int c1 = faceCells(f,1);
217  const VectorT3 en = faceArea[f]/faceAreaMag[f];
218  const VectorT3 dr0(faceCentroid[f]-cellCentroid[c0]);
219 
220  const T dr0_dotn = dot(dr0,en);
221  const VectorT3 dr1 = dr0 - 2.*dr0_dotn*en;
222  cellCentroid[c1] = cellCentroid[c0]+dr0-dr1;
223  }
224 
225  }
226  else
227  {
228  for(int f=0; f<faceCount; f++)
229  {
230  const int c1 = faceCells(f,1);
231  cellCentroid[c1] = faceCentroid[f];
232  }
233  }
234  }
235  }
236 }
CRConnectivity::getCount
int getCount(const int i) const
Definition: CRConnectivity.h:165
Array::zero
virtual void zero()
Definition: Array.h:281
Mesh::isDoubleShell
bool isDoubleShell() const
Definition: Mesh.h:324
FaceGroupPtr
shared_ptr< FaceGroup > FaceGroupPtr
Definition: Mesh.h:46
MeshMetricsCalculator::VectorT3Array
Array< VectorT3 > VectorT3Array
Definition: MeshMetricsCalculator.h:29
StorageSite::getSelfCount
int getSelfCount() const
Definition: StorageSite.h:40
Mesh::getParentFaceGroupSite
const StorageSite & getParentFaceGroupSite() const
Definition: Mesh.h:329
FaceGroup
Definition: Mesh.h:28
Mesh::getAllFaceGroups
const FaceGroupList & getAllFaceGroups() const
Definition: Mesh.h:193
FaceGroup::groupType
string groupType
Definition: Mesh.h:42
Mesh::getAllFaceCells
const CRConnectivity & getAllFaceCells() const
Definition: Mesh.cpp:378
Mesh::getFaces
const StorageSite & getFaces() const
Definition: Mesh.h:108
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
StorageSite::getCountLevel1
int getCountLevel1() const
Definition: StorageSite.h:72
Field::addArray
void addArray(const StorageSite &, shared_ptr< ArrayBase > a)
Definition: Field.cpp:72
Mesh::getFaceCells
const CRConnectivity & getFaceCells(const StorageSite &site) const
Definition: Mesh.cpp:388
Mesh::isShell
bool isShell() const
Definition: Mesh.h:323
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
dot
T dot(const Vector< T, 3 > &a, const Vector< T, 3 > &b)
Definition: Vector.h:253
FaceGroup::site
StorageSite site
Definition: Mesh.h:40
template<class T >
void MeshMetricsCalculator< T >::calculateCellVolumes ( const Mesh &  mesh)
privatevirtual

Definition at line 394 of file MeshMetricsCalculator_impl.h.

References dot(), Mesh::getAllFaceCells(), Mesh::getAllFaceGroups(), Mesh::getCells(), StorageSite::getCount(), StorageSite::getCountLevel1(), Mesh::getDimension(), Mesh::getFaceCells(), Mesh::getFaces(), StorageSite::getSelfCount(), FaceGroup::groupType, Mesh::isShell(), FaceGroup::site, and Array< T >::zero().

395 {
396  const StorageSite& faces = mesh.getFaces();
397  const StorageSite& cells = mesh.getCells();
398 
399  const int cellCount = cells.getCountLevel1();
400  if (cellCount == 0)
401  return;
402 
403  shared_ptr<TArray> vPtr(new TArray(cellCount));
404  TArray& cellVolume = *vPtr;
405 
406  cellVolume.zero();
407  _volumeField.addArray(cells,vPtr);
408 
409  // for shell the volume is zero
410  if (mesh.isShell())
411  return;
412 
413  const VectorT3Array& faceCentroid =
414  dynamic_cast<const VectorT3Array&>(_coordField[faces]);
415  const VectorT3Array& cellCentroid =
416  dynamic_cast<const VectorT3Array&>(_coordField[cells]);
417  const VectorT3Array& faceArea =
418  dynamic_cast<const VectorT3Array&>(_areaField[faces]);
419  const CRConnectivity& faceCells = mesh.getAllFaceCells();
420 
421  const T dim(mesh.getDimension());
422  const int faceCount = faces.getCount();
423 
424 
425 
426  for(int f=0; f<faceCount; f++)
427  {
428  const int c0 = faceCells(f,0);
429  cellVolume[c0] += dot(faceCentroid[f]-cellCentroid[c0],faceArea[f])/dim;
430  const int c1 = faceCells(f,1);
431  cellVolume[c1] -= dot(faceCentroid[f]-cellCentroid[c1],faceArea[f])/dim;
432  }
433 
434  T volumeSum(0);
435  for(int c=0; c<cells.getSelfCount(); c++)
436  volumeSum += cellVolume[c];
437 
438  //cout << "volume sum for Mesh " << mesh.getID() << " = " << volumeSum << endl;
439 
440 
441  // update boundary cells with adjacent interior cells values
442  foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
443  {
444  const FaceGroup& fg = *fgPtr;
445  if ((fg.groupType!="interior") && (fg.groupType!="interface") &&
446  (fg.groupType!="dielectric interface"))
447  {
448  const StorageSite& faces = fg.site;
449  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
450  const int faceCount = faces.getCount();
451 
452  for(int f=0; f<faceCount; f++)
453  {
454  const int c0 = faceCells(f,0);
455  const int c1 = faceCells(f,1);
456  cellVolume[c1] = cellVolume[c0];
457  }
458  }
459  }
460 }
FaceGroupPtr
shared_ptr< FaceGroup > FaceGroupPtr
Definition: Mesh.h:46
StorageSite::getSelfCount
int getSelfCount() const
Definition: StorageSite.h:40
FaceGroup
Definition: Mesh.h:28
Mesh::getAllFaceGroups
const FaceGroupList & getAllFaceGroups() const
Definition: Mesh.h:193
FaceGroup::groupType
string groupType
Definition: Mesh.h:42
Mesh::getAllFaceCells
const CRConnectivity & getAllFaceCells() const
Definition: Mesh.cpp:378
Mesh::getFaces
const StorageSite & getFaces() const
Definition: Mesh.h:108
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
StorageSite::getCountLevel1
int getCountLevel1() const
Definition: StorageSite.h:72
Field::addArray
void addArray(const StorageSite &, shared_ptr< ArrayBase > a)
Definition: Field.cpp:72
Mesh::getFaceCells
const CRConnectivity & getFaceCells(const StorageSite &site) const
Definition: Mesh.cpp:388
Mesh::isShell
bool isShell() const
Definition: Mesh.h:323
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
Mesh::getDimension
int getDimension() const
Definition: Mesh.h:105
dot
T dot(const Vector< T, 3 > &a, const Vector< T, 3 > &b)
Definition: Vector.h:253
FaceGroup::site
StorageSite site
Definition: Mesh.h:40
template<class T >
void MeshMetricsCalculator< T >::calculateFaceAreaMag ( const Mesh &  mesh)
privatevirtual

Definition at line 373 of file MeshMetricsCalculator_impl.h.

References StorageSite::getCount(), Mesh::getFaces(), and mag().

374 {
375  const StorageSite& faces = mesh.getFaces();
376  const VectorT3Array& faceArea =
377  dynamic_cast<const VectorT3Array&>(_areaField[faces]);
378 
379  const int count = faces.getCount();
380  shared_ptr<TArray> famPtr(new TArray(count));
381  TArray& fam = *famPtr;
382 
383  for(int f=0; f<count; f++)
384  {
385  fam[f] = mag(faceArea[f]);
386  }
387 
388  _areaMagField.addArray(faces,famPtr);
389 }
mag
T mag(const Vector< T, 3 > &a)
Definition: Vector.h:260
Mesh::getFaces
const StorageSite & getFaces() const
Definition: Mesh.h:108
Field::addArray
void addArray(const StorageSite &, shared_ptr< ArrayBase > a)
Definition: Field.cpp:72
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
template<class T >
void MeshMetricsCalculator< T >::calculateFaceAreas ( const Mesh &  mesh)
privatevirtual

Definition at line 240 of file MeshMetricsCalculator_impl.h.

References cross(), Mesh::getAllFaceNodes(), StorageSite::getCount(), CRConnectivity::getCount(), Mesh::getFaces(), Mesh::getNodes(), and Array< T >::zero().

241 {
242  const StorageSite& faces = mesh.getFaces();
243  const StorageSite& nodes = mesh.getNodes();
244  const CRConnectivity& faceNodes = mesh.getAllFaceNodes();
245 
246  const int count = faces.getCount();
247  shared_ptr<VectorT3Array> faPtr(new VectorT3Array(count));
248  VectorT3Array& fa = *faPtr;
249 
250  const T half(0.5);
251  const VectorT3Array& nodeCoord =
252  dynamic_cast<const VectorT3Array&>(_coordField[nodes]);
253 
254  for(int f=0; f<count; f++)
255  {
256  const int numNodes = faceNodes.getCount(f);
257 
258  if (numNodes == 2)
259  {
260  const int n0 = faceNodes(f,0);
261  const int n1 = faceNodes(f,1);
262  VectorT3 dr = nodeCoord[n1]-nodeCoord[n0];
263  fa[f][0] = dr[1];
264  fa[f][1] = -dr[0];
265  fa[f][2] = 0.;
266  }
267  else if (numNodes == 3)
268  {
269  const int n0 = faceNodes(f,0);
270  const int n1 = faceNodes(f,1);
271  const int n2 = faceNodes(f,2);
272  VectorT3 dr10 = nodeCoord[n1]-nodeCoord[n0];
273  VectorT3 dr20 = nodeCoord[n2]-nodeCoord[n0];
274  fa[f] = half*cross(dr10,dr20);
275  }
276  else if (numNodes == 4)
277  {
278  const int n0 = faceNodes(f,0);
279  const int n1 = faceNodes(f,1);
280  const int n2 = faceNodes(f,2);
281  const int n3 = faceNodes(f,3);
282  VectorT3 dr20 = nodeCoord[n2]-nodeCoord[n0];
283  VectorT3 dr31 = nodeCoord[n3]-nodeCoord[n1];
284  fa[f] = half*cross(dr20,dr31);
285  }
286  else if (numNodes > 0)
287  {
288  fa[f].zero();
289  for (int nn=0; nn<numNodes; nn++)
290  {
291  const int n0 = faceNodes(f,nn);
292  const int n1 = faceNodes(f,(nn+1)%numNodes);
293  VectorT3 xm = T(0.5)*(nodeCoord[n1]+nodeCoord[n0]);
294  VectorT3 dr = (nodeCoord[n1]-nodeCoord[n0]);
295 
296  fa[f][0] += xm[1]*dr[2];
297  fa[f][1] += xm[2]*dr[0];
298  fa[f][2] += xm[0]*dr[1];
299  }
300  }
301  }
302 
303  _areaField.addArray(faces,faPtr);
304 }
Mesh::getAllFaceNodes
const CRConnectivity & getAllFaceNodes() const
Definition: Mesh.cpp:368
CRConnectivity::getCount
int getCount(const int i) const
Definition: CRConnectivity.h:165
Array::zero
virtual void zero()
Definition: Array.h:281
MeshMetricsCalculator::VectorT3Array
Array< VectorT3 > VectorT3Array
Definition: MeshMetricsCalculator.h:29
Mesh::getNodes
const StorageSite & getNodes() const
Definition: Mesh.h:110
Mesh::getFaces
const StorageSite & getFaces() const
Definition: Mesh.h:108
cross
Vector< T, 3 > cross(const Vector< T, 3 > &a, const Vector< T, 3 > &b)
Definition: Vector.h:242
Field::addArray
void addArray(const StorageSite &, shared_ptr< ArrayBase > a)
Definition: Field.cpp:72
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
template<class T >
void MeshMetricsCalculator< T >::calculateFaceCentroids ( const Mesh &  mesh)
privatevirtual

calculates the face centroids. Needs face area and magnitude for non-planar corrections.

Definition at line 60 of file MeshMetricsCalculator_impl.h.

References cross(), dot(), Mesh::getAllFaceNodes(), StorageSite::getCount(), CRConnectivity::getCount(), Mesh::getFaces(), and Mesh::getNodes().

61 {
62  const StorageSite& faces = mesh.getFaces();
63  const StorageSite& nodes = mesh.getNodes();
64  const CRConnectivity& faceNodes = mesh.getAllFaceNodes();
65 
66  const int count = faces.getCount();
67  shared_ptr<VectorT3Array> fcPtr(new VectorT3Array(count));
68  VectorT3Array& fc = *fcPtr;
69 
70  const VectorT3Array& nodeCoord =
71  dynamic_cast<const VectorT3Array&>(_coordField[nodes]);
72  for(int f=0; f<count; f++)
73  {
74  const int numNodes = faceNodes.getCount(f);
75  if (numNodes > 0)
76  {
77  fc[f] = nodeCoord[faceNodes(f,0)];
78  for (int nf=1; nf<numNodes; nf++)
79  fc[f] += nodeCoord[faceNodes(f,nf)];
80  fc[f] /= T(numNodes);
81  }
82  }
83 
84  const VectorT3Array& faceArea =
85  dynamic_cast<const VectorT3Array&>(_areaField[faces]);
86  const TArray& faceAreaMag =
87  dynamic_cast<const TArray&>(_areaMagField[faces]);
88 
89  // corrections for non-planar quad and polygonal faces
90  const T twoThirds(2./3.);
91  const T half(0.5);
92  for(int f=0; f<count; f++)
93  {
94  const int numNodes = faceNodes.getCount(f);
95  if (numNodes > 3)
96  {
97  const VectorT3 en = faceArea[f]/faceAreaMag[f];
98  T denom(0.0);
99  VectorT3 cfc(NumTypeTraits<VectorT3>::getZero());
100 
101  for (int nn=0; nn<numNodes; nn++)
102  {
103  const int n0 = faceNodes(f,nn);
104  const int n1 = faceNodes(f,(nn+1)%numNodes);
105  const VectorT3 rc0 = nodeCoord[n0] - fc[f];
106  const VectorT3 rc1 = nodeCoord[n1] - fc[f];
107  const VectorT3 triArea = half*cross(rc0,rc1);
108  const T triAreaP = dot(triArea,en);
109  VectorT3 xm = half*(nodeCoord[n0]+nodeCoord[n1]);
110 
111  cfc += twoThirds*(xm-fc[f])*triAreaP;
112  denom += triAreaP;
113  }
114  cfc /= denom;
115  fc[f] += cfc;
116  }
117  }
118 
119  _coordField.addArray(faces,fcPtr);
120 }
Mesh::getAllFaceNodes
const CRConnectivity & getAllFaceNodes() const
Definition: Mesh.cpp:368
CRConnectivity::getCount
int getCount(const int i) const
Definition: CRConnectivity.h:165
MeshMetricsCalculator::VectorT3Array
Array< VectorT3 > VectorT3Array
Definition: MeshMetricsCalculator.h:29
Mesh::getNodes
const StorageSite & getNodes() const
Definition: Mesh.h:110
Mesh::getFaces
const StorageSite & getFaces() const
Definition: Mesh.h:108
cross
Vector< T, 3 > cross(const Vector< T, 3 > &a, const Vector< T, 3 > &b)
Definition: Vector.h:242
Field::addArray
void addArray(const StorageSite &, shared_ptr< ArrayBase > a)
Definition: Field.cpp:72
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
dot
T dot(const Vector< T, 3 > &a, const Vector< T, 3 > &b)
Definition: Vector.h:253
template<class T >
void MeshMetricsCalculator< T >::calculateNodeCoordinates ( const Mesh &  mesh)
privatevirtual

The mesh always has coordinates in double. The coordField that the rest of the code uses to get the coordinates needs to have them in the current AType so this function just creates a copy in that AType.

Definition at line 36 of file MeshMetricsCalculator_impl.h.

References StorageSite::getCount(), Mesh::getNodeCoordinates(), and Mesh::getNodes().

37 {
38  const StorageSite& nodes = mesh.getNodes();
39  const Array<Vector<double,3> >& meshCoords = mesh.getNodeCoordinates();
40 
41  const int count = nodes.getCount();
42  shared_ptr<VectorT3Array> ncPtr(new VectorT3Array(count));
43  VectorT3Array& nc = *ncPtr;
44 
45  for(int i=0; i<count; i++)
46  for(int j=0; j<3; j++)
47  nc[i][j] = T(meshCoords[i][j]);
48 
49  _coordField.addArray(nodes,ncPtr);
50 }
MeshMetricsCalculator::VectorT3Array
Array< VectorT3 > VectorT3Array
Definition: MeshMetricsCalculator.h:29
Mesh::getNodes
const StorageSite & getNodes() const
Definition: Mesh.h:110
Mesh::getNodeCoordinates
const Array< VecD3 > & getNodeCoordinates() const
Definition: Mesh.h:218
Field::addArray
void addArray(const StorageSite &, shared_ptr< ArrayBase > a)
Definition: Field.cpp:72
Array
Definition: Array.h:14
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
template<class T >
void MeshMetricsCalculator< T >::computeGridInterpolationMatrices ( const StorageSite &  grids,
const StorageSite &  faces 
)

Definition at line 2442 of file MeshMetricsCalculator_impl.h.

2443 {
2444  const int numMeshes = _meshes.size();
2445  for (int n=0; n<numMeshes; n++)
2446  {
2447  const Mesh& mesh = *_meshes[n];
2448  computeGridInterpolationMatrices(mesh,g, f);
2449  }
2450 }
MeshMetricsCalculator::computeGridInterpolationMatrices
void computeGridInterpolationMatrices(const StorageSite &grids, const StorageSite &faces)
Definition: MeshMetricsCalculator_impl.h:2442
Mesh
Definition: Mesh.h:49
Model::_meshes
const MeshList _meshes
Definition: Model.h:29
template<class T >
void MeshMetricsCalculator< T >::computeGridInterpolationMatrices ( const Mesh &  mesh,
const StorageSite &  grids,
const StorageSite &  faces 
)
private

Definition at line 2136 of file MeshMetricsCalculator_impl.h.

References dot(), fabs(), and Mesh::getConnectivity().

2139 {
2140 
2141  typedef CRMatrixTranspose<T,T,T> IMatrix;
2142 
2143  const CRConnectivity& faceToGrids = mesh.getConnectivity(faces, grids );
2144 
2145 
2146  shared_ptr<IMatrix> gridToFaces (new IMatrix(faceToGrids));
2147 
2148 #if 0
2149  /* distance weighted */
2150  for(int n=0; n<nFaces; n++)
2151  {
2152 
2153  T wtSum(0);
2154  int nnb(0);
2155 
2156  for(int nc = FGRow[n]; nc < FGRow[n+1]; nc ++)
2157  {
2158 
2159  const int c = FGCol[nc ];
2160  VectorT3 dr(xGrids[c]-xFaces[n]);
2161  T wt = 1.0/dot(dr,dr);
2162  gridToFaceCoeff[nc ] = wt;
2163  wtSum += wt;
2164  nnb++;
2165  }
2166 
2167 
2168 
2169  for(int nc=FGRow[n]; nc<FGRow[n+1]; nc++)
2170  {
2171  gridToFaceCoeff[nc] /= wtSum;
2172  }
2173 
2174  }
2175 #endif
2176 
2177 
2178 #if 0
2179  /* linear interpolation */
2180  // v = a + b*x +c*y
2181  //solve a, b, c
2182  T Q[3][3];
2183  T Qinv[3][3];
2184 
2185 
2186  for(int n=0; n<nFaces; n++)
2187  {
2188  T wt(0);
2189  int nnb(0);
2190 
2191  for(int i=0; i<3; i++){
2192  for(int j=0; j<3; j++){
2193  Q[i][j]=Qinv[i][j]=0;
2194  }
2195  }
2196  // int size = FGRow[n+1]-FGRow[n];
2197 
2198  // if(size == 3){
2199  for(int nc=FGRow[n]; nc<FGRow[n+1]; nc++)
2200  {
2201  const int c = FGCol[nc];
2202  VectorT3 dr(xGrids[c]-xFaces[n]);
2203  //VectorT3 dr(xGrids[c]);
2204  Q[nnb][0]=1.0;
2205  Q[nnb][1]=dr[0];
2206  Q[nnb][2]=dr[1];
2207  nnb++;
2208  }
2209 
2210  matrix<T> matrix;
2211  matrix.Inverse3x3(Q, Qinv);
2212  nnb = 0;
2213  for(int nc=FGRow[n]; nc<FGRow[n+1]; nc++)
2214  {
2215  //wt = Qinv[0][nnb]+xFaces[n][0]*Qinv[1][nnb]+xFaces[n][1]*Qinv[2][nnb];
2216  wt = Qinv[0][nnb];
2217  gridToFaceCoeff[nc] = wt;
2218  nnb++;
2219 
2220  }
2221  //}
2222  }
2223 
2224 
2225 
2226 #endif
2227 
2228 
2229 
2230 #if 0
2231  //bi-linear average
2232  for(int n=0; n<nFaces; n++)
2233  {
2234  const int ncSize = FGRow[n+1] - FGRow[n];
2235  T dX[ncSize];
2236  T dY[ncSize];
2237  T dZ[ncSize];
2238  T wt[ncSize];
2239  int count = 0;
2240  for(int nc=FGRow[n]; nc<FGRow[n+1]; nc++)
2241  {
2242  const int c = FGCol[nc];
2243  VectorT3 dr(xGrids[c]-xFaces[n]);
2244  dX[count] = fabs(dr[0]);
2245  dY[count] = fabs(dr[1]);
2246  dZ[count] = fabs(dr[2]);
2247  count++;
2248  }
2249  const T u = dX[0]/(dX[0]+dX[2]);
2250  const T v = dY[0]/(dY[0]+dY[1]);
2251  wt[3] = u*v;
2252  wt[2] = u*(1.-v);
2253  wt[1] = (1.-u)*v;
2254  wt[0] = (1.-u)*(1.-v);
2255  count = 0;
2256  for(int nc=FGRow[n]; nc<FGRow[n+1]; nc++)
2257  {
2258  gridToFaceCoeff[nc] = wt[count];
2259  count++;
2260  }
2261 
2262  }
2263 
2264 
2265 #endif
2266 
2267  GeomFields::SSPair key1(&faces,&grids);
2268  this->_geomFields._interpolationMatrices[key1] = gridToFaces;
2269 
2270 }
Mesh::getConnectivity
const CRConnectivity & getConnectivity(const StorageSite &from, const StorageSite &to) const
Definition: Mesh.cpp:416
GeomFields::_interpolationMatrices
map< SSPair, shared_ptr< Matrix > > _interpolationMatrices
Definition: GeomFields.h:50
fabs
Tangent fabs(const Tangent &a)
Definition: Tangent.h:312
dot
T dot(const Vector< T, 3 > &a, const Vector< T, 3 > &b)
Definition: Vector.h:253
GeomFields::SSPair
pair< const StorageSite *, const StorageSite * > SSPair
Definition: GeomFields.h:48
template<class T >
void MeshMetricsCalculator< T >::computeIBandSolidInterpolationMatrices ( const StorageSite &  particles)

Definition at line 2430 of file MeshMetricsCalculator_impl.h.

2431 {
2432  const int numMeshes = _meshes.size();
2433  for (int n=0; n<numMeshes; n++)
2434  {
2435  const Mesh& mesh = *_meshes[n];
2436  computeIBandSolidInterpolationMatrices(mesh,p);
2437  }
2438 }
Mesh
Definition: Mesh.h:49
Model::_meshes
const MeshList _meshes
Definition: Model.h:29
MeshMetricsCalculator::computeIBandSolidInterpolationMatrices
void computeIBandSolidInterpolationMatrices(const StorageSite &particles)
Definition: MeshMetricsCalculator_impl.h:2430
template<class T >
void MeshMetricsCalculator< T >::computeIBandSolidInterpolationMatrices ( const Mesh &  mesh,
const StorageSite &  particles 
)
private

Definition at line 2276 of file MeshMetricsCalculator_impl.h.

References Mesh::getCells(), and Mesh::getConnectivity().

2278 {
2279  typedef CRMatrixTranspose<T,T,T> IMatrix;
2280  const StorageSite& cells = mesh.getCells();
2281  const CRConnectivity& cellToParticles
2282  = mesh.getConnectivity(cells,mpmParticles);
2283 
2284  shared_ptr<IMatrix> particlesToCell(new IMatrix(cellToParticles));
2285 
2286  /*
2287  for (int c=0; c<nCells; c++){
2288  const int nP = CPRow[c+1] - CPRow[c];
2289  if (nP == 0){
2290  //fluid cell. do nothing
2291  }
2292 
2293  if (nP == 1){
2294  //only one particle in this cell
2295  particlesToCellCoeff[ CPRow[c] ] = 1.0;
2296  //cout<<"only one particle in cell "<<c<<endl;
2297  }
2298 
2299  if (nP == 2){
2300  for(int nc=CPRow[c]; nc<CPRow[c+1]; nc++){
2301  particlesToCellCoeff[nc] = 0.5;
2302  }
2303  }
2304 
2305  if (nP == 3){
2306  for(int nc=CPRow[c]; nc<CPRow[c+1]; nc++){
2307  particlesToCellCoeff[nc] = 1./3.;
2308  }
2309  }
2310 
2311  if (nP >= 4){
2312 
2313  T wt(0);
2314  int nnb(0);
2315 
2316  T Q[4][4];
2317  T Qinv[4][4];
2318 
2319  for(int i=0; i<4; i++){
2320  for(int j=0; j<4; j++){
2321  Q[i][j]=Qinv[i][j]=0;
2322  }
2323  }
2324 
2325  for(int nc=CPRow[c]; nc<CPRow[c+1]; nc++)
2326  {
2327  const int p = CPCol[nc];
2328  VectorT3 dr(xParticles[p]-xCells[c]);
2329  Q[0][0] += 1.0;
2330  Q[0][1] += dr[0];
2331  Q[0][2] += dr[1];
2332  Q[0][3] += dr[2];
2333  Q[1][1] += dr[0]*dr[0];
2334  Q[1][2] += dr[0]*dr[1];
2335  Q[1][3] += dr[0]*dr[2];
2336  Q[2][2] += dr[1]*dr[1];
2337  Q[2][3] += dr[1]*dr[2];
2338  Q[3][3] += dr[2]*dr[2];
2339  nnb++;
2340  }
2341 
2342  for(int i=0; i<4; i++){
2343  for(int j=0; j<i; j++){
2344  Q[i][j]=Q[j][i];
2345  }
2346  }
2347  matrix<T> matrix;
2348  matrix.Inverse4x4(Q, Qinv);
2349 
2350  for (int nc=CPRow[c]; nc<CPRow[c+1]; nc++)
2351  {
2352  const int p = CPCol[nc];
2353  VectorT3 dr(xParticles[p]-xCells[c]);
2354  wt = Qinv[0][0];
2355  for (int i=1; i<=3; i++){
2356  wt += Qinv[0][i]*dr[i-1];
2357  }
2358  particlesToCellCoeff[nc] = wt;
2359  }
2360  }
2361  }
2362 */
2363  GeomFields::SSPair key2(&cells,&mpmParticles);
2364  this->_geomFields._interpolationMatrices[key2] = particlesToCell;
2365 
2366 }
Mesh::getConnectivity
const CRConnectivity & getConnectivity(const StorageSite &from, const StorageSite &to) const
Definition: Mesh.cpp:416
GeomFields::_interpolationMatrices
map< SSPair, shared_ptr< Matrix > > _interpolationMatrices
Definition: GeomFields.h:50
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
GeomFields::SSPair
pair< const StorageSite *, const StorageSite * > SSPair
Definition: GeomFields.h:48
template<class T >
void MeshMetricsCalculator< T >::computeIBInterpolationMatrices ( const StorageSite &  particles,
const int  option = 0 
)

Definition at line 2371 of file MeshMetricsCalculator_impl.h.

2372 {
2373  const int numMeshes = _meshes.size();
2374  for (int n=0; n<numMeshes; n++)
2375  {
2376  const Mesh& mesh = *_meshes[n];
2377  computeIBInterpolationMatrices(mesh,p, option);
2378  }
2379 }
Mesh
Definition: Mesh.h:49
Model::_meshes
const MeshList _meshes
Definition: Model.h:29
MeshMetricsCalculator::computeIBInterpolationMatrices
void computeIBInterpolationMatrices(const StorageSite &particles, const int option=0)
Definition: MeshMetricsCalculator_impl.h:2371
template<class T >
void MeshMetricsCalculator< T >::computeIBInterpolationMatrices ( const Mesh &  mesh,
const StorageSite &  particles,
const int  option 
)
private

Definition at line 466 of file MeshMetricsCalculator_impl.h.

References determinant(), dot(), Mesh::getAllFaceCells(), Mesh::getCells(), CRConnectivity::getCol(), Mesh::getConnectivity(), StorageSite::getCount(), Mesh::getDimension(), Mesh::getFaces(), Mesh::getIBFaceList(), Mesh::getIBFaces(), Mesh::getLocalToGlobal(), CRConnectivity::getRow(), inverse(), Mesh::isShell(), max(), and min().

469 {
470  if (mesh.isShell() || mesh.getIBFaces().getCount()==0)
471  return;
472 
473  typedef CRMatrixTranspose<T,T,T> IMatrix;
474  typedef map<int,double> IntDoubleMap;
475 
476  const StorageSite& ibFaces = mesh.getIBFaces();
477  const StorageSite& cells = mesh.getCells();
478  const StorageSite& faces = mesh.getFaces();
479  const CRConnectivity& ibFaceToCells = mesh.getConnectivity(ibFaces,cells);
480  const CRConnectivity& ibFaceToParticles
481  = mesh.getConnectivity(ibFaces,mpmParticles);
482 
483  const VectorT3Array& xFaces =
484  dynamic_cast<const VectorT3Array&>(_coordField[faces]);
485 
486  const VectorT3Array& xCells =
487  dynamic_cast<const VectorT3Array&>(_coordField[cells]);
488 
489  const VectorT3Array& xParticles =
490  dynamic_cast<const VectorT3Array&>(_coordField[mpmParticles]);
491 
492  const Array<int>& ibFCRow = ibFaceToCells.getRow();
493  const Array<int>& ibFCCol = ibFaceToCells.getCol();
494 
495  const Array<int>& ibFPRow = ibFaceToParticles.getRow();
496  const Array<int>& ibFPCol = ibFaceToParticles.getCol();
497 
498  const int nIBFaces = ibFaces.getCount();
499 
500  const Array<int>& ibFaceIndices = mesh.getIBFaceList();
501 
502  shared_ptr<IMatrix> cellToIB(new IMatrix(ibFaceToCells));
503  shared_ptr<IMatrix> particlesToIB(new IMatrix(ibFaceToParticles));
504 
505  Array<T>& cellToIBCoeff = cellToIB->getCoeff();
506  Array<T>& particlesToIBCoeff = particlesToIB->getCoeff();
507 
508  const bool is2D = mesh.getDimension() == 2;
509  const bool is3D = mesh.getDimension() == 3;
510 
511  /***********************************************************************/
512  // default is to use linear least square interpolation
513  // if the matrix determinant is too small
514  // then switch to distance weighted interpolation
515  /***********************************************************************/
516 
517  /**********linear least square interpolation*********/
518  // X=x-xf Y=y-yf Z=Z-zf
519  //matrix M=[1 X1 Y1 Z1]
520  // [1 X2 Y2 Z2]
521  // ...........
522  // [1 Xn Yn Zn]
523  //coefficient matrix A=[a, b, c, d]T
524  //velocity element v = M * A
525  //linear relation v = a + b*X + c*Y + d*Z
526  //to make least square
527  //matrix A = (M(T)*M)^(-1)*M(T)*v
528  //so, velocity at face is vface = a + b*Xf + c*Yf + d*Zf = a
529  //which is the first row of matrix A
530  //so, the weight function should be the first row of (M(T)*M)^(-1)*M(T)
531  //note Q = M(T)*M and Qinv = Q^(-1)
532  //the following code is to calculate it
533  //insteading of doing full matrix operation, only nessesary operation on entries are performed
534  //when dealing with coodinates with small numbers, like in micrometers
535  //scaling the coodinates does not change the coefficients but improve the matrix quality
536 
537  // FILE * fp = fopen("/home/lin/work/app-memosa/src/fvm/verification/Structure_Electrostatics_Interaction/2D_beam/test/coeff.dat", "w");
538 
539 #if 1
540 
541 #if 0
542  ofstream debugFileFluid;
543  ofstream debugFileSolid;
544  stringstream ss(stringstream::in | stringstream::out);
545  ss << MPI::COMM_WORLD.Get_rank();
546  string fname1 = "IBinterpolationFluid_proc" + ss.str() + ".dat";
547  string fname2 = "IBinterpolationSolid_proc" + ss.str() + ".dat";
548  debugFileFluid.open( fname1.c_str() );
549  debugFileSolid.open( fname2.c_str() );
550  ss.str("");
551  const Array<int>& localToGlobal = mesh.getLocalToGlobal();
552  const CRConnectivity& faceCells = mesh.getAllFaceCells();
553 #endif
554 
555  for(int n=0; n<nIBFaces; n++)
556  {
557  const int f = ibFaceIndices[n];
558  T wt(0);
559  T wtSum(0);
560  T det(0);
561  int nnb(0);
562  T scale(1.0e6);
563 
564  SquareMatrix<T,4> Q(NumTypeTraits<T>::getZero());
565  SquareMatrix<T,4> Qinv(NumTypeTraits<T>::getZero());
566  SquareMatrix<T,3> QQ(NumTypeTraits<T>::getZero());
567  SquareMatrix<T,3> QQinv(NumTypeTraits<T>::getZero());
568 
569  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++){
570  const int c = ibFCCol[nc];
571  VectorT3 dr((xCells[c]-xFaces[f])*scale);
572  //if (f ==200){
573  // cout << f <<" " << c <<endl;
574  // cout << xCells[c] << endl;
575  // cout << xFaces[f] << endl;
576  // }
577  Q(0,0) += 1.0;
578  Q(0,1) += dr[0];
579  Q(0,2) += dr[1];
580  Q(0,3) += dr[2];
581  Q(1,1) += dr[0]*dr[0];
582  Q(1,2) += dr[0]*dr[1];
583  Q(1,3) += dr[0]*dr[2];
584  Q(2,2) += dr[1]*dr[1];
585  Q(2,3) += dr[1]*dr[2];
586  Q(3,3) += dr[2]*dr[2];
587  nnb++;
588  }
589 
590  for(int np=ibFPRow[n]; np<ibFPRow[n+1]; np++)
591  {
592  const int p = ibFPCol[np];
593  VectorT3 dr((xParticles[p]-xFaces[f])*scale);
594  Q(0,0) += 1.0;
595  Q(0,1) += dr[0];
596  Q(0,2) += dr[1];
597  Q(0,3) += dr[2];
598  Q(1,1) += dr[0]*dr[0];
599  Q(1,2) += dr[0]*dr[1];
600  Q(1,3) += dr[0]*dr[2];
601  Q(2,2) += dr[1]*dr[1];
602  Q(2,3) += dr[1]*dr[2];
603  Q(3,3) += dr[2]*dr[2];
604  nnb++;
605  }
606 
607  //if (nnb < 4)
608  //throw CException("not enough cell or particle neighbors for ib face to interpolate!");
609 
610  //symetric matrix
611  for(int i=0; i<4; i++){
612  for(int j=0; j<i; j++){
613  Q(i,j) = Q(j,i);
614  }
615  }
616 
617  // calculate the determinant of the matrix Q
618  // if 3D mesh, then det(Q)
619  // if 2D mesh, then det(QQ) where QQ is the 3x3 subset of Q
620 
621 
622  if (is2D) {
623  for(int i=0; i<3; i++){
624  for(int j=0; j<3; j++){
625  QQ(i,j)=Q(i,j);
626  }
627  }
628  det = determinant(QQ);
629  }
630 
631  if (is3D) {
632  det = determinant(Q, 4);
633  }
634 
635  // linear least square interpolation if the matrix is not singular
636  //if (nnb >=10){
637  //if (fabs(det) > 1.0 ){
638  if (option == 0) {
639  if(is2D){
640  QQinv = inverse(QQ);
641  for(int i=0; i<3; i++){
642  for(int j=0; j<3; j++){
643  Qinv(i,j)=QQinv(i,j);
644  }
645  }
646  }
647 
648  if(is3D){
649  Qinv = inverse(Q, 4);
650  }
651 
652  //calculate Qinv*M(T) get the first row element, put in coeffMatrix
653  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++) {
654  const int c = ibFCCol[nc];
655  VectorT3 dr((xCells[c]-xFaces[f])*scale);
656  wt = Qinv(0,0);
657  for (int i=1; i<=3; i++){
658  wt += Qinv(0,i)*dr[i-1];
659  }
660  cellToIBCoeff[nc] = wt;
661  wtSum += wt;
662  }
663  for(int np=ibFPRow[n]; np<ibFPRow[n+1]; np++) {
664  const int p = ibFPCol[np];
665  VectorT3 dr((xParticles[p]-xFaces[f])*scale);
666  wt = Qinv(0,0);
667  for (int i=1; i<=3; i++){
668  wt += Qinv(0,i)*dr[i-1];
669  }
670  particlesToIBCoeff[np] = wt;
671  wtSum += wt;
672  }
673 
674 
675  if (wtSum > 1.01 || wtSum < 0.99)
676  cout << "face " << n <<" has wrong wtsum " << wtSum << endl;
677  /*
678  cout << n << endl;
679  cout << "ibface " << xFaces[f][0] << " " << xFaces[f][1] << " " << xFaces[f][2] << " " << endl;
680  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++) {
681  const int c = ibFCCol[nc];
682  cout << "fluid cells " << xCells[c][0] << " " << xCells[c][1] << " " << xCells[c][2] << " " <<cellToIBCoeff[nc] << endl;
683  }
684  for(int np=ibFPRow[n]; np<ibFPRow[n+1]; np++) {
685  const int p = ibFPCol[np];
686  cout << "particles " << xParticles[p][0] << " " << xParticles[p][1] << " " << xParticles[p][2] << " " << particlesToIBCoeff[np] << endl;
687  }
688  */
689 
690 #if 0
691  const int cell0 = localToGlobal[ faceCells(f,0) ];
692  const int cell1 = localToGlobal[ faceCells(f,1) ];
693  debugFileFluid << "ibface = " << n << " " << xFaces[f][0] << " " << xFaces[f][1] << " " << xFaces[f][2] <<
694  " cell0 = " << std::min(cell0,cell1) << " cell1 = " << std::max(cell0,cell1) << endl;
695  map<int, double> cellToValue;
696  map<int, int> globalToLocal;
697 
698  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++){
699  const int localID = ibFCCol[nc];
700  const int c = localToGlobal[ibFCCol[nc]];
701  globalToLocal[c] = localID;
702  cellToValue[c] = cellToIBCoeff[nc];
703 
704  //debugFile << " glblcellID = " << c << ", cellToIBCoeff[" << n << "] = " << cellToIBCoeff[nc] << endl;
705  }
706  foreach( IntDoubleMap::value_type& pos, cellToValue){
707  const int c = pos.first;
708  const double value =pos.second;
709  debugFileFluid << " glblcellID = " << c << " localCellID = " << globalToLocal[c] << ", cellToIBCoeff[" << n << "] = " << value << endl;
710  }
711 
712  debugFileSolid << "ibface = " << n << " " << xFaces[f][0] << " " << xFaces[f][1] << " " << xFaces[f][2] << endl;
713  cellToValue.clear();
714  for(int nc=ibFPRow[n]; nc<ibFPRow[n+1]; nc++){
715  cellToValue[ibFPCol[nc]] = particlesToIBCoeff[nc];
716  }
717  foreach( IntDoubleMap::value_type& pos, cellToValue){
718  const int c = pos.first;
719  const double value =pos.second;
720  debugFileSolid << " GlobalSolidFaceID = " << c << ", solidToIBCoeff[" << n << "] = " << value << endl;
721  }
722 
723 
724 #endif
725 
726  }
727 
728  if (option == 1)
729  { //if matrix is singular, use distance weighted interpolation
730  //cout << "warning: IBM interpolation switched to distance weighted method for face " << f << endl;
731  //cout << xFaces[f][0] << " " << xFaces[f][1] << " " << xFaces[f][2] << " " << endl;
732  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++)
733  {
734  const int c = ibFCCol[nc];
735  VectorT3 dr(xCells[c]-xFaces[f]);
736  T wt = 1.0/dot(dr,dr);
737  cellToIBCoeff[nc] = wt;
738  wtSum += wt;
739  nnb++;
740  //cout << "fluid cells " << xCells[c][0] << " " << xCells[c][1] << " " << xCells[c][2] << " " << endl;
741  }
742  if (nnb == 0)
743  throw CException("no fluid cell for ib face");
744  for(int np=ibFPRow[n]; np<ibFPRow[n+1]; np++)
745  {
746  const int p = ibFPCol[np];
747  VectorT3 dr(xParticles[p]-xFaces[f]);
748  T wt = 1.0/dot(dr,dr);
749  particlesToIBCoeff[np] = wt;
750  wtSum += wt;
751  nnb++;
752  //cout << "particles " << xParticles[p][0] << " " << xParticles[p][1] << " " << xParticles[p][2] << " " << endl;
753  }
754 
755  if (nnb == 0)
756  throw CException("no solid neighbors for ib face");
757 
758  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++)
759  {
760  cellToIBCoeff[nc] /= wtSum;
761  }
762 
763  for(int np=ibFPRow[n]; np<ibFPRow[n+1]; np++)
764  {
765  particlesToIBCoeff[np] /= wtSum;
766  }
767 
768  }
769  }
770 #endif
771 
772 #if 0
773  debugFileFluid.close();
774  debugFileSolid.close();
775 #endif
776 
777 
778  //fclose(fp);
779 #if 0
780 
781  /**********second order least square interpolation*********/
782  // X=x-xf Y=y-yf Z=Z-zf
783  //matrix M=[1 X1 Y1 Z1 X1*X1 Y1*Y1 Z1*Z1 X1*Y1 Y1*Z1 X1*Z1]
784  // [1 X2 Y2 Z2 ...................................]
785  // ...........
786  // [1 Xn Yn Zn Xn*Xn Yn*Yn Zn*Zn Xn*Yn Yn*Zn Xn*Zn]
787  //coefficient matrix A=[a0, a1, a2, a3, a4, a5, a6, a7, a8, a9]T
788  //velocity element v = M * A
789  //linear relation v = a0 + a1*X + a2*Y + a3*Z + a4*X*X + a5*Y*Y + a6*Z*Z + a7*X*Y + a8*Y*Z + a9*X*Z
790  //to make least square
791  //matrix A = (M(T)*M)^(-1)*M(T)*v
792  //so, velocity at face is vface = a0
793  //which is the first row of matrix A
794  //so, the weight function should be the first row of (M(T)*M)^(-1)*M(T)
795  //note Q = M(T)*M and Qinv = Q^(-1)
796  //the following code is to calculate it
797  //insteading of doing full matrix operation, only nessesary operation on entries are performed
798 
799 
800 
801  for(int n=0; n<nIBFaces; n++)
802  {
803  const int f = ibFaceIndices[n];
804  T wt(0);
805  int nnb(0);
806  T scale(1.0e6);
807 
808  if (is2D){
809  const int size = 6;
810  SquareMatrix<T,size> Q(SquareMatrix<T,size>::zero());
811  SquareMatrix<T,size> Qinv(0);
812  //cout << xFaces[f][0] << " " << xFaces[f][1] << " " << xFaces[f][2] << " " << endl;
813  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++) {
814  const int c = ibFCCol[nc];
815  VectorT3 dr((xCells[c]-xFaces[f])*scale);
816  //cout << xCells[c][0] << " " << xCells[c][1] << " " << xCells[c][2] << " " << endl;
817 
818  Q(0,0) += 1.0;
819  Q(0,1) += dr[0];
820  Q(0,2) += dr[1];
821  Q(0,3) += dr[0]*dr[0];
822  Q(0,4) += dr[1]*dr[1];
823  Q(0,5) += dr[0]*dr[1];
824 
825  Q(1,1) += dr[0]*dr[0];
826  Q(1,2) += dr[0]*dr[1];
827  Q(1,3) += dr[0]*dr[0]*dr[0];
828  Q(1,4) += dr[0]*dr[1]*dr[1];
829  Q(1,5) += dr[0]*dr[0]*dr[1];
830 
831  Q(2,2) += dr[1]*dr[1];
832  Q(2,3) += dr[1]*dr[0]*dr[0];
833  Q(2,4) += dr[1]*dr[1]*dr[1];
834  Q(2,5) += dr[1]*dr[0]*dr[1];
835 
836  Q(3,3) += dr[0]*dr[0]*dr[0]*dr[0];
837  Q(3,4) += dr[0]*dr[0]*dr[1]*dr[1];
838  Q(3,5) += dr[0]*dr[0]*dr[0]*dr[1];
839 
840  Q(4,4) += dr[1]*dr[1]*dr[1]*dr[1];
841  Q(4,5) += dr[1]*dr[1]*dr[0]*dr[1];
842 
843  Q(5,5) += dr[0]*dr[1]*dr[0]*dr[1];
844 
845  nnb++;
846  }
847 
848  for(int np=ibFPRow[n]; np<ibFPRow[n+1]; np++) {
849  const int p = ibFPCol[np];
850  VectorT3 dr((xParticles[p]-xFaces[f])*scale);
851  //cout << xParticles[p][0] << " " << xParticles[p][1] << " " << xParticles[p][2] << " " << endl;
852 
853  Q(0,0) += 1.0;
854  Q(0,1) += dr[0];
855  Q(0,2) += dr[1];
856  Q(0,3) += dr[0]*dr[0];
857  Q(0,4) += dr[1]*dr[1];
858  Q(0,5) += dr[0]*dr[1];
859 
860  Q(1,1) += dr[0]*dr[0];
861  Q(1,2) += dr[0]*dr[1];
862  Q(1,3) += dr[0]*dr[0]*dr[0];
863  Q(1,4) += dr[0]*dr[1]*dr[1];
864  Q(1,5) += dr[0]*dr[0]*dr[1];
865 
866  Q(2,2) += dr[1]*dr[1];
867  Q(2,3) += dr[1]*dr[0]*dr[0];
868  Q(2,4) += dr[1]*dr[1]*dr[1];
869  Q(2,5) += dr[1]*dr[0]*dr[1];
870 
871  Q(3,3) += dr[0]*dr[0]*dr[0]*dr[0];
872  Q(3,4) += dr[0]*dr[0]*dr[1]*dr[1];
873  Q(3,5) += dr[0]*dr[0]*dr[0]*dr[1];
874 
875  Q(4,4) += dr[1]*dr[1]*dr[1]*dr[1];
876  Q(4,5) += dr[1]*dr[1]*dr[0]*dr[1];
877 
878  Q(5,5) += dr[0]*dr[1]*dr[0]*dr[1];
879 
880  nnb++;
881  }
882 
883  if (nnb < size)
884  throw CException("not enough cell or particle neighbors for ib face to interpolate!");
885 
886  //symetric matrix
887  for(int i=0; i<size; i++){
888  for(int j=0; j<i; j++){
889  Q(i,j)=Q(j,i);
890  }
891  }
892 
893  //calculate the inverse of Q(6x6)
894  Qinv = inverse(Q, size);
895 
896  //calculate Qinv*M(T) get the first row element, put in coeffMatrix
897  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++) {
898  const int c = ibFCCol[nc];
899  VectorT3 dr((xCells[c]-xFaces[f])*scale);
900  wt = Qinv(0,0);
901  wt += Qinv(0,1)*dr[0];
902  wt += Qinv(0,2)*dr[1];
903  wt += Qinv(0,3)*dr[0]*dr[0];
904  wt += Qinv(0,4)*dr[1]*dr[1];
905  wt += Qinv(0,5)*dr[0]*dr[1];
906  cellToIBCoeff[nc] = wt;
907  //cout<<n<<" cells "<<nc<<" "<<cellToIBCoeff[nc]<<endl;
908  }
909  for(int np=ibFPRow[n]; np<ibFPRow[n+1]; np++) {
910  const int p = ibFPCol[np];
911  VectorT3 dr((xParticles[p]-xFaces[f])*scale);
912  wt = Qinv(0,0);
913  wt += Qinv(0,1)*dr[0];
914  wt += Qinv(0,2)*dr[1];
915  wt += Qinv(0,3)*dr[0]*dr[0];
916  wt += Qinv(0,4)*dr[1]*dr[1];
917  wt += Qinv(0,5)*dr[0]*dr[1];
918  particlesToIBCoeff[np] = wt;
919  // cout<<n<<" particles "<<np<<" "<<particlesToIBCoeff[np]<<endl;
920  }
921  }
922 
923  if (is3D) {
924  const int size = 10;
925  SquareMatrix<T,size> Q(0);
926  SquareMatrix<T,size> Qinv(0);
927 
928  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++)
929  {
930  const int c = ibFCCol[nc];
931  VectorT3 dr((xCells[c]-xFaces[f])*scale);
932  Q(0,0) += 1.0;
933  Q(0,1) += dr[0];
934  Q(0,2) += dr[1];
935  Q(0,3) += dr[2];
936  Q(0,4) += dr[0]*dr[0];
937  Q(0,5) += dr[1]*dr[1];
938  Q(0,6) += dr[2]*dr[2];
939  Q(0,7) += dr[0]*dr[1];
940  Q(0,8) += dr[1]*dr[2];
941  Q(0,9) += dr[0]*dr[2];
942 
943  Q(1,1) += dr[0]*dr[0];
944  Q(1,2) += dr[0]*dr[1];
945  Q(1,3) += dr[0]*dr[2];
946  Q(1,4) += dr[0]*dr[0]*dr[0];
947  Q(1,5) += dr[0]*dr[1]*dr[1];
948  Q(1,6) += dr[0]*dr[2]*dr[2];
949  Q(1,7) += dr[0]*dr[0]*dr[1];
950  Q(1,8) += dr[0]*dr[1]*dr[2];
951  Q(1,9) += dr[0]*dr[0]*dr[2];
952 
953  Q(2,2) += dr[1]*dr[1];
954  Q(2,3) += dr[1]*dr[2];
955  Q(2,4) += dr[1]*dr[0]*dr[0];
956  Q(2,5) += dr[1]*dr[1]*dr[1];
957  Q(2,6) += dr[1]*dr[2]*dr[2];
958  Q(2,7) += dr[1]*dr[0]*dr[1];
959  Q(2,8) += dr[1]*dr[1]*dr[2];
960  Q(2,9) += dr[1]*dr[0]*dr[2];
961 
962  Q(3,3) += dr[2]*dr[2];
963  Q(3,4) += dr[2]*dr[0]*dr[0];
964  Q(3,5) += dr[2]*dr[1]*dr[1];
965  Q(3,6) += dr[2]*dr[2]*dr[2];
966  Q(3,7) += dr[2]*dr[0]*dr[1];
967  Q(3,8) += dr[2]*dr[1]*dr[2];
968  Q(3,8) += dr[2]*dr[0]*dr[2];
969 
970  Q(4,4) += dr[0]*dr[0]*dr[0]*dr[0];
971  Q(4,5) += dr[0]*dr[0]*dr[1]*dr[1];
972  Q(4,6) += dr[0]*dr[0]*dr[2]*dr[2];
973  Q(4,7) += dr[0]*dr[0]*dr[0]*dr[1];
974  Q(4,8) += dr[0]*dr[0]*dr[1]*dr[2];
975  Q(4,9) += dr[0]*dr[0]*dr[0]*dr[2];
976 
977  Q(5,5) += dr[1]*dr[1]*dr[1]*dr[1];
978  Q(5,6) += dr[1]*dr[1]*dr[2]*dr[2];
979  Q(5,7) += dr[1]*dr[1]*dr[0]*dr[1];
980  Q(5,8) += dr[1]*dr[1]*dr[1]*dr[2];
981  Q(5,9) += dr[1]*dr[1]*dr[0]*dr[2];
982 
983  Q(6,6) += dr[2]*dr[2]*dr[2]*dr[2];
984  Q(6,7) += dr[2]*dr[2]*dr[0]*dr[1];
985  Q(6,8) += dr[2]*dr[2]*dr[1]*dr[2];
986  Q(6,9) += dr[2]*dr[2]*dr[0]*dr[2];
987 
988  Q(7,7) += dr[0]*dr[1]*dr[0]*dr[1];
989  Q(7,8) += dr[0]*dr[1]*dr[1]*dr[2];
990  Q(7,9) += dr[0]*dr[1]*dr[0]*dr[2];
991 
992  Q(8,8) += dr[1]*dr[2]*dr[1]*dr[2];
993  Q(8,9) += dr[1]*dr[2]*dr[0]*dr[2];
994 
995  Q(9,9) += dr[0]*dr[2]*dr[0]*dr[2];
996 
997  nnb++;
998  }
999  for(int np=ibFPRow[n]; np<ibFPRow[n+1]; np++)
1000  {
1001  const int p = ibFPCol[np];
1002  VectorT3 dr((xParticles[p]-xFaces[f])*scale);
1003  Q(0,0) += 1.0;
1004  Q(0,1) += dr[0];
1005  Q(0,2) += dr[1];
1006  Q(0,3) += dr[2];
1007  Q(0,4) += dr[0]*dr[0];
1008  Q(0,5) += dr[1]*dr[1];
1009  Q(0,6) += dr[2]*dr[2];
1010  Q(0,7) += dr[0]*dr[1];
1011  Q(0,8) += dr[1]*dr[2];
1012  Q(0,9) += dr[0]*dr[2];
1013 
1014  Q(1,1) += dr[0]*dr[0];
1015  Q(1,2) += dr[0]*dr[1];
1016  Q(1,3) += dr[0]*dr[2];
1017  Q(1,4) += dr[0]*dr[0]*dr[0];
1018  Q(1,5) += dr[0]*dr[1]*dr[1];
1019  Q(1,6) += dr[0]*dr[2]*dr[2];
1020  Q(1,7) += dr[0]*dr[0]*dr[1];
1021  Q(1,8) += dr[0]*dr[1]*dr[2];
1022  Q(1,9) += dr[0]*dr[0]*dr[2];
1023 
1024  Q(2,2) += dr[1]*dr[1];
1025  Q(2,3) += dr[1]*dr[2];
1026  Q(2,4) += dr[1]*dr[0]*dr[0];
1027  Q(2,5) += dr[1]*dr[1]*dr[1];
1028  Q(2,6) += dr[1]*dr[2]*dr[2];
1029  Q(2,7) += dr[1]*dr[0]*dr[1];
1030  Q(2,8) += dr[1]*dr[1]*dr[2];
1031  Q(2,9) += dr[1]*dr[0]*dr[2];
1032 
1033  Q(3,3) += dr[2]*dr[2];
1034  Q(3,4) += dr[2]*dr[0]*dr[0];
1035  Q(3,5) += dr[2]*dr[1]*dr[1];
1036  Q(3,6) += dr[2]*dr[2]*dr[2];
1037  Q(3,7) += dr[2]*dr[0]*dr[1];
1038  Q(3,8) += dr[2]*dr[1]*dr[2];
1039  Q(3,8) += dr[2]*dr[0]*dr[2];
1040 
1041  Q(4,4) += dr[0]*dr[0]*dr[0]*dr[0];
1042  Q(4,5) += dr[0]*dr[0]*dr[1]*dr[1];
1043  Q(4,6) += dr[0]*dr[0]*dr[2]*dr[2];
1044  Q(4,7) += dr[0]*dr[0]*dr[0]*dr[1];
1045  Q(4,8) += dr[0]*dr[0]*dr[1]*dr[2];
1046  Q(4,9) += dr[0]*dr[0]*dr[0]*dr[2];
1047 
1048  Q(5,5) += dr[1]*dr[1]*dr[1]*dr[1];
1049  Q(5,6) += dr[1]*dr[1]*dr[2]*dr[2];
1050  Q(5,7) += dr[1]*dr[1]*dr[0]*dr[1];
1051  Q(5,8) += dr[1]*dr[1]*dr[1]*dr[2];
1052  Q(5,9) += dr[1]*dr[1]*dr[0]*dr[2];
1053 
1054  Q(6,6) += dr[2]*dr[2]*dr[2]*dr[2];
1055  Q(6,7) += dr[2]*dr[2]*dr[0]*dr[1];
1056  Q(6,8) += dr[2]*dr[2]*dr[1]*dr[2];
1057  Q(6,9) += dr[2]*dr[2]*dr[0]*dr[2];
1058 
1059  Q(7,7) += dr[0]*dr[1]*dr[0]*dr[1];
1060  Q(7,8) += dr[0]*dr[1]*dr[1]*dr[2];
1061  Q(7,9) += dr[0]*dr[1]*dr[0]*dr[2];
1062 
1063  Q(8,8) += dr[1]*dr[2]*dr[1]*dr[2];
1064  Q(8,9) += dr[1]*dr[2]*dr[0]*dr[2];
1065 
1066  Q(9,9) += dr[0]*dr[2]*dr[0]*dr[2];
1067 
1068  nnb++;
1069  }
1070 
1071  if (nnb < size)
1072  throw CException("not enough cell or particle neighbors for ib face to interpolate!");
1073 
1074  //symetric matrix
1075  for(int i=0; i<size; i++){
1076  for(int j=0; j<i; j++){
1077  Q(i,j)=Q(j,i);
1078  }
1079  }
1080 
1081  //calculate the inverse of Q(10x10)
1082  Qinv = inverse(Q, size);
1083 
1084 
1085  //calculate Qinv*M(T) get the first row element, put in coeffMatrix
1086  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++)
1087  {
1088  const int c = ibFCCol[nc];
1089  VectorT3 dr((xCells[c]-xFaces[f])*scale);
1090  wt = Qinv(0,0);
1091  wt += Qinv(0,1)*dr[0];
1092  wt += Qinv(0,2)*dr[1];
1093  wt += Qinv(0,3)*dr[2];
1094  wt += Qinv(0,4)*dr[0]*dr[0];
1095  wt += Qinv(0,5)*dr[1]*dr[1];
1096  wt += Qinv(0,6)*dr[2]*dr[2];
1097  wt += Qinv(0,7)*dr[0]*dr[1];
1098  wt += Qinv(0,8)*dr[1]*dr[2];
1099  wt += Qinv(0,9)*dr[0]*dr[2];
1100  cellToIBCoeff[nc] = wt;
1101  //cout<<n<<" cells "<<nc<<" "<<cellToIBCoeff[nc]<<endl;
1102  }
1103  for(int np=ibFPRow[n]; np<ibFPRow[n+1]; np++)
1104  {
1105  const int p = ibFPCol[np];
1106  VectorT3 dr((xParticles[p]-xFaces[f])*scale);
1107  wt = Qinv(0,0);
1108  wt += Qinv(0,1)*dr[0];
1109  wt += Qinv(0,2)*dr[1];
1110  wt += Qinv(0,3)*dr[2];
1111  wt += Qinv(0,4)*dr[0]*dr[0];
1112  wt += Qinv(0,5)*dr[1]*dr[1];
1113  wt += Qinv(0,6)*dr[2]*dr[2];
1114  wt += Qinv(0,7)*dr[0]*dr[1];
1115  wt += Qinv(0,8)*dr[1]*dr[2];
1116  wt += Qinv(0,9)*dr[0]*dr[2];
1117  particlesToIBCoeff[np] = wt;
1118  // cout<<n<<" particles "<<np<<" "<<particlesToIBCoeff[np]<<endl;
1119  }
1120  }
1121  }
1122 #endif
1123  GeomFields::SSPair key1(&ibFaces,&cells);
1124  this->_geomFields._interpolationMatrices[key1] = cellToIB;
1125 
1126  GeomFields::SSPair key2(&ibFaces,&mpmParticles);
1127  this->_geomFields._interpolationMatrices[key2] = particlesToIB;
1128 
1129 }
CRConnectivity::getCol
const Array< int > & getCol() const
Definition: CRConnectivity.h:183
CRConnectivity::getRow
const Array< int > & getRow() const
Definition: CRConnectivity.h:182
Mesh::getIBFaces
const StorageSite & getIBFaces() const
Definition: Mesh.h:111
Mesh::getConnectivity
const CRConnectivity & getConnectivity(const StorageSite &from, const StorageSite &to) const
Definition: Mesh.cpp:416
Mesh::getIBFaceList
const Array< int > & getIBFaceList() const
Definition: Mesh.cpp:686
max
double max(double x, double y)
Definition: Octree.cpp:18
Mesh::getAllFaceCells
const CRConnectivity & getAllFaceCells() const
Definition: Mesh.cpp:378
Mesh::getLocalToGlobal
Array< int > & getLocalToGlobal()
Definition: Mesh.h:240
GeomFields::_interpolationMatrices
map< SSPair, shared_ptr< Matrix > > _interpolationMatrices
Definition: GeomFields.h:50
Mesh::getFaces
const StorageSite & getFaces() const
Definition: Mesh.h:108
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
inverse
SquareMatrix< T, 2 > inverse(SquareMatrix< T, 2 > &a)
Definition: MatrixOperation.h:118
Mesh::isShell
bool isShell() const
Definition: Mesh.h:323
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
min
double min(double x, double y)
Definition: Octree.cpp:23
Mesh::getDimension
int getDimension() const
Definition: Mesh.h:105
dot
T dot(const Vector< T, 3 > &a, const Vector< T, 3 > &b)
Definition: Vector.h:253
GeomFields::SSPair
pair< const StorageSite *, const StorageSite * > SSPair
Definition: GeomFields.h:48
determinant
T determinant(SquareMatrix< T, 2 > &a)
Definition: MatrixOperation.h:67
template<class T >
void MeshMetricsCalculator< T >::computeIBInterpolationMatricesCells ( )

Definition at line 2382 of file MeshMetricsCalculator_impl.h.

2383 {
2384  const int numMeshes = _meshes.size();
2385  for (int n=0; n<numMeshes; n++)
2386  {
2387  const Mesh& mesh = *_meshes[n];
2388  computeIBInterpolationMatricesCells(mesh);
2389  }
2390 }
Mesh
Definition: Mesh.h:49
Model::_meshes
const MeshList _meshes
Definition: Model.h:29
template<class T >
void MeshMetricsCalculator< T >::computeIBInterpolationMatricesCells ( const Mesh &  mesh)
private

Definition at line 1134 of file MeshMetricsCalculator_impl.h.

References determinant(), dot(), fabs(), Mesh::getAllFaceCells(), Mesh::getCells(), CRConnectivity::getCol(), Mesh::getConnectivity(), StorageSite::getCount(), Mesh::getDimension(), Mesh::getFaces(), Mesh::getIBFaceList(), Mesh::getIBFaces(), Mesh::getLocalToGlobal(), CRConnectivity::getRow(), inverse(), and Mesh::isShell().

1135 {
1136  if (mesh.isShell() || mesh.getIBFaces().getCount()==0)
1137  return;
1138 
1139  typedef CRMatrixTranspose<T,T,T> IMatrix;
1140  typedef map<int,double> IntDoubleMap;
1141 
1142  const StorageSite& ibFaces = mesh.getIBFaces();
1143  const StorageSite& cells = mesh.getCells();
1144  const StorageSite& faces = mesh.getFaces();
1145  const CRConnectivity& ibFaceToCells = mesh.getConnectivity(ibFaces,cells);
1146 
1147  const VectorT3Array& xFaces =
1148  dynamic_cast<const VectorT3Array&>(_coordField[faces]);
1149 
1150  const VectorT3Array& xCells =
1151  dynamic_cast<const VectorT3Array&>(_coordField[cells]);
1152 
1153  const Array<int>& ibFCRow = ibFaceToCells.getRow();
1154  const Array<int>& ibFCCol = ibFaceToCells.getCol();
1155 
1156  const int nIBFaces = ibFaces.getCount();
1157 
1158  const Array<int>& ibFaceIndices = mesh.getIBFaceList();
1159 
1160  shared_ptr<IMatrix> cellToIB(new IMatrix(ibFaceToCells));
1161 
1162  Array<T>& cellToIBCoeff = cellToIB->getCoeff();
1163 
1164  const bool is2D = mesh.getDimension() == 2;
1165  const bool is3D = mesh.getDimension() == 3;
1166 
1167  /***********************************************************************/
1168  // default is to use linear least square interpolation
1169  // if the matrix determinant is too small
1170  // then switch to distance weighted interpolation
1171  /***********************************************************************/
1172 
1173  /**********linear least square interpolation*********/
1174  // X=x-xf Y=y-yf Z=Z-zf
1175  //matrix M=[1 X1 Y1 Z1]
1176  // [1 X2 Y2 Z2]
1177  // ...........
1178  // [1 Xn Yn Zn]
1179  //coefficient matrix A=[a, b, c, d]T
1180  //velocity element v = M * A
1181  //linear relation v = a + b*X + c*Y + d*Z
1182  //to make least square
1183  //matrix A = (M(T)*M)^(-1)*M(T)*v
1184  //so, velocity at face is vface = a + b*Xf + c*Yf + d*Zf = a
1185  //which is the first row of matrix A
1186  //so, the weight function should be the first row of (M(T)*M)^(-1)*M(T)
1187  //note Q = M(T)*M and Qinv = Q^(-1)
1188  //the following code is to calculate it
1189  //insteading of doing full matrix operation, only nessesary operation on entries are performed
1190  //when dealing with coodinates with small numbers, like in micrometers
1191  //scaling the coodinates does not change the coefficients but improve the matrix quality
1192 
1193  // FILE * fp = fopen("/home/lin/work/app-memosa/src/fvm/verification/Structure_Electrostatics_Interaction/2D_beam/test/coeff.dat", "w");
1194 
1195 #if 1
1196 
1197 #if 0
1198  ofstream debugFileFluid;
1199  ofstream debugFileSolid;
1200  stringstream ss(stringstream::in | stringstream::out);
1201  ss << MPI::COMM_WORLD.Get_rank();
1202  string fname1 = "IBinterpolationFluid_proc" + ss.str() + ".dat";
1203  string fname2 = "IBinterpolationSolid_proc" + ss.str() + ".dat";
1204  debugFileFluid.open( fname1.c_str() );
1205  debugFileSolid.open( fname2.c_str() );
1206  ss.str("");
1207  const Array<int>& localToGlobal = mesh.getLocalToGlobal();
1208  const CRConnectivity& faceCells = mesh.getAllFaceCells();
1209 #endif
1210 
1211  for(int n=0; n<nIBFaces; n++)
1212  {
1213  const int f = ibFaceIndices[n];
1214  T wt(0);
1215  T wtSum(0);
1216  T det(0);
1217  int nnb(0);
1218  T scale(1.0e6);
1219 
1220  SquareMatrix<T,4> Q(NumTypeTraits<T>::getZero());
1221  SquareMatrix<T,4> Qinv(NumTypeTraits<T>::getZero());
1222  SquareMatrix<T,3> QQ(NumTypeTraits<T>::getZero());
1223  SquareMatrix<T,3> QQinv(NumTypeTraits<T>::getZero());
1224 
1225  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++){
1226  const int c = ibFCCol[nc];
1227  VectorT3 dr((xCells[c]-xFaces[f])*scale);
1228  //if (f ==200){
1229  // cout << f <<" " << c <<endl;
1230  // cout << xCells[c] << endl;
1231  // cout << xFaces[f] << endl;
1232  // }
1233  Q(0,0) += 1.0;
1234  Q(0,1) += dr[0];
1235  Q(0,2) += dr[1];
1236  Q(0,3) += dr[2];
1237  Q(1,1) += dr[0]*dr[0];
1238  Q(1,2) += dr[0]*dr[1];
1239  Q(1,3) += dr[0]*dr[2];
1240  Q(2,2) += dr[1]*dr[1];
1241  Q(2,3) += dr[1]*dr[2];
1242  Q(3,3) += dr[2]*dr[2];
1243  nnb++;
1244  }
1245 
1246  //if (nnb <=4){
1247  // throw CException("not enough cell or particle neighbors for ib face to interpolate!");
1248 
1249  //}
1250  //symetric matrix
1251  for(int i=0; i<4; i++){
1252  for(int j=0; j<i; j++){
1253  Q(i,j) = Q(j,i);
1254  }
1255  }
1256 
1257  // calculate the determinant of the matrix Q
1258  // if 3D mesh, then det(Q)
1259  // if 2D mesh, then det(QQ) where QQ is the 3x3 subset of Q
1260 
1261 
1262  if (is2D) {
1263  for(int i=0; i<3; i++){
1264  for(int j=0; j<3; j++){
1265  QQ(i,j)=Q(i,j);
1266  }
1267  }
1268  det = determinant(QQ);
1269  }
1270 
1271  if (is3D) {
1272  det = determinant(Q, 4);
1273  }
1274 
1275  // linear least square interpolation if the matrix is not singular
1276  //if (nnb >=10){
1277  if (fabs(det) > 1.0 ){
1278  if(is2D){
1279  QQinv = inverse(QQ);
1280  for(int i=0; i<3; i++){
1281  for(int j=0; j<3; j++){
1282  Qinv(i,j)=QQinv(i,j);
1283  }
1284  }
1285  }
1286 
1287  if(is3D){
1288  Qinv = inverse(Q, 4);
1289  }
1290 
1291  //calculate Qinv*M(T) get the first row element, put in coeffMatrix
1292  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++) {
1293  const int c = ibFCCol[nc];
1294  VectorT3 dr((xCells[c]-xFaces[f])*scale);
1295  wt = Qinv(0,0);
1296  for (int i=1; i<=3; i++){
1297  wt += Qinv(0,i)*dr[i-1];
1298  }
1299  cellToIBCoeff[nc] = wt;
1300  wtSum += wt;
1301  }
1302 
1303  if (wtSum > 1.01 || wtSum < 0.99){
1304  cout << "face " << n <<" has wrong wtsum " << wtSum << endl;
1305  cout << n << endl;
1306  cout << "ibface " << xFaces[f][0] << " " << xFaces[f][1] << " " << xFaces[f][2] << " " << endl;
1307  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++) {
1308  const int c = ibFCCol[nc];
1309  cout << "fluid cells " << xCells[c][0] << " " << xCells[c][1] << " " << xCells[c][2] << " " <<cellToIBCoeff[nc] << endl;
1310  }
1311  }
1312  }
1313 
1314 
1315  else { //if matrix is singular, use distance weighted interpolation
1316  //cout << "warning: IBM interpolation switched to distance weighted method for face " << f << endl;
1317  //cout << xFaces[f][0] << " " << xFaces[f][1] << " " << xFaces[f][2] << " " << endl;
1318  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++)
1319  {
1320  const int c = ibFCCol[nc];
1321  VectorT3 dr(xCells[c]-xFaces[f]);
1322  T wt = 1.0/dot(dr,dr);
1323  cellToIBCoeff[nc] = wt;
1324  wtSum += wt;
1325  nnb++;
1326  //cout << "fluid cells " << xCells[c][0] << " " << xCells[c][1] << " " << xCells[c][2] << " " << endl;
1327  }
1328 
1329  if (nnb == 0)
1330  throw CException("no cell or particle neighbors for ib face");
1331 
1332  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++)
1333  {
1334  cellToIBCoeff[nc] /= wtSum;
1335  }
1336 
1337  }
1338  }
1339 #endif
1340 
1341 #if 0
1342  debugFileFluid.close();
1343  debugFileSolid.close();
1344 #endif
1345 
1346 
1347  //fclose(fp);
1348 #if 0
1349 
1350  /**********second order least square interpolation*********/
1351  // X=x-xf Y=y-yf Z=Z-zf
1352  //matrix M=[1 X1 Y1 Z1 X1*X1 Y1*Y1 Z1*Z1 X1*Y1 Y1*Z1 X1*Z1]
1353  // [1 X2 Y2 Z2 ...................................]
1354  // ...........
1355  // [1 Xn Yn Zn Xn*Xn Yn*Yn Zn*Zn Xn*Yn Yn*Zn Xn*Zn]
1356  //coefficient matrix A=[a0, a1, a2, a3, a4, a5, a6, a7, a8, a9]T
1357  //velocity element v = M * A
1358  //linear relation v = a0 + a1*X + a2*Y + a3*Z + a4*X*X + a5*Y*Y + a6*Z*Z + a7*X*Y + a8*Y*Z + a9*X*Z
1359  //to make least square
1360  //matrix A = (M(T)*M)^(-1)*M(T)*v
1361  //so, velocity at face is vface = a0
1362  //which is the first row of matrix A
1363  //so, the weight function should be the first row of (M(T)*M)^(-1)*M(T)
1364  //note Q = M(T)*M and Qinv = Q^(-1)
1365  //the following code is to calculate it
1366  //insteading of doing full matrix operation, only nessesary operation on entries are performed
1367 
1368 
1369 
1370  for(int n=0; n<nIBFaces; n++)
1371  {
1372  const int f = ibFaceIndices[n];
1373  T wt(0);
1374  int nnb(0);
1375  T scale(1.0e6);
1376 
1377  if (is2D){
1378  const int size = 6;
1379  SquareMatrix<T,size> Q(SquareMatrix<T,size>::zero());
1380  SquareMatrix<T,size> Qinv(0);
1381  //cout << xFaces[f][0] << " " << xFaces[f][1] << " " << xFaces[f][2] << " " << endl;
1382  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++) {
1383  const int c = ibFCCol[nc];
1384  VectorT3 dr((xCells[c]-xFaces[f])*scale);
1385  //cout << xCells[c][0] << " " << xCells[c][1] << " " << xCells[c][2] << " " << endl;
1386 
1387  Q(0,0) += 1.0;
1388  Q(0,1) += dr[0];
1389  Q(0,2) += dr[1];
1390  Q(0,3) += dr[0]*dr[0];
1391  Q(0,4) += dr[1]*dr[1];
1392  Q(0,5) += dr[0]*dr[1];
1393 
1394  Q(1,1) += dr[0]*dr[0];
1395  Q(1,2) += dr[0]*dr[1];
1396  Q(1,3) += dr[0]*dr[0]*dr[0];
1397  Q(1,4) += dr[0]*dr[1]*dr[1];
1398  Q(1,5) += dr[0]*dr[0]*dr[1];
1399 
1400  Q(2,2) += dr[1]*dr[1];
1401  Q(2,3) += dr[1]*dr[0]*dr[0];
1402  Q(2,4) += dr[1]*dr[1]*dr[1];
1403  Q(2,5) += dr[1]*dr[0]*dr[1];
1404 
1405  Q(3,3) += dr[0]*dr[0]*dr[0]*dr[0];
1406  Q(3,4) += dr[0]*dr[0]*dr[1]*dr[1];
1407  Q(3,5) += dr[0]*dr[0]*dr[0]*dr[1];
1408 
1409  Q(4,4) += dr[1]*dr[1]*dr[1]*dr[1];
1410  Q(4,5) += dr[1]*dr[1]*dr[0]*dr[1];
1411 
1412  Q(5,5) += dr[0]*dr[1]*dr[0]*dr[1];
1413 
1414  nnb++;
1415  }
1416 
1417  if (nnb < size)
1418  throw CException("not enough cell or particle neighbors for ib face to interpolate!");
1419 
1420  //symetric matrix
1421  for(int i=0; i<size; i++){
1422  for(int j=0; j<i; j++){
1423  Q(i,j)=Q(j,i);
1424  }
1425  }
1426 
1427  //calculate the inverse of Q(6x6)
1428  Qinv = inverse(Q, size);
1429 
1430  //calculate Qinv*M(T) get the first row element, put in coeffMatrix
1431  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++) {
1432  const int c = ibFCCol[nc];
1433  VectorT3 dr((xCells[c]-xFaces[f])*scale);
1434  wt = Qinv(0,0);
1435  wt += Qinv(0,1)*dr[0];
1436  wt += Qinv(0,2)*dr[1];
1437  wt += Qinv(0,3)*dr[0]*dr[0];
1438  wt += Qinv(0,4)*dr[1]*dr[1];
1439  wt += Qinv(0,5)*dr[0]*dr[1];
1440  cellToIBCoeff[nc] = wt;
1441  //cout<<n<<" cells "<<nc<<" "<<cellToIBCoeff[nc]<<endl;
1442  }
1443  }
1444 
1445  if (is3D) {
1446  const int size = 10;
1447  SquareMatrix<T,size> Q(0);
1448  SquareMatrix<T,size> Qinv(0);
1449 
1450  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++)
1451  {
1452  const int c = ibFCCol[nc];
1453  VectorT3 dr((xCells[c]-xFaces[f])*scale);
1454  Q(0,0) += 1.0;
1455  Q(0,1) += dr[0];
1456  Q(0,2) += dr[1];
1457  Q(0,3) += dr[2];
1458  Q(0,4) += dr[0]*dr[0];
1459  Q(0,5) += dr[1]*dr[1];
1460  Q(0,6) += dr[2]*dr[2];
1461  Q(0,7) += dr[0]*dr[1];
1462  Q(0,8) += dr[1]*dr[2];
1463  Q(0,9) += dr[0]*dr[2];
1464 
1465  Q(1,1) += dr[0]*dr[0];
1466  Q(1,2) += dr[0]*dr[1];
1467  Q(1,3) += dr[0]*dr[2];
1468  Q(1,4) += dr[0]*dr[0]*dr[0];
1469  Q(1,5) += dr[0]*dr[1]*dr[1];
1470  Q(1,6) += dr[0]*dr[2]*dr[2];
1471  Q(1,7) += dr[0]*dr[0]*dr[1];
1472  Q(1,8) += dr[0]*dr[1]*dr[2];
1473  Q(1,9) += dr[0]*dr[0]*dr[2];
1474 
1475  Q(2,2) += dr[1]*dr[1];
1476  Q(2,3) += dr[1]*dr[2];
1477  Q(2,4) += dr[1]*dr[0]*dr[0];
1478  Q(2,5) += dr[1]*dr[1]*dr[1];
1479  Q(2,6) += dr[1]*dr[2]*dr[2];
1480  Q(2,7) += dr[1]*dr[0]*dr[1];
1481  Q(2,8) += dr[1]*dr[1]*dr[2];
1482  Q(2,9) += dr[1]*dr[0]*dr[2];
1483 
1484  Q(3,3) += dr[2]*dr[2];
1485  Q(3,4) += dr[2]*dr[0]*dr[0];
1486  Q(3,5) += dr[2]*dr[1]*dr[1];
1487  Q(3,6) += dr[2]*dr[2]*dr[2];
1488  Q(3,7) += dr[2]*dr[0]*dr[1];
1489  Q(3,8) += dr[2]*dr[1]*dr[2];
1490  Q(3,8) += dr[2]*dr[0]*dr[2];
1491 
1492  Q(4,4) += dr[0]*dr[0]*dr[0]*dr[0];
1493  Q(4,5) += dr[0]*dr[0]*dr[1]*dr[1];
1494  Q(4,6) += dr[0]*dr[0]*dr[2]*dr[2];
1495  Q(4,7) += dr[0]*dr[0]*dr[0]*dr[1];
1496  Q(4,8) += dr[0]*dr[0]*dr[1]*dr[2];
1497  Q(4,9) += dr[0]*dr[0]*dr[0]*dr[2];
1498 
1499  Q(5,5) += dr[1]*dr[1]*dr[1]*dr[1];
1500  Q(5,6) += dr[1]*dr[1]*dr[2]*dr[2];
1501  Q(5,7) += dr[1]*dr[1]*dr[0]*dr[1];
1502  Q(5,8) += dr[1]*dr[1]*dr[1]*dr[2];
1503  Q(5,9) += dr[1]*dr[1]*dr[0]*dr[2];
1504 
1505  Q(6,6) += dr[2]*dr[2]*dr[2]*dr[2];
1506  Q(6,7) += dr[2]*dr[2]*dr[0]*dr[1];
1507  Q(6,8) += dr[2]*dr[2]*dr[1]*dr[2];
1508  Q(6,9) += dr[2]*dr[2]*dr[0]*dr[2];
1509 
1510  Q(7,7) += dr[0]*dr[1]*dr[0]*dr[1];
1511  Q(7,8) += dr[0]*dr[1]*dr[1]*dr[2];
1512  Q(7,9) += dr[0]*dr[1]*dr[0]*dr[2];
1513 
1514  Q(8,8) += dr[1]*dr[2]*dr[1]*dr[2];
1515  Q(8,9) += dr[1]*dr[2]*dr[0]*dr[2];
1516 
1517  Q(9,9) += dr[0]*dr[2]*dr[0]*dr[2];
1518 
1519  nnb++;
1520  }
1521  if (nnb < size)
1522  throw CException("not enough cell or particle neighbors for ib face to interpolate!");
1523 
1524  //symetric matrix
1525  for(int i=0; i<size; i++){
1526  for(int j=0; j<i; j++){
1527  Q(i,j)=Q(j,i);
1528  }
1529  }
1530 
1531  //calculate the inverse of Q(10x10)
1532  Qinv = inverse(Q, size);
1533 
1534 
1535  //calculate Qinv*M(T) get the first row element, put in coeffMatrix
1536  for(int nc=ibFCRow[n]; nc<ibFCRow[n+1]; nc++)
1537  {
1538  const int c = ibFCCol[nc];
1539  VectorT3 dr((xCells[c]-xFaces[f])*scale);
1540  wt = Qinv(0,0);
1541  wt += Qinv(0,1)*dr[0];
1542  wt += Qinv(0,2)*dr[1];
1543  wt += Qinv(0,3)*dr[2];
1544  wt += Qinv(0,4)*dr[0]*dr[0];
1545  wt += Qinv(0,5)*dr[1]*dr[1];
1546  wt += Qinv(0,6)*dr[2]*dr[2];
1547  wt += Qinv(0,7)*dr[0]*dr[1];
1548  wt += Qinv(0,8)*dr[1]*dr[2];
1549  wt += Qinv(0,9)*dr[0]*dr[2];
1550  cellToIBCoeff[nc] = wt;
1551  //cout<<n<<" cells "<<nc<<" "<<cellToIBCoeff[nc]<<endl;
1552  }
1553  }
1554  }
1555 #endif
1556  GeomFields::SSPair key1(&faces,&cells);
1557  this->_geomFields._interpolationMatrices[key1] = cellToIB;
1558 }
CRConnectivity::getCol
const Array< int > & getCol() const
Definition: CRConnectivity.h:183
CRConnectivity::getRow
const Array< int > & getRow() const
Definition: CRConnectivity.h:182
Mesh::getIBFaces
const StorageSite & getIBFaces() const
Definition: Mesh.h:111
Mesh::getConnectivity
const CRConnectivity & getConnectivity(const StorageSite &from, const StorageSite &to) const
Definition: Mesh.cpp:416
Mesh::getIBFaceList
const Array< int > & getIBFaceList() const
Definition: Mesh.cpp:686
Mesh::getAllFaceCells
const CRConnectivity & getAllFaceCells() const
Definition: Mesh.cpp:378
Mesh::getLocalToGlobal
Array< int > & getLocalToGlobal()
Definition: Mesh.h:240
GeomFields::_interpolationMatrices
map< SSPair, shared_ptr< Matrix > > _interpolationMatrices
Definition: GeomFields.h:50
Mesh::getFaces
const StorageSite & getFaces() const
Definition: Mesh.h:108
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
inverse
SquareMatrix< T, 2 > inverse(SquareMatrix< T, 2 > &a)
Definition: MatrixOperation.h:118
fabs
Tangent fabs(const Tangent &a)
Definition: Tangent.h:312
Mesh::isShell
bool isShell() const
Definition: Mesh.h:323
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
Mesh::getDimension
int getDimension() const
Definition: Mesh.h:105
dot
T dot(const Vector< T, 3 > &a, const Vector< T, 3 > &b)
Definition: Vector.h:253
GeomFields::SSPair
pair< const StorageSite *, const StorageSite * > SSPair
Definition: GeomFields.h:48
determinant
T determinant(SquareMatrix< T, 2 > &a)
Definition: MatrixOperation.h:67
template<class T >
void MeshMetricsCalculator< T >::computeSolidInterpolationMatrices ( const StorageSite &  particles)

Definition at line 2418 of file MeshMetricsCalculator_impl.h.

2419 {
2420  const int numMeshes = _meshes.size();
2421  for (int n=0; n<numMeshes; n++)
2422  {
2423  const Mesh& mesh = *_meshes[n];
2424  computeSolidInterpolationMatrices(mesh,p);
2425  }
2426 }
Mesh
Definition: Mesh.h:49
MeshMetricsCalculator::computeSolidInterpolationMatrices
void computeSolidInterpolationMatrices(const StorageSite &particles)
Definition: MeshMetricsCalculator_impl.h:2418
Model::_meshes
const MeshList _meshes
Definition: Model.h:29
template<class T >
void MeshMetricsCalculator< T >::computeSolidInterpolationMatrices ( const Mesh &  mesh,
const StorageSite &  particles 
)
private

Definition at line 1565 of file MeshMetricsCalculator_impl.h.

References determinant(), dot(), fabs(), Mesh::getCells(), CRConnectivity::getCol(), Mesh::getConnectivity(), StorageSite::getCount(), Mesh::getDimension(), CRConnectivity::getRow(), inverse(), and Mesh::isShell().

1567 {
1568 
1569  if (mesh.isShell())
1570  return;
1571 
1572  typedef CRMatrixTranspose<T,T,T> IMatrix;
1573 
1574  const StorageSite& cells = mesh.getCells();
1575  const CRConnectivity& solidFacesToCells =
1576  mesh.getConnectivity(solidFaces,cells);
1577 
1578 
1579  const VectorT3Array& xCells =
1580  dynamic_cast<const VectorT3Array&>(_coordField[cells]);
1581 
1582  const VectorT3Array& xFaces =
1583  dynamic_cast<const VectorT3Array&>(_coordField[solidFaces]);
1584 
1585  const Array<int>& sFCRow = solidFacesToCells.getRow();
1586  const Array<int>& sFCCol = solidFacesToCells.getCol();
1587 
1588  const int nSolidFaces = solidFaces.getCount();
1589 
1590  shared_ptr<IMatrix> cellToSolidFaces(new IMatrix(solidFacesToCells));
1591 
1592  Array<T>& cellToSBCoeff = cellToSolidFaces->getCoeff();
1593 
1594  const bool is2D = mesh.getDimension() == 2;
1595  const bool is3D = mesh.getDimension() == 3;
1596 #if 1
1597  for(int f=0; f<nSolidFaces; f++)
1598  {
1599  T wt(0);
1600  T wtSum(0.0);
1601  int nnb(0);
1602  T det(0);
1603  T scale(1.0e6);
1604 
1605  SquareMatrix<T,4> Q(NumTypeTraits<T>::getZero());
1606  SquareMatrix<T,4> Qinv(NumTypeTraits<T>::getZero());
1607  SquareMatrix<T,3> QQ(NumTypeTraits<T>::getZero());
1608  SquareMatrix<T,3> QQinv(NumTypeTraits<T>::getZero());
1609 
1610  for(int nc=sFCRow[f]; nc<sFCRow[f+1]; nc++)
1611  {
1612  const int c = sFCCol[nc];
1613  //if (f ==168){
1614  //cout << f <<" " << c <<endl;
1615  // cout << xCells[c] << endl;
1616  // cout << xFaces[f] << endl;
1617  //}
1618  VectorT3 dr((xCells[c]-xFaces[f])*scale);
1619  Q(0,0) += 1.0;
1620  Q(0,1) += dr[0];
1621  Q(0,2) += dr[1];
1622  Q(0,3) += dr[2];
1623  Q(1,1) += dr[0]*dr[0];
1624  Q(1,2) += dr[0]*dr[1];
1625  Q(1,3) += dr[0]*dr[2];
1626  Q(2,2) += dr[1]*dr[1];
1627  Q(2,3) += dr[1]*dr[2];
1628  Q(3,3) += dr[2]*dr[2];
1629  nnb++;
1630  }
1631  if (nnb == 0)
1632  continue;
1633 
1634  //symetric matrix
1635  for(int i=0; i<4; i++)
1636  for(int j=0; j<i; j++)
1637  Q(i,j)=Q(j,i);
1638 
1639  if (is2D) {
1640  for(int i=0; i<3; i++){
1641  for(int j=0; j<3; j++){
1642  QQ(i,j)=Q(i,j);
1643  }
1644  }
1645  det = determinant(QQ);
1646  }
1647 
1648  if (is3D) {
1649  det = determinant(Q, 4);
1650  }
1651  //cout << "determinant " << f << " " << det <<endl;
1652  //if (nnb>=4){
1653  if (fabs(det) > 1.0){
1654  if(is2D){
1655  QQinv = inverse(QQ);
1656  for(int i=0; i<3; i++){
1657  for(int j=0; j<3; j++){
1658  Qinv(i,j)=QQinv(i,j);
1659  }
1660  }
1661  }
1662 
1663  if(is3D){
1664  Qinv = inverse(Q, 4);
1665  }
1666  //if (f==168)
1667  // cout <<"determinant " << det <<endl;
1668  for(int nc=sFCRow[f]; nc<sFCRow[f+1]; nc++)
1669  {
1670  const int c = sFCCol[nc];
1671  VectorT3 dr((xCells[c]-xFaces[f])*scale);
1672  wt = Qinv(0,0);
1673  for (int i=1; i<=3; i++)
1674  wt += Qinv(0,i)*dr[i-1];
1675 
1676  cellToSBCoeff[nc] = wt;
1677  //if (f==168)
1678  // cout<< " " << f <<" " << cellToSBCoeff[nc] << endl;
1679  }
1680  }
1681 
1682  //distance weighted interpolation
1683  else {
1684  //printf("warning: distance weighted interpolation for solid face is applied %i\n", f);
1685  for(int nc=sFCRow[f]; nc<sFCRow[f+1]; nc++)
1686  {
1687  const int c = sFCCol[nc];
1688  VectorT3 dr(xCells[c]-xFaces[f]);
1689  T wt = 1.0/dot(dr,dr);
1690  cellToSBCoeff[nc] = wt;
1691  wtSum += wt;
1692  nnb++;
1693  //if (f==100){
1694  // cout << "distance weight" <<endl;
1695  // cout << c << endl;
1696  // cout << dr << endl;
1697  // cout << wt << endl;
1698  //}
1699  }
1700 
1701  if (nnb == 0)
1702  throw CException("no cell or particle neighbors for ib face");
1703 
1704  for(int nc=sFCRow[f]; nc<sFCRow[f+1]; nc++)
1705  {
1706  cellToSBCoeff[nc] /= wtSum;
1707  //if (f==168)
1708  // cout <<"distance weight " << cellToSBCoeff[nc] << " " << wtSum << endl;
1709  }
1710  }
1711  }
1712 
1713 #endif
1714 
1715 #if 0
1716 
1717  /**********second order least square interpolation*********/
1718  // X=x-xf Y=y-yf Z=Z-zf
1719  //matrix M=[1 X1 Y1 Z1 X1*X1 Y1*Y1 Z1*Z1 X1*Y1 Y1*Z1 X1*Z1]
1720  // [1 X2 Y2 Z2 ...................................]
1721  // ...........
1722  // [1 Xn Yn Zn Xn*Xn Yn*Yn Zn*Zn Xn*Yn Yn*Zn Xn*Zn]
1723  //coefficient matrix A=[a0, a1, a2, a3, a4, a5, a6, a7, a8, a9]T
1724  //velocity element v = M * A
1725  //linear relation v = a0 + a1*X + a2*Y + a3*Z + a4*X*X + a5*Y*Y + a6*Z*Z + a7*X*Y + a8*Y*Z + a9*X*Z
1726  //to make least square
1727  //matrix A = (M(T)*M)^(-1)*M(T)*v
1728  //so, velocity at face is vface = a0
1729  //which is the first row of matrix A
1730  //so, the weight function should be the first row of (M(T)*M)^(-1)*M(T)
1731  //note Q = M(T)*M and Qinv = Q^(-1)
1732  //the following code is to calculate it
1733  //insteading of doing full matrix operation, only nessesary operation on entries are performed
1734 
1735 
1736 
1737  for(int f=0; f<nSolidFaces; f++)
1738  {
1739  T wt(0);
1740  int nnb(0);
1741  T scale(1.0e6);
1742 
1743  if (mesh.getDimension()== 2) {
1744  const int size = 6;
1745  SquareMatrix<T,size> Q(NumTypeTraits<T>::getZero());
1746  SquareMatrix<T,size> Qinv(NumTypeTraits<T>::getZero());
1747 
1748  for(int nc=sFCRow[f]; nc<sFCRow[f+1]; nc++) {
1749  const int c = sFCCol[nc];
1750  VectorT3 dr((xCells[c]-xFaces[f])*scale);
1751  Q(0,0) += 1.0;
1752  Q(0,1) += dr[0];
1753  Q(0,2) += dr[1];
1754  Q(0,3) += dr[0]*dr[0];
1755  Q(0,4) += dr[1]*dr[1];
1756  Q(0,5) += dr[0]*dr[1];
1757 
1758  Q(1,1) += dr[0]*dr[0];
1759  Q(1,2) += dr[0]*dr[1];
1760  Q(1,3) += dr[0]*dr[0]*dr[0];
1761  Q(1,4) += dr[0]*dr[1]*dr[1];
1762  Q(1,5) += dr[0]*dr[0]*dr[1];
1763 
1764  Q(2,2) += dr[1]*dr[1];
1765  Q(2,3) += dr[1]*dr[0]*dr[0];
1766  Q(2,4) += dr[1]*dr[1]*dr[1];
1767  Q(2,5) += dr[1]*dr[0]*dr[1];
1768 
1769  Q(3,3) += dr[0]*dr[0]*dr[0]*dr[0];
1770  Q(3,4) += dr[0]*dr[0]*dr[1]*dr[1];
1771  Q(3,5) += dr[0]*dr[0]*dr[0]*dr[1];
1772 
1773  Q(4,4) += dr[1]*dr[1]*dr[1]*dr[1];
1774  Q(4,5) += dr[1]*dr[1]*dr[0]*dr[1];
1775 
1776  Q(5,5) += dr[0]*dr[1]*dr[0]*dr[1];
1777 
1778  nnb++;
1779  }
1780 
1781  if (nnb < size){
1782  cout << nnb << endl;
1783  throw CException("not enough fluid neighbors for solid to interpolate!");
1784  }
1785 
1786  //symetric matrix
1787  for(int i=0; i<size; i++){
1788  for(int j=0; j<i; j++){
1789  Q(i,j)=Q(j,i);
1790  }
1791  }
1792 
1793 
1794  //calculate the inverse of Q(6x6)
1795  Qinv = inverse(Q, size);
1796 
1797  //calculate Qinv*M(T) get the first row element, put in coeffMatrix
1798  for(int nc=sFCRow[f]; nc<sFCRow[f+1]; nc++)
1799  {
1800  const int c = sFCCol[nc];
1801  VectorT3 dr((xCells[c]-xFaces[f])*scale);
1802  wt = Qinv(0,0);
1803  wt += Qinv(0,1)*dr[0];
1804  wt += Qinv(0,2)*dr[1];
1805  wt += Qinv(0,3)*dr[0]*dr[0];
1806  wt += Qinv(0,4)*dr[1]*dr[1];
1807  wt += Qinv(0,5)*dr[0]*dr[1];
1808 
1809  cellToSBCoeff[nc] = wt;
1810  //cout<<n<<" cells "<<nc<<" "<<cellToIBCoeff[nc]<<endl;
1811  }
1812  }
1813 
1814  else if (mesh.getDimension()== 3)
1815  {
1816  const int size = 10;
1817  SquareMatrix<T,size> Q(NumTypeTraits<T>::getZero());
1818  SquareMatrix<T,size> Qinv(NumTypeTraits<T>::getZero());
1819 
1820  for(int nc=sFCRow[f]; nc<sFCRow[f+1]; nc++) {
1821  const int c = sFCCol[nc];
1822  VectorT3 dr((xCells[c]-xFaces[f])*scale);
1823  Q(0,0) += 1.0;
1824  Q(0,1) += dr[0];
1825  Q(0,2) += dr[1];
1826  Q(0,3) += dr[2];
1827  Q(0,4) += dr[0]*dr[0];
1828  Q(0,5) += dr[1]*dr[1];
1829  Q(0,6) += dr[2]*dr[2];
1830  Q(0,7) += dr[0]*dr[1];
1831  Q(0,8) += dr[1]*dr[2];
1832  Q(0,9) += dr[0]*dr[2];
1833 
1834  Q(1,1) += dr[0]*dr[0];
1835  Q(1,2) += dr[0]*dr[1];
1836  Q(1,3) += dr[0]*dr[2];
1837  Q(1,4) += dr[0]*dr[0]*dr[0];
1838  Q(1,5) += dr[0]*dr[1]*dr[1];
1839  Q(1,6) += dr[0]*dr[2]*dr[2];
1840  Q(1,7) += dr[0]*dr[0]*dr[1];
1841  Q(1,8) += dr[0]*dr[1]*dr[2];
1842  Q(1,9) += dr[0]*dr[0]*dr[2];
1843 
1844  Q(2,2) += dr[1]*dr[1];
1845  Q(2,3) += dr[1]*dr[2];
1846  Q(2,4) += dr[1]*dr[0]*dr[0];
1847  Q(2,5) += dr[1]*dr[1]*dr[1];
1848  Q(2,6) += dr[1]*dr[2]*dr[2];
1849  Q(2,7) += dr[1]*dr[0]*dr[1];
1850  Q(2,8) += dr[1]*dr[1]*dr[2];
1851  Q(2,9) += dr[1]*dr[0]*dr[2];
1852 
1853  Q(3,3) += dr[2]*dr[2];
1854  Q(3,4) += dr[2]*dr[0]*dr[0];
1855  Q(3,5) += dr[2]*dr[1]*dr[1];
1856  Q(3,6) += dr[2]*dr[2]*dr[2];
1857  Q(3,7) += dr[2]*dr[0]*dr[1];
1858  Q(3,8) += dr[2]*dr[1]*dr[2];
1859  Q(3,8) += dr[2]*dr[0]*dr[2];
1860 
1861  Q(4,4) += dr[0]*dr[0]*dr[0]*dr[0];
1862  Q(4,5) += dr[0]*dr[0]*dr[1]*dr[1];
1863  Q(4,6) += dr[0]*dr[0]*dr[2]*dr[2];
1864  Q(4,7) += dr[0]*dr[0]*dr[0]*dr[1];
1865  Q(4,8) += dr[0]*dr[0]*dr[1]*dr[2];
1866  Q(4,9) += dr[0]*dr[0]*dr[0]*dr[2];
1867 
1868  Q(5,5) += dr[1]*dr[1]*dr[1]*dr[1];
1869  Q(5,6) += dr[1]*dr[1]*dr[2]*dr[2];
1870  Q(5,7) += dr[1]*dr[1]*dr[0]*dr[1];
1871  Q(5,8) += dr[1]*dr[1]*dr[1]*dr[2];
1872  Q(5,9) += dr[1]*dr[1]*dr[0]*dr[2];
1873 
1874  Q(6,6) += dr[2]*dr[2]*dr[2]*dr[2];
1875  Q(6,7) += dr[2]*dr[2]*dr[0]*dr[1];
1876  Q(6,8) += dr[2]*dr[2]*dr[1]*dr[2];
1877  Q(6,9) += dr[2]*dr[2]*dr[0]*dr[2];
1878 
1879  Q(7,7) += dr[0]*dr[1]*dr[0]*dr[1];
1880  Q(7,8) += dr[0]*dr[1]*dr[1]*dr[2];
1881  Q(7,9) += dr[0]*dr[1]*dr[0]*dr[2];
1882 
1883  Q(8,8) += dr[1]*dr[2]*dr[1]*dr[2];
1884  Q(8,9) += dr[1]*dr[2]*dr[0]*dr[2];
1885 
1886  Q(9,9) += dr[0]*dr[2]*dr[0]*dr[2];
1887 
1888  nnb++;
1889  }
1890  if (nnb < size)
1891  throw CException("not enough cell or particle neighbors for solid to interpolate!");
1892 
1893  //symetric matrix
1894  for(int i=0; i<size; i++){
1895  for(int j=0; j<i; j++){
1896  Q(i,j)=Q(j,i);
1897  }
1898  }
1899  //calculate the inverse of Q(10x10)
1900  Qinv = inverse(Q, size);
1901 
1902  //calculate Qinv*M(T) get the first row element, put in coeffMatrix
1903  for(int nc=sFCRow[f]; nc<sFCRow[f+1]; nc++)
1904  {
1905  const int c = sFCCol[nc];
1906  VectorT3 dr((xCells[c]-xFaces[f])*scale);
1907  wt = Qinv(0,0);
1908  wt += Qinv(0,1)*dr[0];
1909  wt += Qinv(0,2)*dr[1];
1910  wt += Qinv(0,3)*dr[2];
1911  wt += Qinv(0,4)*dr[0]*dr[0];
1912  wt += Qinv(0,5)*dr[1]*dr[1];
1913  wt += Qinv(0,6)*dr[2]*dr[2];
1914  wt += Qinv(0,7)*dr[0]*dr[1];
1915  wt += Qinv(0,8)*dr[1]*dr[2];
1916  wt += Qinv(0,9)*dr[0]*dr[2];
1917  cellToSBCoeff[nc] = wt;
1918  //cout<<n<<" cells "<<nc<<" "<<cellToIBCoeff[nc]<<endl;
1919  }
1920  }
1921  }
1922 #endif
1923  GeomFields::SSPair key1(&solidFaces,&cells);
1924  this->_geomFields._interpolationMatrices[key1] = cellToSolidFaces;
1925 }
CRConnectivity::getCol
const Array< int > & getCol() const
Definition: CRConnectivity.h:183
CRConnectivity::getRow
const Array< int > & getRow() const
Definition: CRConnectivity.h:182
Mesh::getConnectivity
const CRConnectivity & getConnectivity(const StorageSite &from, const StorageSite &to) const
Definition: Mesh.cpp:416
GeomFields::_interpolationMatrices
map< SSPair, shared_ptr< Matrix > > _interpolationMatrices
Definition: GeomFields.h:50
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
inverse
SquareMatrix< T, 2 > inverse(SquareMatrix< T, 2 > &a)
Definition: MatrixOperation.h:118
fabs
Tangent fabs(const Tangent &a)
Definition: Tangent.h:312
Mesh::isShell
bool isShell() const
Definition: Mesh.h:323
Mesh::getDimension
int getDimension() const
Definition: Mesh.h:105
dot
T dot(const Vector< T, 3 > &a, const Vector< T, 3 > &b)
Definition: Vector.h:253
GeomFields::SSPair
pair< const StorageSite *, const StorageSite * > SSPair
Definition: GeomFields.h:48
determinant
T determinant(SquareMatrix< T, 2 > &a)
Definition: MatrixOperation.h:67
template<class T >
void MeshMetricsCalculator< T >::createNodeDisplacement ( )
template<class T >
void MeshMetricsCalculator< T >::eraseIBInterpolationMatrices ( const StorageSite &  particles)

Definition at line 2393 of file MeshMetricsCalculator_impl.h.

References Mesh::eraseConnectivity(), Mesh::getCells(), and Mesh::getIBFaces().

2394 {
2395  const int numMeshes = _meshes.size();
2396  for (int n=0; n<numMeshes; n++)
2397  {
2398  const Mesh& mesh = *_meshes[n];
2399  const StorageSite& ibFaces = mesh.getIBFaces();
2400  const StorageSite& cells = mesh.getCells();
2401 
2402  GeomFields::SSPair key1(&ibFaces,&cells);
2403  this->_geomFields._interpolationMatrices.erase(key1);
2404  mesh.eraseConnectivity(ibFaces,cells);
2405 
2406  GeomFields::SSPair key2(&ibFaces,&p);
2407  this->_geomFields._interpolationMatrices.erase(key2);
2408  mesh.eraseConnectivity(ibFaces,p);
2409 
2410  GeomFields::SSPair key3(&p,&cells);
2411  this->_geomFields._interpolationMatrices.erase(key3);
2412  mesh.eraseConnectivity(p,cells);
2413  }
2414 }
Mesh::getIBFaces
const StorageSite & getIBFaces() const
Definition: Mesh.h:111
Mesh::eraseConnectivity
void eraseConnectivity(const StorageSite &rowSite, const StorageSite &colSite) const
Definition: Mesh.cpp:359
Mesh
Definition: Mesh.h:49
GeomFields::_interpolationMatrices
map< SSPair, shared_ptr< Matrix > > _interpolationMatrices
Definition: GeomFields.h:50
Model::_meshes
const MeshList _meshes
Definition: Model.h:29
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
GeomFields::SSPair
pair< const StorageSite *, const StorageSite * > SSPair
Definition: GeomFields.h:48
template<class T >
void MeshMetricsCalculator< T >::init ( )
virtual

Implements Model.

Definition at line 1946 of file MeshMetricsCalculator_impl.h.

References Mesh::getAllFaceCells(), Mesh::getCells(), StorageSite::getCount(), StorageSite::getCountLevel1(), Mesh::getFaces(), Mesh::getPeriodicFacePairs(), Mesh::IBTYPE_FLUID, and Mesh::isShell().

1947 {
1948  const int numMeshes = _meshes.size();
1949  for (int n=0; n<numMeshes; n++)
1950  {
1951  const Mesh& mesh = *_meshes[n];
1952  if (!mesh.isShell())
1953  calculateNodeCoordinates(mesh);
1954  }
1955 
1956  for (int n=0; n<numMeshes; n++)
1957  {
1958  const Mesh& mesh = *_meshes[n];
1959  if (!mesh.isShell())
1960  {
1961  calculateFaceAreas(mesh);
1962  calculateFaceAreaMag(mesh);
1963  calculateFaceCentroids(mesh);
1964  }
1965  }
1966  for (int n=0; n<numMeshes; n++)
1967  {
1968  const Mesh& mesh = *_meshes[n];
1969  calculateCellCentroids(mesh);
1970  }
1971  _coordField.syncLocal();
1972 
1973 
1974  // fix periodic ghost cell coordinates
1975 
1976  for (int n=0; n<numMeshes; n++)
1977  {
1978  const Mesh& mesh = *_meshes[n];
1979  const Mesh::PeriodicFacePairs periodicFacePairs = mesh.getPeriodicFacePairs();
1980 
1981  if (periodicFacePairs.size() == 0)
1982  continue;
1983 
1984  const StorageSite& cells = mesh.getCells();
1985  const StorageSite& faces = mesh.getFaces();
1986 
1987  const CRConnectivity& faceCells = mesh.getAllFaceCells();
1988 
1989  const VectorT3Array& faceCoord =
1990  dynamic_cast<const VectorT3Array&>(_coordField[faces]);
1991 
1992  VectorT3Array& cellCoord =
1993  dynamic_cast<VectorT3Array&>(_coordField[cells]);
1994 
1995 
1996  for(Mesh::PeriodicFacePairs::const_iterator pos = periodicFacePairs.begin();
1997  pos!=periodicFacePairs.end();
1998  ++pos)
1999  {
2000  const int lf = pos->first;
2001  const int rf = pos->second;
2002  const VectorT3 dr = faceCoord[rf] - faceCoord[lf];
2003  cellCoord[faceCells(lf,1)] += dr;
2004  cellCoord[faceCells(rf,1)] -= dr;
2005  }
2006  }
2007 
2008  for (int n=0; n<numMeshes; n++)
2009  {
2010  const Mesh& mesh = *_meshes[n];
2011  calculateCellVolumes(mesh);
2012  }
2013 
2014  for (int n=0; n<numMeshes; n++)
2015  {
2016  const Mesh& mesh = *_meshes[n];
2017  const StorageSite& cells = mesh.getCells();
2018  const StorageSite& faces = mesh.getFaces();
2019 
2020  const int cellCount = cells.getCountLevel1();
2021  if ( (cellCount > 0) ) //&& (!mesh.isShell()) )
2022  {
2023  shared_ptr<IntArray> ibTypePtr(new IntArray(cells.getCountLevel1()));
2024  *ibTypePtr = Mesh::IBTYPE_FLUID;
2025  _geomFields.ibType.addArray(cells,ibTypePtr);
2026 
2027  shared_ptr<IntArray> ibFaceIndexPtr(new IntArray(faces.getCount()));
2028  *ibFaceIndexPtr = -1;
2029  _geomFields.ibFaceIndex.addArray(faces,ibFaceIndexPtr);
2030 
2031  if (_transient)
2032  {
2033  shared_ptr<IntArray> ibTypeN1Ptr(new IntArray(cells.getCountLevel1()));
2034  *ibTypeN1Ptr = Mesh::IBTYPE_FLUID;
2035  _geomFields.ibTypeN1.addArray(cells,ibTypeN1Ptr);
2036  }
2037  }
2038  }
2039 
2040  _volumeField.syncLocal();
2041 }
GeomFields::ibTypeN1
Field ibTypeN1
Definition: GeomFields.h:39
MeshMetricsCalculator::calculateCellVolumes
virtual void calculateCellVolumes(const Mesh &mesh)
Definition: MeshMetricsCalculator_impl.h:394
Mesh::getPeriodicFacePairs
PeriodicFacePairs & getPeriodicFacePairs()
Definition: Mesh.h:337
Mesh
Definition: Mesh.h:49
Mesh::PeriodicFacePairs
map< int, int > PeriodicFacePairs
Definition: Mesh.h:67
GeomFields::ibType
Field ibType
Definition: GeomFields.h:38
MeshMetricsCalculator::calculateFaceCentroids
virtual void calculateFaceCentroids(const Mesh &mesh)
Definition: MeshMetricsCalculator_impl.h:60
Mesh::getAllFaceCells
const CRConnectivity & getAllFaceCells() const
Definition: Mesh.cpp:378
MeshMetricsCalculator::calculateNodeCoordinates
virtual void calculateNodeCoordinates(const Mesh &mesh)
Definition: MeshMetricsCalculator_impl.h:36
Mesh::getFaces
const StorageSite & getFaces() const
Definition: Mesh.h:108
Model::_meshes
const MeshList _meshes
Definition: Model.h:29
GeomFields::ibFaceIndex
Field ibFaceIndex
Definition: GeomFields.h:40
MeshMetricsCalculator::calculateCellCentroids
virtual void calculateCellCentroids(const Mesh &mesh)
Definition: MeshMetricsCalculator_impl.h:130
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
MeshMetricsCalculator::calculateFaceAreas
virtual void calculateFaceAreas(const Mesh &mesh)
Definition: MeshMetricsCalculator_impl.h:240
MeshMetricsCalculator::calculateFaceAreaMag
virtual void calculateFaceAreaMag(const Mesh &mesh)
Definition: MeshMetricsCalculator_impl.h:373
StorageSite::getCountLevel1
int getCountLevel1() const
Definition: StorageSite.h:72
Field::addArray
void addArray(const StorageSite &, shared_ptr< ArrayBase > a)
Definition: Field.cpp:72
Mesh::isShell
bool isShell() const
Definition: Mesh.h:323
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
Field::syncLocal
void syncLocal()
Definition: Field.cpp:334
template<class T >
void MeshMetricsCalculator< T >::recalculate ( )

Definition at line 2046 of file MeshMetricsCalculator_impl.h.

References Mesh::isShell().

2047 {
2048  const int numMeshes = _meshes.size();
2049  for (int n=0; n<numMeshes; n++)
2050  {
2051  const Mesh& mesh = *_meshes[n];
2052  if (!mesh.isShell())
2053  {
2054  calculateFaceAreas(mesh);
2055  calculateFaceAreaMag(mesh);
2056  calculateFaceCentroids(mesh);
2057  calculateCellCentroids(mesh);
2058  }
2059  }
2060 
2061  _volumeField.syncLocal();
2062  _coordField.syncLocal();
2063 }
Mesh
Definition: Mesh.h:49
MeshMetricsCalculator::calculateFaceCentroids
virtual void calculateFaceCentroids(const Mesh &mesh)
Definition: MeshMetricsCalculator_impl.h:60
Model::_meshes
const MeshList _meshes
Definition: Model.h:29
MeshMetricsCalculator::calculateCellCentroids
virtual void calculateCellCentroids(const Mesh &mesh)
Definition: MeshMetricsCalculator_impl.h:130
MeshMetricsCalculator::calculateFaceAreas
virtual void calculateFaceAreas(const Mesh &mesh)
Definition: MeshMetricsCalculator_impl.h:240
MeshMetricsCalculator::calculateFaceAreaMag
virtual void calculateFaceAreaMag(const Mesh &mesh)
Definition: MeshMetricsCalculator_impl.h:373
Mesh::isShell
bool isShell() const
Definition: Mesh.h:323
Field::syncLocal
void syncLocal()
Definition: Field.cpp:334
template<class T >
void MeshMetricsCalculator< T >::recalculate_deform ( )

Definition at line 2098 of file MeshMetricsCalculator_impl.h.

References Mesh::isShell().

2099 {
2100  const int numMeshes = _meshes.size();
2101  for (int n=0; n<numMeshes; n++)
2102  {
2103  const Mesh& mesh = *_meshes[n];
2104  if (!mesh.isShell())
2105  {
2106  calculateFaceAreas(mesh);
2107  calculateFaceAreaMag(mesh);
2108  calculateFaceCentroids(mesh);
2109  }
2110  }
2111 
2112  for (int n=0; n<numMeshes; n++)
2113  {
2114  const Mesh& mesh = *_meshes[n];
2115  if (!mesh.isShell())
2116  calculateCellCentroids(mesh);
2117  }
2118 
2119  _coordField.syncLocal();
2120 
2121  for (int n=0; n<numMeshes; n++)
2122  {
2123  const Mesh& mesh = *_meshes[n];
2124  if (!mesh.isShell())
2125  calculateCellVolumes(mesh);
2126  }
2127 
2128  _volumeField.syncLocal();
2129 }
MeshMetricsCalculator::calculateCellVolumes
virtual void calculateCellVolumes(const Mesh &mesh)
Definition: MeshMetricsCalculator_impl.h:394
Mesh
Definition: Mesh.h:49
MeshMetricsCalculator::calculateFaceCentroids
virtual void calculateFaceCentroids(const Mesh &mesh)
Definition: MeshMetricsCalculator_impl.h:60
Model::_meshes
const MeshList _meshes
Definition: Model.h:29
MeshMetricsCalculator::calculateCellCentroids
virtual void calculateCellCentroids(const Mesh &mesh)
Definition: MeshMetricsCalculator_impl.h:130
MeshMetricsCalculator::calculateFaceAreas
virtual void calculateFaceAreas(const Mesh &mesh)
Definition: MeshMetricsCalculator_impl.h:240
MeshMetricsCalculator::calculateFaceAreaMag
virtual void calculateFaceAreaMag(const Mesh &mesh)
Definition: MeshMetricsCalculator_impl.h:373
Mesh::isShell
bool isShell() const
Definition: Mesh.h:323
Field::syncLocal
void syncLocal()
Definition: Field.cpp:334
template<class T >
void MeshMetricsCalculator< T >::updateTime ( )

Definition at line 2067 of file MeshMetricsCalculator_impl.h.

References Mesh::getCells(), and StorageSite::getCountLevel1().

2068 {
2069  if (_transient)
2070  {
2071  const int numMeshes = _meshes.size();
2072  for (int n=0; n<numMeshes; n++)
2073  {
2074  const Mesh& mesh = *_meshes[n];
2075  const StorageSite& cells = mesh.getCells();
2076  const int cellCount = cells.getCountLevel1();
2077  if (cellCount > 0)
2078  {
2079  IntArray& ibTypeN1 = dynamic_cast<IntArray&>(_geomFields.ibTypeN1[cells]);
2080  const IntArray& ibType =
2081  dynamic_cast<const IntArray& >(_geomFields.ibType[cells]);
2082  for(int c=0; c<cellCount; c++)
2083  ibTypeN1[c] = ibType[c];
2084 
2085  }
2086  }
2087  }
2088  else
2089  throw CException("MeshMetricsCalculator: not transient");
2090 }
GeomFields::ibTypeN1
Field ibTypeN1
Definition: GeomFields.h:39
Mesh
Definition: Mesh.h:49
GeomFields::ibType
Field ibType
Definition: GeomFields.h:38
Model::_meshes
const MeshList _meshes
Definition: Model.h:29
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
StorageSite::getCountLevel1
int getCountLevel1() const
Definition: StorageSite.h:72

Member Data Documentation

template<class T >
Field& MeshMetricsCalculator< T >::_areaField
private

Definition at line 64 of file MeshMetricsCalculator.h.

template<class T >
Field& MeshMetricsCalculator< T >::_areaMagField
private

Definition at line 65 of file MeshMetricsCalculator.h.

template<class T >
Field& MeshMetricsCalculator< T >::_boundaryNodeNormal
private

Definition at line 68 of file MeshMetricsCalculator.h.

template<class T >
Field& MeshMetricsCalculator< T >::_coordField
private

Definition at line 63 of file MeshMetricsCalculator.h.

template<class T >
GeomFields& MeshMetricsCalculator< T >::_geomFields
private

Definition at line 62 of file MeshMetricsCalculator.h.

template<class T >
Field& MeshMetricsCalculator< T >::_nodeDisplacement
private

Definition at line 67 of file MeshMetricsCalculator.h.

template<class T >
bool MeshMetricsCalculator< T >::_transient
private

Definition at line 69 of file MeshMetricsCalculator.h.

template<class T >
Field& MeshMetricsCalculator< T >::_volumeField
private

Definition at line 66 of file MeshMetricsCalculator.h.

The documentation for this class was generated from the following files: