Memosa-FVM  0.2
PlateSourceDiscretization< T, Diag, OffDiag > Class Template Reference

#include <PlateSourceDiscretization.h>

Inheritance diagram for PlateSourceDiscretization< T, Diag, OffDiag >:
Collaboration diagram for PlateSourceDiscretization< T, Diag, OffDiag >:

Public Types

typedef Array< T > TArray
 
typedef Vector< T, 3 > VectorT3
 
typedef Array< VectorT3VectorT3Array
 
typedef Vector< T, 4 > VectorT4
 
typedef Array< VectorT4VectorT4Array
 
typedef Gradient< VectorT3VGradType
 
typedef Array< Gradient
< VectorT3 > > 
VGradArray
 
typedef CRMatrix< Diag,
OffDiag, VectorT3
CCMatrix
 
typedef CCMatrix::DiagArray DiagArray
 
typedef CCMatrix::PairWiseAssembler CCAssembler
 
typedef GradientModel< VectorT3VGradModelType
 
typedef
VGradModelType::GradMatrixType 
VGradMatrix
 

Public Member Functions

 PlateSourceDiscretization (const MeshList &meshes, const GeomFields &geomFields, Field &varField, const Field &ymField, const Field &nuField, const Field &varGradientField, const Field &residualStressField, const Field &thicknessField, const Field &forceField, const T &scf, const Field &devStressField, const Field &VMStressField, Field &plasticStrainField, Field &plasticStrainOutField, Field &plasticStrainN1Field, Field &plasticMomentField, const T &A, const T &B, const T &mm, const T &nn, const T &Sy0, const int &nz, const T &timeStep, const int &creepModel, const bool &creep, bool fullLinearization=true)
 
void discretize (const Mesh &mesh, MultiFieldMatrix &mfmatrix, MultiField &xField, MultiField &rField)
 
void discretizeFaces (const Mesh &mesh, const StorageSite &faces, MultiFieldMatrix &mfmatrix, MultiField &xField, MultiField &rField, const bool isBoundary, const bool isSymmetry)
 
- Public Member Functions inherited from Discretization
 Discretization (const MeshList &meshes)
 
virtual ~Discretization ()
 
 DEFINE_TYPENAME ("Discretization")
 

Private Attributes

const GeomFields_geomFields
 
Field_varField
 
const Field_ymField
 
const Field_nuField
 
const Field_varGradientField
 
const Field_residualStressField
 
const Field_thicknessField
 
const Field_forceField
 
const T & _scf
 
const Field_devStressField
 
const Field_VMStressField
 
Field_plasticStrainField
 
Field_plasticStrainOutField
 
Field_plasticStrainN1Field
 
Field_plasticMomentField
 
const T _A
 
const T _B
 
const T _m
 
const T _n
 
const T _Sy0
 
const int _nz
 
const T _timeStep
 
const int _creepModel
 
const bool _creep
 
const bool _fullLinearization
 

Additional Inherited Members

- Protected Attributes inherited from Discretization
const MeshList_meshes
 

Detailed Description

template<class T, class Diag, class OffDiag>
class PlateSourceDiscretization< T, Diag, OffDiag >

Definition at line 21 of file PlateSourceDiscretization.h.

Member Typedef Documentation

template<class T , class Diag , class OffDiag >
typedef CCMatrix::PairWiseAssembler PlateSourceDiscretization< T, Diag, OffDiag >::CCAssembler

Definition at line 35 of file PlateSourceDiscretization.h.

template<class T , class Diag , class OffDiag >
typedef CRMatrix<Diag,OffDiag,VectorT3> PlateSourceDiscretization< T, Diag, OffDiag >::CCMatrix

Definition at line 33 of file PlateSourceDiscretization.h.

template<class T , class Diag , class OffDiag >
typedef CCMatrix::DiagArray PlateSourceDiscretization< T, Diag, OffDiag >::DiagArray

Definition at line 34 of file PlateSourceDiscretization.h.

template<class T , class Diag , class OffDiag >
typedef Array<T> PlateSourceDiscretization< T, Diag, OffDiag >::TArray

Definition at line 25 of file PlateSourceDiscretization.h.

template<class T , class Diag , class OffDiag >
typedef Vector<T,3> PlateSourceDiscretization< T, Diag, OffDiag >::VectorT3

Definition at line 26 of file PlateSourceDiscretization.h.

template<class T , class Diag , class OffDiag >
typedef Array<VectorT3> PlateSourceDiscretization< T, Diag, OffDiag >::VectorT3Array

Definition at line 27 of file PlateSourceDiscretization.h.

template<class T , class Diag , class OffDiag >
typedef Vector<T,4> PlateSourceDiscretization< T, Diag, OffDiag >::VectorT4

Definition at line 28 of file PlateSourceDiscretization.h.

template<class T , class Diag , class OffDiag >
typedef Array<VectorT4> PlateSourceDiscretization< T, Diag, OffDiag >::VectorT4Array

Definition at line 29 of file PlateSourceDiscretization.h.

template<class T , class Diag , class OffDiag >
typedef Array<Gradient<VectorT3> > PlateSourceDiscretization< T, Diag, OffDiag >::VGradArray

Definition at line 31 of file PlateSourceDiscretization.h.

template<class T , class Diag , class OffDiag >
typedef VGradModelType::GradMatrixType PlateSourceDiscretization< T, Diag, OffDiag >::VGradMatrix

Definition at line 38 of file PlateSourceDiscretization.h.

template<class T , class Diag , class OffDiag >
typedef GradientModel<VectorT3> PlateSourceDiscretization< T, Diag, OffDiag >::VGradModelType

Definition at line 37 of file PlateSourceDiscretization.h.

template<class T , class Diag , class OffDiag >
typedef Gradient<VectorT3> PlateSourceDiscretization< T, Diag, OffDiag >::VGradType

Definition at line 30 of file PlateSourceDiscretization.h.

Constructor & Destructor Documentation

