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

#include <GradientModel.h>

Inheritance diagram for GradientModel< X >:
Inheritance graph
[legend]
Collaboration diagram for GradientModel< X >:
Collaboration graph
[legend]

Public Types

typedef NumTypeTraits< X >
::T_Scalar 		T_Scalar
 
typedef Array< T_Scalar > TArray
 
typedef Array< int > IntArray
 
typedef Vector< T_Scalar, 3 > VectorT3
 
typedef Array< VectorT3 > VectorT3Array
 
typedef Array< X > XArray
 
typedef Gradient< X > GradType
 
typedef Array< GradType > GradArray
 
typedef GradientMatrix< T_Scalar > GradMatrixType
 
typedef
GradMatrixType::PairWiseAssembler 		GradientMatrixAssembler
 

Public Member Functions

 GradientModel (const MeshList &meshes, const Field &varField, Field &gradientField, const GeomFields &geomFields)
 
virtual ~GradientModel ()
 
void init ()
 
void compute ()
 
- Public Member Functions inherited from GradientModelBase
 GradientModelBase (const MeshList &meshes)
 
- 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 ()
 

Static Public Member Functions

static shared_ptr
< GradientMatrixBase > 		getLeastSquaresGradientMatrix3D (const Mesh &mesh, const GeomFields &geomFields)
 
static shared_ptr
< GradientMatrixBase > 		getLeastSquaresGradientMatrix2D (const Mesh &mesh, const GeomFields &geomFields)
 
static GradMatrixType & getGradientMatrix (const Mesh &mesh, const GeomFields &geomFields)
 
- Static Public Member Functions inherited from GradientModelBase
static void clearGradientMatrix (const Mesh &mesh)
 

Private Attributes

const Field & _varField
 
Field & _gradientField
 
const GeomFields & _geomFields
 

Additional Inherited Members

- Protected Attributes inherited from Model
const MeshList _meshes
 
StorageSiteList _varSites
 
StorageSiteList _fluxSites
 
map< string, shared_ptr
< ArrayBase > > 		_persistenceData
 
- Static Protected Attributes inherited from GradientModelBase
static map< const Mesh
*, shared_ptr
< GradientMatrixBase > > 		_gradientMatricesMap
 

Detailed Description

template<class X>
class GradientModel< X >

Definition at line 106 of file GradientModel.h.

Member Typedef Documentation

template<class X>
typedef Array<GradType> GradientModel< X >::GradArray

Definition at line 121 of file GradientModel.h.

template<class X>
typedef GradMatrixType::PairWiseAssembler GradientModel< X >::GradientMatrixAssembler

Definition at line 124 of file GradientModel.h.

template<class X>
typedef GradientMatrix<T_Scalar> GradientModel< X >::GradMatrixType

Definition at line 123 of file GradientModel.h.

template<class X>
typedef Gradient<X> GradientModel< X >::GradType

Definition at line 119 of file GradientModel.h.

template<class X>
typedef Array<int> GradientModel< X >::IntArray

Definition at line 113 of file GradientModel.h.

template<class X>
typedef NumTypeTraits<X>::T_Scalar GradientModel< X >::T_Scalar

Definition at line 109 of file GradientModel.h.

template<class X>
typedef Array<T_Scalar> GradientModel< X >::TArray

Definition at line 111 of file GradientModel.h.

template<class X>
typedef Vector<T_Scalar,3> GradientModel< X >::VectorT3

Definition at line 115 of file GradientModel.h.

template<class X>
typedef Array<VectorT3> GradientModel< X >::VectorT3Array

Definition at line 116 of file GradientModel.h.

template<class X>
typedef Array<X> GradientModel< X >::XArray

Definition at line 118 of file GradientModel.h.

Constructor & Destructor Documentation

template<class X>
GradientModel< X >::GradientModel ( const MeshList &  meshes,
const Field &  varField,
Field &  gradientField,
const GeomFields &  geomFields 
)
inline

Definition at line 439 of file GradientModel.h.

