Ice Sheet System Model  4.18
Code documentation
ThermalAnalysis.cpp
Go to the documentation of this file.
1 #include "./ThermalAnalysis.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 ThermalAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
9 
10  /*Intermediary*/
11  int M,N;
12  int finiteelement;
13  iomodel->FindConstant(&finiteelement,"md.thermal.fe");
14  _assert_(finiteelement==P1Enum);
15 
16  /*Output*/
17  IssmDouble *spcvector = NULL;
18  IssmDouble *spcvectorstatic = NULL;
19  IssmDouble *issurface = NULL;
20 
21  /*Only 3d mesh supported*/
22  if(iomodel->domaintype!=Domain3DEnum) _error_("not supported yet");
23 
24  /*Specific case for PDD, we want the constaints to be updated by the PDD scheme itself*/
25  bool isdynamic = false;
26  if(iomodel->solution_enum==ThermalSolutionEnum){
27  /*No PDD scheme, keep default*/
28  }
29  else if(iomodel->solution_enum==SteadystateSolutionEnum){
30  /*No PDD scheme, keep default*/
31  }
32  else if (iomodel->solution_enum==TransientSolutionEnum){
33  int smb_model;
34  iomodel->FindConstant(&smb_model,"md.smb.model");
35  if(smb_model==SMBpddEnum) isdynamic=true;
36  if(smb_model==SMBd18opddEnum) isdynamic=true;
37  if(smb_model==SMBpddSicopolisEnum) isdynamic=true;
38  }
39  else{
40  _error_("Solution "<<EnumToStringx(iomodel->solution_enum)<<" not supported yet");
41  }
42 
43  /*Fetch data: */
44  iomodel->FetchData(&issurface,&M,&N,"md.mesh.vertexonsurface"); _assert_(N>0); _assert_(M==iomodel->numberofvertices);
45  iomodel->FetchData(&spcvector,&M,&N,"md.thermal.spctemperature");
46  iomodel->FetchData(&spcvectorstatic,&M,&N,"md.thermal.spctemperature");
47 
48  /*Convert spcs from temperatures to enthalpy*/
49  _assert_(N>0); _assert_(M>=iomodel->numberofvertices);
50  for(int i=0;i<iomodel->numberofvertices;i++){
51  for(int j=0;j<N;j++){
52  if (isdynamic){
53  if (issurface[i]==1)spcvectorstatic[i*N+j] = NAN;
54  else spcvector[i*N+j] = NAN;
55  }
56  }
57  }
58 
59  if(isdynamic){
60  IoModelToDynamicConstraintsx(constraints,iomodel,spcvector,iomodel->numberofvertices,1,ThermalAnalysisEnum,finiteelement);
61  IoModelToConstraintsx(constraints,iomodel,spcvectorstatic,M,N,ThermalAnalysisEnum,finiteelement);
62  }
63  else{
64  IoModelToConstraintsx(constraints,iomodel,spcvector,M,N,ThermalAnalysisEnum,finiteelement);
65  }
66 
67  /*Free ressources:*/
68  iomodel->DeleteData(spcvector,"md.thermal.spctemperature");
69  iomodel->DeleteData(spcvectorstatic,"md.thermal.spctemperature");
70  iomodel->DeleteData(issurface,"md.mesh.vertexonsurface");
71 
72 }/*}}}*/
73 void ThermalAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
74 
75  if(iomodel->domaintype==Domain2DhorizontalEnum) _error_("2d meshes not supported yet");
76 
77  /*create penalties for nodes: no node can have a temperature over the melting point*/
78  iomodel->FetchData(1,"md.thermal.spctemperature");
79  CreateSingleNodeToElementConnectivity(iomodel);
80 
81  for(int i=0;i<iomodel->numberofvertices;i++){
82 
83  /*keep only this partition's nodes:*/
84  if(iomodel->my_vertices[i]){
85  if (xIsNan<IssmDouble>(iomodel->Data("md.thermal.spctemperature")[i])){ //No penalty applied on spc nodes!
86  loads->AddObject(new Pengrid(i+1,i,iomodel));
87  }
88  }
89  }
90  iomodel->DeleteData(1,"md.thermal.spctemperature");
91 
92 }/*}}}*/
93 void ThermalAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
94 
95  int finiteelement;
96  iomodel->FindConstant(&finiteelement,"md.thermal.fe");
97  _assert_(finiteelement==P1Enum);
98 
99  if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
100  ::CreateNodes(nodes,iomodel,ThermalAnalysisEnum,finiteelement);
101  iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
102 }/*}}}*/
103 int ThermalAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
104  return 1;
105 }/*}}}*/
106 void ThermalAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
107 
108  int frictionlaw,basalforcing_model,materialstype;
109  int FrictionCoupling;
110  /*Now, is the model 3d? otherwise, do nothing: */
111  if(iomodel->domaintype==Domain2DhorizontalEnum)return;
112 
113  /*Update elements: */
114  int finiteelement;
115  iomodel->FindConstant(&finiteelement,"md.thermal.fe");
116  _assert_(finiteelement==P1Enum);
117  int counter=0;
118  for(int i=0;i<iomodel->numberofelements;i++){
119  if(iomodel->my_elements[i]){
120  Element* element=(Element*)elements->GetObjectByOffset(counter);
121  element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,finiteelement);
122  counter++;
123  }
124  }
125 
126  bool dakota_analysis, ismovingfront;
127  iomodel->FindConstant(&dakota_analysis,"md.qmu.isdakota");
128  iomodel->FindConstant(&ismovingfront,"md.transient.ismovingfront");
129  iomodel->FindConstant(&frictionlaw,"md.friction.law");
130  iomodel->FindConstant(&materialstype,"md.materials.type");
131 
132  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.thickness",ThicknessEnum);
133  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.surface",SurfaceEnum);
134  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.base",BaseEnum);
135  iomodel->FetchDataToInput(inputs2,elements,"md.solidearth.sealevel",SealevelEnum,0);
136  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
137  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
138  if(iomodel->domaintype!=Domain2DhorizontalEnum){
139  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
140  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
141  }
142  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
143  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
144  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.pressure",PressureEnum);
145  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.temperature",TemperatureEnum);
146  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vx",VxEnum);
147  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vy",VyEnum);
148  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vz",VzEnum);
149  InputUpdateFromConstantx(inputs2,elements,0.,VxMeshEnum);
150  InputUpdateFromConstantx(inputs2,elements,0.,VyMeshEnum);
151  InputUpdateFromConstantx(inputs2,elements,0.,VzMeshEnum);
152 
153  /*Rheology type*/
154  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
155  switch(materialstype){
156  case MatenhancediceEnum:
157  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
158  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_E",MaterialsRheologyEEnum);
159  break;
160  case MatdamageiceEnum:
161  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
162  break;
163  case MatestarEnum:
164  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_Ec",MaterialsRheologyEcEnum);
165  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_Es",MaterialsRheologyEsEnum);
166  break;
167  case MaticeEnum:
168  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
169  break;
170  default:
171  _error_("not supported");
172  }
173  if(ismovingfront){
174  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum); // required for updating active nodes
175  }
176  /*Basal forcings variables*/
177  iomodel->FindConstant(&basalforcing_model,"md.basalforcings.model");
178  switch(basalforcing_model){
179  case MantlePlumeGeothermalFluxEnum:
180  break;
181  default:
182  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.geothermalflux",BasalforcingsGeothermalfluxEnum);
183  break;
184  }
185  /*Friction law variables*/
186  switch(frictionlaw){
187  case 1:
188  iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
189  iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficient",FrictionCoefficientEnum);
190  iomodel->FetchDataToInput(inputs2,elements,"md.friction.p",FrictionPEnum);
191  iomodel->FetchDataToInput(inputs2,elements,"md.friction.q",FrictionQEnum);
192  if (FrictionCoupling==3){
193  iomodel->FetchDataToInput(inputs2,elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
194  else if(FrictionCoupling==4){
195  iomodel->FetchDataToInput(inputs2,elements,"md.friction.effective_pressure",EffectivePressureEnum);
196  }
197  break;
198  case 2:
199  iomodel->FetchDataToInput(inputs2,elements,"md.friction.C",FrictionCEnum);
200  iomodel->FetchDataToInput(inputs2,elements,"md.friction.m",FrictionMEnum);
201  break;
202  case 3:
203  iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
204  iomodel->FetchDataToInput(inputs2,elements,"md.friction.C",FrictionCEnum);
205  iomodel->FetchDataToInput(inputs2,elements,"md.friction.As",FrictionAsEnum);
206  iomodel->FetchDataToInput(inputs2,elements,"md.friction.q",FrictionQEnum);
207  if (FrictionCoupling==3){
208  iomodel->FetchDataToInput(inputs2,elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
209  else if(FrictionCoupling==4){
210  iomodel->FetchDataToInput(inputs2,elements,"md.friction.effective_pressure",EffectivePressureEnum);
211  }
212  break;
213  case 4:
214  iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficient",FrictionCoefficientEnum);
215  iomodel->FetchDataToInput(inputs2,elements,"md.friction.p",FrictionPEnum);
216  iomodel->FetchDataToInput(inputs2,elements,"md.friction.q",FrictionQEnum);
217  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.pressure",PressureEnum);
218  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.temperature",TemperatureEnum);
219  iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
220  break;
221  case 5:
222  iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficient",FrictionCoefficientEnum);
223  iomodel->FetchDataToInput(inputs2,elements,"md.friction.p",FrictionPEnum);
224  iomodel->FetchDataToInput(inputs2,elements,"md.friction.q",FrictionQEnum);
225  iomodel->FetchDataToInput(inputs2,elements,"md.friction.water_layer",FrictionWaterLayerEnum);
226  break;
227  case 6:
228  iomodel->FetchDataToInput(inputs2,elements,"md.friction.C",FrictionCEnum);
229  iomodel->FetchDataToInput(inputs2,elements,"md.friction.m",FrictionMEnum);
230  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.pressure",PressureEnum);
231  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.temperature",TemperatureEnum);
232  break;
233  case 7:
234  iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
235  iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficient",FrictionCoefficientEnum);
236  iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficientcoulomb",FrictionCoefficientcoulombEnum);
237  iomodel->FetchDataToInput(inputs2,elements,"md.friction.p",FrictionPEnum);
238  iomodel->FetchDataToInput(inputs2,elements,"md.friction.q",FrictionQEnum);
239  if (FrictionCoupling==3){
240  iomodel->FetchDataToInput(inputs2,elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
241  else if(FrictionCoupling==4){
242  iomodel->FetchDataToInput(inputs2,elements,"md.friction.effective_pressure",EffectivePressureEnum);
243  }
244  break;
245  case 9:
246  iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficient",FrictionCoefficientEnum);
247  iomodel->FetchDataToInput(inputs2,elements,"md.friction.pressure_adjusted_temperature",FrictionPressureAdjustedTemperatureEnum);
248  InputUpdateFromConstantx(inputs2,elements,1.,FrictionPEnum);
249  InputUpdateFromConstantx(inputs2,elements,1.,FrictionQEnum);
250  break;
251  default:
252  _error_("friction law not supported");
253  }
254 }/*}}}*/
255 void ThermalAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
256 
257  int numoutputs;
258  char** requestedoutputs = NULL;
259 
260  parameters->AddObject(iomodel->CopyConstantObject("md.thermal.maxiter",ThermalMaxiterEnum));
261  parameters->AddObject(iomodel->CopyConstantObject("md.thermal.stabilization",ThermalStabilizationEnum));
262  parameters->AddObject(iomodel->CopyConstantObject("md.thermal.penalty_factor",ThermalPenaltyFactorEnum));
263  parameters->AddObject(iomodel->CopyConstantObject("md.thermal.penalty_threshold",ThermalPenaltyThresholdEnum));
264  parameters->AddObject(iomodel->CopyConstantObject("md.thermal.penalty_lock",ThermalPenaltyLockEnum));
265  parameters->AddObject(iomodel->CopyConstantObject("md.thermal.isenthalpy",ThermalIsenthalpyEnum));
266  parameters->AddObject(iomodel->CopyConstantObject("md.thermal.isdynamicbasalspc",ThermalIsdynamicbasalspcEnum));
267  parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum));
268 
269  iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.thermal.requested_outputs");
270  parameters->AddObject(new IntParam(ThermalNumRequestedOutputsEnum,numoutputs));
271  if(numoutputs)parameters->AddObject(new StringArrayParam(ThermalRequestedOutputsEnum,requestedoutputs,numoutputs));
272  iomodel->DeleteData(&requestedoutputs,numoutputs,"md.thermal.requested_outputs");
273 
274  /*Deal with friction parameters*/
275  int frictionlaw;
276  iomodel->FindConstant(&frictionlaw,"md.friction.law");
277  if(frictionlaw==6){
278  parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
279  }
280  if(frictionlaw==4){
281  parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
282  parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
283  parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
284  }
285  if(frictionlaw==1 || frictionlaw==3 || frictionlaw==7){
286  parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
287  parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
288  }
289  if(frictionlaw==9){
290  parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
291  parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
292  parameters->AddObject(new IntParam(FrictionCouplingEnum,0));
293  }
294 
295 }/*}}}*/
296 
297 /*Finite Element Analysis*/
298 void ThermalAnalysis::Core(FemModel* femmodel){/*{{{*/
299  _error_("not implemented");
300 }/*}}}*/
301 ElementVector* ThermalAnalysis::CreateDVector(Element* element){/*{{{*/
302  /*Default, return NULL*/
303  return NULL;
304 }/*}}}*/
305 ElementMatrix* ThermalAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
306 _error_("Not implemented");
307 }/*}}}*/
308 ElementMatrix* ThermalAnalysis::CreateKMatrix(Element* element){/*{{{*/
309 
310  /* Check if ice in element */
311  if(!element->IsIceInElement()) return NULL;
312 
313  /*compute all stiffness matrices for this element*/
314  ElementMatrix* Ke1=CreateKMatrixVolume(element);
315  ElementMatrix* Ke2=CreateKMatrixShelf(element);
316  ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
317 
318  /*clean-up and return*/
319  delete Ke1;
320  delete Ke2;
321  return Ke;
322 }/*}}}*/
323 ElementMatrix* ThermalAnalysis::CreateKMatrixShelf(Element* element){/*{{{*/
324 
325  /* Check if ice in element */
326  if(!element->IsIceInElement()) return NULL;
327 
328  /*Initialize Element matrix and return if necessary*/
329  if(!element->IsOnBase() || !element->IsFloating()) return NULL;
330 
331  IssmDouble dt,Jdet,D;
332  IssmDouble *xyz_list_base = NULL;
333 
334  /*Get basal element*/
335  if(!element->IsOnBase() || !element->IsFloating()) return NULL;
336 
337  /*Fetch number of nodes for this finite element*/
338  int numnodes = element->GetNumberOfNodes();
339 
340  /*Initialize vectors*/
341  ElementMatrix* Ke = element->NewElementMatrix();
342  IssmDouble* basis = xNew<IssmDouble>(numnodes);
343 
344  /*Retrieve all inputs and parameters*/
345  element->GetVerticesCoordinatesBase(&xyz_list_base);
346  element->FindParam(&dt,TimesteppingTimeStepEnum);
347  IssmDouble gravity = element->FindParam(ConstantsGEnum);
348  IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
349  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
350  IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum);
351  IssmDouble mixed_layer_capacity= element->FindParam(MaterialsMixedLayerCapacityEnum);
352  IssmDouble thermal_exchange_vel= element->FindParam(MaterialsThermalExchangeVelocityEnum);
353 
354  /* Start looping on the number of gaussian points: */
355  Gauss* gauss=element->NewGaussBase(4);
356  for(int ig=gauss->begin();ig<gauss->end();ig++){
357  gauss->GaussPoint(ig);
358 
359  element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
360  element->NodalFunctions(basis,gauss);
361 
362  D=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel/(heatcapacity*rho_ice);
363  if(reCast<bool,IssmDouble>(dt)) D=dt*D;
364  TripleMultiply(basis,numnodes,1,0,
365  &D,1,1,0,
366  basis,1,numnodes,0,
367  &Ke->values[0],1);
368 
369  }
370 
371  /*Clean up and return*/
372  delete gauss;
373  xDelete<IssmDouble>(basis);
374  xDelete<IssmDouble>(xyz_list_base);
375  return Ke;
376 }/*}}}*/
377 ElementMatrix* ThermalAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
378 
379  /* Check if ice in element */
380  if(!element->IsIceInElement()) return NULL;
381 
382  /*Intermediaries */
383  int stabilization;
384  IssmDouble Jdet,dt,u,v,w,um,vm,wm,vel;
385  IssmDouble h,hx,hy,hz,vx,vy,vz,D_scalar;
386  IssmDouble tau_parameter,diameter;
387  IssmDouble tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver;
388  IssmDouble* xyz_list = NULL;
389 
390  /*Fetch number of nodes and dof for this finite element*/
391  int numnodes = element->GetNumberOfNodes();
392 
393  /*Initialize Element vector and other vectors*/
394  ElementMatrix* Ke = element->NewElementMatrix();
395  IssmDouble* basis = xNew<IssmDouble>(numnodes);
396  IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
397  IssmDouble* Bprime = xNew<IssmDouble>(3*numnodes);
398  IssmDouble K[3][3];
399 
400  /*Retrieve all inputs and parameters*/
401  element->GetVerticesCoordinates(&xyz_list);
402  element->FindParam(&dt,TimesteppingTimeStepEnum);
403  element->FindParam(&stabilization,ThermalStabilizationEnum);
404  IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
405  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
406  IssmDouble gravity = element->FindParam(ConstantsGEnum);
407  IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum);
408  IssmDouble thermalconductivity = element->FindParam(MaterialsThermalconductivityEnum);
409  IssmDouble kappa = thermalconductivity/(rho_ice*heatcapacity);
410  Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
411  Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
412  Input2* vz_input = element->GetInput2(VzEnum); _assert_(vz_input);
413  Input2* vxm_input = element->GetInput2(VxMeshEnum); _assert_(vxm_input);
414  Input2* vym_input = element->GetInput2(VyMeshEnum); _assert_(vym_input);
415  Input2* vzm_input = element->GetInput2(VzMeshEnum); _assert_(vzm_input);
416 
417  /* Start looping on the number of gaussian points: */
418  Gauss* gauss=element->NewGauss(4);
419  for(int ig=gauss->begin();ig<gauss->end();ig++){
420  gauss->GaussPoint(ig);
421 
422  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
423  element->NodalFunctions(basis,gauss);
424  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
425 
426  D_scalar=gauss->weight*Jdet;
427  if(dt!=0.) D_scalar=D_scalar*dt;
428 
429  /*Conduction: */
430  for(int i=0;i<numnodes;i++){
431  for(int j=0;j<numnodes;j++){
432  Ke->values[i*numnodes+j] += D_scalar*kappa*(
433  dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
434  );
435  }
436  }
437 
438  /*Advection: */
439  vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
440  vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
441  vz_input->GetInputValue(&w,gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
442  for(int i=0;i<numnodes;i++){
443  for(int j=0;j<numnodes;j++){
444  Ke->values[i*numnodes+j] += D_scalar*(
445  vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i] +vz*dbasis[2*numnodes+j]*basis[i]
446  );
447  }
448  }
449 
450  /*Transient: */
451  if(dt!=0.){
452  D_scalar=gauss->weight*Jdet;
453  for(int i=0;i<numnodes;i++){
454  for(int j=0;j<numnodes;j++){
455  Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
456  }
457  }
458  D_scalar=D_scalar*dt;
459  }
460 
461  /*Artifficial diffusivity*/
462  if(stabilization==1){
463  element->ElementSizes(&hx,&hy,&hz);
464  vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
465  h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
466  K[0][0]=h/(2.*vel)*fabs(vx*vx); K[0][1]=h/(2.*vel)*fabs(vx*vy); K[0][2]=h/(2.*vel)*fabs(vx*vz);
467  K[1][0]=h/(2.*vel)*fabs(vy*vx); K[1][1]=h/(2.*vel)*fabs(vy*vy); K[1][2]=h/(2.*vel)*fabs(vy*vz);
468  K[2][0]=h/(2.*vel)*fabs(vz*vx); K[2][1]=h/(2.*vel)*fabs(vz*vy); K[2][2]=h/(2.*vel)*fabs(vz*vz);
469  for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j];
470  GetBAdvecprime(Bprime,element,xyz_list,gauss);
471 
472  TripleMultiply(Bprime,3,numnodes,1,
473  &K[0][0],3,3,0,
474  Bprime,3,numnodes,0,
475  &Ke->values[0],1);
476  }
477  else if(stabilization==2){
478  diameter=element->MinEdgeLength(xyz_list);
479  tau_parameter=element->StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);
480  for(int i=0;i<numnodes;i++){
481  for(int j=0;j<numnodes;j++){
482  Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*
483  ((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i])*
484  ((u-um)*dbasis[0*numnodes+j]+(v-vm)*dbasis[1*numnodes+j]+(w-wm)*dbasis[2*numnodes+j]);
485  }
486  }
487  if(dt!=0.){
488  D_scalar=gauss->weight*Jdet;
489  for(int i=0;i<numnodes;i++){
490  for(int j=0;j<numnodes;j++){
491  Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*basis[j]*((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i]);
492  }
493  }
494  }
495  }
496  /*anisotropic SUPG*/
497  else if(stabilization==3){
498  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
499  element->ElementSizes(&hx,&hy,&hz);
500  element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u-um,v-vm,w-wm,hx,hy,hz,kappa);
501  tau_parameter_hor=tau_parameter_anisotropic[0];
502  tau_parameter_ver=tau_parameter_anisotropic[1];
503  for(int i=0;i<numnodes;i++){
504  for(int j=0;j<numnodes;j++){
505  Ke->values[i*numnodes+j]+=D_scalar*
506  (sqrt(tau_parameter_hor)*(u-um)*dbasis[0*numnodes+i]+sqrt(tau_parameter_hor)*(v-vm)*dbasis[1*numnodes+i]+sqrt(tau_parameter_ver)*(w-wm)*dbasis[2*numnodes+i])*
507  (sqrt(tau_parameter_hor)*(u-um)*dbasis[0*numnodes+j]+sqrt(tau_parameter_hor)*(v-vm)*dbasis[1*numnodes+j]+sqrt(tau_parameter_ver)*(w-wm)*dbasis[2*numnodes+j]);
508  }
509  }
510  }
511  }
512 
513  /*Clean up and return*/
514  xDelete<IssmDouble>(xyz_list);
515  xDelete<IssmDouble>(Bprime);
516  xDelete<IssmDouble>(basis);
517  xDelete<IssmDouble>(dbasis);
518  delete gauss;
519  return Ke;
520 }/*}}}*/
521 ElementVector* ThermalAnalysis::CreatePVector(Element* element){/*{{{*/
522 
523  /* Check if ice in element */
524  if(!element->IsIceInElement()) return NULL;
525 
526  /*compute all load vectors for this element*/
527  ElementVector* pe1=CreatePVectorVolume(element);
528  ElementVector* pe2=CreatePVectorSheet(element);
529  ElementVector* pe3=CreatePVectorShelf(element);
530  ElementVector* pe =new ElementVector(pe1,pe2,pe3);
531 
532  /*clean-up and return*/
533  delete pe1;
534  delete pe2;
535  delete pe3;
536  return pe;
537 }/*}}}*/
538 ElementVector* ThermalAnalysis::CreatePVectorSheet(Element* element){/*{{{*/
539 
540  /* Check if ice in element */
541  if(!element->IsIceInElement()) return NULL;
542 
543  /* Geothermal flux on ice sheet base and basal friction */
544  if(!element->IsOnBase() || element->IsFloating()) return NULL;
545 
546  IssmDouble dt,Jdet,geothermalflux,vx,vy,vz;
547  IssmDouble alpha2,scalar,basalfriction,heatflux;
548  IssmDouble *xyz_list_base = NULL;
549 
550  /*Fetch number of nodes for this finite element*/
551  int numnodes = element->GetNumberOfNodes();
552 
553  /*Initialize vectors*/
554  ElementVector* pe = element->NewElementVector();
555  IssmDouble* basis = xNew<IssmDouble>(numnodes);
556 
557  /*Retrieve all inputs and parameters*/
558  element->GetVerticesCoordinatesBase(&xyz_list_base);
559  element->FindParam(&dt,TimesteppingTimeStepEnum);
560  Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
561  Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
562  Input2* vz_input = element->GetInput2(VzEnum); _assert_(vz_input);
563  Input2* geothermalflux_input = element->GetInput2(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
564  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
565  IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum);
566 
567  /*Build friction element, needed later: */
568  Friction* friction=new Friction(element,3);
569 
570  /* Start looping on the number of gaussian points: */
571  Gauss* gauss = element->NewGaussBase(4);
572  for(int ig=gauss->begin();ig<gauss->end();ig++){
573  gauss->GaussPoint(ig);
574 
575  element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
576  element->NodalFunctions(basis,gauss);
577 
578  geothermalflux_input->GetInputValue(&geothermalflux,gauss);
579  friction->GetAlpha2(&alpha2,gauss);
580  vx_input->GetInputValue(&vx,gauss);
581  vy_input->GetInputValue(&vy,gauss);
582  vz_input->GetInputValue(&vz,gauss);
583  vz = 0.;//FIXME
584  basalfriction = alpha2*(vx*vx + vy*vy + vz*vz);
585  heatflux = (basalfriction+geothermalflux)/(rho_ice*heatcapacity);
586 
587  scalar = gauss->weight*Jdet*heatflux;
588  if(dt!=0.) scalar=dt*scalar;
589 
590  for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
591  }
592 
593  /*Clean up and return*/
594  delete gauss;
595  delete friction;
596  xDelete<IssmDouble>(basis);
597  xDelete<IssmDouble>(xyz_list_base);
598  return pe;
599 }/*}}}*/
600 ElementVector* ThermalAnalysis::CreatePVectorShelf(Element* element){/*{{{*/
601 
602  /* Check if ice in element */
603  if(!element->IsIceInElement()) return NULL;
604 
605  IssmDouble t_pmp,dt,Jdet,scalar_ocean,pressure;
606  IssmDouble *xyz_list_base = NULL;
607 
608  /*Get basal element*/
609  if(!element->IsOnBase() || !element->IsFloating()) return NULL;
610 
611  /*Fetch number of nodes for this finite element*/
612  int numnodes = element->GetNumberOfNodes();
613 
614  /*Initialize vectors*/
615  ElementVector* pe = element->NewElementVector();
616  IssmDouble* basis = xNew<IssmDouble>(numnodes);
617 
618  /*Retrieve all inputs and parameters*/
619  element->GetVerticesCoordinatesBase(&xyz_list_base);
620  element->FindParam(&dt,TimesteppingTimeStepEnum);
621  Input2* pressure_input=element->GetInput2(PressureEnum); _assert_(pressure_input);
622  IssmDouble gravity = element->FindParam(ConstantsGEnum);
623  IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
624  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
625  IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum);
626  IssmDouble mixed_layer_capacity= element->FindParam(MaterialsMixedLayerCapacityEnum);
627  IssmDouble thermal_exchange_vel= element->FindParam(MaterialsThermalExchangeVelocityEnum);
628 
629  /* Start looping on the number of gaussian points: */
630  Gauss* gauss=element->NewGaussBase(4);
631  for(int ig=gauss->begin();ig<gauss->end();ig++){
632  gauss->GaussPoint(ig);
633 
634  element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
635  element->NodalFunctions(basis,gauss);
636 
637  pressure_input->GetInputValue(&pressure,gauss);
638  t_pmp=element->TMeltingPoint(pressure);
639 
640  scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel*(t_pmp)/(heatcapacity*rho_ice);
641  if(reCast<bool,IssmDouble>(dt)) scalar_ocean=dt*scalar_ocean;
642 
643  for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_ocean*basis[i];
644  }
645 
646  /*Clean up and return*/
647  delete gauss;
648  xDelete<IssmDouble>(basis);
649  xDelete<IssmDouble>(xyz_list_base);
650  return pe;
651 }/*}}}*/
652 ElementVector* ThermalAnalysis::CreatePVectorVolume(Element* element){/*{{{*/
653 
654  /* Check if ice in element */
655  if(!element->IsIceInElement()) return NULL;
656 
657  /*Intermediaries*/
658  int stabilization;
659  IssmDouble Jdet,phi,dt;
660  IssmDouble temperature;
661  IssmDouble tau_parameter,diameter,hx,hy,hz;
662  IssmDouble tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver;
663  IssmDouble u,v,w;
664  IssmDouble scalar_def,scalar_transient;
665  IssmDouble* xyz_list = NULL;
666 
667  /*Fetch number of nodes and dof for this finite element*/
668  int numnodes = element->GetNumberOfNodes();
669 
670  /*Initialize Element vector*/
671  ElementVector* pe = element->NewElementVector();
672  IssmDouble* basis = xNew<IssmDouble>(numnodes);
673  IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
674 
675  /*Retrieve all inputs and parameters*/
676  element->GetVerticesCoordinates(&xyz_list);
677  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
678  IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum);
679  IssmDouble thermalconductivity = element->FindParam(MaterialsThermalconductivityEnum);
680  IssmDouble kappa = thermalconductivity/(rho_ice*heatcapacity);
681  element->FindParam(&dt,TimesteppingTimeStepEnum);
682  element->FindParam(&stabilization,ThermalStabilizationEnum);
683  Input2* vx_input=element->GetInput2(VxEnum); _assert_(vx_input);
684  Input2* vy_input=element->GetInput2(VyEnum); _assert_(vy_input);
685  Input2* vz_input=element->GetInput2(VzEnum); _assert_(vz_input);
686  Input2* temperature_input = NULL;
687  if(reCast<bool,IssmDouble>(dt)){temperature_input = element->GetInput2(TemperatureEnum); _assert_(temperature_input);}
688 
689  /* Start looping on the number of gaussian points: */
690  Gauss* gauss=element->NewGauss(4);
691  for(int ig=gauss->begin();ig<gauss->end();ig++){
692  gauss->GaussPoint(ig);
693 
694  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
695  element->NodalFunctions(basis,gauss);
696  element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input);
697 
698  scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss->weight;
699  if(reCast<bool,IssmDouble>(dt)) scalar_def=scalar_def*dt;
700 
701  for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_def*basis[i];
702 
703  /* Build transient now */
704  if(reCast<bool,IssmDouble>(dt)){
705  temperature_input->GetInputValue(&temperature, gauss);
706  scalar_transient=temperature*Jdet*gauss->weight;
707  for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_transient*basis[i];
708  }
709 
710  if(stabilization==2){
711  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
712  diameter=element->MinEdgeLength(xyz_list);
713  vx_input->GetInputValue(&u,gauss);
714  vy_input->GetInputValue(&v,gauss);
715  vz_input->GetInputValue(&w,gauss);
716 
717  tau_parameter=element->StabilizationParameter(u,v,w,diameter,kappa);
718 
719  for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
720  if(reCast<bool,IssmDouble>(dt)){
721  for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
722  }
723  }
724  /* anisotropic SUPG */
725  else if(stabilization==3){
726  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
727  element->ElementSizes(&hx,&hy,&hz);
728  vx_input->GetInputValue(&u,gauss);
729  vy_input->GetInputValue(&v,gauss);
730  vz_input->GetInputValue(&w,gauss);
731  element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u,v,w,hx,hy,hz,kappa);
732  tau_parameter_hor=tau_parameter_anisotropic[0];
733  tau_parameter_ver=tau_parameter_anisotropic[1];
734 
735  for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_def*(tau_parameter_hor*u*dbasis[0*numnodes+i]+tau_parameter_hor*v*dbasis[1*numnodes+i]+tau_parameter_ver*w*dbasis[2*numnodes+i]);
736  }
737  }
738 
739  /*Clean up and return*/
740  xDelete<IssmDouble>(basis);
741  xDelete<IssmDouble>(dbasis);
742  xDelete<IssmDouble>(xyz_list);
743  delete gauss;
744  return pe;
745 
746 }/*}}}*/
747 void ThermalAnalysis::GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
748  /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*1.
749  * For node i, Bi' can be expressed in the actual coordinate system
750  * by:
751  * Bi_advec =[ h ]
752  * [ h ]
753  * [ h ]
754  * where h is the interpolation function for node i.
755  *
756  * We assume B has been allocated already, of size: 3x(1*NUMNODESP1)
757  */
758 
759  /*Fetch number of nodes for this finite element*/
760  int numnodes = element->GetNumberOfNodes();
761 
762  /*Get nodal functions*/
763  IssmDouble* basis=xNew<IssmDouble>(numnodes);
764  element->NodalFunctions(basis,gauss);
765 
766  /*Build B: */
767  for(int i=0;i<numnodes;i++){
768  B[numnodes*0+i] = basis[i];
769  B[numnodes*1+i] = basis[i];
770  B[numnodes*2+i] = basis[i];
771  }
772 
773  /*Clean-up*/
774  xDelete<IssmDouble>(basis);
775 }/*}}}*/
776 void ThermalAnalysis::GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
777  /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*1.
778  * For node i, Bi' can be expressed in the actual coordinate system
779  * by:
780  * Biprime_advec=[ dh/dx ]
781  * [ dh/dy ]
782  * [ dh/dz ]
783  * where h is the interpolation function for node i.
784  *
785  * We assume B has been allocated already, of size: 3x(1*numnodes)
786  */
787 
788  /*Fetch number of nodes for this finite element*/
789  int numnodes = element->GetNumberOfNodes();
790 
791  /*Get nodal functions derivatives*/
792  IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
793  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
794 
795  /*Build B: */
796  for(int i=0;i<numnodes;i++){
797  B[numnodes*0+i] = dbasis[0*numnodes+i];
798  B[numnodes*1+i] = dbasis[1*numnodes+i];
799  B[numnodes*2+i] = dbasis[2*numnodes+i];
800  }
801 
802  /*Clean-up*/
803  xDelete<IssmDouble>(dbasis);
804 }/*}}}*/
805 void ThermalAnalysis::GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
806  /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*1.
807  * For node i, Bi' can be expressed in the actual coordinate system
808  * by:
809  * Bi_conduct=[ dh/dx ]
810  * [ dh/dy ]
811  * [ dh/dz ]
812  * where h is the interpolation function for node i.
813  *
814  * We assume B has been allocated already, of size: 3x(1*numnodes)
815  */
816 
817  /*Fetch number of nodes for this finite element*/
818  int numnodes = element->GetNumberOfNodes();
819 
820  /*Get nodal functions derivatives*/
821  IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
822  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
823 
824  /*Build B: */
825  for(int i=0;i<numnodes;i++){
826  B[numnodes*0+i] = dbasis[0*numnodes+i];
827  B[numnodes*1+i] = dbasis[1*numnodes+i];
828  B[numnodes*2+i] = dbasis[2*numnodes+i];
829  }
830 
831  /*Clean-up*/
832  xDelete<IssmDouble>(dbasis);
833 }/*}}}*/
834 void ThermalAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
835  element->GetSolutionFromInputsOneDof(solution,TemperatureEnum);
836 }/*}}}*/
837 void ThermalAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
838  _error_("Not implemented yet");
839 }/*}}}*/
840 void ThermalAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
841 
842  bool converged;
843  int i,rheology_law;
844  int *doflist = NULL;
845  IssmDouble *xyz_list = NULL;
846  bool hack = false;
847 
848  /*Fetch number of nodes and dof for this finite element*/
849  int numnodes = element->GetNumberOfNodes();
850 
851  /*Fetch dof list and allocate solution vector*/
852  element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
853  IssmDouble* values = xNew<IssmDouble>(numnodes);
854  IssmDouble* surface = xNew<IssmDouble>(numnodes);
855  IssmDouble* B = xNew<IssmDouble>(numnodes);
856 
857  /*Use the dof list to index into the solution vector: */
858  for(i=0;i<numnodes;i++){
859  values[i]=solution[doflist[i]];
860 
861  /*Check solution*/
862  if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
863  if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
864  //if(values[i]<0) _printf_("temperature < 0°K found in solution vector\n");
865  //if(values[i]>275) _printf_("temperature > 275°K found in solution vector (Paterson's rheology associated is negative)\n");
866  }
867 
868  /*Force temperature between [Tpmp-50 Tpmp] to disable penalties*/
869  if(hack){
870  IssmDouble* pressure = xNew<IssmDouble>(numnodes);
871  element->GetInputListOnNodes(&pressure[0],PressureEnum);
872  for(i=0;i<numnodes;i++){
873  if(values[i]>element->TMeltingPoint(pressure[i])) values[i]=element->TMeltingPoint(pressure[i]);
874  if(values[i]<element->TMeltingPoint(pressure[i])-50.) values[i]=element->TMeltingPoint(pressure[i])-50.;
875  }
876  xDelete<IssmDouble>(pressure);
877  }
878 
879  /*Get all inputs and parameters*/
880  element->GetInputValue(&converged,ConvergedEnum);
881  if(converged){
882  element->AddInput2(TemperatureEnum,values,element->GetElementType());
883 
884  IssmDouble* n = xNew<IssmDouble>(numnodes);
885  if(element->material->ObjectEnum()==MatestarEnum){
886  for(i=0;i<numnodes;i++) n[i]=3.;
887  }
888  else{
889  element->GetInputListOnNodes(&n[0],MaterialsRheologyNEnum);
890  }
891 
892  /*Update Rheology only if converged (we must make sure that the temperature is below melting point
893  * otherwise the rheology could be negative*/
894  element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
895  element->GetInputListOnNodes(&surface[0],SurfaceEnum);
896 
897  switch(rheology_law){
898  case NoneEnum:
899  /*Do nothing: B is not temperature dependent*/
900  break;
901  case BuddJackaEnum:
902  for(i=0;i<numnodes;i++) B[i]=BuddJacka(values[i]);
903  element->AddInput2(MaterialsRheologyBEnum,&B[0],element->GetElementType());
904  break;
905  case CuffeyEnum:
906  for(i=0;i<numnodes;i++) B[i]=Cuffey(values[i]);
907  element->AddInput2(MaterialsRheologyBEnum,&B[0],element->GetElementType());
908  break;
909  case PatersonEnum:
910  for(i=0;i<numnodes;i++) B[i]=Paterson(values[i]);
911  element->AddInput2(MaterialsRheologyBEnum,&B[0],element->GetElementType());
912  break;
913  case NyeH2OEnum:
914  for(i=0;i<numnodes;i++) B[i]=NyeH2O(values[i]);
915  element->AddInput2(MaterialsRheologyBEnum,&B[0],element->GetElementType());
916  break;
917  case NyeCO2Enum:
918  for(i=0;i<numnodes;i++) B[i]=NyeCO2(values[i]);
919  element->AddInput2(MaterialsRheologyBEnum,&B[0],element->GetElementType());
920  break;
921  case ArrheniusEnum:{
922  element->GetVerticesCoordinates(&xyz_list);
923  for(i=0;i<numnodes;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i*3+2],n[i]);
924  element->AddInput2(MaterialsRheologyBEnum,&B[0],element->GetElementType());
925  break;
926  }
927  default:
928  _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
929  }
930  xDelete<IssmDouble>(n);
931  }
932  else{
933  element->AddInput2(TemperaturePicardEnum,values,element->GetElementType());
934  }
935 
936  /*Free ressources:*/
937  xDelete<IssmDouble>(values);
938  xDelete<IssmDouble>(surface);
939  xDelete<IssmDouble>(B);
940  xDelete<IssmDouble>(xyz_list);
941  xDelete<int>(doflist);
942 }/*}}}*/
943 void ThermalAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
944  SetActiveNodesLSMx(femmodel);
945 }/*}}}*/
MaterialsThermalconductivityEnum
@ MaterialsThermalconductivityEnum
Definition: EnumDefinitions.h:268
FrictionCoefficientEnum
@ FrictionCoefficientEnum
Definition: EnumDefinitions.h:574
ThermalPenaltyFactorEnum
@ ThermalPenaltyFactorEnum
Definition: EnumDefinitions.h:420
IoModel::solution_enum
int solution_enum
Definition: IoModel.h:63
BaseEnum
@ BaseEnum
Definition: EnumDefinitions.h:495
Element::GetElementType
virtual int GetElementType(void)=0
Cuffey
IssmDouble Cuffey(IssmDouble temperature)
Definition: Cuffey.cpp:11
Paterson
IssmDouble Paterson(IssmDouble temperature)
Definition: Paterson.cpp:11
SMBd18opddEnum
@ SMBd18opddEnum
Definition: EnumDefinitions.h:1243
Element::GetInputListOnNodes
void GetInputListOnNodes(IssmDouble *pvalue, int enumtype)
Definition: Element.cpp:1106
ThermalAnalysis::GetSolutionFromInputs
void GetSolutionFromInputs(Vector< IssmDouble > *solution, Element *element)
Definition: ThermalAnalysis.cpp:834
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
ThermalAnalysis::CreateKMatrixVolume
ElementMatrix * CreateKMatrixVolume(Element *element)
Definition: ThermalAnalysis.cpp:377
ThermalAnalysis::InputUpdateFromSolution
void InputUpdateFromSolution(IssmDouble *solution, Element *element)
Definition: ThermalAnalysis.cpp:840
ThermalStabilizationEnum
@ ThermalStabilizationEnum
Definition: EnumDefinitions.h:425
IssmDouble
double IssmDouble
Definition: types.h:37
Element::IsOnBase
bool IsOnBase()
Definition: Element.cpp:1984
Nodes
Declaration of Nodes class.
Definition: Nodes.h:19
Element::GetNumberOfNodes
virtual int GetNumberOfNodes(void)=0
FrictionPEnum
@ FrictionPEnum
Definition: EnumDefinitions.h:578
ThermalRequestedOutputsEnum
@ ThermalRequestedOutputsEnum
Definition: EnumDefinitions.h:424
FrictionEffectivePressureEnum
@ FrictionEffectivePressureEnum
Definition: EnumDefinitions.h:576
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
ThermalAnalysis::GetBAdvecprime
void GetBAdvecprime(IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss)
Definition: ThermalAnalysis.cpp:776
Element::NewGaussBase
virtual Gauss * NewGaussBase(int order)=0
DataSet::AddObject
int AddObject(Object *object)
Definition: DataSet.cpp:252
MaterialsRheologyEcEnum
@ MaterialsRheologyEcEnum
Definition: EnumDefinitions.h:647
MaskOceanLevelsetEnum
@ MaskOceanLevelsetEnum
Definition: EnumDefinitions.h:640
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
TransientSolutionEnum
@ TransientSolutionEnum
Definition: EnumDefinitions.h:1317
MaskIceLevelsetEnum
@ MaskIceLevelsetEnum
Definition: EnumDefinitions.h:641
ThermalAnalysis::CreatePVector
ElementVector * CreatePVector(Element *element)
Definition: ThermalAnalysis.cpp:521
TimesteppingTimeStepEnum
@ TimesteppingTimeStepEnum
Definition: EnumDefinitions.h:433
ThermalNumRequestedOutputsEnum
@ ThermalNumRequestedOutputsEnum
Definition: EnumDefinitions.h:419
Constraints
Declaration of Constraints class.
Definition: Constraints.h:13
ConvergedEnum
@ ConvergedEnum
Definition: EnumDefinitions.h:514
FrictionGammaEnum
@ FrictionGammaEnum
Definition: EnumDefinitions.h:147
ThermalIsdynamicbasalspcEnum
@ ThermalIsdynamicbasalspcEnum
Definition: EnumDefinitions.h:416
MeshVertexonbaseEnum
@ MeshVertexonbaseEnum
Definition: EnumDefinitions.h:653
Elements
Declaration of Elements class.
Definition: Elements.h:17
ThermalAnalysis::CreateJacobianMatrix
ElementMatrix * CreateJacobianMatrix(Element *element)
Definition: ThermalAnalysis.cpp:305
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
MaterialsRhoIceEnum
@ MaterialsRhoIceEnum
Definition: EnumDefinitions.h:264
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
BuddJackaEnum
@ BuddJackaEnum
Definition: EnumDefinitions.h:1000
FrictionQEnum
@ FrictionQEnum
Definition: EnumDefinitions.h:580
ThermalPenaltyThresholdEnum
@ ThermalPenaltyThresholdEnum
Definition: EnumDefinitions.h:422
Element::JacobianDeterminantBase
virtual void JacobianDeterminantBase(IssmDouble *Jdet, IssmDouble *xyz_list_base, Gauss *gauss)=0
MaterialsRheologyNEnum
@ MaterialsRheologyNEnum
Definition: EnumDefinitions.h:651
Element::IsFloating
bool IsFloating()
Definition: Element.cpp:1987
ThermalAnalysis::GetBConduct
void GetBConduct(IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss)
Definition: ThermalAnalysis.cpp:805
FrictionCouplingEnum
@ FrictionCouplingEnum
Definition: EnumDefinitions.h:143
ThermalAnalysis::CreateKMatrixShelf
ElementMatrix * CreateKMatrixShelf(Element *element)
Definition: ThermalAnalysis.cpp:323
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
NyeH2OEnum
@ NyeH2OEnum
Definition: EnumDefinitions.h:1205
Parameters::AddObject
void AddObject(Param *newparam)
Definition: Parameters.cpp:67
ThermalAnalysis::CreatePVectorShelf
ElementVector * CreatePVectorShelf(Element *element)
Definition: ThermalAnalysis.cpp:600
IoModel::my_vertices
bool * my_vertices
Definition: IoModel.h:72
FrictionAsEnum
@ FrictionAsEnum
Definition: EnumDefinitions.h:571
PatersonEnum
@ PatersonEnum
Definition: EnumDefinitions.h:1228
ThermalAnalysis::CreatePVectorVolume
ElementVector * CreatePVectorVolume(Element *element)
Definition: ThermalAnalysis.cpp:652
IoModelToDynamicConstraintsx
void IoModelToDynamicConstraintsx(Constraints *constraints, IoModel *iomodel, const char *spc_name, int analysis_type, int finite_element, int dof)
Definition: IoModelToConstraintsx.cpp:31
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
ThermalAnalysis::UpdateConstraints
void UpdateConstraints(FemModel *femmodel)
Definition: ThermalAnalysis.cpp:943
Element
Definition: Element.h:41
Element::NodalFunctions
virtual void NodalFunctions(IssmDouble *basis, Gauss *gauss)=0
IoModel::numberofvertices
int numberofvertices
Definition: IoModel.h:99
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
VzMeshEnum
@ VzMeshEnum
Definition: EnumDefinitions.h:856
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
NyeCO2
IssmDouble NyeCO2(IssmDouble temperature)
Definition: NyeCO2.cpp:11
VyMeshEnum
@ VyMeshEnum
Definition: EnumDefinitions.h:851
Element::ElementSizes
virtual void ElementSizes(IssmDouble *phx, IssmDouble *phy, IssmDouble *phz)=0
ThermalAnalysis::GetBAdvec
void GetBAdvec(IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss)
Definition: ThermalAnalysis.cpp:747
ThermalSolutionEnum
@ ThermalSolutionEnum
Definition: EnumDefinitions.h:1303
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
FrictionLawEnum
@ FrictionLawEnum
Definition: EnumDefinitions.h:148
ThermalAnalysis.h
: header file for generic external result object
IoModel::CopyConstantObject
Param * CopyConstantObject(const char *constant_name, int param_enum)
Definition: IoModel.cpp:418
ThermalAnalysis::UpdateElements
void UpdateElements(Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
Definition: ThermalAnalysis.cpp:106
ConstantsGEnum
@ ConstantsGEnum
Definition: EnumDefinitions.h:102
BuddJacka
IssmDouble BuddJacka(IssmDouble temperature)
Definition: BuddJacka.cpp:10
NoneApproximationEnum
@ NoneApproximationEnum
Definition: EnumDefinitions.h:1201
ThermalAnalysis::CreateNodes
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
Definition: ThermalAnalysis.cpp:93
ThermalAnalysis::CreateDVector
ElementVector * CreateDVector(Element *element)
Definition: ThermalAnalysis.cpp:301
Element::AddInput2
virtual void AddInput2(int input_enum, IssmDouble *values, int interpolation_enum)
Definition: Element.h:216
MatestarEnum
@ MatestarEnum
Definition: EnumDefinitions.h:1168
Element::NewGauss
virtual Gauss * NewGauss(void)=0
FrictionPressureAdjustedTemperatureEnum
@ FrictionPressureAdjustedTemperatureEnum
Definition: EnumDefinitions.h:579
VzEnum
@ VzEnum
Definition: EnumDefinitions.h:853
MaterialsThermalExchangeVelocityEnum
@ MaterialsThermalExchangeVelocityEnum
Definition: EnumDefinitions.h:267
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
Pengrid
Definition: Pengrid.h:21
ThermalIsenthalpyEnum
@ ThermalIsenthalpyEnum
Definition: EnumDefinitions.h:417
MantlePlumeGeothermalFluxEnum
@ MantlePlumeGeothermalFluxEnum
Definition: EnumDefinitions.h:1158
GsetEnum
@ GsetEnum
Definition: EnumDefinitions.h:1093
MatdamageiceEnum
@ MatdamageiceEnum
Definition: EnumDefinitions.h:1165
IoModel::FindConstant
void FindConstant(bool *pvalue, const char *constant_name)
Definition: IoModel.cpp:2362
StringArrayParam
Definition: StringArrayParam.h:20
NyeH2O
IssmDouble NyeH2O(IssmDouble temperature)
Definition: NyeH2O.cpp:11
Element::StabilizationParameterAnisotropic
virtual void StabilizationParameterAnisotropic(IssmDouble *tau_parameter_anisotropic, IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble hx, IssmDouble hy, IssmDouble hz, IssmDouble kappa)=0
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
MaterialsHeatcapacityEnum
@ MaterialsHeatcapacityEnum
Definition: EnumDefinitions.h:253
ThermalAnalysis::DofsPerNode
int DofsPerNode(int **doflist, int domaintype, int approximation)
Definition: ThermalAnalysis.cpp:103
IntParam
Definition: IntParam.h:20
CuffeyEnum
@ CuffeyEnum
Definition: EnumDefinitions.h:1024
SurfaceEnum
@ SurfaceEnum
Definition: EnumDefinitions.h:823
MaterialsRheologyEEnum
@ MaterialsRheologyEEnum
Definition: EnumDefinitions.h:645
MaterialsRheologyLawEnum
@ MaterialsRheologyLawEnum
Definition: EnumDefinitions.h:262
IoModel::FetchData
void FetchData(bool *pboolean, const char *data_name)
Definition: IoModel.cpp:933
Inputs2
Declaration of Inputs class.
Definition: Inputs2.h:23
InputUpdateFromConstantx
void InputUpdateFromConstantx(FemModel *femmodel, bool constant, int name)
Definition: InputUpdateFromConstantx.cpp:10
MaterialsRheologyEsEnum
@ MaterialsRheologyEsEnum
Definition: EnumDefinitions.h:649
ThermalAnalysis::Core
void Core(FemModel *femmodel)
Definition: ThermalAnalysis.cpp:298
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
MaterialsRhoSeawaterEnum
@ MaterialsRhoSeawaterEnum
Definition: EnumDefinitions.h:265
ThermalPenaltyLockEnum
@ ThermalPenaltyLockEnum
Definition: EnumDefinitions.h:421
Friction
Definition: Friction.h:13
CreateSingleNodeToElementConnectivity
void CreateSingleNodeToElementConnectivity(IoModel *iomodel)
Definition: CreateSingleNodeToElementConnectivity.cpp:10
Object::ObjectEnum
virtual int ObjectEnum()=0
Input2
Definition: Input2.h:18
Element::GetDofListLocal
void GetDofListLocal(int **pdoflist, int approximation_enum, int setenum)
Definition: Element.cpp:984
FemModel
Definition: FemModel.h:31
ThermalAnalysis::CreateLoads
void CreateLoads(Loads *loads, IoModel *iomodel)
Definition: ThermalAnalysis.cpp:73
Element::TMeltingPoint
IssmDouble TMeltingPoint(IssmDouble pressure)
Definition: Element.cpp:4583
ThermalAnalysisEnum
@ ThermalAnalysisEnum
Definition: EnumDefinitions.h:1302
Element::StabilizationParameter
virtual IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa)=0
TemperatureEnum
@ TemperatureEnum
Definition: EnumDefinitions.h:831
SealevelEnum
@ SealevelEnum
Definition: EnumDefinitions.h:675
Element::IsIceInElement
bool IsIceInElement()
Definition: Element.cpp:2021
IoModel::Data
IssmDouble * Data(const char *data_name)
Definition: IoModel.cpp:437
ThermalAnalysis::CreatePVectorSheet
ElementVector * CreatePVectorSheet(Element *element)
Definition: ThermalAnalysis.cpp:538
Loads
Declaration of Loads class.
Definition: Loads.h:16
Element::GetSolutionFromInputsOneDof
void GetSolutionFromInputsOneDof(Vector< IssmDouble > *solution, int solutionenum)
Definition: Element.cpp:1281
NoneEnum
@ NoneEnum
Definition: EnumDefinitions.h:1202
Element::GetVerticesCoordinatesBase
virtual void GetVerticesCoordinatesBase(IssmDouble **xyz_list)=0
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
SetActiveNodesLSMx
void SetActiveNodesLSMx(FemModel *femmodel)
Definition: SetActiveNodesLSMx.cpp:12
MaterialsRheologyBEnum
@ MaterialsRheologyBEnum
Definition: EnumDefinitions.h:643
Gauss::begin
virtual int begin(void)=0
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
MeshVertexonsurfaceEnum
@ MeshVertexonsurfaceEnum
Definition: EnumDefinitions.h:655
Arrhenius
IssmDouble Arrhenius(IssmDouble temperature, IssmDouble depth, IssmDouble n)
Definition: Arrhenius.cpp:9
NyeCO2Enum
@ NyeCO2Enum
Definition: EnumDefinitions.h:1204
Element::JacobianDeterminant
virtual void JacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)=0
FrictionEffectivePressureLimitEnum
@ FrictionEffectivePressureLimitEnum
Definition: EnumDefinitions.h:145
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
ThermalAnalysis::CreateKMatrix
ElementMatrix * CreateKMatrix(Element *element)
Definition: ThermalAnalysis.cpp:308
Element::GetInputValue
void GetInputValue(bool *pvalue, int enum_type)
Definition: Element.cpp:1177
ThicknessEnum
@ ThicknessEnum
Definition: EnumDefinitions.h:840
Friction::GetAlpha2
void GetAlpha2(IssmDouble *palpha2, Gauss *gauss)
Definition: Friction.cpp:164
IoModel::FetchDataToInput
void FetchDataToInput(Inputs2 *inputs2, Elements *elements, const char *vector_name, int input_enum)
Definition: IoModel.cpp:1651
FrictionCoefficientcoulombEnum
@ FrictionCoefficientcoulombEnum
Definition: EnumDefinitions.h:575
MatenhancediceEnum
@ MatenhancediceEnum
Definition: EnumDefinitions.h:1166
ThermalAnalysis::GradientJ
void GradientJ(Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
Definition: ThermalAnalysis.cpp:837
IoModelToConstraintsx
void IoModelToConstraintsx(Constraints *constraints, IoModel *iomodel, const char *spc_name, int analysis_type, int finite_element, int dof)
Definition: IoModelToConstraintsx.cpp:10
ElementVector
Definition: ElementVector.h:20
ThermalMaxiterEnum
@ ThermalMaxiterEnum
Definition: EnumDefinitions.h:418
SteadystateSolutionEnum
@ SteadystateSolutionEnum
Definition: EnumDefinitions.h:1283
MaticeEnum
@ MaticeEnum
Definition: EnumDefinitions.h:1169
Input2::GetInputValue
virtual void GetInputValue(IssmDouble *pvalue, Gauss *gauss)
Definition: Input2.h:38
SMBpddSicopolisEnum
@ SMBpddSicopolisEnum
Definition: EnumDefinitions.h:1253
BasalforcingsGeothermalfluxEnum
@ BasalforcingsGeothermalfluxEnum
Definition: EnumDefinitions.h:477
Gauss::weight
IssmDouble weight
Definition: Gauss.h:11
Element::ViscousHeating
virtual void ViscousHeating(IssmDouble *pphi, IssmDouble *xyz_list, Gauss *gauss, Input2 *vx_input, Input2 *vy_input, Input2 *vz_input)
Definition: Element.h:365
FrictionMEnum
@ FrictionMEnum
Definition: EnumDefinitions.h:577
EffectivePressureEnum
@ EffectivePressureEnum
Definition: EnumDefinitions.h:548
IoModel
Definition: IoModel.h:48
IoModel::domaintype
int domaintype
Definition: IoModel.h:78
FrictionWaterLayerEnum
@ FrictionWaterLayerEnum
Definition: EnumDefinitions.h:583
ElementMatrix
Definition: ElementMatrix.h:19
FrictionCEnum
@ FrictionCEnum
Definition: EnumDefinitions.h:572
Vector< IssmDouble >
ThermalAnalysis::UpdateParameters
void UpdateParameters(Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
Definition: ThermalAnalysis.cpp:255
Gauss
Definition: Gauss.h:8
Element::material
Material * material
Definition: Element.h:50
VxMeshEnum
@ VxMeshEnum
Definition: EnumDefinitions.h:847
ArrheniusEnum
@ ArrheniusEnum
Definition: EnumDefinitions.h:977
ThermalAnalysis::CreateConstraints
void CreateConstraints(Constraints *constraints, IoModel *iomodel)
Definition: ThermalAnalysis.cpp:8
SMBpddEnum
@ SMBpddEnum
Definition: EnumDefinitions.h:1252
Element::NodalFunctionsDerivatives
virtual void NodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
MaterialsMixedLayerCapacityEnum
@ MaterialsMixedLayerCapacityEnum
Definition: EnumDefinitions.h:260
TemperaturePicardEnum
@ TemperaturePicardEnum
Definition: EnumDefinitions.h:833
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
Element::MinEdgeLength
virtual IssmDouble MinEdgeLength(IssmDouble *xyz_list)=0
Element::NewElementMatrix
ElementMatrix * NewElementMatrix(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2497
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:16 for Ice Sheet System Model by   doxygen 1.8.19