template<class T , class Diag , class OffDiag >
PlateSourceDiscretization< T, Diag, OffDiag >::PlateSourceDiscretization ( const MeshList meshes,
const GeomFields geomFields,
Field varField,
const Field ymField,
const Field nuField,
const Field varGradientField,
const Field residualStressField,
const Field thicknessField,
const Field forceField,
const T &  scf,
const Field devStressField,
const Field VMStressField,
Field plasticStrainField,
Field plasticStrainOutField,
Field plasticStrainN1Field,
Field plasticMomentField,
const T &  A,
const T &  B,
const T &  mm,
const T &  nn,
const T &  Sy0,
const int &  nz,
const T &  timeStep,
const int &  creepModel,
const bool &  creep,
bool  fullLinearization = true 
)
inline

Definition at line 40 of file PlateSourceDiscretization.h.

65  :
66  Discretization(meshes),
67  _geomFields(geomFields),
68  _varField(varField),
69  _ymField(ymField),
70  _nuField(nuField),
71  _varGradientField(varGradientField),
72  _residualStressField(residualStressField),
73  _thicknessField(thicknessField),
74  _forceField(forceField),
75  _scf(scf),
76  _devStressField(devStressField),
77  _VMStressField(VMStressField),
78  _plasticStrainField(plasticStrainField),
79  _plasticStrainOutField(plasticStrainOutField),
80  _plasticStrainN1Field(plasticStrainN1Field),
81  _plasticMomentField(plasticMomentField),
82  _A(A),
83  _B(B),
84  _m(mm),
85  _n(nn),
86  _Sy0(Sy0),
87  _nz(nz),
88  _timeStep(timeStep),
89  _creepModel(creepModel),
90  _creep(creep),
91  _fullLinearization(fullLinearization)
92  {}
Discretization(const MeshList &meshes)

Member Function Documentation

template<class T , class Diag , class OffDiag >
void PlateSourceDiscretization< T, Diag, OffDiag >::discretize ( const Mesh mesh,
MultiFieldMatrix mfmatrix,
MultiField xField,
MultiField rField 
)
inlinevirtual

Implements Discretization.

Definition at line 94 of file PlateSourceDiscretization.h.

References PlateSourceDiscretization< T, Diag, OffDiag >::_A, PlateSourceDiscretization< T, Diag, OffDiag >::_B, PlateSourceDiscretization< T, Diag, OffDiag >::_creep, PlateSourceDiscretization< T, Diag, OffDiag >::_creepModel, PlateSourceDiscretization< T, Diag, OffDiag >::_devStressField, PlateSourceDiscretization< T, Diag, OffDiag >::_forceField, PlateSourceDiscretization< T, Diag, OffDiag >::_geomFields, PlateSourceDiscretization< T, Diag, OffDiag >::_m, PlateSourceDiscretization< T, Diag, OffDiag >::_n, PlateSourceDiscretization< T, Diag, OffDiag >::_nuField, PlateSourceDiscretization< T, Diag, OffDiag >::_nz, PlateSourceDiscretization< T, Diag, OffDiag >::_plasticMomentField, PlateSourceDiscretization< T, Diag, OffDiag >::_plasticStrainField, PlateSourceDiscretization< T, Diag, OffDiag >::_plasticStrainN1Field, PlateSourceDiscretization< T, Diag, OffDiag >::_plasticStrainOutField, PlateSourceDiscretization< T, Diag, OffDiag >::_residualStressField, PlateSourceDiscretization< T, Diag, OffDiag >::_Sy0, PlateSourceDiscretization< T, Diag, OffDiag >::_thicknessField, PlateSourceDiscretization< T, Diag, OffDiag >::_timeStep, PlateSourceDiscretization< T, Diag, OffDiag >::_varField, PlateSourceDiscretization< T, Diag, OffDiag >::_VMStressField, PlateSourceDiscretization< T, Diag, OffDiag >::_ymField, PlateSourceDiscretization< T, Diag, OffDiag >::discretizeFaces(), Mesh::getAllFaceGroups(), Mesh::getCells(), StorageSite::getCountLevel1(), CRMatrix< T_Diag, T_OffDiag, X >::getDiag(), Mesh::getInteriorFaceGroup(), MultiFieldMatrix::getMatrix(), StorageSite::getSelfCount(), FaceGroup::groupType, FaceGroup::site, sqrt(), and GeomFields::volume.