441  :
442  GradientModelBase(meshes),
443  _varField(varField),
444  _gradientField(gradientField),
445  _geomFields(geomFields)
446  {
447  logCtor();
448  }
GradientModel::_geomFields
const GeomFields & _geomFields
Definition: GradientModel.h:607
logCtor
#define logCtor()
Definition: RLogInterface.h:26
GradientModel::_gradientField
Field & _gradientField
Definition: GradientModel.h:606
GradientModel::_varField
const Field & _varField
Definition: GradientModel.h:605
GradientModelBase::GradientModelBase
GradientModelBase(const MeshList &meshes)
Definition: GradientModel.h:95
template<class X>
virtual GradientModel< X >::~GradientModel ( )
inlinevirtual

Definition at line 450 of file GradientModel.h.

451  {}

Member Function Documentation

template<class X>
void GradientModel< X >::compute ( )
inline

Definition at line 472 of file GradientModel.h.

Referenced by SpeciesModel< T >::Impl::linearize(), ElectricModel< T >::Impl::linearizeElectroStatics(), BatteryModel< T >::Impl::linearizePC(), BatteryModel< T >::Impl::linearizePC_Thermal(), BatteryModel< T >::Impl::linearizePotential(), BatteryModel< T >::Impl::linearizeSpecies(), and BatteryModel< T >::Impl::linearizeThermal().

