Ice Sheet System Model  4.18
Code documentation
DamageEvolutionAnalysis.cpp
Go to the documentation of this file.
1 #include "./DamageEvolutionAnalysis.h"
2 #include "../toolkits/toolkits.h"
3 #include "../classes/classes.h"
4 #include "../shared/shared.h"
5 #include "../modules/modules.h"
6 
7 /*Model processing*/
8 void DamageEvolutionAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
9 
10  int finiteelement;
11  iomodel->FindConstant(&finiteelement,"md.damage.elementinterp");
12 
13  /*Fetch parameters: */
14  int stabilization;
15  iomodel->FindConstant(&stabilization,"md.damage.stabilization");
16 
17  /*Do not add constraints in DG, they are weakly imposed*/
18  if(stabilization!=3){
19  IoModelToConstraintsx(constraints,iomodel,"md.damage.spcdamage",DamageEvolutionAnalysisEnum,finiteelement);
20  }
21 
22  /*FCT, constraints are imposed using penalties*/
23  if(stabilization==4){
24  constraints->ActivatePenaltyMethod(DamageEvolutionAnalysisEnum);
25  }
26 }/*}}}*/
27 void DamageEvolutionAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
28 
29  /*Nothing for now*/
30 
31 }/*}}}*/
32 void DamageEvolutionAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
33 
34  int finiteelement;
35 
36  iomodel->FindConstant(&finiteelement,"md.damage.elementinterp");
37  ::CreateNodes(nodes,iomodel,DamageEvolutionAnalysisEnum,finiteelement);
38 }/*}}}*/
39 int DamageEvolutionAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
40  return 1;
41 }/*}}}*/
42 void DamageEvolutionAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
43 
44  int finiteelement;
45  bool ismovingfront;
46 
47  iomodel->FindConstant(&finiteelement,"md.damage.elementinterp");
48  iomodel->FindConstant(&ismovingfront,"md.transient.ismovingfront");
49 
50  /*Update elements: */
51  iomodel->FetchData(1,"md.flowequation.element_equation");
52  int counter=0;
53  for(int i=0;i<iomodel->numberofelements;i++){
54  if(iomodel->my_elements[i]){
55  Element* element=(Element*)elements->GetObjectByOffset(counter);
56  element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,finiteelement);
57 
58  /*Need to know the type of approximation for this element*/
59  if(iomodel->Data("md.flowequation.element_equation")){
60  inputs2->SetInput(ApproximationEnum,counter,IoCodeToEnumElementEquation(reCast<int>(iomodel->Data("md.flowequation.element_equation")[i])));
61  }
62  counter++;
63  }
64  }
65  iomodel->DeleteData(1,"md.flowequation.element_equation");
66 
67  /*First, reset all F to 0 */
68  for(int i=0;i<elements->Size();i++){
69  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
70  element->SetElementInput(inputs2,DamageFEnum,0.);
71  }
72 
73 
74  /*What input do I need to run my damage evolution model?*/
75  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vx",VxEnum);
76  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vy",VyEnum);
77  if(iomodel->domaintype==Domain3DEnum) iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vz",VzEnum);
78  iomodel->FetchDataToInput(inputs2,elements,"md.damage.D",DamageDEnum);
79  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
80  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.pressure",PressureEnum);
81 
82 }/*}}}*/
83 void DamageEvolutionAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
84 
85  /*Intermediaries*/
86  int numoutputs;
87  char** requestedoutputs = NULL;
88 
89  /*retrieve some parameters: */
90  parameters->AddObject(iomodel->CopyConstantObject("md.damage.law",DamageLawEnum));
91  parameters->AddObject(iomodel->CopyConstantObject("md.damage.stabilization",DamageStabilizationEnum));
92  parameters->AddObject(iomodel->CopyConstantObject("md.damage.max_damage",DamageMaxDamageEnum));
93 
94  /*Requested outputs*/
95  iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.damage.requested_outputs");
96  parameters->AddObject(new IntParam(DamageEvolutionNumRequestedOutputsEnum,numoutputs));
97  if(numoutputs)parameters->AddObject(new StringArrayParam(DamageEvolutionRequestedOutputsEnum,requestedoutputs,numoutputs));
98  iomodel->DeleteData(&requestedoutputs,numoutputs,"md.damage.requested_outputs");
99 
100  /*Retrieve law dependent parameters: */
101  int law;
102  iomodel->FindConstant(&law,"md.damage.law");
103  if (law==0){
104  parameters->AddObject(iomodel->CopyConstantObject("md.damage.stress_threshold",DamageStressThresholdEnum));
105  parameters->AddObject(iomodel->CopyConstantObject("md.damage.kappa",DamageKappaEnum));
106  }
107  else if (law>0){
108  parameters->AddObject(iomodel->CopyConstantObject("md.damage.c1",DamageC1Enum));
109  parameters->AddObject(iomodel->CopyConstantObject("md.damage.c2",DamageC2Enum));
110  parameters->AddObject(iomodel->CopyConstantObject("md.damage.c3",DamageC3Enum));
111  parameters->AddObject(iomodel->CopyConstantObject("md.damage.c4",DamageC4Enum));
112  parameters->AddObject(iomodel->CopyConstantObject("md.damage.stress_threshold",DamageStressThresholdEnum));
113  parameters->AddObject(iomodel->CopyConstantObject("md.damage.stress_ubound",DamageStressUBoundEnum));
114  parameters->AddObject(iomodel->CopyConstantObject("md.damage.kappa",DamageKappaEnum));
115  parameters->AddObject(iomodel->CopyConstantObject("md.damage.healing",DamageHealingEnum));
116  parameters->AddObject(iomodel->CopyConstantObject("md.damage.equiv_stress",DamageEquivStressEnum));
117  }
118 
119 }/*}}}*/
120 
121 /*Finite Element Analysis*/
122 void DamageEvolutionAnalysis::Core(FemModel* femmodel){/*{{{*/
123  _error_("not implemented");
124 }/*}}}*/
125 void DamageEvolutionAnalysis::CreateDamageFInput(Element* element){/*{{{*/
126 
127  /*Fetch number of vertices and allocate output*/
128  int numnodes = element->GetNumberOfNodes();
129  IssmDouble* f = xNew<IssmDouble>(numnodes);
130 
131  /*Calculate damage evolution source term: */
132  for (int i=0;i<numnodes;i++){
133 
134  /* healing could be handled here */
135 
136  /* no source term; damage handled in stress balance */
137  f[i]=0.;
138  }
139 
140  /*Add input*/
141  element->AddInput2(DamageFEnum,f,element->GetElementType());
142 
143  /*Clean up and return*/
144  xDelete<IssmDouble>(f);
145 }/*}}}*/
146 void DamageEvolutionAnalysis::CreateDamageFInputArctan(Element* element){/*{{{*/
147  IssmDouble c1, c2, stress_threshold, stress_ubound;
148  IssmDouble damage;
149  IssmDouble yts;
150  IssmDouble principalDevStress1, principalDevStress2;
151  IssmDouble tensileStress, compressiveStress;
152 
153  int equivstress, domaintype, dim;
154 
155  /*Fetch number of vertices and allocate output*/
156  int numnodes = element->GetNumberOfNodes();
157  IssmDouble* f = xNew<IssmDouble>(numnodes);
158 
159  /*retrieve parameters:*/
160  element->FindParam(&c1,DamageC1Enum);
161  element->FindParam(&c2,DamageC2Enum);
162  element->FindParam(&yts,ConstantsYtsEnum);
163  element->FindParam(&stress_threshold,DamageStressThresholdEnum);
164  element->FindParam(&stress_ubound,DamageStressUBoundEnum);
165  element->FindParam(&domaintype,DomainTypeEnum);
166 
167  /*Get problem dimension*/
168  switch(domaintype){
169  case Domain2DhorizontalEnum: dim = 2; break;
170  case Domain3DEnum: dim = 3; break;
171  default: _error_("not implemented");
172  }
173  /*Compute stress tensor and Stress Max Principal: */
174  element->ComputeDeviatoricStressTensor();
175 
176  Input2* principalDevStress1_input = element->GetInput2(DeviatoricStress1Enum); _assert_(principalDevStress1_input);
177  Input2* principalDevStress2_input = element->GetInput2(DeviatoricStress2Enum); _assert_(principalDevStress2_input);
178 
179  Input2* damage_input = NULL;
180  if(domaintype==Domain2DhorizontalEnum){
181  damage_input = element->GetInput2(DamageDbarEnum); _assert_(damage_input);
182  }
183  else{
184  damage_input = element->GetInput2(DamageDEnum); _assert_(damage_input);
185  }
186 
187  /*Calculate damage evolution source term */
188  Gauss* gauss=element->NewGauss();
189 
190  /* To keep arctan output (bounded by -pi/2 and pi/2) within the specified boundaries */
191  c1 /= (PI/2);
192  c2 /= (PI/2);
193  /* To have per second output with per annum parameters */
194  c1 /= yts;
195  c2 /= yts;
196 
197  for (int i=0;i<numnodes;i++){
198  f[i] = 0;
199 
200  gauss->GaussNode(element->GetElementType(),i);
201 
202  damage_input->GetInputValue(&damage,gauss);
203  principalDevStress1_input->GetInputValue(&principalDevStress1,gauss);
204  principalDevStress2_input->GetInputValue(&principalDevStress2,gauss);
205 
206  tensileStress = sqrt(1.5*(pow(max(principalDevStress1, 0.), 2) + pow(max(principalDevStress2, 0.), 2)));
207  compressiveStress = sqrt(1.5*(pow(min(principalDevStress1, 0.), 2) + pow(min(principalDevStress2, 0.), 2)));
208 
209  /* Calculate principal effective stresses */
210  if(dim==2){
211  f[i] = 0;
212  if(tensileStress > stress_threshold)
213  f[i] += c1*atan((tensileStress/stress_threshold - 1)/(1-damage));
214 
215  if(compressiveStress < stress_ubound)
216  f[i] += c2*atan((compressiveStress/stress_ubound - 1)/(1-damage));
217  }
218  else{
219  _error_("Only 2D is implemented.");
220  }
221  }
222 
223  /*Add input*/
224  element->AddInput2(DamageFEnum,f,element->GetElementType());
225 
226  /*Clean up and return*/
227  xDelete<IssmDouble>(f);
228  delete gauss;
229 }/*}}}*/
230 
231 void DamageEvolutionAnalysis::CreateDamageFInputExp(Element* element){/*{{{*/
232 
233  /*Intermediaries */
234  IssmDouble epsf,stress_threshold,eps0;
235  IssmDouble damage,B,n,epseff;
236  IssmDouble eps_xx,eps_yy,eps_xy,eps1,eps2,epstmp;
237  int domaintype;
238 
239  /*Fetch number of vertices and allocate output*/
240  int numnodes = element->GetNumberOfNodes();
241  IssmDouble* f = xNew<IssmDouble>(numnodes);
242 
243  /*retrieve parameters:*/
244  element->FindParam(&epsf,DamageC1Enum);
245  element->FindParam(&stress_threshold,DamageStressThresholdEnum);
246  element->FindParam(&domaintype,DomainTypeEnum);
247 
248  /*Compute stress tensor: */
249  element->ComputeStrainRate();
250 
251  /*retrieve what we need: */
252  Input2* eps_xx_input = element->GetInput2(StrainRatexxEnum); _assert_(eps_xx_input);
253  Input2* eps_xy_input = element->GetInput2(StrainRatexyEnum); _assert_(eps_xy_input);
254  Input2* eps_yy_input = element->GetInput2(StrainRateyyEnum); _assert_(eps_yy_input);
255  Input2* n_input=element->GetInput2(MaterialsRheologyNEnum); _assert_(n_input);
256  Input2* damage_input = NULL;
257  Input2* B_input = NULL;
258  if(domaintype==Domain2DhorizontalEnum){
259  damage_input = element->GetInput2(DamageDbarEnum); _assert_(damage_input);
260  B_input=element->GetInput2(MaterialsRheologyBbarEnum); _assert_(B_input);
261  }
262  else{
263  damage_input = element->GetInput2(DamageDEnum); _assert_(damage_input);
264  B_input=element->GetInput2(MaterialsRheologyBEnum); _assert_(B_input);
265  }
266 
267  /*Calculate damage evolution source term: */
268  Gauss* gauss=element->NewGauss();
269  for (int i=0;i<numnodes;i++){
270  gauss->GaussNode(element->GetElementType(),i);
271 
272  eps_xx_input->GetInputValue(&eps_xx,gauss);
273  eps_xy_input->GetInputValue(&eps_xy,gauss);
274  eps_yy_input->GetInputValue(&eps_yy,gauss);
275  B_input->GetInputValue(&B,gauss);
276  n_input->GetInputValue(&n,gauss);
277  damage_input->GetInputValue(&damage,gauss);
278 
279  /*Calculate principal effective strain rates*/
280  eps1=(eps_xx+eps_yy)/2.+sqrt(pow((eps_xx-eps_yy)/2.,2)+pow(eps_xy,2));
281  eps2=(eps_xx+eps_yy)/2.-sqrt(pow((eps_xx-eps_yy)/2.,2)+pow(eps_xy,2));
282  if(fabs(eps2)>fabs(eps1)){epstmp=eps2; eps2=eps1; eps1=epstmp;}
283 
284  /*Calculate effective strain rate and threshold strain rate*/
285  epseff=1./sqrt(2.)*sqrt(eps1*eps1-eps1*eps2+eps2*eps2);
286  eps0=pow(stress_threshold/B,n);
287 
288  if(epseff>eps0){
289  f[i]=1.-pow(eps0/epseff,1./n)*exp(-(epseff-eps0)/(epsf-eps0))-damage;
290  }
291  else f[i]=0;
292 
293  /*Edits from MM*/
294  if(f[i]>10.) f[i]=10.;
295  if(f[i]<-10.) f[i]=-10.;
296  }
297 
298  /*Add input*/
299  element->AddInput2(DamageFEnum,f,P1DGEnum);
300 
301  /*Clean up and return*/
302  xDelete<IssmDouble>(f);
303  delete gauss;
304 }/*}}}*/
305 void DamageEvolutionAnalysis::CreateDamageFInputPralong(Element* element){/*{{{*/
306 
307  /*Intermediaries */
308  IssmDouble c1,c2,c3,healing,stress_threshold;
309  IssmDouble s_xx,s_xy,s_xz,s_yy,s_yz,s_zz,s1,s2,s3,stmp;
310  IssmDouble J2s,Chi,Psi,PosPsi,NegPsi;
311  IssmDouble damage,tau_xx,tau_xy,tau_xz,tau_yy,tau_yz,tau_zz,stressMaxPrincipal;
312  int equivstress,domaintype,dim;
313 
314  /*Fetch number of vertices and allocate output*/
315  int numnodes = element->GetNumberOfNodes();
316  IssmDouble* f = xNew<IssmDouble>(numnodes);
317 
318  /*retrieve parameters:*/
319  element->FindParam(&c1,DamageC1Enum);
320  element->FindParam(&c2,DamageC2Enum);
321  element->FindParam(&c3,DamageC3Enum);
322  element->FindParam(&healing,DamageHealingEnum);
323  element->FindParam(&stress_threshold,DamageStressThresholdEnum);
324  element->FindParam(&domaintype,DomainTypeEnum);
325 
326  /*Get problem dimension*/
327  switch(domaintype){
328  case Domain2DhorizontalEnum: dim = 2; break;
329  case Domain3DEnum: dim = 3; break;
330  default: _error_("not implemented");
331  }
332  /*Compute stress tensor and Stress Max Principal: */
333  element->ComputeDeviatoricStressTensor();
334  if(dim==3){
335  /*Only works in 3d because the pressure is defined*/
336  element->StressMaxPrincipalCreateInput();
337  }
338  /*retrieve what we need: */
339  Input2* tau_xx_input = element->GetInput2(DeviatoricStressxxEnum); _assert_(tau_xx_input);
340  Input2* tau_xy_input = element->GetInput2(DeviatoricStressxyEnum); _assert_(tau_xy_input);
341  Input2* tau_yy_input = element->GetInput2(DeviatoricStressyyEnum); _assert_(tau_yy_input);
342  Input2* tau_xz_input = NULL;
343  Input2* tau_yz_input = NULL;
344  Input2* tau_zz_input = NULL;
345  Input2* stressMaxPrincipal_input = NULL;
346  if(dim==3){
347  tau_xz_input = element->GetInput2(DeviatoricStressxzEnum); _assert_(tau_xz_input);
348  tau_yz_input = element->GetInput2(DeviatoricStressyzEnum); _assert_(tau_yz_input);
349  tau_zz_input = element->GetInput2(DeviatoricStresszzEnum); _assert_(tau_zz_input);
350  stressMaxPrincipal_input = element->GetInput2(StressMaxPrincipalEnum); _assert_(stressMaxPrincipal_input);
351  }
352  Input2* damage_input = NULL;
353  if(domaintype==Domain2DhorizontalEnum){
354  damage_input = element->GetInput2(DamageDbarEnum); _assert_(damage_input);
355  }
356  else{
357  damage_input = element->GetInput2(DamageDEnum); _assert_(damage_input);
358  }
359 
360  /*retrieve the desired type of equivalent stress*/
361  element->FindParam(&equivstress,DamageEquivStressEnum);
362 
363  /*Calculate damage evolution source term: */
364  Gauss* gauss=element->NewGauss();
365  for (int i=0;i<numnodes;i++){
366  gauss->GaussNode(element->GetElementType(),i);
367 
368  damage_input->GetInputValue(&damage,gauss);
369  tau_xx_input->GetInputValue(&tau_xx,gauss);
370  tau_xy_input->GetInputValue(&tau_xy,gauss);
371  tau_yy_input->GetInputValue(&tau_yy,gauss);
372  if(dim==3){
373  tau_xz_input->GetInputValue(&tau_xz,gauss);
374  tau_yz_input->GetInputValue(&tau_yz,gauss);
375  tau_zz_input->GetInputValue(&tau_zz,gauss);
376  }
377  /*Calculate effective stress components*/
378  s_xx=tau_xx/(1.-damage);
379  s_xy=tau_xy/(1.-damage);
380  s_yy=tau_yy/(1.-damage);
381  if(dim==3){
382  s_xz=tau_xz/(1.-damage);
383  s_yz=tau_yz/(1.-damage);
384  s_zz=tau_zz/(1.-damage);
385  }
386  /*Calculate principal effective stresses*/
387  if(dim==2){
388  s1=(s_xx+s_yy)/2.+sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
389  s2=(s_xx+s_yy)/2.-sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
390  if(fabs(s2)>fabs(s1)){stmp=s2; s2=s1; s1=stmp;}
391 
392  if(equivstress==0){ /* von Mises */
393  Chi=sqrt(s1*s1-s1*s2+s2*s2);
394  }
395  else if(equivstress==1){ /* max principal stress */
396  Chi=s1;
397  }
398  Psi=Chi-stress_threshold;
399  NegPsi=max(-Chi,0.); /* healing only for compressive stresses */
400  PosPsi=max(Psi,0.);
401  f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
402  }
403  else{
404  if(equivstress==1){/* max principal stress */
405  stressMaxPrincipal_input->GetInputValue(&stressMaxPrincipal,gauss);
406  Chi=stressMaxPrincipal/(1.-damage);
407  }
408  else if(equivstress==0){/* von Mises */
409  Chi=sqrt(((s_xx-s_yy)*(s_xx-s_yy)+(s_yy-s_zz)*(s_yy-s_zz)+(s_zz-s_xx)*(s_zz-s_xx)+6.*(s_xy*s_xy+s_yz*s_yz+s_xz*s_xz))/2.);
410  }
411  Psi=Chi-stress_threshold;
412  NegPsi=max(-Chi,0.); /* healing only for compressive stresses */
413  PosPsi=max(Psi,0.);
414  f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
415  }
416  }
417  /*Add input*/
418  element->AddInput2(DamageFEnum,f,P1DGEnum);
419 
420  /*Clean up and return*/
421  xDelete<IssmDouble>(f);
422  delete gauss;
423 }/*}}}*/
424 ElementVector* DamageEvolutionAnalysis::CreateDVector(Element* element){/*{{{*/
425  /*Default, return NULL*/
426  return NULL;
427 }/*}}}*/
428 ElementMatrix* DamageEvolutionAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
429 _error_("Not implemented");
430 }/*}}}*/
431 ElementMatrix* DamageEvolutionAnalysis::CreateKMatrix(Element* element){/*{{{*/
432 
433  /* Check if ice in element */
434  if(!element->IsIceInElement()) return NULL;
435  /*Intermediaries*/
436  int domaintype,dim;
437  int stabilization;
438  IssmDouble Jdet,dt,D_scalar,h,hx,hy,hz;
439  IssmDouble vel,vx,vy,vz,dvxdx,dvydy,dvzdz,dvx[3],dvy[3],dvz[3];
440  IssmDouble *xyz_list = NULL;
441 
442  /*Get problem dimension*/
443  element->FindParam(&domaintype,DomainTypeEnum);
444  switch(domaintype){
445  case Domain2DhorizontalEnum: dim = 2; break;
446  case Domain3DEnum: dim = 3; break;
447  default: _error_("Not implemented yet");
448  }
449 
450  /*Fetch number of nodes and dof for this finite element*/
451  int numnodes = element->GetNumberOfNodes();
452 
453  /*Initialize Element vector*/
454  ElementMatrix* Ke = element->NewElementMatrix();
455  IssmDouble* basis = xNew<IssmDouble>(numnodes);
456  IssmDouble* B = xNew<IssmDouble>(dim*numnodes);
457  IssmDouble* Bprime = xNew<IssmDouble>(dim*numnodes);
458  IssmDouble* D = xNewZeroInit<IssmDouble>(dim*dim);
459 
460  /*Retrieve all inputs and parameters*/
461  element->GetVerticesCoordinates(&xyz_list);
462  element->FindParam(&dt,TimesteppingTimeStepEnum);
463  //printf("dt %f\n", dt);
464  element->FindParam(&stabilization,DamageStabilizationEnum);
465  Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
466  Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
467  Input2* vz_input = NULL;
468  if(dim==3){
469  vz_input=element->GetInput2(VzEnum); _assert_(vz_input);
470  }
471 
472  if(dim==2) h=element->CharacteristicLength();
473 
474  /* Start looping on the number of gaussian points: */
475  Gauss* gauss=element->NewGauss(2);
476  for(int ig=gauss->begin();ig<gauss->end();ig++){
477  gauss->GaussPoint(ig);
478 
479  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
480  element->NodalFunctions(basis,gauss);
481 
482  vx_input->GetInputValue(&vx,gauss);
483  vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
484  vy_input->GetInputValue(&vy,gauss);
485  vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
486 
487  if(dim==3){
488  vz_input->GetInputValue(&vz,gauss);
489  vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
490  }
491 
492  D_scalar=gauss->weight*Jdet;
493  TripleMultiply(basis,1,numnodes,1,
494  &D_scalar,1,1,0,
495  basis,1,numnodes,0,
496  &Ke->values[0],1);
497 
498  GetB(B,element,dim,xyz_list,gauss);
499  GetBprime(Bprime,element,dim,xyz_list,gauss);
500 
501  dvxdx=dvx[0];
502  dvydy=dvy[1];
503  if(dim==3) dvzdz=dvz[2];
504  D_scalar=dt*gauss->weight*Jdet;
505 
506  D[0*dim+0]=D_scalar*dvxdx;
507  D[1*dim+1]=D_scalar*dvydy;
508  if(dim==3) D[2*dim+2]=D_scalar*dvzdz;
509 
510  TripleMultiply(B,dim,numnodes,1,
511  D,dim,dim,0,
512  B,dim,numnodes,0,
513  &Ke->values[0],1);
514 
515  D[0*dim+0]=D_scalar*vx;
516  D[1*dim+1]=D_scalar*vy;
517  if(dim==3) D[2*dim+2]=D_scalar*vz;
518 
519  TripleMultiply(B,dim,numnodes,1,
520  D,dim,dim,0,
521  Bprime,dim,numnodes,0,
522  &Ke->values[0],1);
523 
524  if(stabilization==2){
525  if(dim==3){
526  vel=sqrt(vx*vx+vy*vy+vz*vz)+1.e-8;
527  D[0*dim+0]=h/(2.0*vel)*vx*vx;
528  D[1*dim+0]=h/(2.0*vel)*vy*vx;
529  D[2*dim+0]=h/(2.0*vel)*vz*vx;
530  D[0*dim+1]=h/(2.0*vel)*vx*vy;
531  D[1*dim+1]=h/(2.0*vel)*vy*vy;
532  D[2*dim+1]=h/(2.0*vel)*vy*vz;
533  D[0*dim+2]=h/(2.0*vel)*vx*vz;
534  D[1*dim+2]=h/(2.0*vel)*vy*vz;
535  D[2*dim+2]=h/(2.0*vel)*vz*vz;
536  }
537  else{
538  /*Streamline upwinding*/
539  vel=sqrt(vx*vx+vy*vy)+1.e-8;
540  D[0*dim+0]=h/(2.0*vel)*vx*vx;
541  D[1*dim+0]=h/(2.0*vel)*vy*vx;
542  D[0*dim+1]=h/(2.0*vel)*vx*vy;
543  D[1*dim+1]=h/(2.0*vel)*vy*vy;
544  }
545  }
546  else if(stabilization==1){
547  if(dim==2){
548  vx_input->GetInputAverage(&vx);
549  vy_input->GetInputAverage(&vy);
550  D[0*dim+0]=h/2.0*fabs(vx);
551  D[1*dim+1]=h/2.0*fabs(vy);
552  }
553  else if(dim==3){
554  element->ElementSizes(&hx,&hy,&hz);
555  vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
556  h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
557  D[0*dim+0]=h/(2.*vel)*fabs(vx*vx); D[0*dim+1]=h/(2.*vel)*fabs(vx*vy); D[0*dim+2]=h/(2.*vel)*fabs(vx*vz);
558  D[1*dim+0]=h/(2.*vel)*fabs(vy*vx); D[1*dim+1]=h/(2.*vel)*fabs(vy*vy); D[1*dim+2]=h/(2.*vel)*fabs(vy*vz);
559  D[2*dim+0]=h/(2.*vel)*fabs(vz*vx); D[2*dim+1]=h/(2.*vel)*fabs(vz*vy); D[2*dim+2]=h/(2.*vel)*fabs(vz*vz);
560  }
561  }
562  if(stabilization==1 || stabilization==2){
563  if(dim==2){
564  D[0*dim+0]=D_scalar*D[0*dim+0];
565  D[1*dim+0]=D_scalar*D[1*dim+0];
566  D[0*dim+1]=D_scalar*D[0*dim+1];
567  D[1*dim+1]=D_scalar*D[1*dim+1];
568  }
569  else if(dim==3){
570  D[0*dim+0]=D_scalar*D[0*dim+0];
571  D[1*dim+0]=D_scalar*D[1*dim+0];
572  D[2*dim+0]=D_scalar*D[2*dim+0];
573  D[0*dim+1]=D_scalar*D[0*dim+1];
574  D[1*dim+1]=D_scalar*D[1*dim+1];
575  D[2*dim+1]=D_scalar*D[2*dim+1];
576  D[0*dim+2]=D_scalar*D[0*dim+2];
577  D[1*dim+2]=D_scalar*D[1*dim+2];
578  D[2*dim+2]=D_scalar*D[2*dim+2];
579  }
580  TripleMultiply(Bprime,dim,numnodes,1,
581  D,dim,dim,0,
582  Bprime,dim,numnodes,0,
583  &Ke->values[0],1);
584  }
585 
586  }
587 
588  /*Clean up and return*/
589  xDelete<IssmDouble>(xyz_list);
590  xDelete<IssmDouble>(basis);
591  xDelete<IssmDouble>(B);
592  xDelete<IssmDouble>(Bprime);
593  xDelete<IssmDouble>(D);
594  delete gauss;
595  return Ke;
596 }/*}}}*/
597 ElementVector* DamageEvolutionAnalysis::CreatePVector(Element* element){/*{{{*/
598 
599  /* Check if ice in element */
600  if(!element->IsIceInElement()) return NULL;
601 
602  /*Intermediaries*/
603  int domaintype,damagelaw;
604  IssmDouble Jdet,dt;
605  IssmDouble f,damage;
606  IssmDouble* xyz_list = NULL;
607  /*Get element*/
608  element->FindParam(&domaintype,DomainTypeEnum);
609 
610  /*Fetch number of nodes and dof for this finite element*/
611  int numnodes = element->GetNumberOfNodes();
612 
613  /*Initialize Element vector and other vectors*/
614  ElementVector* pe = element->NewElementVector();
615  IssmDouble* basis = xNew<IssmDouble>(numnodes);
616 
617  /*Retrieve all inputs and parameters*/
618  element->GetVerticesCoordinates(&xyz_list);
619  element->FindParam(&dt,TimesteppingTimeStepEnum);
620  element->FindParam(&damagelaw,DamageLawEnum);
621  switch(damagelaw){
622  case 0:
623  this->CreateDamageFInput(element);
624  break;
625  case 1:
626  this->CreateDamageFInputPralong(element);
627  break;
628  case 2:
629  this->CreateDamageFInputExp(element);
630  break;
631  case 3:
632  this->CreateDamageFInputArctan(element);
633  break;
634  default:
635  _error_("not implemented yet");
636  }
637 
638  Input2* damaged_input = NULL;
639  Input2* damagef_input = element->GetInput2(DamageFEnum); _assert_(damagef_input);
640  if(domaintype==Domain2DhorizontalEnum){
641  damaged_input = element->GetInput2(DamageDbarEnum); _assert_(damaged_input);
642  }
643  else{
644  damaged_input = element->GetInput2(DamageDEnum); _assert_(damaged_input);
645  }
646 
647  /* Start looping on the number of gaussian points: */
648  Gauss* gauss=element->NewGauss(2);
649  for(int ig=gauss->begin();ig<gauss->end();ig++){
650  gauss->GaussPoint(ig);
651 
652  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
653  element->NodalFunctions(basis,gauss);
654 
655  damaged_input->GetInputValue(&damage,gauss);
656  damagef_input->GetInputValue(&f,gauss);
657 
658  for(int i=0;i<numnodes;i++){
659  pe->values[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i];
660  }
661  }
662  /*Clean up and return*/
663  xDelete<IssmDouble>(xyz_list);
664  xDelete<IssmDouble>(basis);
665  delete gauss;
666  return pe;
667 }/*}}}*/
668 void DamageEvolutionAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
669  /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*2.
670  * For node i, Bi can be expressed in the actual coordinate system
671  * by:
672  * Bi=[ N ]
673  * [ N ]
674  * where N is the finiteelement function for node i.
675  *
676  * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
677  */
678 
679  /*Fetch number of nodes for this finite element*/
680  int numnodes = element->GetNumberOfNodes();
681 
682  /*Get nodal functions*/
683  IssmDouble* basis=xNew<IssmDouble>(numnodes);
684  element->NodalFunctions(basis,gauss);
685 
686  /*Build B: */
687  for(int i=0;i<numnodes;i++){
688  for(int j=0;j<dim;j++){
689  B[numnodes*j+i] = basis[i];
690  }
691  }
692 
693  /*Clean-up*/
694  xDelete<IssmDouble>(basis);
695 }/*}}}*/
696 void DamageEvolutionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
697  /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2.
698  * For node i, Bi' can be expressed in the actual coordinate system
699  * by:
700  * Bi_prime=[ dN/dx ]
701  * [ dN/dy ]
702  * where N is the finiteelement function for node i.
703  *
704  * We assume B' has been allocated already, of size: 3x(2*numnodes)
705  */
706 
707  /*Fetch number of nodes for this finite element*/
708  int numnodes = element->GetNumberOfNodes();
709 
710  /*Get nodal functions derivatives*/
711  IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
712  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
713 
714  /*Build B': */
715  for(int i=0;i<numnodes;i++){
716  for(int j=0;j<dim;j++){
717  Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
718  }
719  }
720 
721  /*Clean-up*/
722  xDelete<IssmDouble>(dbasis);
723 
724 }/*}}}*/
725 void DamageEvolutionAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
726  element->GetSolutionFromInputsOneDof(solution,DamageDbarEnum);
727 }/*}}}*/
728 void DamageEvolutionAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
729  _error_("Not implemented yet");
730 }/*}}}*/
731 void DamageEvolutionAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
732 
733  /*Fetch number of nodes and dof for this finite element*/
734  int numnodes = element->GetNumberOfNodes();
735 
736  /*Fetch dof list and allocate solution vector*/
737  int *doflist = NULL;
738  element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
739  IssmDouble* newdamage = xNew<IssmDouble>(numnodes);
740 
741  /*Get user-supplied max_damage: */
742  IssmDouble max_damage = element->FindParam(DamageMaxDamageEnum);
743 
744  /*Use the dof list to index into the solution vector: */
745  for(int i=0;i<numnodes;i++){
746  newdamage[i]=solution[doflist[i]];
747  /*Check solution*/
748  if(xIsNan<IssmDouble>(newdamage[i])) _error_("NaN found in solution vector");
749  if(xIsInf<IssmDouble>(newdamage[i])) _error_("Inf found in solution vector");
750  /*Enforce D < max_damage and D > 0 */
751  if(newdamage[i]>max_damage) newdamage[i]=max_damage;
752  else if(newdamage[i]<0.) newdamage[i]=0.;
753  }
754 
755  /*Get all inputs and parameters*/
756  int domaintype;
757  element->FindParam(&domaintype,DomainTypeEnum);
758  if(domaintype==Domain2DhorizontalEnum){
759  element->AddInput2(DamageDbarEnum,newdamage,element->GetElementType());
760  }
761  else{
762  element->AddInput2(DamageDEnum,newdamage,element->GetElementType());
763  }
764 
765  /*Free ressources:*/
766  xDelete<IssmDouble>(newdamage);
767  xDelete<int>(doflist);
768 }/*}}}*/
769 void DamageEvolutionAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
770  SetActiveNodesLSMx(femmodel);
771 }/*}}}*/
772 
773 /*Flux Correction Transport*/
774 ElementMatrix* DamageEvolutionAnalysis::CreateFctKMatrix(Element* element){/*{{{*/
775 
776  /* Check if ice in element */
777  if(!element->IsIceInElement()) return NULL;
778 
779  /*Intermediaries */
780  IssmDouble Jdet;
781  IssmDouble vx,vy;
782  IssmDouble* xyz_list = NULL;
783 
784  /*Fetch number of nodes and dof for this finite element*/
785  int numnodes = element->GetNumberOfNodes();
786  int dim = 2;
787 
788  /*Initialize Element vector and other vectors*/
789  ElementMatrix* Ke = element->NewElementMatrix();
790  IssmDouble* B = xNew<IssmDouble>(dim*numnodes);
791  IssmDouble* Bprime = xNew<IssmDouble>(dim*numnodes);
792  IssmDouble* D = xNewZeroInit<IssmDouble>(dim*dim);
793 
794  /*Retrieve all inputs and parameters*/
795  element->GetVerticesCoordinates(&xyz_list);
796  Input2* vxaverage_input=element->GetInput2(VxEnum); _assert_(vxaverage_input);
797  Input2* vyaverage_input=element->GetInput2(VyEnum); _assert_(vyaverage_input);
798 
799  /* Start looping on the number of gaussian points: */
800  Gauss* gauss=element->NewGauss(2);
801  for(int ig=gauss->begin();ig<gauss->end();ig++){
802  gauss->GaussPoint(ig);
803 
804  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
805  GetB(B,element,dim,xyz_list,gauss);
806  GetBprime(Bprime,element,dim,xyz_list,gauss);
807  vxaverage_input->GetInputValue(&vx,gauss);
808  vyaverage_input->GetInputValue(&vy,gauss);
809 
810  D[0*dim+0] = -gauss->weight*vx*Jdet;
811  D[1*dim+1] = -gauss->weight*vy*Jdet;
812 
813  TripleMultiply(B,dim,numnodes,1,
814  D,dim,dim,0,
815  Bprime,dim,numnodes,0,
816  &Ke->values[0],1);
817 
818  }
819 
820  /*Clean up and return*/
821  xDelete<IssmDouble>(xyz_list);
822  xDelete<IssmDouble>(B);
823  xDelete<IssmDouble>(Bprime);
824  xDelete<IssmDouble>(D);
825  delete gauss;
826  return Ke;
827 }/*}}}*/
828 ElementMatrix* DamageEvolutionAnalysis::CreateMassMatrix(Element* element){/*{{{*/
829 
830  /* Check if ice in element */
831  if(!element->IsIceInElement()) return NULL;
832 
833  /*Intermediaries*/
834  IssmDouble D,Jdet;
835  IssmDouble* xyz_list = NULL;
836 
837  /*Fetch number of nodes and dof for this finite element*/
838  int numnodes = element->GetNumberOfNodes();
839 
840  /*Initialize Element vector and other vectors*/
841  ElementMatrix* Me = element->NewElementMatrix();
842  IssmDouble* basis = xNew<IssmDouble>(numnodes);
843 
844  /*Retrieve all inputs and parameters*/
845  element->GetVerticesCoordinates(&xyz_list);
846 
847  /* Start looping on the number of gaussian points: */
848  Gauss* gauss=element->NewGauss(2);
849  for(int ig=gauss->begin();ig<gauss->end();ig++){
850  gauss->GaussPoint(ig);
851 
852  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
853  element->NodalFunctions(basis,gauss);
854 
855  D=gauss->weight*Jdet;
856  TripleMultiply(basis,1,numnodes,1,
857  &D,1,1,0,
858  basis,1,numnodes,0,
859  &Me->values[0],1);
860  }
861 
862  /*Clean up and return*/
863  xDelete<IssmDouble>(xyz_list);
864  xDelete<IssmDouble>(basis);
865  delete gauss;
866  return Me;
867 }/*}}}*/
868 void DamageEvolutionAnalysis::FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel){/*{{{*/
869 
870  /*Output*/
871  Matrix<IssmDouble>* Kff = NULL;
872  Matrix<IssmDouble>* Kfs = NULL;
873 
874  /*Initialize Jacobian Matrix*/
875  AllocateSystemMatricesx(&Kff,&Kfs,NULL,NULL,femmodel);
876 
877  /*Create and assemble matrix*/
878  for(int i=0;i<femmodel->elements->Size();i++){
879  Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
880  ElementMatrix* Ke = this->CreateFctKMatrix(element);
881  if(Ke) Ke->AddToGlobal(Kff,Kfs);
882  delete Ke;
883  }
884  Kff->Assemble();
885  Kfs->Assemble();
886 
887  /*Assign output pointer*/
888  *pKff=Kff;
889  if(pKfs){
890  *pKfs=Kfs;
891  }
892  else{
893  delete Kfs;
894  }
895 }/*}}}*/
896 void DamageEvolutionAnalysis::MassMatrix(Matrix<IssmDouble>** pMff,FemModel* femmodel){/*{{{*/
897 
898  /*Initialize Mass matrix*/
899  Matrix<IssmDouble> *Mff = NULL;
900  AllocateSystemMatricesx(&Mff,NULL,NULL,NULL,femmodel);
901 
902  /*Create and assemble matrix*/
903  for(int i=0;i<femmodel->elements->Size();i++){
904  Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
905  ElementMatrix* MLe = this->CreateMassMatrix(element);
906  if(MLe){
907  MLe->AddToGlobal(Mff);
908  }
909  delete MLe;
910  }
911  Mff->Assemble();
912 
913  /*Assign output pointer*/
914  *pMff=Mff;
915 }/*}}}*/
DataSet::Size
int Size()
Definition: DataSet.cpp:399
Matrix< IssmDouble >
Element::GetElementType
virtual int GetElementType(void)=0
Element::ComputeStrainRate
void ComputeStrainRate()
Definition: Element.cpp:244
DamageC3Enum
@ DamageC3Enum
Definition: EnumDefinitions.h:109
DeviatoricStressxzEnum
@ DeviatoricStressxzEnum
Definition: EnumDefinitions.h:526
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
Nodes
Declaration of Nodes class.
Definition: Nodes.h:19
Element::GetNumberOfNodes
virtual int GetNumberOfNodes(void)=0
DamageEvolutionAnalysis::UpdateParameters
void UpdateParameters(Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
Definition: DamageEvolutionAnalysis.cpp:83
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
DeviatoricStressyzEnum
@ DeviatoricStressyzEnum
Definition: EnumDefinitions.h:528
DamageEvolutionAnalysis::Core
void Core(FemModel *femmodel)
Definition: DamageEvolutionAnalysis.cpp:122
DamageEvolutionAnalysis::CreateDamageFInput
void CreateDamageFInput(Element *element)
Definition: DamageEvolutionAnalysis.cpp:125
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
DeviatoricStress2Enum
@ DeviatoricStress2Enum
Definition: EnumDefinitions.h:531
MaskIceLevelsetEnum
@ MaskIceLevelsetEnum
Definition: EnumDefinitions.h:641
TimesteppingTimeStepEnum
@ TimesteppingTimeStepEnum
Definition: EnumDefinitions.h:433
DeviatoricStressxxEnum
@ DeviatoricStressxxEnum
Definition: EnumDefinitions.h:524
Constraints
Declaration of Constraints class.
Definition: Constraints.h:13
DamageEvolutionAnalysisEnum
@ DamageEvolutionAnalysisEnum
Definition: EnumDefinitions.h:1026
IoCodeToEnumElementEquation
int IoCodeToEnumElementEquation(int enum_in)
Definition: IoCodeConversions.cpp:311
Elements
Declaration of Elements class.
Definition: Elements.h:17
Matrix::Assemble
void Assemble(void)
Definition: Matrix.h:185
TripleMultiply
int TripleMultiply(IssmDouble *a, int nrowa, int ncola, int itrna, IssmDouble *b, int nrowb, int ncolb, int itrnb, IssmDouble *c, int nrowc, int ncolc, int itrnc, IssmDouble *d, int iaddd)
Definition: MatrixUtils.cpp:20
PressureEnum
@ PressureEnum
Definition: EnumDefinitions.h:664
P1DGEnum
@ P1DGEnum
Definition: EnumDefinitions.h:1215
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
AllocateSystemMatricesx
void AllocateSystemMatricesx(Matrix< IssmDouble > **pKff, Matrix< IssmDouble > **pKfs, Vector< IssmDouble > **pdf, Vector< IssmDouble > **ppf, FemModel *femmodel)
Definition: AllocateSystemMatricesx.cpp:9
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
DamageC2Enum
@ DamageC2Enum
Definition: EnumDefinitions.h:108
Gauss::GaussNode
virtual void GaussNode(int finitelement, int iv)=0
DamageEvolutionAnalysis::UpdateElements
void UpdateElements(Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
Definition: DamageEvolutionAnalysis.cpp:42
ConstantsYtsEnum
@ ConstantsYtsEnum
Definition: EnumDefinitions.h:104
MaterialsRheologyNEnum
@ MaterialsRheologyNEnum
Definition: EnumDefinitions.h:651
DeviatoricStressxyEnum
@ DeviatoricStressxyEnum
Definition: EnumDefinitions.h:525
DamageC1Enum
@ DamageC1Enum
Definition: EnumDefinitions.h:107
DamageStabilizationEnum
@ DamageStabilizationEnum
Definition: EnumDefinitions.h:119
DamageEvolutionNumRequestedOutputsEnum
@ DamageEvolutionNumRequestedOutputsEnum
Definition: EnumDefinitions.h:113
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
DamageEvolutionAnalysis::CreateDVector
ElementVector * CreateDVector(Element *element)
Definition: DamageEvolutionAnalysis.cpp:424
DamageEvolutionAnalysis::CreateKMatrix
ElementMatrix * CreateKMatrix(Element *element)
Definition: DamageEvolutionAnalysis.cpp:431
Parameters::AddObject
void AddObject(Param *newparam)
Definition: Parameters.cpp:67
Input2::GetInputDerivativeValue
virtual void GetInputDerivativeValue(IssmDouble *derivativevalues, IssmDouble *xyz_list, Gauss *gauss)
Definition: Input2.h:37
DamageEvolutionAnalysis::GetSolutionFromInputs
void GetSolutionFromInputs(Vector< IssmDouble > *solution, Element *element)
Definition: DamageEvolutionAnalysis.cpp:725
MaterialsRheologyBbarEnum
@ MaterialsRheologyBbarEnum
Definition: EnumDefinitions.h:644
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
Element
Definition: Element.h:41
Element::NodalFunctions
virtual void NodalFunctions(IssmDouble *basis, Gauss *gauss)=0
DeviatoricStresszzEnum
@ DeviatoricStresszzEnum
Definition: EnumDefinitions.h:529
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
DamageEvolutionAnalysis::CreateDamageFInputArctan
void CreateDamageFInputArctan(Element *element)
Definition: DamageEvolutionAnalysis.cpp:146
Element::ElementSizes
virtual void ElementSizes(IssmDouble *phx, IssmDouble *phy, IssmDouble *phz)=0
DamageEvolutionAnalysis::CreateMassMatrix
ElementMatrix * CreateMassMatrix(Element *element)
Definition: DamageEvolutionAnalysis.cpp:828
Element::NewElementVector
ElementVector * NewElementVector(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2505
IoModel::DeleteData
void DeleteData(int num,...)
Definition: IoModel.cpp:500
IoModel::numberofelements
int numberofelements
Definition: IoModel.h:96
Constraints::ActivatePenaltyMethod
void ActivatePenaltyMethod(int in_analysis)
Definition: Constraints.cpp:22
IoModel::CopyConstantObject
Param * CopyConstantObject(const char *constant_name, int param_enum)
Definition: IoModel.cpp:418
NoneApproximationEnum
@ NoneApproximationEnum
Definition: EnumDefinitions.h:1201
DomainTypeEnum
@ DomainTypeEnum
Definition: EnumDefinitions.h:124
Element::AddInput2
virtual void AddInput2(int input_enum, IssmDouble *values, int interpolation_enum)
Definition: Element.h:216
DamageEvolutionAnalysis::CreateFctKMatrix
ElementMatrix * CreateFctKMatrix(Element *element)
Definition: DamageEvolutionAnalysis.cpp:774
Element::NewGauss
virtual Gauss * NewGauss(void)=0
ApproximationEnum
@ ApproximationEnum
Definition: EnumDefinitions.h:470
DamageEquivStressEnum
@ DamageEquivStressEnum
Definition: EnumDefinitions.h:112
DeviatoricStress1Enum
@ DeviatoricStress1Enum
Definition: EnumDefinitions.h:530
VzEnum
@ VzEnum
Definition: EnumDefinitions.h:853
DeviatoricStressyyEnum
@ DeviatoricStressyyEnum
Definition: EnumDefinitions.h:527
DamageEvolutionAnalysis.h
: header file for generic external result object
GsetEnum
@ GsetEnum
Definition: EnumDefinitions.h:1093
DamageEvolutionAnalysis::GradientJ
void GradientJ(Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
Definition: DamageEvolutionAnalysis.cpp:728
IoModel::FindConstant
void FindConstant(bool *pvalue, const char *constant_name)
Definition: IoModel.cpp:2362
DamageEvolutionAnalysis::InputUpdateFromSolution
void InputUpdateFromSolution(IssmDouble *solution, Element *element)
Definition: DamageEvolutionAnalysis.cpp:731
DamageDEnum
@ DamageDEnum
Definition: EnumDefinitions.h:516
StringArrayParam
Definition: StringArrayParam.h:20
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
DamageLawEnum
@ DamageLawEnum
Definition: EnumDefinitions.h:117
DamageEvolutionAnalysis::CreateDamageFInputPralong
void CreateDamageFInputPralong(Element *element)
Definition: DamageEvolutionAnalysis.cpp:305
DamageFEnum
@ DamageFEnum
Definition: EnumDefinitions.h:520
IntParam
Definition: IntParam.h:20
IoModel::FetchData
void FetchData(bool *pboolean, const char *data_name)
Definition: IoModel.cpp:933
Inputs2
Declaration of Inputs class.
Definition: Inputs2.h:23
DamageEvolutionAnalysis::FctKMatrix
void FctKMatrix(Matrix< IssmDouble > **pKff, Matrix< IssmDouble > **pKfs, FemModel *femmodel)
Definition: DamageEvolutionAnalysis.cpp:868
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
DamageEvolutionAnalysis::CreateDamageFInputExp
void CreateDamageFInputExp(Element *element)
Definition: DamageEvolutionAnalysis.cpp:231
PI
const double PI
Definition: constants.h:11
FemModel::elements
Elements * elements
Definition: FemModel.h:44
DamageDbarEnum
@ DamageDbarEnum
Definition: EnumDefinitions.h:518
Input2
Definition: Input2.h:18
Element::GetDofListLocal
void GetDofListLocal(int **pdoflist, int approximation_enum, int setenum)
Definition: Element.cpp:984
DamageMaxDamageEnum
@ DamageMaxDamageEnum
Definition: EnumDefinitions.h:118
FemModel
Definition: FemModel.h:31
DamageStressThresholdEnum
@ DamageStressThresholdEnum
Definition: EnumDefinitions.h:120
DamageEvolutionAnalysis::DofsPerNode
int DofsPerNode(int **doflist, int domaintype, int approximation)
Definition: DamageEvolutionAnalysis.cpp:39
Element::SetElementInput
virtual void SetElementInput(int enum_in, IssmDouble values)
Definition: Element.h:333
DamageEvolutionAnalysis::GetBprime
void GetBprime(IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss)
Definition: DamageEvolutionAnalysis.cpp:696
DamageC4Enum
@ DamageC4Enum
Definition: EnumDefinitions.h:110
StrainRatexxEnum
@ StrainRatexxEnum
Definition: EnumDefinitions.h:804
Element::IsIceInElement
bool IsIceInElement()
Definition: Element.cpp:2021
IoModel::Data
IssmDouble * Data(const char *data_name)
Definition: IoModel.cpp:437
Loads
Declaration of Loads class.
Definition: Loads.h:16
Inputs2::SetInput
void SetInput(int enum_in, int index, bool value)
Definition: Inputs2.cpp:572
Element::GetSolutionFromInputsOneDof
void GetSolutionFromInputsOneDof(Vector< IssmDouble > *solution, int solutionenum)
Definition: Element.cpp:1281
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
SetActiveNodesLSMx
void SetActiveNodesLSMx(FemModel *femmodel)
Definition: SetActiveNodesLSMx.cpp:12
DamageEvolutionAnalysis::CreatePVector
ElementVector * CreatePVector(Element *element)
Definition: DamageEvolutionAnalysis.cpp:597
MaterialsRheologyBEnum
@ MaterialsRheologyBEnum
Definition: EnumDefinitions.h:643
StrainRateyyEnum
@ StrainRateyyEnum
Definition: EnumDefinitions.h:807
Gauss::begin
virtual int begin(void)=0
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
Element::CharacteristicLength
virtual IssmDouble CharacteristicLength(void)=0
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
min
IssmDouble min(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:14
DamageEvolutionAnalysis::CreateLoads
void CreateLoads(Loads *loads, IoModel *iomodel)
Definition: DamageEvolutionAnalysis.cpp:27
DamageEvolutionRequestedOutputsEnum
@ DamageEvolutionRequestedOutputsEnum
Definition: EnumDefinitions.h:114
StressMaxPrincipalEnum
@ StressMaxPrincipalEnum
Definition: EnumDefinitions.h:810
Element::JacobianDeterminant
virtual void JacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)=0
Gauss::GaussPoint
virtual void GaussPoint(int ig)=0
Element::Update
virtual void Update(Inputs2 *inputs2, int index, IoModel *iomodel, int analysis_counter, int analysis_type, int finite_element)=0
IoModel::FetchDataToInput
void FetchDataToInput(Inputs2 *inputs2, Elements *elements, const char *vector_name, int input_enum)
Definition: IoModel.cpp:1651
IoModelToConstraintsx
void IoModelToConstraintsx(Constraints *constraints, IoModel *iomodel, const char *spc_name, int analysis_type, int finite_element, int dof)
Definition: IoModelToConstraintsx.cpp:10
ElementMatrix::AddToGlobal
void AddToGlobal(Matrix< IssmDouble > *Kff, Matrix< IssmDouble > *Kfs)
Definition: ElementMatrix.cpp:271
ElementVector
Definition: ElementVector.h:20
Input2::GetInputValue
virtual void GetInputValue(IssmDouble *pvalue, Gauss *gauss)
Definition: Input2.h:38
StrainRatexyEnum
@ StrainRatexyEnum
Definition: EnumDefinitions.h:805
DamageEvolutionAnalysis::MassMatrix
void MassMatrix(Matrix< IssmDouble > **pMff, FemModel *femmodel)
Definition: DamageEvolutionAnalysis.cpp:896
Gauss::weight
IssmDouble weight
Definition: Gauss.h:11
IoModel
Definition: IoModel.h:48
DamageEvolutionAnalysis::CreateConstraints
void CreateConstraints(Constraints *constraints, IoModel *iomodel)
Definition: DamageEvolutionAnalysis.cpp:8
DamageEvolutionAnalysis::GetB
void GetB(IssmDouble *B, Element *element, int dim, IssmDouble *xyz_list, Gauss *gauss)
Definition: DamageEvolutionAnalysis.cpp:668
max
IssmDouble max(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:24
DamageEvolutionAnalysis::CreateNodes
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
Definition: DamageEvolutionAnalysis.cpp:32
IoModel::domaintype
int domaintype
Definition: IoModel.h:78
ElementMatrix
Definition: ElementMatrix.h:19
Vector< IssmDouble >
Input2::GetInputAverage
virtual void GetInputAverage(IssmDouble *pvalue)
Definition: Input2.h:33
DamageEvolutionAnalysis::CreateJacobianMatrix
ElementMatrix * CreateJacobianMatrix(Element *element)
Definition: DamageEvolutionAnalysis.cpp:428
DamageEvolutionAnalysis::UpdateConstraints
void UpdateConstraints(FemModel *femmodel)
Definition: DamageEvolutionAnalysis.cpp:769
Gauss
Definition: Gauss.h:8
DamageHealingEnum
@ DamageHealingEnum
Definition: EnumDefinitions.h:115
Element::NodalFunctionsDerivatives
virtual void NodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
DamageKappaEnum
@ DamageKappaEnum
Definition: EnumDefinitions.h:116
DamageStressUBoundEnum
@ DamageStressUBoundEnum
Definition: EnumDefinitions.h:121
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
Element::ComputeDeviatoricStressTensor
virtual void ComputeDeviatoricStressTensor(void)=0
Element::NewElementMatrix
ElementMatrix * NewElementMatrix(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2497
Element::StressMaxPrincipalCreateInput
void StressMaxPrincipalCreateInput(void)
Definition: Element.cpp:4172
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:14 for Ice Sheet System Model by   doxygen 1.8.19