96  {
97  const StorageSite& cells = mesh.getCells();
98 
99  const MultiField::ArrayIndex cVarIndex(&_varField,&cells);
100  CCMatrix& matrix =
101  dynamic_cast<CCMatrix&>(mfmatrix.getMatrix(cVarIndex,cVarIndex));
102 
103  DiagArray& diag = matrix.getDiag();
104 
105  const TArray& cellVolume =
106  dynamic_cast<const TArray&>(_geomFields.volume[cells]);
107 
108  const TArray& forceCell =
109  dynamic_cast<const TArray&>(_forceField[cells]);
110 
111  VectorT3Array& rCell = dynamic_cast<VectorT3Array&>(rField[cVarIndex]);
112  const VectorT3Array& xCell = dynamic_cast<const VectorT3Array&>(xField[cVarIndex]);
113 
114  const VGradArray& residualCell =
115  dynamic_cast<const VGradArray&>(_residualStressField[cells]);
116 
117  const TArray& ymCell =
118  dynamic_cast<const TArray&>(_ymField[cells]);
119 
120  const TArray& nuCell =
121  dynamic_cast<const TArray&>(_nuField[cells]);
122 
123  const TArray& tCell =
124  dynamic_cast<const TArray&>(_thicknessField[cells]);
125 
126  const VectorT4Array& devStress =
127  dynamic_cast<const VectorT4Array&>(_devStressField[cells]);
128 
129  const TArray& VMStress =
130  dynamic_cast<const TArray&>(_VMStressField[cells]);
131 
132  VectorT4Array& plasticStrain =
133  dynamic_cast<VectorT4Array&>(_plasticStrainField[cells]);
134 
135  VectorT3Array& plasticStrainOut =
136  dynamic_cast<VectorT3Array&>(_plasticStrainOutField[cells]);
137 
138  VectorT4Array& plasticStrainN1 =
139  dynamic_cast<VectorT4Array&>(_plasticStrainN1Field[cells]);
140 
141  VectorT3Array& plasticMoment =
142  dynamic_cast<VectorT3Array&>(_plasticMomentField[cells]);
143 
144  for(int c=0; c<cells.getSelfCount(); c++)
145  {
146  rCell[c][2]-=forceCell[c]*cellVolume[c];
147  }
148 
149  //residual stress
150  for(int c=0; c<cells.getSelfCount(); c++)
151  {
152  const T HSigmazzV = tCell[c]*(residualCell[c])[2][2]*cellVolume[c];
153  rCell[c][0] += HSigmazzV*xCell[c][0];
154  (diag[c])(0,0) += HSigmazzV;
155  rCell[c][1] += HSigmazzV*xCell[c][1];
156  (diag[c])(1,1) += HSigmazzV;
157  }
158  //
159 
160  const int nCells = cells.getCountLevel1();
161  const T zero(0.0);
162  const T half(0.5);
163  const T one(1.0);
164  const T two(2.0);
165  const T three(3.0);
166  const T four(4.0);
167  const T six(6.0);
168 
169  if(_creep)
170  {
171  if (_creepModel==1)
172  {
173  for(int n=0; n<nCells; n++)
174  {
175  int nn = n*(_nz+1);
176  for(int k=0; k<=_nz; k++)
177  {
178 
179  T VMPlasticStrain = sqrt(half*
180  (pow(plasticStrain[nn+k][0]-plasticStrain[nn+k][1],2.0) +
181  pow(plasticStrain[nn+k][1]-plasticStrain[nn+k][2],2.0) +
182  pow(plasticStrain[nn+k][2]-plasticStrain[nn+k][0],2.0) +
183  six*(pow(plasticStrain[nn+k][3],2.0))));
184 
185  //T VMPlasticStrain = sqrt(pow(plasticStrain[nn+k][0],2.0));
186  T Sy = _Sy0*(one + _B*pow(VMPlasticStrain,_n));
187  T mult = _A*(pow((VMStress[nn+k]/Sy),_m))/VMStress[nn+k];
188  if(k==_nz/2)
189  mult = zero;
190  for(int i=0;i<4;i++)
191  plasticStrain[nn+k][i] = plasticStrainN1[nn+k][i]+mult*devStress[nn+k][i]*_timeStep;
192  }
193  plasticStrainOut[n][0] = plasticStrain[nn+_nz][0];
194  plasticStrainOut[n][1] = plasticStrain[nn+_nz][1];
195  plasticStrainOut[n][2] = plasticStrain[nn+_nz][3];
196  }
197  }
198 
199  for(int n=0; n<nCells; n++)
200  {
201  T var1 = ymCell[n]/(one - pow(nuCell[n],two));
202  T var2 = one - nuCell[n];
203  T var3 = (tCell[n]/T(_nz))/three;
204  T tHalf(0.0);
205  T tempxx(0.0);
206  T tempyy(0.0);
207  T tempxy(0.0);
208  tHalf = tCell[n]/two;
209  int nn = n*(_nz+1);
210 
211  tempxx += (-tHalf)*(plasticStrain[nn][0] + nuCell[n]*plasticStrain[nn][1]);
212  tempyy += (-tHalf)*(plasticStrain[nn][1] + nuCell[n]*plasticStrain[nn][0]);
213  tempxy += (-tHalf)*var2*plasticStrain[nn][3];
214  for(int k=1; k<_nz; k++)
215  {
216  T n1 = four;
217  if((k%2)==0)n1 = two;
218  T zz(tCell[n]*(T(k)-(T(_nz)/T(2)))/T(_nz));
219  tempxx += T(n1)*zz*(plasticStrain[nn+k][0] + nuCell[n]*plasticStrain[nn+k][1]);
220  tempyy += T(n1)*zz*(plasticStrain[nn+k][1] + nuCell[n]*plasticStrain[nn+k][0]);
221  tempxy += T(n1)*zz*var2*plasticStrain[nn+k][3];
222  }
223  tempxx += (tHalf)*(plasticStrain[nn+_nz][0] + nuCell[n]*plasticStrain[nn+_nz][1]);
224  tempyy += (tHalf)*(plasticStrain[nn+_nz][1] + nuCell[n]*plasticStrain[nn+_nz][0]);
225  tempxy += (tHalf)*var2*plasticStrain[nn+_nz][3];
226 
227  plasticMoment[n][0] = var1*var3*tempxx;
228  plasticMoment[n][1] = var1*var3*tempyy;
229  plasticMoment[n][2] = var1*var3*tempxy;
230 
231  }
232  }
233 
234  const StorageSite& iFaces = mesh.getInteriorFaceGroup().site;
235 
236  discretizeFaces(mesh, iFaces, mfmatrix, xField, rField, false, false);
237 
238  /*
239  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
240  {
241  const FaceGroup& fg = *fgPtr;
242  const StorageSite& faces = fg.site;
243  discretizeFaces(mesh, faces, mfmatrix, xField, rField, false, false);
244  }
245  */
246 
247  // boundaries
248  foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
249  {
250  const FaceGroup& fg = *fgPtr;
251  const StorageSite& faces = fg.site;
252  if (fg.groupType!="interior")
253  {
254  discretizeFaces(mesh, faces, mfmatrix, xField, rField,
255  fg.groupType!="interface",
256  fg.groupType=="symmetry");
257  }
258  }
259  }
const FaceGroup & getInteriorFaceGroup() const
Definition: Mesh.h:181
Matrix & getMatrix(const Index &rowIndex, const Index &colIndex)
shared_ptr< FaceGroup > FaceGroupPtr
Definition: Mesh.h:46
int getSelfCount() const
Definition: StorageSite.h:40
Definition: Mesh.h:28
const FaceGroupList & getAllFaceGroups() const
Definition: Mesh.h:193
CRMatrix< Diag, OffDiag, VectorT3 > CCMatrix
void discretizeFaces(const Mesh &mesh, const StorageSite &faces, MultiFieldMatrix &mfmatrix, MultiField &xField, MultiField &rField, const bool isBoundary, const bool isSymmetry)
string groupType
Definition: Mesh.h:42
Tangent sqrt(const Tangent &a)
Definition: Tangent.h:317
Array< Gradient< VectorT3 > > VGradArray
pair< const Field *, const StorageSite * > ArrayIndex
Definition: MultiField.h:21
const StorageSite & getCells() const
Definition: Mesh.h:109
Field volume
Definition: GeomFields.h:26
int getCountLevel1() const
Definition: StorageSite.h:72
StorageSite site
Definition: Mesh.h:40
template<class T , class Diag , class OffDiag >
void PlateSourceDiscretization< T, Diag, OffDiag >::discretizeFaces ( const Mesh mesh,
const StorageSite faces,
MultiFieldMatrix mfmatrix,
MultiField xField,
MultiField rField,
const bool  isBoundary,
const bool  isSymmetry 
)
inline

