Memosa-FVM  0.2
StructureModel_impl.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 #ifdef FVM_PARALLEL
6 #include <mpi.h>
7 #endif
8 
9 #include "Mesh.h"
10 
11 #include "NumType.h"
12 #include "Array.h"
13 #include "Field.h"
14 #include "CRConnectivity.h"
15 #include "LinearSystem.h"
16 #include "StorageSite.h"
17 #include "MultiFieldMatrix.h"
18 #include "CRMatrix.h"
19 #include "FluxJacobianMatrix.h"
20 #include "DiagonalMatrix.h"
21 #include "GenericBCS.h"
22 #include "Vector.h"
23 #include "VectorTranspose.h"
24 #include "DiffusionDiscretization.h"
25 #include "TimeDerivativeStructureDiscretization.h"
26 #include "StructureSourceDiscretization.h"
27 #include "StructurePlasticDiscretization.h"
28 #include "SourceDiscretization.h"
29 #include "Underrelaxer.h"
30 #include "AMG.h"
31 #include "Linearizer.h"
32 #include "GradientModel.h"
33 #include "GenericIBDiscretization.h"
34 #include "StressTensor.h"
35 #include "SquareTensor.h"
36 
37 
38 
39 template<class X, class Diag, class OffDiag>
40 class StructureBCS
41 {
42 public:
43 
44  typedef typename NumTypeTraits<X>::T_Scalar T_Scalar;
45 
46  typedef Array<T_Scalar> TArray;
47  typedef Vector<T_Scalar,3> VectorT3;
48 
49  typedef CRMatrix<Diag,OffDiag,X> CCMatrix;
50  typedef typename CCMatrix::PairWiseAssembler CCAssembler;
51 
52  typedef FluxJacobianMatrix<Diag,X> FMatrix;
53  typedef DiagonalMatrix<Diag,X> BBMatrix;
54 
55  typedef Array<Diag> DiagArray;
56  typedef Array<OffDiag> OffDiagArray;
57 
58  typedef Array<X> XArray;
59  typedef Array<VectorT3> VectorT3Array;
60 
61 
62  StructureBCS(const StorageSite& faces,
63  const Mesh& mesh,
64  const GeomFields& geomFields,
65  Field& varField,
66  MultiFieldMatrix& matrix,
67  MultiField& xField,
68  MultiField& rField,
69  const bool explicitMode
70  ) :
71  _faces(faces),
72  _cells(mesh.getCells()),
73  _faceCells(mesh.getFaceCells(_faces)),
74  _varField(varField),
75  _xIndex(&_varField,&_cells),
76  _dRdX(dynamic_cast<CCMatrix&>(matrix.getMatrix(_xIndex,_xIndex))),
77  _assembler(_dRdX.getPairWiseAssembler(_faceCells)),
78  _dRdXDiag(_dRdX.getDiag()),
79  _x(dynamic_cast<XArray&>(xField[_xIndex])),
80  _r(dynamic_cast<XArray&>(rField[_xIndex])),
81  _areaMagField(geomFields.areaMag),
82  _faceAreaMag(dynamic_cast<const TArray&>(_areaMagField[_faces])),
83  _areaField(geomFields.area),
84  _faceArea(dynamic_cast<const VectorT3Array&>(_areaField[_faces])),
85  _explicitMode(explicitMode)
86  {}
87 
88  X applyDirichletBC(int f, const X& bValue) const
89  {
90  const int c1 = _faceCells(f,1);
91 
92  const X fluxB = -_r[c1];
93 
94  const X dXC1 = bValue - _x[c1];
95 
96  _dRdX.eliminateDirichlet(c1,_r,dXC1, _explicitMode);
97  _x[c1] = bValue;
98  _r[c1] = NumTypeTraits<X>::getZero();
99  if (!_explicitMode)
100  _dRdX.setDirichlet(c1);
101  return fluxB;
102  }
103 
104  X applyDirichletBC(const X& bValue) const
105  {
106  X fluxB(NumTypeTraits<X>::getZero());
107  for(int i=0; i<_faces.getCount(); i++)
108  fluxB += applyDirichletBC(i,bValue);
109  return fluxB;
110  }
111 
112  X applyDirichletBC(const FloatValEvaluator<X>& bValue) const
113  {
114  X fluxB(NumTypeTraits<X>::getZero());
115  for(int i=0; i<_faces.getCount(); i++)
116  fluxB += applyDirichletBC(i,bValue[i]);
117  return fluxB;
118  }
119 
120  X applyNeumannBC(const int f,
121  const X& specifiedFlux) const
122  {
123  //const int c0 = _faceCells(f,0);
124  const int c1 = _faceCells(f,1);
125  const X fluxB = -_r[c1];
126 
127  // the current value of flux and its Jacobians
128  const X dFlux = specifiedFlux*_faceAreaMag[f] - fluxB;
129 
130  // setup the equation for the boundary value; the coefficients
131  // are already computed so just need to set the rhs
132  _r[c1] = dFlux;
133 
134  // _dRdX.eliminate(c1,_r);
135  // mark this row as a "boundary" row so that we will update it
136  // after the overall system is solved
137  _dRdX.setBoundary(c1);
138  return fluxB;
139  }
140 
141 
142  X applyNeumannBC(const X& bFlux) const
143  {
144  X fluxB(NumTypeTraits<X>::getZero());
145 
146  for(int i=0; i<_faces.getCount(); i++)
147  fluxB += applyNeumannBC(i,bFlux);
148  return fluxB;
149  }
150 
151  X applyZeroDerivativeBC() const
152  {
153  X fluxTotal(NumTypeTraits<X>::getZero());
154  for(int f=0; f<_faces.getCount(); f++)
155  {
156  const int c0 = _faceCells(f,0);
157  const int c1 = _faceCells(f,1);
158  const X fluxB = -_r[c1];
159 
160  // setup the equation for the boundary value; we want xb to be equal to xc
161 
162  _r[c1] = _x[c0] - _x[c1];
163  _dRdXDiag[c1] = -NumTypeTraits<Diag>::getUnity();
164  this->_assembler.getCoeff10(f) = NumTypeTraits<OffDiag>::getUnity();
165 
166 
167  // mark this row as a "boundary" row so that we will update it
168  // after the overall system is solved
169  _dRdX.setBoundary(c1);
170  fluxTotal += fluxB;
171  }
172  return fluxTotal;
173  }
174 
175  X applyNeumannBC(const FloatValEvaluator<X>& bFlux) const
176  {
177  X fluxB(NumTypeTraits<X>::getZero());
178  for(int i=0; i<_faces.getCount(); i++)
179  fluxB += applyNeumannBC(i,bFlux[i]);
180  return fluxB;
181  }
182 
183  void applyInterfaceBC(const int f) const
184  {
185  // the boundary cell could be either c0 or c1 at an interface
186  int cb = _faceCells(f,1);
187  T_Scalar sign(NumTypeTraits<T_Scalar>::getUnity());
188  if (cb < _cells.getSelfCount())
189  {
190  cb = _faceCells(f,0);
191  sign *= -1.0;
192  }
193 
194 
195  _r[cb] = T_Scalar(0);
196 
197  if (sign>0)
198  _assembler.getCoeff10(f) = NumTypeTraits<OffDiag>::getZero();
199  else
200  _assembler.getCoeff01(f) = NumTypeTraits<OffDiag>::getZero();
201 
202  }
203 
204  void applyInterfaceBC() const
205  {
206  for(int i=0; i<_faces.getCount(); i++)
207  applyInterfaceBC(i);
208  }
209 
210  void applySymmetryBC() const
211  {
212  for(int f=0; f<this->_faces.getCount(); f++)
213  {
214  const int c0 = this->_faceCells(f,0);
215  const int c1 = this->_faceCells(f,1);
216 
217 
218  const VectorT3 en = this->_faceArea[f]/this->_faceAreaMag[f];
219  const T_Scalar xC0_dotn = dot(this->_x[c0],en);
220  const X xB = this->_x[c0] - 2.*xC0_dotn * en;
221 
222  Diag dxBdxC0(Diag::getZero());
223  dxBdxC0(0,0) = 1.0 - 2.*en[0]*en[0];
224  dxBdxC0(0,1) = - 2.*en[0]*en[1];
225  dxBdxC0(0,2) = - 2.*en[0]*en[2];
226 
227  dxBdxC0(1,0) = - 2.*en[1]*en[0];
228  dxBdxC0(1,1) = 1.0 - 2.*en[1]*en[1];
229  dxBdxC0(1,2) = - 2.*en[1]*en[2];
230 
231  dxBdxC0(2,0) = - 2.*en[2]*en[0];
232  dxBdxC0(2,1) = - 2.*en[2]*en[1];
233  dxBdxC0(2,2) = 1.0 - 2.*en[2]*en[2];
234 
235 
236  const X xc1mxB = xB-this->_x[c1];
237 
238  // boundary value equation
239  // set all neighbour coeffs to zero first and ap to -1
240  this->_dRdX.setDirichlet(c1);
241 
242  // dependance on c0
243  this->_assembler.getCoeff10(f) = dxBdxC0;
244  this->_r[c1] = xc1mxB;
245 
246  }
247  }
248 
249 
250 
251 protected:
252  const StorageSite& _faces;
253  const StorageSite& _cells;
254  const CRConnectivity& _faceCells;
255  const Field& _varField;
256  const MultiField::ArrayIndex _xIndex;
257  CCMatrix& _dRdX;
258  CCAssembler& _assembler;
259  DiagArray& _dRdXDiag;
260  XArray& _x;
261  XArray& _r;
262  const Field& _areaMagField;
263  const TArray& _faceAreaMag;
264  const Field& _areaField;
265  const VectorT3Array& _faceArea;
266  const bool _explicitMode;
267 };
268 
269 template<class T>
270 class StructureModel<T>::Impl
271 {
272 public:
273  typedef Array<T> TArray;
274  typedef Vector<T,3> VectorT3;
275  typedef VectorTranspose<T,3> VectorT3T;
276 
277  typedef Array<VectorT3> VectorT3Array;
278  //typedef DiagonalTensor<T,3> DiagTensorT3;
279  //typedef SquareTensor<T,3> SquareTensorT3;
280  typedef SquareTensor<T,3> DiagTensorT3;
281 
282  typedef CRMatrix<DiagTensorT3,DiagTensorT3,VectorT3> VVMatrix;
283  typedef typename VVMatrix::DiagArray VVDiagArray;
284 
285 
286  typedef Gradient<VectorT3> VGradType;
287  typedef Array<Gradient<VectorT3> > VGradArray;
288 
289  Impl(const GeomFields& geomFields,
290  StructureFields& structureFields,
291  const MeshList& meshes) :
292  _meshes(meshes),
293  _geomFields(geomFields),
294  _structureFields(structureFields),
295  _deformationGradientModel(_meshes,_structureFields.deformation,
296  _structureFields.deformationGradient,_geomFields),
297  _initialDeformationNorm(),
298  _niters(0)
299  {
300  const int numMeshes = _meshes.size();
301  for (int n=0; n<numMeshes; n++)
302  {
303  const Mesh& mesh = *_meshes[n];
304 
305  StructureVC<T> *vc(new StructureVC<T>());
306  vc->vcType = "structure";
307  _vcMap[mesh.getID()] = vc;
308  foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
309  {
310  const FaceGroup& fg = *fgPtr;
311  if ((_bcMap.find(fg.id) == _bcMap.end())&&(fg.groupType != "interior"))
312  {
313  StructureBC<T> *bc(new StructureBC<T>());
314 
315  _bcMap[fg.id] = bc;
316  if (fg.groupType == "wall")
317  {
318  bc->bcType = "SpecifiedTraction";
319  }
320  else if (fg.groupType == "interface")
321  {
322  bc->bcType = "Interface";
323  }
324  else if (fg.groupType == "symmetry")
325  {
326  bc->bcType = "Symmetry";
327  }
328  else if ((fg.groupType == "velocity-inlet") ||
329  (fg.groupType == "pressure-inlet"))
330  {
331  bc->bcType = "SpecifiedDeformation";
332  }
333  else if (fg.groupType == "pressure-outlet")
334  {
335  bc->bcType = "SpecifiedForce";
336  }
337  else
338  throw CException("StructuralModel: unknown face group type "
339  + fg.groupType);
340  }
341  }
342  }
343  }
344 
345  void init()
346  {
347  const int numMeshes = _meshes.size();
348  for (int n=0; n<numMeshes; n++)
349  {
350  const Mesh& mesh = *_meshes[n];
351 
352  const StructureVC<T>& vc = *_vcMap[mesh.getID()];
353 
354  const StorageSite& cells = mesh.getCells();
355  // const StorageSite& faces = mesh.getFaces();
356 
357  shared_ptr<VectorT3Array> sCell(new VectorT3Array(cells.getCountLevel1()));
358 
359  VectorT3 initialDeformation;
360  initialDeformation[0] = _options["initialXDeformation"];
361  initialDeformation[1] = _options["initialYDeformation"];
362  initialDeformation[2] = _options["initialZDeformation"];
363  *sCell = initialDeformation;
364 
365  _structureFields.deformation.addArray(cells,sCell);
366 
367  if (_options.transient)
368  {
369  _structureFields.volume0.addArray(cells,
370  dynamic_pointer_cast<ArrayBase>
371  (_geomFields.volume[cells].newCopy()));
372  _structureFields.deformationN1.addArray(cells,
373  dynamic_pointer_cast<ArrayBase>(sCell->newCopy()));
374  _structureFields.deformationN2.addArray(cells,
375  dynamic_pointer_cast<ArrayBase>(sCell->newCopy()));
376  if (_options.timeDiscretizationOrder > 1)
377  _structureFields.deformationN3.addArray(cells,
378  dynamic_pointer_cast<ArrayBase>(sCell->newCopy()));
379  if(_options.variableTimeStep)
380  {
381  _options.timeStepN1 = _options["timeStep"];
382  _options.timeStepN2 = _options["timeStep"];
383  }
384  }
385 
386  if (_options.creep)
387  {
388  _structureFields.elasticDeformation.addArray(cells,
389  dynamic_pointer_cast<ArrayBase>(sCell->newCopy()));
390  shared_ptr<VGradArray> devStressField(new VGradArray(cells.getCountLevel1()));
391  devStressField->zero();
392  _structureFields.devStress.addArray(cells,devStressField);
393  shared_ptr<TArray> VMStressField(new TArray(cells.getCountLevel1()));
394  VMStressField->zero();
395  _structureFields.VMStress.addArray(cells,VMStressField);
396  shared_ptr<VGradArray> plasticStrainField(new VGradArray(cells.getCountLevel1()));
397  plasticStrainField->zero();
398  _structureFields.plasticStrain.addArray(cells,plasticStrainField);
399  shared_ptr<TArray> acCell(new TArray(cells.getCountLevel1()));
400  *acCell = _options.A;
401  _structureFields.creepConstant.addArray(cells, acCell);
402  }
403 
404  shared_ptr<TArray> rhoCell(new TArray(cells.getCountLevel1()));
405  *rhoCell = vc["density"];
406  _structureFields.density.addArray(cells,rhoCell);
407 
408  shared_ptr<TArray> etaCell(new TArray(cells.getCountLevel1()));
409  *etaCell = vc["eta"];
410  _structureFields.eta.addArray(cells,etaCell);
411 
412  shared_ptr<TArray> eta1Cell(new TArray(cells.getCountLevel1()));
413  *eta1Cell = vc["eta1"];
414  _structureFields.eta1.addArray(cells,eta1Cell);
415 
416  shared_ptr<TArray> alphaCell(new TArray(cells.getCountLevel1()));
417  *alphaCell = vc["alpha"];
418  _structureFields.alpha.addArray(cells,alphaCell);
419 
420  shared_ptr<TArray> tCell(new TArray(cells.getCountLevel1()));
421  *tCell = _options["operatingTemperature"];
422  _structureFields.temperature.addArray(cells,tCell);
423 
424  // compute values of deformation flux
425 
426  // store deformation flux at interfaces
427  /*
428  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
429  {
430  const FaceGroup& fg = *fgPtr;
431  const StorageSite& faces = fg.site;
432  shared_ptr<VectorT3Array> deformationFlux(new VectorT3Array(faces.getCount()));
433  deformationFlux->zero();
434  _structureFields.deformationFlux.addArray(faces,deformationFlux);
435  }
436  */
437 
438  // store deformation flux at boundary faces
439  foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
440  {
441  const FaceGroup& fg = *fgPtr;
442  const StorageSite& faces = fg.site;
443  shared_ptr<VectorT3Array> deformationFlux(new VectorT3Array(faces.getCount()));
444  deformationFlux->zero();
445  if (fg.groupType != "interior")
446  {
447  _structureFields.deformationFlux.addArray(faces,deformationFlux);
448  }
449  }
450 
451  }
452  _structureFields.eta.syncLocal();
453  _structureFields.eta1.syncLocal();
454  _structureFields.alpha.syncLocal();
455  _structureFields.density.syncLocal();
456 
457 
458 
459  _niters =0;
460  _initialDeformationNorm = MFRPtr();
461  }
462 
463  StructureBCMap& getBCMap() {return _bcMap;}
464  StructureVCMap& getVCMap() {return _vcMap;}
465  StructureModelOptions<T>& getOptions() {return _options;}
466 
467  void updateTime()
468  {
469  const int numMeshes = _meshes.size();
470  for (int n=0; n<numMeshes; n++)
471  {
472  const Mesh& mesh = *_meshes[n];
473 
474  const StorageSite& cells = mesh.getCells();
475  VectorT3Array& w =
476  dynamic_cast<VectorT3Array&>(_structureFields.deformation[cells]);
477  VectorT3Array& wN1 =
478  dynamic_cast<VectorT3Array&>(_structureFields.deformationN1[cells]);
479  VectorT3Array& wN2 =
480  dynamic_cast<VectorT3Array&>(_structureFields.deformationN2[cells]);
481  if (_options.timeDiscretizationOrder > 1)
482  {
483  VectorT3Array& wN3 =
484  dynamic_cast<VectorT3Array&>(_structureFields.deformationN3[cells]);
485  wN3 = wN2;
486  }
487  wN2 = wN1;
488  wN1 = w;
489  if(_options.variableTimeStep)
490  {
491  if (_options.timeDiscretizationOrder > 1)
492  {
493  _options.timeStepN2 = _options.timeStepN1;
494  }
495  _options.timeStepN1 = _options["timeStep"];
496  }
497 
498  }
499  }
500 
501  void creepInit()
502  {
503  const int numMeshes = _meshes.size();
504  for (int n=0; n<numMeshes; n++)
505  {
506  const Mesh& mesh = *_meshes[n];
507 
508  const StorageSite& cells = mesh.getCells();
509  VectorT3Array& w =
510  dynamic_cast<VectorT3Array&>(_structureFields.deformation[cells]);
511  VectorT3Array& ew =
512  dynamic_cast<VectorT3Array&>(_structureFields.elasticDeformation[cells]);
513  ew = w;
514 
515  _deformationGradientModel.compute();
516  const VGradArray& wGradCell =
517  dynamic_cast<const VGradArray&>(_structureFields.deformationGradient[cells]);
518 
519  const TArray& muCell =
520  dynamic_cast<const TArray&>(_structureFields.eta[cells]);
521 
522  const TArray& lambdaCell =
523  dynamic_cast<const TArray&>(_structureFields.eta1[cells]);
524 
525 
526  VGradArray& devStress =
527  dynamic_cast<VGradArray&>(_structureFields.devStress[cells]);
528  TArray& VMStress =
529  dynamic_cast<TArray&>(_structureFields.VMStress[cells]);
530 
531  const int nCells = cells.getCount();
532  const T onethird(1.0/3.0);
533  const T half(0.5);
534  const T twothirds(2.0/3.0);
535  const T six(6.0);
536  const T zero(0.0);
537 
538  const VectorT3Array& xc =
539  dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[cells]);
540  T xm(1.0);
541  T xv(1.0);
542  T xI(1.0);
543  const T one(1.0);
544  const T two(2.0);
545  const T three(3.0);
546  const T four(4.0);
547  const T eight(8.0);
548  const T twelve(12.0);
549  const T xl(400.e-6);
550  const T mid(200.e-6);
551  const T th(4.e-6);
552  const T midh(2.e-6);
553  xI = pow(th,three)/twelve;
554 
555  for(int k=0; k<nCells; k++)
556  {
557  const VGradType& wg = wGradCell[k];
558  VGradType wgPlusTranspose = wGradCell[k];
559  VGradType stress = wGradCell[k];
560 
561  for(int i=0;i<3;i++)
562  for(int j=0;j<3;j++)
563  wgPlusTranspose[i][j] += wg[j][i];
564 
565  stress[0][0] = wgPlusTranspose[0][0]*muCell[k]+
566  (wg[0][0]+wg[1][1]+wg[2][2])*lambdaCell[k];
567  stress[0][1] = wgPlusTranspose[0][1]*muCell[k];
568  stress[0][2] = wgPlusTranspose[0][2]*muCell[k];
569 
570  stress[1][0] = wgPlusTranspose[1][0]*muCell[k];
571  stress[1][1] = wgPlusTranspose[1][1]*muCell[k]+
572  (wg[0][0]+wg[1][1]+wg[2][2])*lambdaCell[k];
573  stress[1][2] = wgPlusTranspose[1][2]*muCell[k];
574 
575  stress[2][0] = wgPlusTranspose[2][0]*muCell[k];
576  stress[2][1] = wgPlusTranspose[2][1]*muCell[k];
577  stress[2][2] = wgPlusTranspose[2][2]*muCell[k]+
578  (wg[0][0]+wg[1][1]+wg[2][2])*lambdaCell[k];
579 
580  if (mesh.getDimension() == 2)
581  {
582  stress[2][0] = zero;
583  stress[2][1] = zero;
584  stress[2][2] = zero;
585  }
586 
587  /*
588  if(xc[k][0]<=mid)
589  {
590  xm = -xl*two/eight+two*xc[k][0]/two;
591  xv = one;
592  }
593  else
594  {
595  xm = (two/eight)*(three*xl-four*xc[k][0]);
596  xv = -one;
597  }
598 
599  stress[0][0] = -xm*(xc[k][1]-midh)/xI;
600  stress[0][1] = (three/two)*(xv/th)*(pow(midh,two)-
601  pow((xc[k][1]-midh),two))/pow(midh,two);
602  stress[0][2] = zero;
603  stress[1][0] = (three/two)*(xv/th)*(pow(midh,two)-
604  pow((xc[k][1]-midh),two))/pow(midh,two);
605  stress[1][1] = zero;
606  stress[1][2] = zero;
607  */
608 
609  devStress[k] = stress;
610  const T trace = stress[0][0]+stress[1][1]+stress[2][2];
611  (devStress[k])[0][0] = (devStress[k])[0][0] - onethird*trace;
612  (devStress[k])[1][1] = (devStress[k])[1][1] - onethird*trace;
613  (devStress[k])[2][2] = (devStress[k])[2][2] - onethird*trace;
614 
615  VMStress[k] = sqrt(half*(pow((stress[0][0]-stress[1][1]),2.0) +
616  pow((stress[1][1]-stress[2][2]),2.0) +
617  pow((stress[2][2]-stress[0][0]),2.0) +
618  six*(pow((stress[0][1]),2.0) +
619  pow((stress[1][2]),2.0) +
620  pow((stress[2][0]),2.0))));
621  }
622 
623  }
624  }
625 
626  void computeVMStress()
627  {
628  const int numMeshes = _meshes.size();
629  for (int n=0; n<numMeshes; n++)
630  {
631  const Mesh& mesh = *_meshes[n];
632 
633  const StorageSite& cells = mesh.getCells();
634  VectorT3Array& w =
635  dynamic_cast<VectorT3Array&>(_structureFields.deformation[cells]);
636  _deformationGradientModel.compute();
637  const VGradArray& wGradCell =
638  dynamic_cast<const VGradArray&>(_structureFields.deformationGradient[cells]);
639 
640  const TArray& muCell =
641  dynamic_cast<const TArray&>(_structureFields.eta[cells]);
642 
643  const TArray& lambdaCell =
644  dynamic_cast<const TArray&>(_structureFields.eta1[cells]);
645 
646  VGradArray& devStress =
647  dynamic_cast<VGradArray&>(_structureFields.devStress[cells]);
648  const VGradArray& plasticStrainCell =
649  dynamic_cast<VGradArray&>(_structureFields.plasticStrain[cells]);
650  TArray& VMStress =
651  dynamic_cast<TArray&>(_structureFields.VMStress[cells]);
652 
653  const int nCells = cells.getCount();
654  const T onethird(1.0/3.0);
655  const T half(0.5);
656  const T twothirds(2.0/3.0);
657  const T two(2.0);
658  const T three(3.0);
659  const T six(6.0);
660  const T zero(0.0);
661 
662  if (mesh.getDimension() == 3)
663  {
664  for(int k=0; k<nCells; k++)
665  {
666  const VGradType& wg = wGradCell[k];
667  VGradType wgPlusTranspose = wGradCell[k];
668  VGradType stress = wGradCell[k];
669  const VGradType pS = plasticStrainCell[k];
670 
671  for(int i=0;i<3;i++)
672  for(int j=0;j<3;j++)
673  wgPlusTranspose[i][j] += wg[j][i];
674 
675  stress[0][0] = wgPlusTranspose[0][0]*muCell[k]-
676  two*pS[0][0]*muCell[k]+
677  (wg[0][0]+wg[1][1]+wg[2][2])*lambdaCell[k];
678  stress[0][1] = wgPlusTranspose[0][1]*muCell[k]-
679  two*pS[0][1]*muCell[k];
680  stress[0][2] = wgPlusTranspose[0][2]*muCell[k]-
681  two*pS[0][2]*muCell[k];
682 
683  stress[1][0] = wgPlusTranspose[1][0]*muCell[k]-
684  two*pS[1][0]*muCell[k];
685  stress[1][1] = wgPlusTranspose[1][1]*muCell[k]-
686  two*pS[1][1]*muCell[k]+
687  (wg[0][0]+wg[1][1]+wg[2][2])*lambdaCell[k];
688  stress[1][2] = wgPlusTranspose[1][2]*muCell[k]-
689  two*pS[1][2]*muCell[k];
690 
691  stress[2][0] = wgPlusTranspose[2][0]*muCell[k]-
692  two*pS[2][0]*muCell[k];
693  stress[2][1] = wgPlusTranspose[2][1]*muCell[k]-
694  two*pS[2][1]*muCell[k];
695  stress[2][2] = wgPlusTranspose[2][2]*muCell[k]-
696  two*pS[2][2]*muCell[k]+
697  (wg[0][0]+wg[1][1]+wg[2][2])*lambdaCell[k];
698 
699  if(_options.creepModel!=3)
700  {
701  stress[0][0]-=(pS[0][0]+pS[1][1]+pS[2][2])*lambdaCell[k];
702  stress[1][1]-=(pS[0][0]+pS[1][1]+pS[2][2])*lambdaCell[k];
703  stress[2][2]-=(pS[0][0]+pS[1][1]+pS[2][2])*lambdaCell[k];
704  }
705 
706  devStress[k] = stress;
707  const T trace = stress[0][0]+stress[1][1]+stress[2][2];
708  (devStress[k])[0][0] = (devStress[k])[0][0] - onethird*trace;
709  (devStress[k])[1][1] = (devStress[k])[1][1] - onethird*trace;
710  (devStress[k])[2][2] = (devStress[k])[2][2] - onethird*trace;
711 
712  VMStress[k] = sqrt(half*(pow((stress[0][0]-stress[1][1]),2.0) +
713  pow((stress[1][1]-stress[2][2]),2.0) +
714  pow((stress[2][2]-stress[0][0]),2.0) +
715  six*(pow((stress[0][1]),2.0) +
716  pow((stress[1][2]),2.0) +
717  pow((stress[2][0]),2.0))));
718  }
719  }
720  else
721  {
722  for(int k=0; k<nCells; k++)
723  {
724  const VGradType& wg = wGradCell[k];
725  VGradType wgPlusTranspose = wGradCell[k];
726  VGradType stress = wGradCell[k];
727  const VGradType pS = plasticStrainCell[k];
728 
729  for(int i=0;i<3;i++)
730  for(int j=0;j<3;j++)
731  wgPlusTranspose[i][j] += wg[j][i];
732 
733  stress[0][0] = wgPlusTranspose[0][0]*muCell[k]-
734  two*pS[0][0]*muCell[k]+
735  (wg[0][0]+wg[1][1]+wg[2][2])*lambdaCell[k];
736  stress[0][1] = wgPlusTranspose[0][1]*muCell[k]-
737  two*pS[0][1]*muCell[k];
738  stress[0][2] = wgPlusTranspose[0][2]*muCell[k]-
739  two*pS[0][2]*muCell[k];
740 
741  stress[1][0] = wgPlusTranspose[1][0]*muCell[k]-
742  two*pS[1][0]*muCell[k];
743  stress[1][1] = wgPlusTranspose[1][1]*muCell[k]-
744  two*pS[1][1]*muCell[k]+
745  (wg[0][0]+wg[1][1]+wg[2][2])*lambdaCell[k];
746  stress[1][2] = wgPlusTranspose[1][2]*muCell[k]-
747  two*pS[1][2]*muCell[k];
748 
749 
750  if(_options.creepModel!=3)
751  {
752  stress[0][0]+=(pS[2][2])*lambdaCell[k];
753  stress[1][1]+=(pS[2][2])*lambdaCell[k];
754  }
755 
756  stress[2][0] = zero;
757  stress[2][1] = zero;
758  stress[2][2] = zero;
759 
760  devStress[k] = stress;
761  const T trace = stress[0][0]+stress[1][1]+stress[2][2];
762  (devStress[k])[0][0] = (devStress[k])[0][0] - onethird*trace;
763  (devStress[k])[1][1] = (devStress[k])[1][1] - onethird*trace;
764  (devStress[k])[2][2] = (devStress[k])[2][2] - onethird*trace;
765  //devStress[k][2][2] = zero;
766  /*
767  VMStress[k] = sqrt((three/two)*(pow(devStress[k][0][0],2.0) +
768  pow(devStress[k][1][1],2.0) +
769  pow(devStress[k][2][2],2.0) +
770  two*(pow((devStress[k][0][1]),2.0) +
771  pow((devStress[k][1][2]),2.0) +
772  pow((devStress[k][2][0]),2.0))));
773  */
774  VMStress[k] = sqrt(half*(pow((stress[0][0]-stress[1][1]),2.0) +
775  pow((stress[1][1]-stress[2][2]),2.0) +
776  pow((stress[2][2]-stress[0][0]),2.0) +
777  six*(pow((stress[0][1]),2.0) +
778  pow((stress[1][2]),2.0) +
779  pow((stress[2][0]),2.0))));
780 
781  }
782  }
783  }
784  }
785 
786  void initDeformationLinearization(LinearSystem& ls)
787  {
788 
789  const int numMeshes = _meshes.size();
790  for (int n=0; n<numMeshes; n++)
791  {
792  const Mesh& mesh = *_meshes[n];
793 
794  const StorageSite& cells = mesh.getCells();
795  MultiField::ArrayIndex wIndex(&_structureFields.deformation,&cells);
796 
797  ls.getX().addArray(wIndex,_structureFields.deformation.getArrayPtr(cells));
798 
799  const CRConnectivity& cellCells2 = mesh.getCellCells2();
800 
801  shared_ptr<Matrix> m(new CRMatrix<DiagTensorT3,DiagTensorT3,VectorT3>(cellCells2));
802 
803  ls.getMatrix().addMatrix(wIndex,wIndex,m);
804 
805  }
806  }
807 
808  void applyBC(LinearSystem& ls, bool explicitMode)
809  {
810  bool allNeumann = true;
811 
812  const int numMeshes = _meshes.size();
813  for (int n=0; n<numMeshes; n++)
814  {
815  const Mesh& mesh = *_meshes[n];
816 
817  foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
818  {
819  const FaceGroup& fg = *fgPtr;
820  if (fg.groupType != "interior")
821  {
822  const StorageSite& faces = fg.site;
823  const int nFaces = faces.getCount();
824  const TArray& faceAreaMag =
825  dynamic_cast<const TArray&>(_geomFields.areaMag[faces]);
826  const StructureBC<T>& bc = *_bcMap[fg.id];
827 
828 
829  StructureBCS<VectorT3,DiagTensorT3,DiagTensorT3> gbc(faces,mesh,
830  _geomFields,
831  _structureFields.deformation,
832  ls.getMatrix(), ls.getX(), ls.getB(),
833  explicitMode);
834 
835  VectorT3 fluxB(NumTypeTraits<VectorT3>::getZero());
836 
837  if (bc.bcType == "SpecifiedDeformation")
838  {
839  FloatValEvaluator<VectorT3>
840  bDeformation(bc.getVal("specifiedXDeformation"),
841  bc.getVal("specifiedYDeformation"),
842  bc.getVal("specifiedZDeformation"),
843  faces);
844  fluxB = gbc.applyDirichletBC(bDeformation);
845 
846  allNeumann = false;
847  }
848  else if (bc.bcType == "Symmetry")
849  {
850  gbc.applySymmetryBC();
851  allNeumann = false;
852  }
853  else if (bc.bcType == "Interface")
854  {
855  gbc.applyInterfaceBC();
856  }
857  else if (bc.bcType == "ZeroDerivative")
858  {
859  gbc.applyZeroDerivativeBC();
860  }
861  else if (bc.bcType == "SpecifiedTraction")
862  {
863  TractionValEvaluator<T>
864  bTraction(bc.getVal("specifiedXXTraction"),
865  bc.getVal("specifiedXYTraction"),
866  bc.getVal("specifiedXZTraction"),
867  bc.getVal("specifiedYXTraction"),
868  bc.getVal("specifiedYYTraction"),
869  bc.getVal("specifiedYZTraction"),
870  bc.getVal("specifiedZXTraction"),
871  bc.getVal("specifiedZYTraction"),
872  bc.getVal("specifiedZZTraction"),
873  _geomFields,
874  faces);
875  for(int f=0; f<nFaces; f++)
876  {
877 
878  fluxB += gbc.applyNeumannBC(f,bTraction[f]);
879 
880  }
881  }
882  else if (bc.bcType == "SpecifiedForce")
883  {
884  FloatValEvaluator<VectorT3>
885  bForce(bc.getVal("specifiedXForce"),
886  bc.getVal("specifiedYForce"),
887  bc.getVal("specifiedZForce"),
888  faces);
889  for(int f=0; f<nFaces; f++)
890  {
891 
892  fluxB += gbc.applyNeumannBC(f,bForce[f]/faceAreaMag[f]);
893 
894  }
895  }
896  else if (bc.bcType == "SpecifiedDistForce")
897  {
898  FloatValEvaluator<VectorT3>
899  bDistForce(bc.getVal("specifiedXDistForce"),
900  bc.getVal("specifiedYDistForce"),
901  bc.getVal("specifiedZDistForce"),
902  faces);
903  for(int f=0; f<nFaces; f++)
904  {
905  fluxB += gbc.applyNeumannBC(f,bDistForce[f]);
906 
907  }
908  }
909  else
910  throw CException(bc.bcType + " not implemented for StructureModel");
911  //cout << "force sum for " << fg.id << " = " << fluxB << endl;
912  }
913  }
914  /*
915  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
916  {
917  const FaceGroup& fg = *fgPtr;
918  const StorageSite& faces = fg.site;
919  GenericBCS<VectorT3,DiagTensorT3,DiagTensorT3> gbc(faces,mesh,
920  _geomFields,
921  _structureFields.deformation,
922  _structureFields.deformationFlux,
923  ls.getMatrix(), ls.getX(), ls.getB());
924 
925  gbc.applyInterfaceBC();
926  }
927  */
928  }
929 
930 
931 #ifdef FVM_PARALLEL
932  int count = 1;
933  int allNeumannInt = int( allNeumann);
934  MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &allNeumannInt, count, MPI::INT, MPI::PROD);
935  allNeumann = bool(allNeumannInt);
936 #endif
937 
938 
939  if(allNeumann && !explicitMode)
940  {
941  const Mesh& mesh = *_meshes[0];
942  const StorageSite& cells = mesh.getCells();
943  MultiField& b = ls.getB();
944  MultiFieldMatrix& matrix = ls.getMatrix();
945  MultiField::ArrayIndex wIndex(&_structureFields.deformation,&cells);
946  VectorT3Array& rCell = dynamic_cast<VectorT3Array&>(b[wIndex]);
947  VectorT3Array& w = dynamic_cast<VectorT3Array&>
948  (_structureFields.deformation[cells]);
949  VVMatrix& vvMatrix =
950  dynamic_cast<VVMatrix&>(matrix.getMatrix(wIndex,wIndex));
951  rCell[0] = T(0);
952  w[0] = T(0);
953  vvMatrix.setDirichlet(0);
954  }
955 
956  }
957 
958 
959  void linearizeDeformation(LinearSystem& ls, bool explicitMode)
960  {
961  _deformationGradientModel.compute();
962  DiscrList discretizations;
963 
964  if(!_options.creep)
965  {
966  shared_ptr<Discretization>
967  sd(new StructureSourceDiscretization<T,DiagTensorT3,DiagTensorT3>
968  (_meshes,_geomFields,
969  _structureFields.deformation,
970  _structureFields.eta,
971  _structureFields.eta1,
972  _structureFields.alpha,
973  _structureFields.deformationGradient,
974  _structureFields.temperature,
975  _options["operatingTemperature"],
976  _options["residualXXStress"],
977  _options["residualYYStress"],
978  _options["residualZZStress"],
979  _options.thermo,
980  _options.residualStress));
981 
982  // shared_ptr<Discretization>
983  // bfd(new SourceDiscretization<VectorT3>
984  // (_meshes,_geomFields,
985  // _structureFields.deformation,
986  // _structureFields.bodyForce));
987 
988  discretizations.push_back(sd);
989  //discretizations.push_back(bfd);
990  }
991 
992  if (_options.creep)
993  {
994  shared_ptr<Discretization>
995  scd(new StructurePlasticDiscretization<T,DiagTensorT3,DiagTensorT3>
996  (_meshes,_geomFields,
997  _structureFields.deformation,
998  _structureFields.eta,
999  _structureFields.eta1,
1000  _structureFields.alpha,
1001  _structureFields.deformationGradient,
1002  _structureFields.temperature,
1003  _options["operatingTemperature"],
1004  _options["residualXXStress"],
1005  _options["residualYYStress"],
1006  _options["residualZZStress"],
1007  _options.thermo,
1008  _options.residualStress,
1009  _structureFields.devStress,
1010  _structureFields.VMStress,
1011  _structureFields.plasticStrain,
1012  _structureFields.creepConstant,
1013  _options.A,
1014  _options.B,
1015  _options.m,
1016  _options.n,
1017  _options.Sy0,
1018  _options["timeStep"],
1019  _options.creepModel));
1020 
1021  discretizations.push_back(scd);
1022  }
1023 
1024  if (_options.transient && !explicitMode)
1025  {
1026  shared_ptr<Discretization>
1027  td(new TimeDerivativeStructureDiscretization
1028  <VectorT3,DiagTensorT3,DiagTensorT3>
1029  (_meshes,_geomFields,
1030  _structureFields.deformation,
1031  _structureFields.deformationN1,
1032  _structureFields.deformationN2,
1033  _structureFields.deformationN3,
1034  _structureFields.density,
1035  _structureFields.volume0,
1036  _options.variableTimeStep,
1037  _options["timeStep"],
1038  _options.timeStepN1,
1039  _options.timeStepN2));
1040 
1041  discretizations.push_back(td);
1042  }
1043 
1044  /*
1045  shared_ptr<Discretization>
1046  ibm(new GenericIBDiscretization<VectorT3,DiagTensorT3,T>
1047  (_meshes,_geomFields,_structureFields.deformation));
1048 
1049  discretizations.push_back(ibm);
1050  */
1051  Linearizer linearizer;
1052 
1053  linearizer.linearize(discretizations,_meshes,ls.getMatrix(),
1054  ls.getX(), ls.getB());
1055 
1056  if (!explicitMode)
1057  applyBC(ls,explicitMode);
1058 
1059 #if 0
1060  shared_ptr<Discretization>
1061  ud(new Underrelaxer<VectorT3,DiagTensorT3,DiagTensorT3>
1062  (_meshes,_structureFields.deformation,
1063  _options["deformationURF"]));
1064 
1065  DiscrList discretizations2;
1066  discretizations2.push_back(ud);
1067 
1068  linearizer.linearize(discretizations2,_meshes,ls.getMatrix(),
1069  ls.getX(), ls.getB());
1070 #endif
1071 
1072  }
1073 
1074 
1075  MFRPtr solveDeformation()
1076  {
1077  LinearSystem ls;
1078  initDeformationLinearization(ls);
1079  ls.initAssembly();
1080  linearizeDeformation(ls,false);
1081  ls.initSolve();
1082  //AMG solver(ls);
1083  MFRPtr rNorm = _options.getDeformationLinearSolver().solve(ls);
1084  if (!_initialDeformationNorm) _initialDeformationNorm = rNorm;
1085  _options.getDeformationLinearSolver().cleanup();
1086  ls.postSolve();
1087  ls.updateSolution();
1088  //postStructureSolve(ls);
1089  return rNorm;
1090  }
1091 
1092  void initExplicitAdvance()
1093  {
1094  this->_lsK = shared_ptr<LinearSystem>(new LinearSystem());
1095  LinearSystem& ls = *(this->_lsK);
1096  initDeformationLinearization(ls);
1097  ls.initAssembly();
1098  linearizeDeformation(ls,true);
1099  // add delta explicitly since we won't be calling initSolve
1100  ls.replaceDelta(dynamic_pointer_cast<MultiField>(ls.getX().newClone()));
1101  ls.replaceResidual(dynamic_pointer_cast<MultiField>(ls.getX().newClone()));
1102 
1103  }
1104 
1105  void finishExplicitAdvance()
1106  {
1107  this->_lsK = shared_ptr<LinearSystem>();
1108  }
1109 
1110  void advanceExplicit(const int nSteps, const double deltaT)
1111  {
1112  const int nMeshes = this->_meshes.size();
1113  LinearSystem& ls = *(this->_lsK);
1114 
1115  TimeDerivativeStructureDiscretization
1116  <VectorT3,DiagTensorT3,DiagTensorT3>
1117  td(_meshes,_geomFields,
1118  _structureFields.deformation,
1119  _structureFields.deformationN1,
1120  _structureFields.deformationN2,
1121  _structureFields.deformationN3,
1122  _structureFields.density,
1123  _structureFields.volume0,
1124  _options.variableTimeStep,
1125  deltaT,
1126  _options.timeStepN1,
1127  _options.timeStepN2);
1128 
1129 
1130  for(int n=0; n<nSteps; n++)
1131  {
1132 
1133  ls.getMatrix().multiply(ls.getB(), ls.getX());
1134 
1135 
1136  applyBC(ls, true);
1137 
1138  ls.getDelta().zero();
1139 
1140  for(int n=0; n<nMeshes; n++)
1141  {
1142  const Mesh& mesh = *(this->_meshes[n]);
1143  td.explicitAdvance(mesh,ls.getX(),ls.getB(),ls.getDelta());
1144  }
1145 
1146  // this will update the boundary delta's
1147  ls.postSolve();
1148 
1149  ls.updateSolution();
1150 
1151  updateTime();
1152  }
1153  }
1154 
1155 #if 0
1156 
1157  void postStructureSolve(LinearSystem& ls)
1158  {
1159  MultiField& sField = ls.getDelta();
1160 
1161  const int numMeshes = _meshes.size();
1162  for(int n=0;n<numMeshes;n++)
1163  {
1164  const Mesh& mesh = *_meshes[n];
1165 
1166  const StorageSite& cells = mesh.getCells();
1167 
1168  MultiField::ArrayIndex sIndex(&_structureFields.deformation,&cells);
1169  VectorT3Array& w = dynamic_cast<VectorT3Array&>
1170  (_structureFields.deformation[cells]);
1171  const VectorT3Array& ww = dynamic_cast<const VectorT3Array&>
1172  (sField[sIndex]);
1173  //const T deformationURF(_options["deformationURF"]);
1174 
1175  const int nCells = cells.getCountLevel1();
1176  for(int c=0;c<nCells;c++)
1177  {
1178  w[c] += ww[c];
1179  }
1180  }
1181  }
1182 
1183 #endif
1184 
1185  bool advance(const int niter)
1186  {
1187 
1188  for(int n=0; n<niter; n++)
1189  {
1190  MFRPtr dNorm = solveDeformation();
1191 
1192  if (_niters < 5)
1193  {
1194  _initialDeformationNorm->setMax(*dNorm);
1195  }
1196 
1197  MFRPtr dNormRatio(dNorm->normalize(*_initialDeformationNorm));
1198 
1199 
1200 #ifdef FVM_PARALLEL
1201  if ( MPI::COMM_WORLD.Get_rank() == 0 ){ //only root process
1202  if (_options.printNormalizedResiduals)
1203  cout << _niters << ": " << *dNormRatio << endl;
1204  else
1205  cout << _niters << ": " << *dNorm << endl;
1206  }
1207 #endif
1208 
1209 #ifndef FVM_PARALLEL
1210  if (_options.printNormalizedResiduals)
1211  cout << _niters << ": " << *dNormRatio << endl;
1212  else
1213  cout << _niters << ": " << *dNorm << endl;
1214 #endif
1215 
1216 
1217  _niters++;
1218  if (*dNormRatio < _options.deformationTolerance)
1219  return true;
1220  }
1221  return false;
1222  }
1223 
1224 #if 0
1225  void printBCs()
1226  {
1227  foreach(typename StructureBCMap::value_type& pos, _bcMap)
1228  {
1229  cout << "Face Group " << pos.first << ":" << endl;
1230  cout << " bc type " << pos.second->bcType << endl;
1231  foreach(typename StructureBC<T>::value_type& vp, *pos.second)
1232  {
1233  cout << " " << vp.first << " " << vp.second.constant << endl;
1234  }
1235  }
1236  }
1237 
1238 
1239  VectorT3 getDeformationFluxIntegral(const Mesh& mesh, const int faceGroupId)
1240  {
1241  VectorT3 r(VectorT3::getZero());
1242  bool found = false;
1243  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
1244  {
1245  const FaceGroup& fg = *fgPtr;
1246  if (fg.id == faceGroupId)
1247  {
1248  const StorageSite& faces = fg.site;
1249  const int nFaces = faces.getCount();
1250  const VectorT3Array& deformationFlux =
1251  dynamic_cast<const VectorT3Array&>(_structureFields.deformationFlux[faces]);
1252  for(int f=0; f<nFaces; f++)
1253  r += deformationFlux[f];
1254  found=true;
1255  }
1256  }
1257  if (!found)
1258  throw CException("getDeformationFluxIntegral: invalid faceGroupID");
1259  return r;
1260  }
1261 
1262  VectorT3 getDeformationDerivativeIntegral(const Mesh& mesh)
1263  {
1264  VectorT3 r(VectorT3::getZero());
1265  const StorageSite& cells = mesh.getCells();
1266  const int nCells = cells.getSelfCount();
1267 
1268  const TArray& density =
1269  dynamic_cast<const TArray&>(_structureFields.density[cells]);
1270  const VectorT3Array& w =
1271  dynamic_cast<const VectorT3Array&>(_structureFields.deformation[cells]);
1272  const VectorT3Array& wN1 =
1273  dynamic_cast<const VectorT3Array&>(_structureFields.deformationN1[cells]);
1274  const VectorT3Array& wN2 =
1275  dynamic_cast<const VectorT3Array&>(_structureFields.deformationN2[cells]);
1276 
1277  const TArray& volume =
1278  dynamic_cast<const TArray&>(_geomFields.volume[cells]);
1279 
1280  const T deltaT = _options["timeStep"];
1281 
1282  T onePointFive(1.5);
1283  T two(2.0);
1284  T pointFive(0.5);
1285 
1286  for(int c=0; c<nCells; c++)
1287  {
1288  const T rhoVbydT = density[c]*volume[c]/deltaT;
1289  r += rhoVbydT*(onePointFive*w[c]- two*wN1[c]
1290  + pointFive*wN2[c]);
1291  }
1292  return r;
1293  }
1294 
1295  boost::shared_ptr<ArrayBase> getStressTensor(const Mesh& mesh, const ArrayBase& gcellIds)
1296  {
1297  typedef Array<StressTensor<T> > StressTensorArray;
1298 
1299  const StorageSite& cells = mesh.getCells();
1300 
1301  const Array<int>& cellIds = dynamic_cast<const Array<int> &>(gcellIds);
1302  const int nCells = cellIds.getLength();
1303 
1304  _deformationGradientModel.compute();
1305 
1306  const VGradArray& wGrad =
1307  dynamic_cast<const VGradArray&>(_structureFields.deformationGradient[cells]);
1308 
1309  /* const TArray& pCell =
1310  dynamic_cast<const TArray&>(_flowFields.pressure[cells]);
1311  */
1312  const TArray& eta = dynamic_cast<const TArray&>(_structureFields.eta[cells]);
1313  const TArray& eta1 = dynamic_cast<const TArray&>(_structureFields.eta1[cells]);
1314 
1315  boost::shared_ptr<StressTensorArray> stressTensorPtr( new StressTensorArray(nCells));
1316  StressTensorArray& stressTensor = *stressTensorPtr;
1317 
1318  for(int n=0; n<nCells; n++)
1319  {
1320  const int c = cellIds[n];
1321  const VGradType& wg = wGrad[c];
1322  VGradType wgPlusTranspose = wGrad[c];
1323 
1324  for(int i=0;i<3;i++)
1325  for(int j=0;j<3;j++)
1326  wgPlusTranspose[i][j] += wg[j][i];
1327 
1328  stressTensor[n][0] = wgPlusTranspose[0][0]*eta[c]+
1329  (wg[0][0]+wg[1][1]+wg[2][2])*eta1[c];
1330  stressTensor[n][1] = wgPlusTranspose[1][1]*eta[c]+
1331  (wg[0][0]+wg[1][1]+wg[2][2])*eta1[c];
1332  stressTensor[n][2] = wgPlusTranspose[2][2]*eta[c]+
1333  (wg[0][0]+wg[1][1]+wg[2][2])*eta1[c];
1334  stressTensor[n][3] = wgPlusTranspose[0][1]*eta[c];
1335  stressTensor[n][4] = wgPlusTranspose[1][2]*eta[c];
1336  stressTensor[n][5] = wgPlusTranspose[2][0]*eta[c];
1337  }
1338 
1339  return stressTensorPtr;
1340  }
1341 
1342 #endif
1343 
1344  void getTraction(const Mesh& mesh)
1345  {
1346  const StorageSite& cells = mesh.getCells();
1347 
1348  const int nCells = cells.getCount();
1349 
1350  shared_ptr<VectorT3Array> tractionXPtr(new VectorT3Array(nCells));
1351  tractionXPtr->zero();
1352  _structureFields.tractionX.addArray(cells,tractionXPtr);
1353  VectorT3Array& tractionX = *tractionXPtr;
1354 
1355  shared_ptr<VectorT3Array> tractionYPtr(new VectorT3Array(nCells));
1356  tractionYPtr->zero();
1357  _structureFields.tractionY.addArray(cells,tractionYPtr);
1358  VectorT3Array& tractionY = *tractionYPtr;
1359 
1360  shared_ptr<VectorT3Array> tractionZPtr(new VectorT3Array(nCells));
1361  tractionZPtr->zero();
1362  _structureFields.tractionZ.addArray(cells,tractionZPtr);
1363  VectorT3Array& tractionZ = *tractionZPtr;
1364 
1365  _deformationGradientModel.compute();
1366 
1367  const VGradArray& wGrad =
1368  dynamic_cast<const VGradArray&>(_structureFields.deformationGradient[cells]);
1369 
1370  const VGradArray& plasticStrainCell =
1371  dynamic_cast<VGradArray&>(_structureFields.plasticStrain[cells]);
1372 
1373  const TArray& eta = dynamic_cast<const TArray&>(_structureFields.eta[cells]);
1374  const TArray& eta1 = dynamic_cast<const TArray&>(_structureFields.eta1[cells]);
1375  const TArray& alpha = dynamic_cast<const TArray&>(_structureFields.alpha[cells]);
1376 
1377  const TArray& temperature = dynamic_cast<const TArray&>(_structureFields.temperature[cells]);
1378 
1379  const T two(2.0);
1380  const T three(3.0);
1381  const T zero(0.0);
1382 
1383  for(int n=0; n<nCells; n++)
1384  {
1385  const VGradType& wg = wGrad[n];
1386  VGradType wgPlusTranspose = wGrad[n];
1387  const VGradType& pS = plasticStrainCell[n];
1388  for(int i=0;i<3;i++)
1389  for(int j=0;j<3;j++)
1390  wgPlusTranspose[i][j] += wg[j][i];
1391 
1392  tractionX[n][0] = wgPlusTranspose[0][0]*eta[n]+
1393  (wg[0][0]+wg[1][1]+wg[2][2])*eta1[n];
1394  tractionX[n][1] = wgPlusTranspose[0][1]*eta[n];
1395  tractionX[n][2] = wgPlusTranspose[0][2]*eta[n];
1396 
1397  tractionY[n][0] = wgPlusTranspose[1][0]*eta[n];
1398  tractionY[n][1] = wgPlusTranspose[1][1]*eta[n]+
1399  (wg[0][0]+wg[1][1]+wg[2][2])*eta1[n];
1400  tractionY[n][2] = wgPlusTranspose[1][2]*eta[n];
1401 
1402  tractionZ[n][0] = wgPlusTranspose[2][0]*eta[n];
1403  tractionZ[n][1] = wgPlusTranspose[2][1]*eta[n];
1404  tractionZ[n][2] = wgPlusTranspose[2][2]*eta[n]+
1405  (wg[0][0]+wg[1][1]+wg[2][2])*eta1[n];
1406 
1407  if(_options.residualStress)
1408  {
1409  tractionX[n][0] += _options["residualXXStress"];
1410  tractionY[n][1] += _options["residualYYStress"];
1411  tractionZ[n][2] += _options["residualZZStress"];
1412  }
1413 
1414  if(_options.thermo)
1415  {
1416  if(mesh.getDimension()==2)
1417  {
1418  tractionX[n][0] -= (two*eta1[n]+two*eta[n])*alpha[n]*(temperature[n]-_options["operatingTemperature"]);
1419  tractionY[n][1] -= (two*eta1[n]+two*eta[n])*alpha[n]*(temperature[n]-_options["operatingTemperature"]);
1420  }
1421  else
1422  {
1423  tractionX[n][0] -= (three*eta1[n]+two*eta[n])*alpha[n]*(temperature[n]-_options["operatingTemperature"]);
1424  tractionY[n][1] -= (three*eta1[n]+two*eta[n])*alpha[n]*(temperature[n]-_options["operatingTemperature"]);
1425  tractionZ[n][2] -= (three*eta1[n]+two*eta[n])*alpha[n]*(temperature[n]-_options["operatingTemperature"]);
1426  }
1427  }
1428 
1429  if(_options.creep)
1430  {
1431  tractionX[n][0]-=(two*eta[n]*pS[0][0] +
1432  (pS[0][0] + pS[1][1])*eta1[n]);
1433  tractionX[n][1]-=(two*eta[n]*pS[0][1]);
1434  tractionX[n][2]-=(two*eta[n]*pS[0][2]);
1435 
1436  tractionY[n][0]-=(two*eta[n]*pS[1][0]);
1437  tractionY[n][1]-=(two*eta[n]*pS[1][1] +
1438  (pS[0][0] + pS[1][1])*eta1[n]);
1439  tractionY[n][2]-=(two*eta[n]*pS[1][2]);
1440  }
1441 
1442  if (mesh.getDimension() == 2)
1443  {
1444  tractionZ[n] = zero;
1445  }
1446  }
1447  }
1448 
1449  void getStrain(const Mesh& mesh)
1450  {
1451  const StorageSite& cells = mesh.getCells();
1452 
1453  const int nCells = cells.getCount();
1454 
1455  shared_ptr<VectorT3Array> strainXPtr(new VectorT3Array(nCells));
1456  strainXPtr->zero();
1457  _structureFields.strainX.addArray(cells,strainXPtr);
1458  VectorT3Array& strainX = *strainXPtr;
1459 
1460  shared_ptr<VectorT3Array> strainYPtr(new VectorT3Array(nCells));
1461  strainYPtr->zero();
1462  _structureFields.strainY.addArray(cells,strainYPtr);
1463  VectorT3Array& strainY = *strainYPtr;
1464 
1465  shared_ptr<VectorT3Array> strainZPtr(new VectorT3Array(nCells));
1466  strainZPtr->zero();
1467  _structureFields.strainZ.addArray(cells,strainZPtr);
1468  VectorT3Array& strainZ = *strainZPtr;
1469 
1470  _deformationGradientModel.compute();
1471 
1472  const VGradArray& wGrad =
1473  dynamic_cast<const VGradArray&>(_structureFields.deformationGradient[cells]);
1474 
1475  const T half(1./2.);
1476 
1477  for(int n=0; n<nCells; n++)
1478  {
1479  const VGradType& wg = wGrad[n];
1480  VGradType wgPlusTranspose = wGrad[n];
1481 
1482  for(int i=0;i<3;i++)
1483  for(int j=0;j<3;j++)
1484  wgPlusTranspose[i][j] += wg[j][i];
1485 
1486  strainX[n][0] = half*wgPlusTranspose[0][0];
1487  strainX[n][1] = half*wgPlusTranspose[0][1];
1488  strainX[n][2] = half*wgPlusTranspose[0][2];
1489 
1490  strainY[n][0] = half*wgPlusTranspose[1][0];
1491  strainY[n][1] = half*wgPlusTranspose[1][1];
1492  strainY[n][2] = half*wgPlusTranspose[1][2];
1493 
1494  strainZ[n][0] = half*wgPlusTranspose[2][0];
1495  strainZ[n][1] = half*wgPlusTranspose[2][1];
1496  strainZ[n][2] = half*wgPlusTranspose[2][2];
1497  }
1498  }
1499 
1500  void getPlasticDiagStrain(const Mesh& mesh)
1501  {
1502  const StorageSite& cells = mesh.getCells();
1503 
1504  const int nCells = cells.getSelfCount();
1505 
1506  shared_ptr<VectorT3Array> plasticDiagStrainPtr(new VectorT3Array(nCells));
1507  plasticDiagStrainPtr->zero();
1508  _structureFields.plasticDiagStrain.addArray(cells,plasticDiagStrainPtr);
1509  VectorT3Array& plasticDiagStrain = *plasticDiagStrainPtr;
1510 
1511  const VGradArray& plasticStrain =
1512  dynamic_cast<const VGradArray&>(_structureFields.plasticStrain[cells]);
1513 
1514  for(int n=0; n<nCells; n++)
1515  {
1516  plasticDiagStrain[n][0]=plasticStrain[n][0][0];
1517  plasticDiagStrain[n][1]=plasticStrain[n][1][1];
1518  plasticDiagStrain[n][2]=plasticStrain[n][0][1];
1519  }
1520  }
1521 
1522 void updateForceOnBoundary(const StorageSite& faceSite, const ArrayBase& bforceA, const map<int,int>& commonFacesMap,
1523  ArrayBase& fxA, ArrayBase& fyA, ArrayBase& fzA)
1524 {
1525  //bforce came from fluid+elec side
1526  const VectorT3Array& bforce = dynamic_cast<const VectorT3Array&>(bforceA);
1527  //following will be updated
1528  TArray& fx = dynamic_cast<TArray&> (fxA);
1529  TArray& fy = dynamic_cast<TArray&> (fyA);
1530  TArray& fz = dynamic_cast<TArray&> (fzA);
1531  const int offset = faceSite.getOffset();
1532  for (int i = 0; i < faceSite.getCount(); i++ ){
1533  const int faceID = i + offset; //localface ID
1534  //commonFacesMap will get right index in bforce
1535  const int indx = commonFacesMap.find(faceID)->second;
1536  fx[i] = bforce[indx][0];
1537  fy[i] = bforce[indx][1];
1538  fz[i] = bforce[indx][2];
1539  }
1540 
1541 }
1542 
1543  /*
1544  VectorT3 getMomentumFluxIntegralonIBFaces(const Mesh& mesh)
1545  {
1546  VectorT3 r(VectorT3::getZero());
1547  const StorageSite& ibFaces = mesh.getIBFaces();
1548  const StorageSite& faces = mesh.getFaces();
1549  const Array<int>& ibFaceIndices = mesh.getIBFaceList();
1550  const int nibf = ibFaces.getCount();
1551  const VectorT3Array& deformationFlux =
1552  dynamic_cast<const VectorT3Array&>(_structureFields.deformationFlux[faces]);
1553  for ( int f = 0; f < nibf; f ++){
1554  const int ibFaceIndex = ibFaceIndices[f];
1555  r += deformationFlux[ibFaceIndex];
1556  }
1557  if (nibf == 0)
1558  throw CException("getDeformationFluxIntegralonIBFaces: no ibFaces found!");
1559  return r;
1560  }
1561 
1562 
1563  void printPressureIntegrals()
1564  {
1565  const int numMeshes = _meshes.size();
1566  for (int n=0; n<numMeshes; n++)
1567  {
1568  const Mesh& mesh = *_meshes[n];
1569  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
1570  {
1571  const FaceGroup& fg = *fgPtr;
1572 
1573  VectorT3 r(VectorT3::getZero());
1574 
1575  const StorageSite& faces = fg.site;
1576  const VectorT3Array& faceArea =
1577  dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
1578  const int nFaces = faces.getCount();
1579  const TArray& facePressure = dynamic_cast<const TArray&>(_flowFields.pressure[faces]);
1580  for(int f=0; f<nFaces; f++)
1581  r += faceArea[f]*facePressure[f];
1582 
1583  cout << "Mesh " << mesh.getID() << " faceGroup " << fg.id << " : " << r << endl;
1584  }
1585  }
1586  }
1587 
1588  */
1589 
1590 #if 0
1591 
1592  void printDeformationFluxIntegrals()
1593  {
1594  const int numMeshes = _meshes.size();
1595  for (int n=0; n<numMeshes; n++)
1596  {
1597  const Mesh& mesh = *_meshes[n];
1598  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
1599  {
1600  const FaceGroup& fg = *fgPtr;
1601 
1602  VectorT3 r(VectorT3::getZero());
1603 
1604  const StorageSite& faces = fg.site;
1605  const int nFaces = faces.getCount();
1606  const VectorT3Array& deformationFlux =
1607  dynamic_cast<const VectorT3Array&>(_structureFields.deformationFlux[faces]);
1608  for(int f=0; f<nFaces; f++)
1609  r += deformationFlux[f];
1610 
1611  cout << "Mesh " << mesh.getID() << " faceGroup " << fg.id << " : " << r << endl;
1612  }
1613  }
1614  }
1615 
1616 #endif
1617 
1618  /*
1619  void printMassFluxIntegrals()
1620  {
1621  const int numMeshes = _meshes.size();
1622  for (int n=0; n<numMeshes; n++)
1623  {
1624  const Mesh& mesh = *_meshes[n];
1625  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
1626  {
1627  const FaceGroup& fg = *fgPtr;
1628 
1629  T r(0.);
1630 
1631  const StorageSite& faces = fg.site;
1632  const int nFaces = faces.getCount();
1633  const TArray& massFlux = dynamic_cast<const TArray&>(_flowFields.massFlux[faces]);
1634  for(int f=0; f<nFaces; f++)
1635  r += massFlux[f];
1636 
1637  cout << "Mesh " << mesh.getID() << " faceGroup " << fg.id << " : " << r << endl;
1638  }
1639  }
1640  }
1641 
1642 #include "FlowModelPressureBC.h"
1643 #include "FlowModelVelocityBC.h"
1644 #include "FlowModelInterior.h"
1645 #include "FlowModelSlipJump.h"
1646 
1647  */
1648 
1649 
1650 private:
1651  const MeshList _meshes;
1652  const GeomFields& _geomFields;
1653  StructureFields& _structureFields;
1654 
1655  StructureBCMap _bcMap;
1656  StructureVCMap _vcMap;
1657 
1658  StructureModelOptions<T> _options;
1659  GradientModel<VectorT3> _deformationGradientModel;
1660  // GradientModel<T> _pressureGradientModel;
1661 
1662  MFRPtr _initialDeformationNorm;
1663 // MFRPtr _initialCoupledNorm;
1664  int _niters;
1665 
1666  shared_ptr<LinearSystem> _lsK;
1667  // shared_ptr<Field> _previousVelocity;
1668  // shared_ptr<Field> _momApField;
1669 
1670  // bool _useReferencePressure;
1671  // T _referencePP;
1672  //AMG _momSolver;
1673  //AMG _continuitySolver;
1674 };
1675 
1676 template<class T>
1677 StructureModel<T>::StructureModel(const GeomFields& geomFields,
1678  StructureFields& structureFields,
1679  const MeshList& meshes) :
1680  Model(meshes),
1681  _impl(new Impl(geomFields,structureFields,meshes))
1682 {
1683  logCtor();
1684 }
1685 
1686 
1687 template<class T>
1688 StructureModel<T>::~StructureModel()
1689 {
1690  logDtor();
1691 }
1692 
1693 template<class T>
1694 void
1695 StructureModel<T>::init()
1696 {
1697  _impl->init();
1698 }
1699 
1700 template<class T>
1701 typename StructureModel<T>::StructureBCMap&
1702 StructureModel<T>::getBCMap() {return _impl->getBCMap();}
1703 
1704 template<class T>
1705 typename StructureModel<T>::StructureVCMap&
1706 StructureModel<T>::getVCMap() {return _impl->getVCMap();}
1707 
1708 template<class T>
1709 StructureModelOptions<T>&
1710 StructureModel<T>::getOptions() {return _impl->getOptions();}
1711 
1712 #if 0
1713 
1714 template<class T>
1715 void
1716 StructureModel<T>::printBCs()
1717 {
1718  _impl->printBCs();
1719 }
1720 
1721 #endif
1722 
1723 template<class T>
1724 bool
1725 StructureModel<T>::advance(const int niter)
1726 {
1727  return _impl->advance(niter);
1728 }
1729 
1730 template<class T>
1731 void
1732 StructureModel<T>::advanceExplicit(const int nsteps, const double deltaT)
1733 {
1734  _impl->advanceExplicit(nsteps,deltaT);
1735 }
1736 
1737 template<class T>
1738 void
1739 StructureModel<T>::initExplicitAdvance()
1740 {
1741  return _impl->initExplicitAdvance();
1742 }
1743 
1744 
1745 template<class T>
1746 void
1747 StructureModel<T>::updateForceOnBoundary(const StorageSite& faceSite, const ArrayBase& bforceA, const map<int,int>& commonFacesMap,
1748  ArrayBase& fxA, ArrayBase& fyA, ArrayBase& fzA)
1749 
1750 {
1751  _impl->updateForceOnBoundary(faceSite, bforceA, commonFacesMap, fxA, fyA, fzA);
1752 ;
1753 }
1754 
1755 template<class T>
1756 void
1757 StructureModel<T>::finishExplicitAdvance()
1758 {
1759  return _impl->finishExplicitAdvance();
1760 }
1761 
1762 template<class T>
1763 void
1764 StructureModel<T>::creepInit()
1765 {
1766  _impl->creepInit();
1767 }
1768 
1769 template<class T>
1770 void
1771 StructureModel<T>::computeVMStress()
1772 {
1773  _impl->computeVMStress();
1774 }
1775 
1776 template<class T>
1777 void
1778 StructureModel<T>::getStrain(const Mesh& mesh)
1779 {
1780  return _impl->getStrain(mesh);
1781 }
1782 
1783 template<class T>
1784 void
1785 StructureModel<T>::getPlasticDiagStrain(const Mesh& mesh)
1786 {
1787  return _impl->getPlasticDiagStrain(mesh);
1788 }
1789 
1790 template<class T>
1791 void
1792 StructureModel<T>::updateTime()
1793 {
1794  _impl->updateTime();
1795 }
1796 
1797 #if 0
1798 
1799 template<class T>
1800 void
1801 StructureModel<T>::printDeformationFluxIntegrals()
1802 {
1803  _impl->printDeformationFluxIntegrals();
1804 }
1805 
1806 
1807 template<class T>
1808 Vector<T,3>
1809 StructureModel<T>::getDeformationFluxIntegral(const Mesh& mesh, const int faceGroupId)
1810 {
1811  return _impl->getDeformationFluxIntegral(mesh,faceGroupId);
1812 }
1813 
1814 template<class T>
1815 Vector<T,3>
1816 StructureModel<T>::getDeformationDerivativeIntegral(const Mesh& mesh)
1817 {
1818  return _impl->getDeformationDerivativeIntegral(mesh);
1819 }
1820 
1821 #endif
1822 
1823 template<class T>
1824 void
1825 StructureModel<T>::getTraction(const Mesh& mesh)
1826 {
1827  return _impl->getTraction(mesh);
1828 }
1829 
1830 #if 0
1831 
1832 template<class T>
1833 boost::shared_ptr<ArrayBase>
1834 StructureModel<T>::getStressTensor(const Mesh& mesh, const ArrayBase& cellIds)
1835 {
1836  return _impl->getStressTensor(mesh, cellIds);
1837 }
1838 
1839 #endif
StructureBCS::DiagArray
Array< Diag > DiagArray
Definition: StructureModel_impl.h:55
Mesh::getBoundaryFaceGroups
const FaceGroupList & getBoundaryFaceGroups() const
Definition: Mesh.h:187
StructureModel::advanceExplicit
void advanceExplicit(const int nsteps, const double deltaT)
Definition: StructureModel_impl.h:1732
StructureBCS::TArray
Array< T_Scalar > TArray
Definition: StructureModel_impl.h:46
StructureModel::Impl::getTraction
void getTraction(const Mesh &mesh)
Definition: StructureModel_impl.h:1344
StructureBCS::applyNeumannBC
X applyNeumannBC(const int f, const X &specifiedFlux) const
Definition: StructureModel_impl.h:120
StructureModel::Impl::VVMatrix
CRMatrix< DiagTensorT3, DiagTensorT3, VectorT3 > VVMatrix
Definition: StructureModel_impl.h:282
LinearSystem::initAssembly
void initAssembly()
Definition: LinearSystem.cpp:32
StructureModel::Impl::linearizeDeformation
void linearizeDeformation(LinearSystem &ls, bool explicitMode)
Definition: StructureModel_impl.h:959
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
StructureBCS::_explicitMode
const bool _explicitMode
Definition: StructureModel_impl.h:266
StructureBCS::_assembler
CCAssembler & _assembler
Definition: StructureModel_impl.h:258
StructureBCS::StructureBCS
StructureBCS(const StorageSite &faces, const Mesh &mesh, const GeomFields &geomFields, Field &varField, MultiFieldMatrix &matrix, MultiField &xField, MultiField &rField, const bool explicitMode)
Definition: StructureModel_impl.h:62
StructureBCS::_faceAreaMag
const TArray & _faceAreaMag
Definition: StructureModel_impl.h:263
VectorT3
Vector< double, 3 > VectorT3
Definition: CellMark.cpp:6
FaceGroup
Definition: Mesh.h:28
StructureModel::getOptions
StructureModelOptions< T > & getOptions()
Definition: StructureModel_impl.h:1710
StructureModel::StructureVCMap
std::map< int, StructureVC< T > * > StructureVCMap
Definition: StructureModel.h:24
Mesh::getAllFaceGroups
const FaceGroupList & getAllFaceGroups() const
Definition: Mesh.h:193
StructureModel::Impl::VVDiagArray
VVMatrix::DiagArray VVDiagArray
Definition: StructureModel_impl.h:283
Field
Definition: Field.h:14
StructureModel::Impl::VectorT3T
VectorTranspose< T, 3 > VectorT3T
Definition: StructureModel_impl.h:275
StructureBCS::_faceArea
const VectorT3Array & _faceArea
Definition: StructureModel_impl.h:265
StructureModel::Impl::DiagTensorT3
SquareTensor< T, 3 > DiagTensorT3
Definition: StructureModel_impl.h:280
StructureBCS::_varField
const Field & _varField
Definition: StructureModel_impl.h:255
Mesh
Definition: Mesh.h:49
CRMatrix::PairWiseAssembler::getCoeff10
OffDiag & getCoeff10(const int np)
Definition: CRMatrix.h:131
StructureModel::StructureBCMap
std::map< int, StructureBC< T > * > StructureBCMap
Definition: StructureModel.h:23
StructureBCS::VectorT3Array
Array< VectorT3 > VectorT3Array
Definition: StructureModel_impl.h:59
MultiField::zero
virtual void zero()
Definition: MultiField.cpp:115
logCtor
#define logCtor()
Definition: RLogInterface.h:26
StructureModel::init
virtual void init()
Definition: StructureModel_impl.h:1695
StructureBCS::applyDirichletBC
X applyDirichletBC(const FloatValEvaluator< X > &bValue) const
Definition: StructureModel_impl.h:112
StructureModel::Impl::_deformationGradientModel
GradientModel< VectorT3 > _deformationGradientModel
Definition: StructureModel_impl.h:1659
StructureModel::getVCMap
StructureVCMap & getVCMap()
Definition: StructureModel_impl.h:1706
FaceGroup::groupType
string groupType
Definition: Mesh.h:42
LinearSystem::getDelta
MultiField & getDelta()
Definition: LinearSystem.h:34
StructureModel::StructureModel
StructureModel(const GeomFields &geomFields, StructureFields &structureFields, const MeshList &meshes)
Definition: StructureModel_impl.h:1677
FloatVarDict::getVal
FloatVal< T > getVal(const string varName) const
Definition: FloatVarDict.h:85
StructureBCS::CCAssembler
CCMatrix::PairWiseAssembler CCAssembler
Definition: StructureModel_impl.h:50
StructureBCS::applySymmetryBC
void applySymmetryBC() const
Definition: StructureModel_impl.h:210
StructureBCS::FMatrix
FluxJacobianMatrix< Diag, X > FMatrix
Definition: StructureModel_impl.h:52
StructureModel::Impl::initDeformationLinearization
void initDeformationLinearization(LinearSystem &ls)
Definition: StructureModel_impl.h:786
VectorT3Array
Array< Vector< double, 3 > > VectorT3Array
Definition: CellMark.cpp:7
sqrt
Tangent sqrt(const Tangent &a)
Definition: Tangent.h:317
StructureModel::advance
bool advance(const int niter)
Definition: StructureModel_impl.h:1725
StructureModel::getStrain
void getStrain(const Mesh &mesh)
Definition: StructureModel_impl.h:1778
StructureModel::Impl::_structureFields
StructureFields & _structureFields
Definition: StructureModel_impl.h:1653
StructureBCS::BBMatrix
DiagonalMatrix< Diag, X > BBMatrix
Definition: StructureModel_impl.h:53
FaceGroup::id
const int id
Definition: Mesh.h:41
LinearSystem::updateSolution
void updateSolution()
Definition: LinearSystem.cpp:265
StructureBCS::applyDirichletBC
X applyDirichletBC(const X &bValue) const
Definition: StructureModel_impl.h:104
StructureBCS::applyInterfaceBC
void applyInterfaceBC(const int f) const
Definition: StructureModel_impl.h:183
StructureModel::Impl::getPlasticDiagStrain
void getPlasticDiagStrain(const Mesh &mesh)
Definition: StructureModel_impl.h:1500
StructureModel::Impl::getBCMap
StructureBCMap & getBCMap()
Definition: StructureModel_impl.h:463
MultiField::ArrayIndex
pair< const Field *, const StorageSite * > ArrayIndex
Definition: MultiField.h:21
StructureModel::updateForceOnBoundary
void updateForceOnBoundary(const StorageSite &faceSite, const ArrayBase &bforceA, const map< int, int > &commonFacesMap, ArrayBase &fxA, ArrayBase &fyA, ArrayBase &fzA)
Definition: StructureModel_impl.h:1747
StructureBCS::_cells
const StorageSite & _cells
Definition: StructureModel_impl.h:253
LinearSystem::initSolve
void initSolve()
Definition: LinearSystem.cpp:40
StructureBCS::VectorT3
Vector< T_Scalar, 3 > VectorT3
Definition: StructureModel_impl.h:47
StructureBCS::_areaField
const Field & _areaField
Definition: StructureModel_impl.h:264
StructureModel::getBCMap
StructureBCMap & getBCMap()
Definition: StructureModel_impl.h:1702
StructureModel::Impl::getStrain
void getStrain(const Mesh &mesh)
Definition: StructureModel_impl.h:1449
DiscrList
vector< shared_ptr< Discretization > > DiscrList
Definition: Discretization.h:31
StructureModel::getTraction
void getTraction(const Mesh &mesh)
Definition: StructureModel_impl.h:1825
StructureBC::bcType
string bcType
Definition: StructureBC.h:34
Model::_meshes
const MeshList _meshes
Definition: Model.h:29
Model
Definition: Model.h:13
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
StructureModel::Impl::applyBC
void applyBC(LinearSystem &ls, bool explicitMode)
Definition: StructureModel_impl.h:808
StructureBCS::applyInterfaceBC
void applyInterfaceBC() const
Definition: StructureModel_impl.h:204
logDtor
#define logDtor()
Definition: RLogInterface.h:33
StorageSite::getCountLevel1
int getCountLevel1() const
Definition: StorageSite.h:72
CRMatrix::PairWiseAssembler::getCoeff01
OffDiag & getCoeff01(const int np)
Definition: CRMatrix.h:126
StructureModel::Impl::_geomFields
const GeomFields & _geomFields
Definition: StructureModel_impl.h:1652
StructureBCS::_dRdXDiag
DiagArray & _dRdXDiag
Definition: StructureModel_impl.h:259
StructureBCS::_faces
const StorageSite & _faces
Definition: StructureModel_impl.h:252
StructureModel::Impl::_lsK
shared_ptr< LinearSystem > _lsK
Definition: StructureModel_impl.h:1666
Array
Definition: Array.h:14
StructureBCS::applyZeroDerivativeBC
X applyZeroDerivativeBC() const
Definition: StructureModel_impl.h:151
MultiFieldMatrix::addMatrix
void addMatrix(const Index &rowI, const Index &colI, shared_ptr< Matrix > m)
Definition: MultiFieldMatrix.h:137
StorageSite::getOffset
int getOffset() const
Definition: StorageSite.h:87
CRMatrix::setBoundary
void setBoundary(const int nr)
Definition: CRMatrix.h:1056
StructureModel::Impl::VGradArray
Array< Gradient< VectorT3 > > VGradArray
Definition: StructureModel_impl.h:287
StructureModel::Impl::getOptions
StructureModelOptions< T > & getOptions()
Definition: StructureModel_impl.h:465
StructureModel::Impl::VGradType
Gradient< VectorT3 > VGradType
Definition: StructureModel_impl.h:286
StructureModel::Impl::VectorT3Array
Array< VectorT3 > VectorT3Array
Definition: StructureModel_impl.h:277
StructureBCS::_xIndex
const MultiField::ArrayIndex _xIndex
Definition: StructureModel_impl.h:256
Mesh::getCellCells2
const CRConnectivity & getCellCells2() const
Definition: Mesh.cpp:495
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
StructureBCS::applyNeumannBC
X applyNeumannBC(const X &bFlux) const
Definition: StructureModel_impl.h:142
StructureVC::vcType
string vcType
Definition: StructureBC.h:47
LinearSystem::getB
MultiField & getB()
Definition: LinearSystem.h:33
StructureModel::Impl::_options
StructureModelOptions< T > _options
Definition: StructureModel_impl.h:1658
StructureBCS::applyDirichletBC
X applyDirichletBC(int f, const X &bValue) const
Definition: StructureModel_impl.h:88
StructureModel::Impl::updateForceOnBoundary
void updateForceOnBoundary(const StorageSite &faceSite, const ArrayBase &bforceA, const map< int, int > &commonFacesMap, ArrayBase &fxA, ArrayBase &fyA, ArrayBase &fzA)
Definition: StructureModel_impl.h:1522
MFRPtr
shared_ptr< MultiFieldReduction > MFRPtr
Definition: MultiFieldReduction.h:56
MultiField::addArray
void addArray(const ArrayIndex &aIndex, shared_ptr< ArrayBase > a)
Definition: MultiField.cpp:270
LinearSystem::getX
MultiField & getX()
Definition: LinearSystem.h:32
CRMatrix::setDirichlet
void setDirichlet(const int nr)
Definition: CRMatrix.h:887
StructureBCS::T_Scalar
NumTypeTraits< X >::T_Scalar T_Scalar
Definition: StructureModel_impl.h:44
StructureBCS::CCMatrix
CRMatrix< Diag, OffDiag, X > CCMatrix
Definition: StructureModel_impl.h:49
Mesh::getDimension
int getDimension() const
Definition: Mesh.h:105
MultiFieldMatrix::multiply
virtual void multiply(IContainer &yB, const IContainer &xB) const
Definition: MultiFieldMatrix.cpp:65
StructureModel::Impl::Impl
Impl(const GeomFields &geomFields, StructureFields &structureFields, const MeshList &meshes)
Definition: StructureModel_impl.h:289
StructureModel::Impl::getVCMap
StructureVCMap & getVCMap()
Definition: StructureModel_impl.h:464
dot
T dot(const Vector< T, 3 > &a, const Vector< T, 3 > &b)
Definition: Vector.h:253
StructureBCS::applyNeumannBC
X applyNeumannBC(const FloatValEvaluator< X > &bFlux) const
Definition: StructureModel_impl.h:175
StructureModel::Impl::advanceExplicit
void advanceExplicit(const int nSteps, const double deltaT)
Definition: StructureModel_impl.h:1110
Vector< double, 3 >::getZero
static Vector getZero()
Definition: Vector.h:182
Mesh::getID
int getID() const
Definition: Mesh.h:106
StructureBCS::OffDiagArray
Array< OffDiag > OffDiagArray
Definition: StructureModel_impl.h:56
CRMatrix::eliminateDirichlet
void eliminateDirichlet(const int j, Array< X > &b, const X &delta_j, const bool explicitMode=false)
Definition: CRMatrix.h:1041
StructureBCS::_areaMagField
const Field & _areaMagField
Definition: StructureModel_impl.h:262
StructureModel::getPlasticDiagStrain
void getPlasticDiagStrain(const Mesh &mesh)
Definition: StructureModel_impl.h:1785
Linearizer::linearize
virtual void linearize(DiscrList &discretizations, const MeshList &meshes, MultiFieldMatrix &matrix, MultiField &x, MultiField &b)
Definition: Linearizer.cpp:17
StructureBCS::_faceCells
const CRConnectivity & _faceCells
Definition: StructureModel_impl.h:254
LinearSystem::getMatrix
MultiFieldMatrix & getMatrix()
Definition: LinearSystem.h:37
MeshList
vector< Mesh * > MeshList
Definition: Mesh.h:439
FaceGroup::site
StorageSite site
Definition: Mesh.h:40
Array::getLength
int getLength() const
Definition: Array.h:87
StructureModel::Impl::advance
bool advance(const int niter)
Definition: StructureModel_impl.h:1185