Memosa-FVM  0.2
PlateSourceDiscretization.h
Go to the documentation of this file.
1 // This file os part of FVM
2 // Copyright (c) 2012 FVM Authors
3 // See LICENSE file for terms.
4 
5 #ifndef _PLATESOURCEDISCRETIZATION_H_
6 #define _PLATESOURCEDISCRETIZATION_H_
7 
8 #include "Field.h"
9 #include "MultiField.h"
10 #include "MultiFieldMatrix.h"
11 #include "Mesh.h"
12 #include "Discretization.h"
13 #include "StorageSite.h"
14 #include "GeomFields.h"
15 #include "CRConnectivity.h"
16 #include "CRMatrixRect.h"
17 #include "Vector.h"
18 #include "GradientModel.h"
19 
20 template<class T, class Diag, class OffDiag>
21 class PlateSourceDiscretization : public Discretization
22 {
23 public:
24 
25  typedef Array<T> TArray;
26  typedef Vector<T,3> VectorT3;
27  typedef Array<VectorT3> VectorT3Array;
28  typedef Vector<T,4> VectorT4;
29  typedef Array<VectorT4> VectorT4Array;
30  typedef Gradient<VectorT3> VGradType;
31  typedef Array<Gradient<VectorT3> > VGradArray;
32 
33  typedef CRMatrix<Diag,OffDiag,VectorT3> CCMatrix;
34  typedef typename CCMatrix::DiagArray DiagArray;
35  typedef typename CCMatrix::PairWiseAssembler CCAssembler;
36 
37  typedef GradientModel<VectorT3> VGradModelType;
38  typedef typename VGradModelType::GradMatrixType VGradMatrix;
39 
40  PlateSourceDiscretization(const MeshList& meshes,
41  const GeomFields& geomFields,
42  Field& varField,
43  const Field& ymField,
44  const Field& nuField,
45  const Field& varGradientField,
46  const Field& residualStressField,
47  const Field& thicknessField,
48  const Field& forceField,
49  const T& scf,
50  const Field& devStressField,
51  const Field& VMStressField,
52  Field& plasticStrainField,
53  Field& plasticStrainOutField,
54  Field& plasticStrainN1Field,
55  Field& plasticMomentField,
56  const T& A,
57  const T& B,
58  const T& mm,
59  const T& nn,
60  const T& Sy0,
61  const int& nz,
62  const T& timeStep,
63  const int& creepModel,
64  const bool& creep,
65  bool fullLinearization=true) :
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  {}
93 
94  void discretize(const Mesh& mesh, MultiFieldMatrix& mfmatrix,
95  MultiField& xField, MultiField& rField)
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  }
260 
261 
262  void discretizeFaces(const Mesh& mesh, const StorageSite& faces,
263  MultiFieldMatrix& mfmatrix,
264  MultiField& xField, MultiField& rField,
265  const bool isBoundary, const bool isSymmetry)
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 
316  const VGradMatrix& vgMatrix = VGradModelType::getGradientMatrix(mesh,_geomFields);
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 
442  VectorT3 source(NumTypeTraits<VectorT3>::getZero());
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 
578  VectorT3 s0(NumTypeTraits<VectorT3>::getZero());
579  VectorT3 s1(NumTypeTraits<VectorT3>::getZero());
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  }
801 
802 
803 private:
804  const GeomFields& _geomFields;
805  Field& _varField;
806  const Field& _ymField;
807  const Field& _nuField;
808  const Field& _varGradientField;
809  const Field& _residualStressField;
810  const Field& _thicknessField;
811  const Field& _forceField;
812  const T& _scf;
813  const Field& _devStressField;
814  const Field& _VMStressField;
815  Field& _plasticStrainField;
816  Field& _plasticStrainOutField;
817  Field& _plasticStrainN1Field;
818  Field& _plasticMomentField;
819  const T _A;
820  const T _B;
821  const T _m;
822  const T _n;
823  const T _Sy0;
824  const int _nz;
825  const T _timeStep;
826  const int _creepModel;
827  const bool _creep;
828  const bool _fullLinearization;
829 };
830 
831 #endif
CRConnectivity::getCol
const Array< int > & getCol() const
Definition: CRConnectivity.h:183
CRConnectivity::getRow
const Array< int > & getRow() const
Definition: CRConnectivity.h:182
Mesh::getInteriorFaceGroup
const FaceGroup & getInteriorFaceGroup() const
Definition: Mesh.h:181
MultiFieldMatrix::getMatrix
Matrix & getMatrix(const Index &rowIndex, const Index &colIndex)
Definition: MultiFieldMatrix.cpp:42
FaceGroupPtr
shared_ptr< FaceGroup > FaceGroupPtr
Definition: Mesh.h:46
StorageSite::getSelfCount
int getSelfCount() const
Definition: StorageSite.h:40
GeomFields::coordinate
Field coordinate
Definition: GeomFields.h:19
harmonicAverage
T harmonicAverage(const T &x0, const T &x1)
Definition: DiffusionDiscretization.h:20
FaceGroup
Definition: Mesh.h:28
Mesh::getAllFaceGroups
const FaceGroupList & getAllFaceGroups() const
Definition: Mesh.h:193
Field
Definition: Field.h:14
PlateSourceDiscretization::CCMatrix
CRMatrix< Diag, OffDiag, VectorT3 > CCMatrix
Definition: PlateSourceDiscretization.h:33
CRMatrix::getCoeff
OffDiag & getCoeff(const int i, const int j)
Definition: CRMatrix.h:836
Mesh
Definition: Mesh.h:49
CRMatrix::PairWiseAssembler::getCoeff10
OffDiag & getCoeff10(const int np)
Definition: CRMatrix.h:131
GradientMatrix::getCoeff
Coord & getCoeff(const int i, const int j)
Definition: GradientMatrix.h:129
CRMatrix::getDiag
Array< Diag > & getDiag()
Definition: CRMatrix.h:856
PlateSourceDiscretization::discretizeFaces
void discretizeFaces(const Mesh &mesh, const StorageSite &faces, MultiFieldMatrix &mfmatrix, MultiField &xField, MultiField &rField, const bool isBoundary, const bool isSymmetry)
Definition: PlateSourceDiscretization.h:262
PlateSourceDiscretization::VGradMatrix
VGradModelType::GradMatrixType VGradMatrix
Definition: PlateSourceDiscretization.h:38
FaceGroup::groupType
string groupType
Definition: Mesh.h:42
PlateSourceDiscretization::CCAssembler
CCMatrix::PairWiseAssembler CCAssembler
Definition: PlateSourceDiscretization.h:35
sqrt
Tangent sqrt(const Tangent &a)
Definition: Tangent.h:317
PlateSourceDiscretization::VGradArray
Array< Gradient< VectorT3 > > VGradArray
Definition: PlateSourceDiscretization.h:31
MultiField::ArrayIndex
pair< const Field *, const StorageSite * > ArrayIndex
Definition: MultiField.h:21
Mesh::getCellCells
const CRConnectivity & getCellCells() const
Definition: Mesh.cpp:480
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
GeomFields::volume
Field volume
Definition: GeomFields.h:26
PlateSourceDiscretization::discretize
void discretize(const Mesh &mesh, MultiFieldMatrix &mfmatrix, MultiField &xField, MultiField &rField)
Definition: PlateSourceDiscretization.h:94
StorageSite::getCountLevel1
int getCountLevel1() const
Definition: StorageSite.h:72
CRMatrix::PairWiseAssembler::getCoeff01
OffDiag & getCoeff01(const int np)
Definition: CRMatrix.h:126
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
Array
Definition: Array.h:14
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
PlateSourceDiscretization::VGradModelType
GradientModel< VectorT3 > VGradModelType
Definition: PlateSourceDiscretization.h:37
GeomFields::area
Field area
Definition: GeomFields.h:23
PlateSourceDiscretization::PlateSourceDiscretization
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)
Definition: PlateSourceDiscretization.h:40
CRMatrix::getPairWiseAssembler
PairWiseAssembler & getPairWiseAssembler(const CRConnectivity &pairs)
Definition: CRMatrix.h:867
GeomFields::areaMag
Field areaMag
Definition: GeomFields.h:25
dot
T dot(const Vector< T, 3 > &a, const Vector< T, 3 > &b)
Definition: Vector.h:253
MeshList
vector< Mesh * > MeshList
Definition: Mesh.h:439
FaceGroup::site
StorageSite site
Definition: Mesh.h:40