Definition at line 262 of file PlateSourceDiscretization.h.

References PlateSourceDiscretization< T, Diag, OffDiag >::_creep, PlateSourceDiscretization< T, Diag, OffDiag >::_fullLinearization, PlateSourceDiscretization< T, Diag, OffDiag >::_geomFields, PlateSourceDiscretization< T, Diag, OffDiag >::_nuField, PlateSourceDiscretization< T, Diag, OffDiag >::_plasticMomentField, PlateSourceDiscretization< T, Diag, OffDiag >::_residualStressField, PlateSourceDiscretization< T, Diag, OffDiag >::_scf, PlateSourceDiscretization< T, Diag, OffDiag >::_thicknessField, PlateSourceDiscretization< T, Diag, OffDiag >::_varField, PlateSourceDiscretization< T, Diag, OffDiag >::_varGradientField, PlateSourceDiscretization< T, Diag, OffDiag >::_ymField, GeomFields::area, GeomFields::areaMag, GeomFields::coordinate, dot(), Mesh::getCellCells(), Mesh::getCells(), GradientMatrix< T_Scalar >::getCoeff(), CRMatrix< T_Diag, T_OffDiag, X >::getCoeff(), CRMatrix< T_Diag, T_OffDiag, X >::PairWiseAssembler::getCoeff01(), CRMatrix< T_Diag, T_OffDiag, X >::PairWiseAssembler::getCoeff10(), CRConnectivity::getCol(), StorageSite::getCount(), CRMatrix< T_Diag, T_OffDiag, X >::getDiag(), Mesh::getFaceCells(), GradientModel< X >::getGradientMatrix(), MultiFieldMatrix::getMatrix(), CRMatrix< T_Diag, T_OffDiag, X >::getPairWiseAssembler(), CRConnectivity::getRow(), StorageSite::getSelfCount(), harmonicAverage(), and GeomFields::volume.

Referenced by PlateSourceDiscretization< T, Diag, OffDiag >::discretize().