473  {
474 
475  const int numMeshes = _meshes.size();
476 
477  for (int n=0; n<numMeshes; n++)
478  {
479 
480  const Mesh& mesh = *_meshes[n];
481 
482  if (mesh.isShell() == false){
483  const StorageSite& cells = mesh.getCells();
484 
485  GradMatrixType& gradMatrix = getGradientMatrix(mesh,_geomFields);
486 
487  const XArray& var = dynamic_cast<const XArray&>(_varField[cells]);
488  shared_ptr<GradArray> gradPtr = gradMatrix.getGradient(var);
489 
490  _gradientField.addArray(cells,gradPtr);
491 
492  // fix values in cells adjacent to IB Faces
493 
494  const StorageSite& ibFaces = mesh.getIBFaces();
495  const int nIBFaces = ibFaces.getCount();
496  if (nIBFaces > 0)
497  {
498  const Array<int>& ibFaceList = mesh.getIBFaceList();
499 
500  const CRConnectivity& faceCells = mesh.getAllFaceCells();
501 
502  const IntArray& ibType =
503  dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
504 
505  GradientMatrixAssembler& assembler =
506  gradMatrix.getPairWiseAssembler(faceCells);
507 
508  const XArray& varIB =
509  dynamic_cast<const XArray&>(_varField[ibFaces]);
510 
511  for (int nf=0; nf<nIBFaces; nf++)
512  {
513  const int f = ibFaceList[nf];
514  const int c0 = faceCells(f,0);
515  const int c1 = faceCells(f,1);
516 
517  if (ibType[c0]==Mesh::IBTYPE_FLUID)
518  {
519  (*gradPtr)[c0].accumulate(assembler.getCoeff01(f),
520  varIB[nf]-var[c1]);
521  }
522  else
523  {
524  (*gradPtr)[c1].accumulate(assembler.getCoeff10(f),
525  varIB[nf]-var[c0]);
526  }
527  }
528  }
529 
530  // copy boundary values from adjacent cells
531 
532  foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
533  {
534  const FaceGroup& fg = *fgPtr;
535  const StorageSite& faces = fg.site;
536  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
537  const int faceCount = faces.getCount();
538  if ((fg.groupType!="interior") && (fg.groupType!="interface")
539  && (fg.groupType!="dielectric interface"))
540  {
541  if (fg.groupType == "symmetry")
542  {
543  const VectorT3Array& faceArea =
544  dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
545  const TArray& faceAreaMag =
546  dynamic_cast<const TArray&>(_geomFields.areaMag[faces]);
547  for(int f=0; f<faceCount; f++)
548  {
549  const int c0 = faceCells(f,0);
550  const int c1 = faceCells(f,1);
551  const VectorT3 en = faceArea[f]/faceAreaMag[f];
552  reflectGradient((*gradPtr)[c1], (*gradPtr)[c0], en);
553  }
554  }
555  else
556  {
557  for(int f=0; f<faceCount; f++)
558  {
559  const int c0 = faceCells(f,0);
560  const int c1 = faceCells(f,1);
561 
562  (*gradPtr)[c1] = (*gradPtr)[c0];
563  }
564  }
565  }
566  }
567  }
568  }
569 
570 
571  //create values
572  for (int n=0; n<numMeshes; n++ ){
573  const Mesh& mesh = *_meshes[n];
574  if (mesh.isShell() == false){
575  GradMatrixType& gradMatrix = getGradientMatrix(mesh,_geomFields);
576  gradMatrix.createScatterGatherValuesBuffer();
577  //partitioner interfaces
578  gradMatrix.syncValues();
579  }
580  }
581 
582 //skipping for multiple meshes, we have to fixe it,
583 #if 0
584  for (int n=0; n<numMeshes; n++ ){
585  const Mesh& mesh = *_meshes[n];
586  if (mesh.isShell() == false){
587  GradMatrixType& gradMatrix = getGradientMatrix(mesh,_geomFields);
588  //mesh interfaces
589  gradMatrix.recvScatterGatherValuesBufferLocal();
590  }
591  }
592 
593 #endif
594 
595  typedef map<const Mesh*, shared_ptr<GradientMatrixBase> > gradType;
596  foreach( const gradType::value_type& mpos, _gradientMatricesMap){
597  GradientMatrixBase& gMtrxBase = *mpos.second;
598  gMtrxBase.syncLocal();
599  }
600  _gradientField.syncLocal();
601 
602  }
FaceGroupPtr
shared_ptr< FaceGroup > FaceGroupPtr
Definition: Mesh.h:46
Mesh::getIBFaces
const StorageSite & getIBFaces() const
Definition: Mesh.h:111
GradientModel::TArray
Array< T_Scalar > TArray
Definition: GradientModel.h:111
FaceGroup
Definition: Mesh.h:28
Mesh::getAllFaceGroups
const FaceGroupList & getAllFaceGroups() const
Definition: Mesh.h:193
GradientMatrixBase::syncLocal
virtual void syncLocal()
Definition: GradientMatrix.h:26
reflectGradient
void reflectGradient(Gradient< T > &gr, const Gradient< T > &g0, const Vector< T, 3 > &en)
Definition: GradientModel.h:23
GradientModel::_geomFields
const GeomFields & _geomFields
Definition: GradientModel.h:607
Mesh::getIBFaceList
const Array< int > & getIBFaceList() const
Definition: Mesh.cpp:686
Mesh
Definition: Mesh.h:49
FaceGroup::groupType
string groupType
Definition: Mesh.h:42
GeomFields::ibType
Field ibType
Definition: GeomFields.h:38
Mesh::getAllFaceCells
const CRConnectivity & getAllFaceCells() const
Definition: Mesh.cpp:378
GradientModel::GradientMatrixAssembler
GradMatrixType::PairWiseAssembler GradientMatrixAssembler
Definition: GradientModel.h:124
GradientModel::_gradientField
Field & _gradientField
Definition: GradientModel.h:606
Model::_meshes
const MeshList _meshes
Definition: Model.h:29
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
GradientModelBase::_gradientMatricesMap
static map< const Mesh *, shared_ptr< GradientMatrixBase > > _gradientMatricesMap
Definition: GradientModel.h:101
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
GradientModel::getGradientMatrix
static GradMatrixType & getGradientMatrix(const Mesh &mesh, const GeomFields &geomFields)
Definition: GradientModel.h:454
Mesh::isShell
bool isShell() const
Definition: Mesh.h:323
GradientModel::XArray
Array< X > XArray
Definition: GradientModel.h:118
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
GradientModel::_varField
const Field & _varField
Definition: GradientModel.h:605
GeomFields::area
Field area
Definition: GeomFields.h:23
GradientModel::GradMatrixType
GradientMatrix< T_Scalar > GradMatrixType
Definition: GradientModel.h:123
GeomFields::areaMag
Field areaMag
Definition: GeomFields.h:25
Field::syncLocal
void syncLocal()
Definition: Field.cpp:334
GradientModel::IntArray
Array< int > IntArray
Definition: GradientModel.h:113
FaceGroup::site
StorageSite site
Definition: Mesh.h:40
template<class X>
static GradMatrixType& GradientModel< X >::getGradientMatrix ( const Mesh &  mesh,
const GeomFields &  geomFields 
)
inlinestatic