266  {
267  const StorageSite& cells = mesh.getCells();
268 
269  const MultiField::ArrayIndex cVarIndex(&_varField,&cells);
270 
271  const VectorT3Array& faceArea =
272  dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
273 
274  const TArray& faceAreaMag =
275  dynamic_cast<const TArray&>(_geomFields.areaMag[faces]);
276 
277  const VectorT3Array& cellCentroid =
278  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[cells]);
279 
280  const VectorT3Array& faceCentroid =
281  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[faces]);
282 
283  const TArray& cellVolume =
284  dynamic_cast<const TArray&>(_geomFields.volume[cells]);
285 
286  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
287 
288  VectorT3Array& rCell = dynamic_cast<VectorT3Array&>(rField[cVarIndex]);
289  const VectorT3Array& xCell = dynamic_cast<const VectorT3Array&>(xField[cVarIndex]);
290 
291  const VGradArray& vGradCell =
292  dynamic_cast<const VGradArray&>(_varGradientField[cells]);
293 
294  const VGradArray& residualCell =
295  dynamic_cast<const VGradArray&>(_residualStressField[cells]);
296 
297  const TArray& ymCell =
298  dynamic_cast<const TArray&>(_ymField[cells]);
299 
300  const TArray& nuCell =
301  dynamic_cast<const TArray&>(_nuField[cells]);
302 
303  const TArray& tCell =
304  dynamic_cast<const TArray&>(_thicknessField[cells]);
305 
306  VectorT3Array& plasticMoment =
307  dynamic_cast<VectorT3Array&>(_plasticMomentField[cells]);
308 
309  CCMatrix& matrix = dynamic_cast<CCMatrix&>(mfmatrix.getMatrix(cVarIndex,
310  cVarIndex));
311  CCAssembler& assembler = matrix.getPairWiseAssembler(faceCells);
312  DiagArray& diag = matrix.getDiag();
313 
314  const int nFaces = faces.getCount();
315 
317  const CRConnectivity& cellCells = mesh.getCellCells();
318  const Array<int>& ccRow = cellCells.getRow();
319  const Array<int>& ccCol = cellCells.getCol();
320 
321  const int nInteriorCells = cells.getSelfCount();
322 
323  const T zero(0.0);
324  const T one(1.0);
325  const T two(2.0);
326  const T three(3.0);
327  const T twelve(12.0);
328 
329  for(int f=0; f<nFaces; f++)
330  {
331  const int c0 = faceCells(f,0);
332  const int c1 = faceCells(f,1);
333 
334  const VectorT3& Af = faceArea[f];
335  const VectorT3 en = Af/faceAreaMag[f];
336 
337  VectorT3 ds=cellCentroid[c1]-cellCentroid[c0];
338  VectorT3 dzeta0 = faceCentroid[f]-cellCentroid[c0];
339  VectorT3 dzeta1 = faceCentroid[f]-cellCentroid[c1];
340 
341  const T diffMetric = faceAreaMag[f]*faceAreaMag[f]/dot(faceArea[f],ds);
342  const VectorT3 secondaryCoeff = (faceArea[f]-ds*diffMetric);
343 
344  const T dfx0(dzeta0[0]);
345  const T dfy0(dzeta0[1]);
346 
347  const T dfx1(dzeta1[0]);
348  const T dfy1(dzeta1[1]);
349 
350  T vol0 = cellVolume[c0];
351  T vol1 = cellVolume[c1];
352 
353  T wt0 = vol0/(vol0+vol1);
354  T wt1 = vol1/(vol0+vol1);
355 
356  T bwt0(wt0);
357  T bwt1(wt1);
358  if (isBoundary && !isSymmetry)
359  {
360  wt0 = T(1.0);
361  wt1 = T(0.);
362  //bwt0 = T(0.);
363  // bwt1 = T(1.);
364  }
365 
366  T faceD(1.0);
367  T D0(1.0);
368  T D1(1.0);
369 
370  T faceG(1.0);
371  T G0(1.0);
372  T G1(1.0);
373 
374  T faceB0(1.0);
375  T faceB1(1.0);
376 
377  T faceNu(1.0);
378 
379  T faceM0(1.0);
380  T faceM1(1.0);
381  T faceM2(1.0);
382 
383  T faceT(1.0);
384 
385  D0 = ymCell[c0]*pow(tCell[c0],three)/(twelve*(one - nuCell[c0]*nuCell[c0]));
386  D1 = ymCell[c1]*pow(tCell[c1],three)/(twelve*(one - nuCell[c1]*nuCell[c1]));
387 
388  G0 = _scf*ymCell[c0]*tCell[c0]/(two*(one+nuCell[c0]));
389  G1 = _scf*ymCell[c1]*tCell[c1]/(two*(one+nuCell[c1]));
390 
391  Diag& a00 = diag[c0];
392  Diag& a11 = diag[c1];
393  OffDiag& a01 = assembler.getCoeff01(f);
394  OffDiag& a10 = assembler.getCoeff10(f);
395 
396  if (vol0 == 0.)
397  {
398  faceD = D1;
399  faceG = G1;
400  faceNu = nuCell[c1];
401  }
402  else if (vol1 == 0.)
403  {
404  faceD = D0;
405  faceG = G0;
406  faceNu = nuCell[c0];
407  }
408  else
409  {
410  faceD = harmonicAverage(D0,D1);
411  faceG = harmonicAverage(G0,G1);
412  faceNu = harmonicAverage(nuCell[c0],nuCell[c1]);
413  }
414 
415  faceD = D0*wt0 + D1*wt1;
416  faceG = G0*wt0 + G1*wt1;
417  faceNu = nuCell[c0]*wt0 + nuCell[c1]*wt1;
418 
419  faceB0 = xCell[c0][0]*bwt0 + xCell[c1][0]*bwt1;
420  faceB1 = xCell[c0][1]*bwt0 + xCell[c1][1]*bwt1;
421 
422  faceT = tCell[c0]*bwt0 + tCell[c1]*bwt1;
423 
424  if(_creep)
425  {
426  faceM0 = plasticMoment[c0][0]*bwt0 + plasticMoment[c1][0]*bwt1;
427  faceM1 = plasticMoment[c0][1]*bwt0 + plasticMoment[c1][1]*bwt1;
428  faceM2 = plasticMoment[c0][2]*bwt0 + plasticMoment[c1][2]*bwt1;
429  }
430 
431  //residual stress contribution//
432  const VGradType residualStressF = (residualCell[c0]*wt0 + residualCell[c1]*wt1);
433  const T diffRMetric = dot(faceArea[f],(residualStressF*faceArea[f]))/dot(faceArea[f],ds);
434  VectorT3 secondaryRCoeff = (residualStressF*faceArea[f]-ds*diffRMetric);
435  secondaryRCoeff[0] *= pow(faceT,three)/twelve;
436  secondaryRCoeff[1] *= pow(faceT,three)/twelve;
437  secondaryRCoeff[2] *= faceT;
438  //
439 
440  const VGradType gradF = (vGradCell[c0]*wt0 + vGradCell[c1]*wt1);
441 
443 
444  // primary w flux, contribution from grad w
445 
446  T wFlux = faceG*diffMetric*(xCell[c1][2]-xCell[c0][2]);
447 
448  // secondary w flux
449  wFlux += faceG*(gradF*secondaryCoeff)[2];
450 
451  // contribution from beta term
452 
453  wFlux += faceG*(faceB0*Af[0]+faceB1*Af[1]);
454 
455  // primary MxFlux
456  T MxFlux = -faceD*diffMetric*(xCell[c1][0]-xCell[c0][0]);
457 
458  // secondary MxFlux
459  MxFlux -= faceD*(gradF*secondaryCoeff)[0];
460 
461  // primary MyFlux
462  T MyFlux = -faceD*diffMetric*(xCell[c1][1]-xCell[c0][1]);
463 
464  // secondary MyFlux
465  MyFlux -= faceD*(gradF*secondaryCoeff)[1];
466 
467  // source for cell 0, uses dfx0 and dfy0
468  source[0] = -faceD*((faceNu*gradF[1][1])*Af[0]+
469  ((one-faceNu)/two)*(gradF[0][1])*Af[1]-
470  ((one+faceNu)/two)*(gradF[1][0])*Af[1])+
471  + dfx0*wFlux + MxFlux;
472 
473  source[1] = -faceD*(((one-faceNu)/two)*(gradF[1][0])*Af[0]-
474  ((one+faceNu)/two)*(gradF[0][1])*Af[0]+
475  (faceNu*gradF[0][0])*Af[1])+
476  + dfy0*wFlux + MyFlux ;
477 
478  source[2] = -wFlux;
479 
480  if(_creep)
481  {
482  source[0] += faceM0*Af[0] + faceM2*Af[1];
483  source[1] += faceM2*Af[0] + faceM1*Af[1];
484  }
485 
486  // add flux to the residual of c0
487  rCell[c0] += source;
488 
489 
490  // source for cell 1, uses dfx1 and dfy1
491  source[0] = -faceD*((faceNu*gradF[1][1])*Af[0]+
492  ((one-faceNu)/two)*(gradF[0][1])*Af[1]-
493  ((one+faceNu)/two)*(gradF[1][0])*Af[1])+
494  + dfx1*wFlux + MxFlux;
495 
496  source[1] = -faceD*(((one-faceNu)/two)*(gradF[1][0])*Af[0]-
497  ((one+faceNu)/two)*(gradF[0][1])*Af[0]+
498  (faceNu*gradF[0][0])*Af[1])+
499  + dfy1*wFlux + MyFlux ;
500 
501  source[2] = -wFlux;
502 
503  if(_creep)
504  {
505  source[0] += faceM0*Af[0] + faceM2*Af[1];
506  source[1] += faceM2*Af[0] + faceM1*Af[1];
507  }
508 
509  rCell[c1] -= source;
510 
511 
512  // linearization of wflux, primary part
513 
514  a00(0,2) += -diffMetric*faceG*dfx0;
515  a00(1,2) += -diffMetric*faceG*dfy0;
516  a00(2,2) += diffMetric*faceG;
517 
518  a01(0,2) += diffMetric*faceG*dfx0;
519  a01(1,2) += diffMetric*faceG*dfy0;
520  a01(2,2) += -diffMetric*faceG;
521 
522  a11(0,2) += -diffMetric*faceG*dfx1;
523  a11(1,2) += -diffMetric*faceG*dfy1;
524  a11(2,2) += diffMetric*faceG;
525 
526  a10(0,2) += diffMetric*faceG*dfx1;
527  a10(1,2) += diffMetric*faceG*dfy1;
528  a10(2,2) += -diffMetric*faceG;
529 
530  // linearization of wflux, primary part
531 
532  a00(0,0) += diffMetric*faceD;
533  a00(1,1) += diffMetric*faceD;
534 
535  a01(0,0) += -diffMetric*faceD;
536  a01(1,1) += -diffMetric*faceD;
537 
538  a11(0,0) += diffMetric*faceD;
539  a11(1,1) += diffMetric*faceD;
540 
541  a10(0,0) += -diffMetric*faceD;
542  a10(1,1) += -diffMetric*faceD;
543 
544  // linearization of beta term
545  Diag coeffPair;
546 
547  coeffPair(0,0)=faceG*dfx0*Af[0];
548  coeffPair(0,1)=faceG*dfx0*Af[1];
549  coeffPair(0,2)=zero;
550 
551  coeffPair(1,0)=faceG*dfy0*Af[0];
552  coeffPair(1,1)=faceG*dfy0*Af[1];
553  coeffPair(1,2)=zero;
554 
555  coeffPair(2,0)=-faceG*Af[0];
556  coeffPair(2,1)=-faceG*Af[1];
557  coeffPair(2,2)=zero;
558 
559  a00 += bwt0*coeffPair;
560  a01 += bwt1*coeffPair;
561 
562  coeffPair(0,0)=faceG*dfx1*Af[0];
563  coeffPair(0,1)=faceG*dfx1*Af[1];
564  coeffPair(0,2)=zero;
565 
566  coeffPair(1,0)=faceG*dfy1*Af[0];
567  coeffPair(1,1)=faceG*dfy1*Af[1];
568  coeffPair(1,2)=zero;
569 
570  coeffPair(2,0)=-faceG*Af[0];
571  coeffPair(2,1)=-faceG*Af[1];
572  coeffPair(2,2)=zero;
573 
574 
575  a10 -= bwt0*coeffPair;
576  a11 -= bwt1*coeffPair;
577 
580  if (_fullLinearization)
581  {
582  // loop over level 0 neighbors
583  for(int nnb = ccRow[c0]; nnb<ccRow[c0+1]; nnb++)
584  {
585  const int nb = ccCol[nnb];
586 
587  // get the coefficient from the gradient matrix that uses level 2 connectivity
588  // getting a copy rather than reference since we might modify it
589  VectorT3 g_nb = vgMatrix.getCoeff(c0,nb);
590 #if 1
591  if (isSymmetry)
592  {
593  const T gnb_dot_nx2 = T(2.0)*dot(en,g_nb);
594  g_nb = gnb_dot_nx2*en;
595  }
596 #endif
597  Diag coeff;
598 
599  coeff(0,0)=-wt0*faceD*(-((one+faceNu)/two)*Af[1]*g_nb[1]);
600  coeff(0,1)=-wt0*faceD*(((one-faceNu)/two)*Af[1]*g_nb[0]+faceNu*Af[0]*g_nb[1]);
601  coeff(0,2)=wt0*faceG*dfx0*(g_nb[0]*secondaryCoeff[0]+g_nb[1]*secondaryCoeff[1]);
602 
603  coeff(1,0)=-wt0*faceD*(((one-faceNu)/two)*Af[0]*g_nb[1]+faceNu*Af[1]*g_nb[0]);
604  coeff(1,1)=-wt0*faceD*(-((one+faceNu)/two)*Af[0]*g_nb[0]);
605  coeff(1,2)=wt0*faceG*dfy0*(g_nb[0]*secondaryCoeff[0]+g_nb[1]*secondaryCoeff[1]);
606 
607  coeff(2,0)=zero;
608  coeff(2,1)=zero;
609  coeff(2,2)=-wt0*faceG*(g_nb[0]*secondaryCoeff[0]+g_nb[1]*secondaryCoeff[1]);
610 
611  //secondary coefficient
612  for(int k=0; k<3; k++)
613  {
614  coeff(0,0) -= wt0*faceD*secondaryCoeff[k]*g_nb[k];
615  coeff(1,1) -= wt0*faceD*secondaryCoeff[k]*g_nb[k];
616  }
617 
618  //residual stress
619  for(int i=0; i<3; i++)
620  {
621  for(int k=0; k<3; k++)
622  coeff(i,i) -= wt0*secondaryRCoeff[k]*g_nb[k];
623  }
624  //
625 
626  OffDiag& a0_nb = matrix.getCoeff(c0,nb);
627 
628  a0_nb += coeff;
629  a00 -= coeff;
630 
631 
632  coeff(0,0)=-wt0*faceD*(-((one+faceNu)/two)*Af[1]*g_nb[1]);
633  coeff(0,1)=-wt0*faceD*(((one-faceNu)/two)*Af[1]*g_nb[0]+faceNu*Af[0]*g_nb[1]);
634  coeff(0,2)=wt0*faceG*dfx1*(g_nb[0]*secondaryCoeff[0]+g_nb[1]*secondaryCoeff[1]);
635 
636  coeff(1,0)=-wt0*faceD*(((one-faceNu)/two)*Af[0]*g_nb[1]+faceNu*Af[1]*g_nb[0]);
637  coeff(1,1)=-wt0*faceD*(-((one+faceNu)/two)*Af[0]*g_nb[0]);
638  coeff(1,2)=wt0*faceG*dfy1*(g_nb[0]*secondaryCoeff[0]+g_nb[1]*secondaryCoeff[1]);
639 
640  coeff(2,0)=zero;
641  coeff(2,1)=zero;
642  coeff(2,2)=-wt0*faceG*(g_nb[0]*secondaryCoeff[0]+g_nb[1]*secondaryCoeff[1]);
643 
644  //secondary coefficient
645  for(int k=0; k<3; k++)
646  {
647  coeff(0,0) -= wt0*faceD*secondaryCoeff[k]*g_nb[k];
648  coeff(1,1) -= wt0*faceD*secondaryCoeff[k]*g_nb[k];
649  }
650 
651  //residual stress
652  for(int i=0; i<3; i++)
653  {
654  for(int k=0; k<3; k++)
655  coeff(i,i) -= wt0*secondaryRCoeff[k]*g_nb[k];
656  }
657  //
658 
659  a10 += coeff;
660 
661  if (c1 != nb)
662  {
663  OffDiag& a1_nb = matrix.getCoeff(c1,nb);
664  a1_nb -= coeff;
665  }
666  else
667  {
668  a11 -= coeff;
669  }
670  }
671 
672 
673  if (!isBoundary)
674  {
675  for(int nnb = ccRow[c1]; nnb<ccRow[c1+1]; nnb++)
676  {
677  const int nb = ccCol[nnb];
678  const VectorT3& g_nb = vgMatrix.getCoeff(c1,nb);
679 
680  Diag coeff;
681 
682  coeff(0,0)=-wt1*faceD*(-((one+faceNu)/two)*Af[1]*g_nb[1]);
683  coeff(0,1)=-wt1*faceD*(((one-faceNu)/two)*Af[1]*g_nb[0]+faceNu*Af[0]*g_nb[1]);
684  coeff(0,2)=wt1*faceG*dfx1*(g_nb[0]*secondaryCoeff[0]+g_nb[1]*secondaryCoeff[1]);
685 
686  coeff(1,0)=-wt1*faceD*(((one-faceNu)/two)*Af[0]*g_nb[1]+faceNu*Af[1]*g_nb[0]);
687  coeff(1,1)=-wt1*faceD*(-((one+faceNu)/two)*Af[0]*g_nb[0]);
688  coeff(1,2)=wt1*faceG*dfy1*(g_nb[0]*secondaryCoeff[0]+g_nb[1]*secondaryCoeff[1]);
689 
690  coeff(2,0)=zero;
691  coeff(2,1)=zero;
692  coeff(2,2)=-wt1*faceG*(g_nb[0]*secondaryCoeff[0]+g_nb[1]*secondaryCoeff[1]);
693 
694  //secondary coefficient
695  for(int k=0; k<3; k++)
696  {
697  coeff(0,0) -= wt1*faceD*secondaryCoeff[k]*g_nb[k];
698  coeff(1,1) -= wt1*faceD*secondaryCoeff[k]*g_nb[k];
699  }
700 
701  //residual stress
702  for(int i=0; i<3; i++)
703  {
704  for(int k=0; k<3; k++)
705  coeff(i,i) -= wt1*secondaryRCoeff[k]*g_nb[k];
706  }
707  //
708 
709  OffDiag& a1_nb = matrix.getCoeff(c1,nb);
710  a1_nb -= coeff;
711  a11 += coeff;
712 
713 
714  coeff(0,0)=-wt1*faceD*(-((one+faceNu)/two)*Af[1]*g_nb[1]);
715  coeff(0,1)=-wt1*faceD*(((one-faceNu)/two)*Af[1]*g_nb[0]+faceNu*Af[0]*g_nb[1]);
716  coeff(0,2)=wt1*faceG*dfx0*(g_nb[0]*secondaryCoeff[0]+g_nb[1]*secondaryCoeff[1]);
717 
718  coeff(1,0)=-wt1*faceD*(((one-faceNu)/two)*Af[0]*g_nb[1]+faceNu*Af[1]*g_nb[0]);
719  coeff(1,1)=-wt1*faceD*(-((one+faceNu)/two)*Af[0]*g_nb[0]);
720  coeff(1,2)=wt1*faceG*dfy0*(g_nb[0]*secondaryCoeff[0]+g_nb[1]*secondaryCoeff[1]);
721 
722  coeff(2,0)=zero;
723  coeff(2,1)=zero;
724  coeff(2,2)=-wt1*faceG*(g_nb[0]*secondaryCoeff[0]+g_nb[1]*secondaryCoeff[1]);
725 
726  //secondary coefficient
727  for(int k=0; k<3; k++)
728  {
729  coeff(0,0) -= wt1*faceD*secondaryCoeff[k]*g_nb[k];
730  coeff(1,1) -= wt1*faceD*secondaryCoeff[k]*g_nb[k];
731  }
732 
733  //residual stress
734  for(int i=0; i<3; i++)
735  {
736  for(int k=0; k<3; k++)
737  coeff(i,i) -= wt1*secondaryRCoeff[k]*g_nb[k];
738  }
739  //
740 
741  a01 -= coeff;
742 
743  if (c0 != nb)
744  {
745  OffDiag& a0_nb = matrix.getCoeff(c0,nb);
746 
747  a0_nb += coeff;
748  }
749  else
750  a00 += coeff;
751  }
752  }
753  }
754  //contribution due to residual stress
755 
756  //primary part
757 
758  VectorT3 residualSource(NumTypeTraits<VectorT3>::getZero());
759  residualSource -= diffRMetric*(xCell[c1]-xCell[c0]);
760  residualSource[0] *= pow(faceT,three)/twelve;
761  residualSource[1] *= pow(faceT,three)/twelve;
762  residualSource[2] *= faceT;
763 
764  // secondary part
765 
766  residualSource -= gradF*secondaryRCoeff;
767 
768  // add flux due to residual Stress to the residual of c0 and c1
769  rCell[c0] += residualSource;
770  rCell[c1] -= residualSource;
771 
772  const T diffRCoeff = -diffRMetric;
773 
774 
775  a01(0,0) +=(pow(faceT,three)/twelve)*diffRCoeff;
776  a10(0,0) +=(pow(faceT,three)/twelve)*diffRCoeff;
777 
778 
779  a00(0,0) -=(pow(faceT,three)/twelve)*diffRCoeff;
780  a11(0,0) -=(pow(faceT,three)/twelve)*diffRCoeff;
781 
782  a01(1,1) +=(pow(faceT,three)/twelve)*diffRCoeff;
783  a10(1,1) +=(pow(faceT,three)/twelve)*diffRCoeff;
784 
785 
786  a00(1,1) -=(pow(faceT,three)/twelve)*diffRCoeff;
787  a11(1,1) -=(pow(faceT,three)/twelve)*diffRCoeff;
788 
789 
790  a01(2,2) +=faceT*diffRCoeff;
791  a10(2,2) +=faceT*diffRCoeff;
792 
793 
794  a00(2,2) -=faceT*diffRCoeff;
795  a11(2,2) -=faceT*diffRCoeff;
796 
797  //residual stress
798 
799  }
800  }
const Array< int > & getCol() const
const Array< int > & getRow() const
Matrix & getMatrix(const Index &rowIndex, const Index &colIndex)
int getSelfCount() const
Definition: StorageSite.h:40
Field coordinate
Definition: GeomFields.h:19
T harmonicAverage(const T &x0, const T &x1)
CRMatrix< Diag, OffDiag, VectorT3 > CCMatrix
VGradModelType::GradMatrixType VGradMatrix
CCMatrix::PairWiseAssembler CCAssembler
Array< Gradient< VectorT3 > > VGradArray
pair< const Field *, const StorageSite * > ArrayIndex
Definition: MultiField.h:21
const CRConnectivity & getCellCells() const
Definition: Mesh.cpp:480
const StorageSite & getCells() const
Definition: Mesh.h:109
Field volume
Definition: GeomFields.h:26
const CRConnectivity & getFaceCells(const StorageSite &site) const
Definition: Mesh.cpp:388
static GradMatrixType & getGradientMatrix(const Mesh &mesh, const GeomFields &geomFields)
int getCount() const
Definition: StorageSite.h:39
Field area
Definition: GeomFields.h:23
Field areaMag
Definition: GeomFields.h:25
T dot(const Vector< T, 3 > &a, const Vector< T, 3 > &b)
Definition: Vector.h:253

Member Data Documentation

template<class T , class Diag , class OffDiag >
const T PlateSourceDiscretization< T, Diag, OffDiag >::_A
private
template<class T , class Diag , class OffDiag >
const T PlateSourceDiscretization< T, Diag, OffDiag >::_B
private
template<class T , class Diag , class OffDiag >
const bool PlateSourceDiscretization< T, Diag, OffDiag >::_creep
private
template<class T , class Diag , class OffDiag >
const int PlateSourceDiscretization< T, Diag, OffDiag >::_creepModel
private
template<class T , class Diag , class OffDiag >
const Field& PlateSourceDiscretization< T, Diag, OffDiag >::_devStressField
private
template<class T , class Diag , class OffDiag >
const Field& PlateSourceDiscretization< T, Diag, OffDiag >::_forceField
private
template<class T , class Diag , class OffDiag >
const bool PlateSourceDiscretization< T, Diag, OffDiag >::_fullLinearization
private
template<class T , class Diag , class OffDiag >
const GeomFields& PlateSourceDiscretization< T, Diag, OffDiag >::_geomFields
private
template<class T , class Diag , class OffDiag >
const T PlateSourceDiscretization< T, Diag, OffDiag >::_m
private
template<class T , class Diag , class OffDiag >
const T PlateSourceDiscretization< T, Diag, OffDiag >::_n
private
template<class T , class Diag , class OffDiag >
const Field& PlateSourceDiscretization< T, Diag, OffDiag >::_nuField
private
template<class T , class Diag , class OffDiag >
const int PlateSourceDiscretization< T, Diag, OffDiag >::_nz
private
template<class T , class Diag , class OffDiag >
Field& PlateSourceDiscretization< T, Diag, OffDiag >::_plasticMomentField
private
template<class T , class Diag , class OffDiag >
Field& PlateSourceDiscretization< T, Diag, OffDiag >::_plasticStrainField
private
template<class T , class Diag , class OffDiag >
Field& PlateSourceDiscretization< T, Diag, OffDiag >::_plasticStrainN1Field
private
template<class T , class Diag , class OffDiag >
Field& PlateSourceDiscretization< T, Diag, OffDiag >::_plasticStrainOutField
private
template<class T , class Diag , class OffDiag >
const Field& PlateSourceDiscretization< T, Diag, OffDiag >::_residualStressField
private
template<class T , class Diag , class OffDiag >
const T& PlateSourceDiscretization< T, Diag, OffDiag >::_scf
private
template<class T , class Diag , class OffDiag >
const T PlateSourceDiscretization< T, Diag, OffDiag >::_Sy0
private
template<class T , class Diag , class OffDiag >
const Field& PlateSourceDiscretization< T, Diag, OffDiag >::_thicknessField
private
template<class T , class Diag , class OffDiag >
const T PlateSourceDiscretization< T, Diag, OffDiag >::_timeStep
private
template<class T , class Diag , class OffDiag >
Field& PlateSourceDiscretization< T, Diag, OffDiag >::_varField
private
template<class T , class Diag , class OffDiag >
const Field& PlateSourceDiscretization< T, Diag, OffDiag >::_varGradientField
private
template<class T , class Diag , class OffDiag >
const Field& PlateSourceDiscretization< T, Diag, OffDiag >::_VMStressField
private
template<class T , class Diag , class OffDiag >
const Field& PlateSourceDiscretization< T, Diag, OffDiag >::_ymField
private

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