Definition at line 454 of file GradientModel.h.

Referenced by COMETBoundaryConditions< X, Diag, OffDiag >::applyPressureInletBC(), COMETBoundaryConditions< X, Diag, OffDiag >::applyRealWallBC(), COMETBoundaryConditions< X, Diag, OffDiag >::applyZeroGradientBC(), GradientModel< Vector< T, 3 > >::compute(), StructureSourceDiscretization< T, Diag, OffDiag >::discretizeFaces(), StructurePlasticDiscretization< T, Diag, OffDiag >::discretizeFaces(), and PlateSourceDiscretization< T, Diag, OffDiag >::discretizeFaces().

455  {
456  if (_gradientMatricesMap.find(&mesh) == _gradientMatricesMap.end())
457  {
458  if (mesh.getDimension() == 2)
459  {
460  _gradientMatricesMap[&mesh] = getLeastSquaresGradientMatrix2D(mesh,geomFields);
461  }
462  else
463  {
464  _gradientMatricesMap[&mesh] = getLeastSquaresGradientMatrix3D(mesh,geomFields);
465  }
466  }
467  return dynamic_cast<GradMatrixType&>(*_gradientMatricesMap[&mesh]);
468  }
GradientModel::getLeastSquaresGradientMatrix2D
static shared_ptr< GradientMatrixBase > getLeastSquaresGradientMatrix2D(const Mesh &mesh, const GeomFields &geomFields)
Definition: GradientModel.h:286
GradientModelBase::_gradientMatricesMap
static map< const Mesh *, shared_ptr< GradientMatrixBase > > _gradientMatricesMap
Definition: GradientModel.h:101
GradientModel::getLeastSquaresGradientMatrix3D
static shared_ptr< GradientMatrixBase > getLeastSquaresGradientMatrix3D(const Mesh &mesh, const GeomFields &geomFields)
Definition: GradientModel.h:128
GradientModel::GradMatrixType
GradientMatrix< T_Scalar > GradMatrixType
Definition: GradientModel.h:123
Mesh::getDimension
int getDimension() const
Definition: Mesh.h:105
template<class X>
static shared_ptr<GradientMatrixBase> GradientModel< X >::getLeastSquaresGradientMatrix2D ( const Mesh &  mesh,
const GeomFields &  geomFields 
)
inlinestatic

Definition at line 286 of file GradientModel.h.

Referenced by GradientModel< Vector< T, 3 > >::getGradientMatrix().

287  {
288  const StorageSite& cells = mesh.getCells();
289  const StorageSite& faces = mesh.getFaces();
290 
291  const CRConnectivity& faceCells = mesh.getAllFaceCells();
292 
293 
294  const int cellCount = cells.getSelfCount();
295  const int faceCount = faces.getSelfCount();
296 
297  GradMatrixType* gMPtr(new GradMatrixType(mesh));
298  GradMatrixType& gM = *gMPtr;
299  GradientMatrixAssembler& assembler = gM.getPairWiseAssembler(faceCells);
300 
301  const CRConnectivity& cellCells = gM.getConnectivity();
302 
303  VectorT3Array& coeffs = gM.getCoeffs();
304 
305  const VectorT3Array& cellCentroid =
306  dynamic_cast<const VectorT3Array&>(geomFields.coordinate[cells]);
307 
308  const VectorT3Array& faceCentroid =
309  dynamic_cast<const VectorT3Array&>(geomFields.coordinate[faces]);
310 
311  const VectorT3Array& faceArea =
312  dynamic_cast<const VectorT3Array& >(geomFields.area[faces]);
313 
314  const TArray& cellVolume =
315  dynamic_cast<const TArray&>(geomFields.volume[cells]);
316 
317  const IntArray& ibType = dynamic_cast<const IntArray&>(geomFields.ibType[cells]);
318 
319  const T_Scalar epsilon(1e-26);
320 
321  coeffs.zero();
322 
323  Array<bool> isDegenerate(cells.getCount());
324  isDegenerate = false;
325 
326  for(int f=0; f<faceCount; f++)
327  {
328  const int c0 = faceCells(f,0);
329  const int c1 = faceCells(f,1);
330  VectorT3 ds = cellCentroid[c1]-cellCentroid[c0];
331 
332  //fix the distance for ibfaces
333  if (((ibType[c0] == Mesh::IBTYPE_FLUID)
334  && (ibType[c1] == Mesh::IBTYPE_BOUNDARY)) ||
335  ((ibType[c1] == Mesh::IBTYPE_FLUID)
336  && (ibType[c0] == Mesh::IBTYPE_BOUNDARY)))
337  {
338  if (ibType[c0] == Mesh::IBTYPE_FLUID)
339  {
340  ds = faceCentroid[f]-cellCentroid[c0];
341  }
342  else
343  {
344  ds = cellCentroid[c1]-faceCentroid[f];
345  }
346  }
347  const T_Scalar dsMag = mag(ds);
348  assembler.getCoeff01(f)=ds/dsMag;
349  assembler.getCoeff10(f)=-ds/dsMag;
350  }
351 
352  const Array<int>& row = cellCells.getRow();
353  //const Array<int>& col = cellCells.getCol();
354 
355  for(int nc=0; nc<cellCount; nc++)
356  {
357  T_Scalar Ixx(0), Iyy(0);
358  T_Scalar Ixy(0);
359 
360  for(int inb=row[nc]; inb<row[nc+1]; inb++)
361  {
362  const VectorT3& ds = coeffs[inb];
363  Ixx += ds[0]*ds[0];
364  Iyy += ds[1]*ds[1];
365  Ixy += ds[0]*ds[1];
366  }
367 
368  const T_Scalar det = Ixx*Iyy-Ixy*Ixy;
369  //T_Scalar det = 0;
370  if (det > epsilon)
371  {
372  const T_Scalar Kxx = Iyy/det;
373  const T_Scalar Kxy = -Ixy/det;
374  const T_Scalar Kyy = Ixx/det;
375  for(int inb=row[nc]; inb<row[nc+1]; inb++)
376  {
377  // make a copy
378  VectorT3 ds(coeffs[inb]);
379 
380  //const int j = col[inb];
381  //T_Scalar dsMag = mag(cellCentroid[j]-cellCentroid[nc]);
382 
383  coeffs[inb][0] = (Kxx*ds[0] + Kxy*ds[1]);
384  coeffs[inb][1] = (Kxy*ds[0] + Kyy*ds[1]);
385  coeffs[inb][2] = 0;
386  }
387  }
388  else
389  {
390  isDegenerate[nc] = true;
391  }
392  }
393 
394  for(int f=0; f<faceCount; f++)
395  {
396  const int c0 = faceCells(f,0);
397  const int c1 = faceCells(f,1);
398  VectorT3 ds = cellCentroid[c1]-cellCentroid[c0];
399 
400  //fix the distance for ibfaces
401  if (((ibType[c0] == Mesh::IBTYPE_FLUID)
402  && (ibType[c1] == Mesh::IBTYPE_BOUNDARY)) ||
403  ((ibType[c1] == Mesh::IBTYPE_FLUID)
404  && (ibType[c0] == Mesh::IBTYPE_BOUNDARY)))
405  {
406  if (ibType[c0] == Mesh::IBTYPE_FLUID)
407  {
408  ds = faceCentroid[f]-cellCentroid[c0];
409  }
410  else
411  {
412  ds = cellCentroid[c1]-faceCentroid[f];
413  }
414  }
415  const T_Scalar dsMag = mag(ds);
416  assembler.getCoeff01(f) /= dsMag;
417  assembler.getCoeff10(f) /= dsMag;
418  }
419 
420  for(int f=0; f<faceCount; f++)
421  {
422  const int c0 = faceCells(f,0);
423  const int c1 = faceCells(f,1);
424  if (isDegenerate[c0])
425  {
426  assembler.getCoeff01(f)= T_Scalar(0.5)*faceArea[f]/cellVolume[c0];
427  }
428  if (isDegenerate[c1])
429  {
430  assembler.getCoeff10(f)= T_Scalar(-0.5)*faceArea[f]/cellVolume[c1];
431  }
432  }
433 
434 
435  return shared_ptr<GradientMatrixBase>(gMPtr);
436  }
epsilon
#define epsilon
Definition: Mesh.cpp:17
Array::zero
virtual void zero()
Definition: Array.h:281
CRConnectivity::getRow
const Array< int > & getRow() const
Definition: CRConnectivity.h:182
StorageSite::getSelfCount
int getSelfCount() const
Definition: StorageSite.h:40
GeomFields::coordinate
Field coordinate
Definition: GeomFields.h:19
GradientModel::TArray
Array< T_Scalar > TArray
Definition: GradientModel.h:111
mag
T mag(const Vector< T, 3 > &a)
Definition: Vector.h:260
GeomFields::ibType
Field ibType
Definition: GeomFields.h:38
Mesh::getAllFaceCells
const CRConnectivity & getAllFaceCells() const
Definition: Mesh.cpp:378
GradientModel::GradientMatrixAssembler
GradMatrixType::PairWiseAssembler GradientMatrixAssembler
Definition: GradientModel.h:124
Mesh::getFaces
const StorageSite & getFaces() const
Definition: Mesh.h:108
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
GeomFields::volume
Field volume
Definition: GeomFields.h:26
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
GeomFields::area
Field area
Definition: GeomFields.h:23
GradientModel::T_Scalar
NumTypeTraits< X >::T_Scalar T_Scalar
Definition: GradientModel.h:109
GradientModel::GradMatrixType
GradientMatrix< T_Scalar > GradMatrixType
Definition: GradientModel.h:123
GradientModel::IntArray
Array< int > IntArray
Definition: GradientModel.h:113
template<class X>
static shared_ptr<GradientMatrixBase> GradientModel< X >::getLeastSquaresGradientMatrix3D ( const Mesh &  mesh,
const GeomFields &  geomFields 
)
inlinestatic

Definition at line 128 of file GradientModel.h.

Referenced by GradientModel< Vector< T, 3 > >::getGradientMatrix().

129  {
130  const StorageSite& cells = mesh.getCells();
131  const StorageSite& faces = mesh.getFaces();
132 
133  const CRConnectivity& faceCells = mesh.getAllFaceCells();
134 
135  const int cellCount = cells.getSelfCount();
136  const int faceCount = faces.getSelfCount();
137  GradMatrixType* gMPtr(new GradMatrixType(mesh));
138  GradMatrixType& gM = *gMPtr;
139  GradientMatrixAssembler& assembler = gM.getPairWiseAssembler(faceCells);
140 
141  const CRConnectivity& cellCells = gM.getConnectivity();
142 
143  VectorT3Array& coeffs = gM.getCoeffs();
144 
145  const VectorT3Array& cellCentroid =
146  dynamic_cast<const VectorT3Array&>(geomFields.coordinate[cells]);
147 
148  const VectorT3Array& faceCentroid =
149  dynamic_cast<const VectorT3Array&>(geomFields.coordinate[faces]);
150 
151  const VectorT3Array& faceArea =
152  dynamic_cast<const VectorT3Array& >(geomFields.area[faces]);
153 
154  const TArray& cellVolume =
155  dynamic_cast<const TArray&>(geomFields.volume[cells]);
156 
157  const T_Scalar epsilon(1e-6);
158 
159  const IntArray& ibType = dynamic_cast<const IntArray&>(geomFields.ibType[cells]);
160 
161  coeffs.zero();
162 
163  Array<bool> isDegenerate(cells.getCount());
164  isDegenerate = false;
165 
166  for(int f=0; f<faceCount; f++)
167  {
168  const int c0 = faceCells(f,0);
169  const int c1 = faceCells(f,1);
170  VectorT3 ds = cellCentroid[c1]-cellCentroid[c0];
171 
172  //fix the distance for ibfaces
173  if (((ibType[c0] == Mesh::IBTYPE_FLUID)
174  && (ibType[c1] == Mesh::IBTYPE_BOUNDARY)) ||
175  ((ibType[c1] == Mesh::IBTYPE_FLUID)
176  && (ibType[c0] == Mesh::IBTYPE_BOUNDARY)))
177  {
178  if (ibType[c0] == Mesh::IBTYPE_FLUID)
179  {
180  ds = faceCentroid[f]-cellCentroid[c0];
181  }
182  else
183  {
184  ds = cellCentroid[c1]-faceCentroid[f];
185  }
186  }
187  const T_Scalar dsMag = mag(ds);
188  assembler.getCoeff01(f)=ds/dsMag;
189  assembler.getCoeff10(f)=-ds/dsMag;
190  }
191 
192  const Array<int>& row = cellCells.getRow();
193  //const Array<int>& col = cellCells.getCol();
194 
195  for(int nc=0; nc<cellCount; nc++)
196  {
197  T_Scalar Ixx(0), Iyy(0), Izz(0);
198  T_Scalar Ixy(0), Ixz(0), Iyz(0);
199 
200  for(int inb=row[nc]; inb<row[nc+1]; inb++)
201  {
202  const VectorT3& ds = coeffs[inb];
203  Ixx += ds[0]*ds[0];
204  Iyy += ds[1]*ds[1];
205  Izz += ds[2]*ds[2];
206  Ixy += ds[0]*ds[1];
207  Ixz += ds[0]*ds[2];
208  Iyz += ds[1]*ds[2];
209  }
210 
211 
212  const T_Scalar det = Ixx*(Iyy*Izz-Iyz*Iyz) -Ixy*(Ixy*Izz-Iyz*Ixz)
213  + Ixz*(Ixy*Iyz-Iyy*Ixz);
214  //T_Scalar det = NumTypeTraits<T_Scalar>::getZero();
215  if (det > epsilon)
216  {
217  const T_Scalar Kxx = (Iyy*Izz-Iyz*Iyz)/det;
218  const T_Scalar Kxy = -(Ixy*Izz-Iyz*Ixz)/det;
219  const T_Scalar Kxz = (Ixy*Iyz-Iyy*Ixz)/det;
220  const T_Scalar Kyy = (Ixx*Izz-Ixz*Ixz)/det;
221  const T_Scalar Kyz = -(Ixx*Iyz-Ixy*Ixz)/det;
222  const T_Scalar Kzz = (Ixx*Iyy-Ixy*Ixy)/det;
223  for(int inb=row[nc]; inb<row[nc+1]; inb++)
224  {
225  // make a copy
226  VectorT3 ds(coeffs[inb]);
227  //const int j = col[inb];
228  //T_Scalar dsMag = mag(cellCentroid[j]-cellCentroid[nc]);
229  coeffs[inb][0] = (Kxx*ds[0] + Kxy*ds[1] + Kxz*ds[2]);
230  coeffs[inb][1] = (Kxy*ds[0] + Kyy*ds[1] + Kyz*ds[2]);
231  coeffs[inb][2] = (Kxz*ds[0] + Kyz*ds[1] + Kzz*ds[2]);
232  }
233  }
234  else
235  {
236  isDegenerate[nc] = true;
237  }
238  }
239 
240  for(int f=0; f<faceCount; f++)
241  {
242  const int c0 = faceCells(f,0);
243  const int c1 = faceCells(f,1);
244  VectorT3 ds = cellCentroid[c1]-cellCentroid[c0];
245 
246  //fix the distance for ibfaces
247  if (((ibType[c0] == Mesh::IBTYPE_FLUID)
248  && (ibType[c1] == Mesh::IBTYPE_BOUNDARY)) ||
249  ((ibType[c1] == Mesh::IBTYPE_FLUID)
250  && (ibType[c0] == Mesh::IBTYPE_BOUNDARY)))
251  {
252  if (ibType[c0] == Mesh::IBTYPE_FLUID)
253  {
254  ds = faceCentroid[f]-cellCentroid[c0];
255  }
256  else
257  {
258  ds = cellCentroid[c1]-faceCentroid[f];
259  }
260  }
261  const T_Scalar dsMag = mag(ds);
262  assembler.getCoeff01(f) /= dsMag;
263  assembler.getCoeff10(f) /= dsMag;
264  }
265 
266 
267  for(int f=0; f<faceCount; f++)
268  {
269  const int c0 = faceCells(f,0);
270  const int c1 = faceCells(f,1);
271  if (isDegenerate[c0])
272  {
273  assembler.getCoeff01(f)= T_Scalar(0.5)*faceArea[f]/cellVolume[c0];
274  }
275  if (isDegenerate[c1])
276  {
277  assembler.getCoeff10(f)= T_Scalar(-0.5)*faceArea[f]/cellVolume[c1];
278  }
279  }
280 
281  return shared_ptr<GradientMatrixBase>(gMPtr);
282  }
epsilon
#define epsilon
Definition: Mesh.cpp:17
Array::zero
virtual void zero()
Definition: Array.h:281
CRConnectivity::getRow
const Array< int > & getRow() const
Definition: CRConnectivity.h:182
StorageSite::getSelfCount
int getSelfCount() const
Definition: StorageSite.h:40
GeomFields::coordinate
Field coordinate
Definition: GeomFields.h:19
GradientModel::TArray
Array< T_Scalar > TArray
Definition: GradientModel.h:111
mag
T mag(const Vector< T, 3 > &a)
Definition: Vector.h:260
GeomFields::ibType
Field ibType
Definition: GeomFields.h:38
Mesh::getAllFaceCells
const CRConnectivity & getAllFaceCells() const
Definition: Mesh.cpp:378
GradientModel::GradientMatrixAssembler
GradMatrixType::PairWiseAssembler GradientMatrixAssembler
Definition: GradientModel.h:124
Mesh::getFaces
const StorageSite & getFaces() const
Definition: Mesh.h:108
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
GeomFields::volume
Field volume
Definition: GeomFields.h:26
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
GeomFields::area
Field area
Definition: GeomFields.h:23
GradientModel::T_Scalar
NumTypeTraits< X >::T_Scalar T_Scalar
Definition: GradientModel.h:109
GradientModel::GradMatrixType
GradientMatrix< T_Scalar > GradMatrixType
Definition: GradientModel.h:123
GradientModel::IntArray
Array< int > IntArray
Definition: GradientModel.h:113
template<class X>
void GradientModel< X >::init ( )
inlinevirtual

Implements Model.

Definition at line 470 of file GradientModel.h.

470 {}

Member Data Documentation

template<class X>
const GeomFields& GradientModel< X >::_geomFields
private

Definition at line 607 of file GradientModel.h.

Referenced by GradientModel< Vector< T, 3 > >::compute().

template<class X>
Field& GradientModel< X >::_gradientField
private

Definition at line 606 of file GradientModel.h.

Referenced by GradientModel< Vector< T, 3 > >::compute().

template<class X>
const Field& GradientModel< X >::_varField
private

Definition at line 605 of file GradientModel.h.

Referenced by GradientModel< Vector< T, 3 > >::compute().

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