Ice Sheet System Model  4.18
Code documentation
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
StressbalanceVerticalAnalysis.cpp
Go to the documentation of this file.
1 #include "./StressbalanceVerticalAnalysis.h"
2 #include "../toolkits/toolkits.h"
3 #include "../classes/classes.h"
4 #include "../shared/shared.h"
5 #include "../modules/modules.h"
6 #include "../solutionsequences/solutionsequences.h"
7 
8 /*Model processing*/
9 void StressbalanceVerticalAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
10 
11  /*Intermediary*/
12  bool isSIA,isSSA,isL1L2,isHO,isFS,iscoupling;
13  int Mz,Nz;
14  IssmDouble *spcvz = NULL;
15 
16  /*return if not 3d mesh*/
17  if(iomodel->domaintype!=Domain3DEnum) return;
18 
19  /*Fetch parameters: */
20  iomodel->FindConstant(&isSIA,"md.flowequation.isSIA");
21  iomodel->FindConstant(&isSSA,"md.flowequation.isSSA");
22  iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2");
23  iomodel->FindConstant(&isHO,"md.flowequation.isHO");
24  iomodel->FindConstant(&isFS,"md.flowequation.isFS");
25 
26  /*Do we have coupling*/
27  if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
28  iscoupling = true;
29  else
30  iscoupling = false;
31 
32  /*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/
33  if(!iscoupling){
34  IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvz",StressbalanceVerticalAnalysisEnum,P1Enum,0);
35  }
36  else{
37  /*Fetch data: */
38  iomodel->FetchData(1,"md.flowequation.borderFS");
39  /*Fetch Spc*/
40  iomodel->FetchData(&spcvz,&Mz,&Nz,"md.stressbalance.spcvz");
41  if(Nz>1) _error_("not supported yet (needs to be coded)");
42 
43  /*Initialize counter*/
44  int count=0;
45 
46  /*Create spcs from x,y,z, as well as the spc values on those spcs: */
47  for(int i=0;i<iomodel->numberofvertices;i++){
48 
49  /*keep only this partition's nodes:*/
50  if(iomodel->my_vertices[i]){
51 
52  if (reCast<int,IssmDouble>(iomodel->Data("md.flowequation.borderFS")[i])){
53  constraints->AddObject(new SpcStatic(count+1,i+1,0,0,StressbalanceVerticalAnalysisEnum)); //spc to zero as vertical velocity is done in Horiz for FS
54  count++;
55  }
56  else if (!xIsNan<IssmDouble>(spcvz[i])){
57  constraints->AddObject(new SpcStatic(count+1,i+1,0,
58  spcvz[i],StressbalanceVerticalAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
59  count++;
60 
61  }
62  }
63  }
64 
65  /*Free data: */
66  iomodel->DeleteData(1,"md.flowequation.borderFS");
67  iomodel->DeleteData(spcvz,"md.stressbalance.spcvz");
68  }
69 
70 }/*}}}*/
71 void StressbalanceVerticalAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
72 
73  /*No loads*/
74 
75 }/*}}}*/
76 void StressbalanceVerticalAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
77 
78  /*return if not 3d mesh*/
79  if(iomodel->domaintype!=Domain3DEnum) return;
80 
81  iomodel->FetchData(3,"md.mesh.vertexonbase","md.mesh.vertexonsurface","md.flowequation.vertex_equation");
82  ::CreateNodes(nodes,iomodel,StressbalanceVerticalAnalysisEnum,P1Enum);
83  iomodel->DeleteData(3,"md.mesh.vertexonbase","md.mesh.vertexonsurface","md.flowequation.vertex_equation");
84 }/*}}}*/
85 int StressbalanceVerticalAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
86  return 1;
87 }/*}}}*/
88 void StressbalanceVerticalAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
89 
90  /*return if not 3d mesh*/
91  if(iomodel->domaintype!=Domain3DEnum) return;
92 
93  /*Update elements: */
94  int counter=0;
95  for(int i=0;i<iomodel->numberofelements;i++){
96  if(iomodel->my_elements[i]){
97  Element* element=(Element*)elements->GetObjectByOffset(counter);
98  element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,P1Enum);
99  counter++;
100  }
101  }
102 
103  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.thickness",ThicknessEnum);
104  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.surface",SurfaceEnum);
105  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.base",BaseEnum);
106  iomodel->FetchDataToInput(inputs2,elements,"md.solidearth.sealevel",SealevelEnum,0);
107  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
108  if(iomodel->domaintype!=Domain2DhorizontalEnum){
109  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
110  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
111  }
112  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
113  //iomodel->FetchDataToInput(inputs2,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
114 
115 
116  /*Add basal forcings to compute melt rate*/
117  int basalforcing_model;
118  iomodel->FindConstant(&basalforcing_model,"md.basalforcings.model");
119  switch(basalforcing_model){
120  case FloatingMeltRateEnum:
121  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.floatingice_melting_rate",BasalforcingsFloatingiceMeltingRateEnum);
122  break;
123  case LinearFloatingMeltRateEnum:
124  break;
125  case MismipFloatingMeltRateEnum:
126  break;
127  case MantlePlumeGeothermalFluxEnum:
128  break;
129  case SpatialLinearFloatingMeltRateEnum:
130  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum);
131  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum);
132  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum);
133  break;
134  case BasalforcingsPicoEnum:
135  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.basin_id",BasalforcingsPicoBasinIdEnum);
136  break;
137  case BasalforcingsIsmip6Enum:
138  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum);
139  break;
140  case BeckmannGoosseFloatingMeltRateEnum:
141  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
142  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
143  break;
144  default:
145  _error_("Basal forcing model "<<EnumToStringx(basalforcing_model)<<" not supported yet");
146  }
147  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vx",VxEnum,0.);
148  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vy",VyEnum,0.);
149 }/*}}}*/
150 void StressbalanceVerticalAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
151 
152  /*No specific parameters*/
153 
154 }/*}}}*/
155 
156 /*Finite Element Analysis*/
157 void StressbalanceVerticalAnalysis::Core(FemModel* femmodel){/*{{{*/
158 
159  if(VerboseSolution()) _printf0_(" computing vertical velocities\n");
160  femmodel->SetCurrentConfiguration(StressbalanceVerticalAnalysisEnum);
161  solutionsequence_linear(femmodel);
162 }/*}}}*/
163 ElementVector* StressbalanceVerticalAnalysis::CreateDVector(Element* element){/*{{{*/
164  /*Default, return NULL*/
165  return NULL;
166 }/*}}}*/
167 ElementMatrix* StressbalanceVerticalAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
168 _error_("Not implemented");
169 }/*}}}*/
170 ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrix(Element* element){/*{{{*/
171 
172  /* Check if ice in element */
173  if(!element->IsIceInElement()) return NULL;
174 
175  bool hack = false;
176 
177  /*compute all stiffness matrices for this element*/
178  ElementMatrix* Ke1=CreateKMatrixVolume(element);
179  ElementMatrix* Ke2=NULL;
180  if(hack) Ke2=CreateKMatrixBase(element);
181  else Ke2=CreateKMatrixSurface(element);
182  ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
183 
184  /*clean-up and return*/
185  delete Ke1;
186  delete Ke2;
187  return Ke;
188 
189 }/*}}}*/
190 ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixBase(Element* element){/*{{{*/
191 
192  if(!element->IsOnBase()) return NULL;
193 
194  /*Intermediaries*/
195  IssmDouble D,Jdet,normal[3];
196  IssmDouble *xyz_list = NULL;
197 
198  /*Fetch number of nodes and dof for this finite element*/
199  int numnodes = element->GetNumberOfNodes();
200 
201  /*Initialize Element matrix and vectors*/
202  ElementMatrix* Ke = element->NewElementMatrix(NoneApproximationEnum);
203  IssmDouble* basis = xNew<IssmDouble>(numnodes);
204 
205  /*Retrieve all inputs and parameters*/
206  element->GetVerticesCoordinatesBase(&xyz_list);
207 
208  /* Start looping on the number of gaussian points: */
209  Gauss* gauss = element->NewGaussBase(2);
210  element->NormalBase(&normal[0],xyz_list);
211  for(int ig=gauss->begin();ig<gauss->end();ig++){
212  gauss->GaussPoint(ig);
213 
214  element->JacobianDeterminantBase(&Jdet,xyz_list,gauss);
215  element->NodalFunctions(basis,gauss);
216  D = -gauss->weight*Jdet*normal[2];
217 
218  for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D*basis[i]*basis[j];
219  }
220 
221  /*Clean up and return*/
222  delete gauss;
223  xDelete<IssmDouble>(xyz_list);
224  xDelete<IssmDouble>(basis);
225  return Ke;
226 }/*}}}*/
227 ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixSurface(Element* element){/*{{{*/
228 
229  if(!element->IsOnSurface()) return NULL;
230 
231  /*Intermediaries*/
232  IssmDouble D,Jdet,normal[3];
233  IssmDouble *xyz_list = NULL;
234 
235  /*Fetch number of nodes and dof for this finite element*/
236  int numnodes = element->GetNumberOfNodes();
237 
238  /*Initialize Element matrix and vectors*/
239  ElementMatrix* Ke = element->NewElementMatrix(NoneApproximationEnum);
240  IssmDouble* basis = xNew<IssmDouble>(numnodes);
241 
242  /*Retrieve all inputs and parameters*/
243  element->GetVerticesCoordinatesTop(&xyz_list);
244 
245  /* Start looping on the number of gaussian points: */
246  Gauss* gauss = element->NewGaussTop(2);
247  element->NormalTop(&normal[0],xyz_list);
248  for(int ig=gauss->begin();ig<gauss->end();ig++){
249  gauss->GaussPoint(ig);
250 
251  element->JacobianDeterminantTop(&Jdet,xyz_list,gauss);
252  element->NodalFunctions(basis,gauss);
253  D = -gauss->weight*Jdet*normal[2];
254 
255  for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D*basis[i]*basis[j];
256  }
257 
258  /*Clean up and return*/
259  delete gauss;
260  xDelete<IssmDouble>(xyz_list);
261  xDelete<IssmDouble>(basis);
262  return Ke;
263 }/*}}}*/
264 ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
265 
266  /*Intermediaries*/
267  IssmDouble Jdet;
268  IssmDouble *xyz_list = NULL;
269 
270  /*Fetch number of nodes and dof for this finite element*/
271  int numnodes = element->GetNumberOfNodes();
272 
273  /*Initialize Element matrix and vectors*/
274  ElementMatrix* Ke = element->NewElementMatrix(NoneApproximationEnum);
275  IssmDouble* basis = xNew<IssmDouble>(numnodes);
276  IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
277 
278  /*Retrieve all inputs and parameters*/
279  element->GetVerticesCoordinates(&xyz_list);
280 
281  /* Start looping on the number of gaussian points: */
282  Gauss* gauss = element->NewGauss(2);
283  for(int ig=gauss->begin();ig<gauss->end();ig++){
284  gauss->GaussPoint(ig);
285 
286  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
287  element->NodalFunctions(basis,gauss);
288  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
289 
290  for(int i=0;i<numnodes;i++){
291  for(int j=0;j<numnodes;j++){
292  Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
293  basis[j]*dbasis[2*numnodes+i]
294  );
295  }
296  }
297  }
298 
299  /*Clean up and return*/
300  delete gauss;
301  xDelete<IssmDouble>(xyz_list);
302  xDelete<IssmDouble>(dbasis);
303  xDelete<IssmDouble>(basis);
304  return Ke;
305 
306 }/*}}}*/
307 ElementVector* StressbalanceVerticalAnalysis::CreatePVector(Element* element){/*{{{*/
308 
309  /* Check if ice in element */
310  if(!element->IsIceInElement()) return NULL;
311 
312  bool hack = false;
313 
314  /*compute all load vectors for this element*/
315  ElementVector* pe1=CreatePVectorVolume(element);
316  ElementVector* pe2=NULL;
317  if(hack) pe2=CreatePVectorSurface(element);
318  else pe2=CreatePVectorBase(element);
319  ElementVector* pe =new ElementVector(pe1,pe2);
320 
321  /*clean-up and return*/
322  delete pe1;
323  delete pe2;
324  return pe;
325 }/*}}}*/
326 ElementVector* StressbalanceVerticalAnalysis::CreatePVectorBase(Element* element){/*{{{*/
327 
328  /*Intermediaries */
329  int approximation;
330  IssmDouble *xyz_list = NULL;
331  IssmDouble *xyz_list_base = NULL;
332  IssmDouble Jdet,slope[3];
333  IssmDouble vx,vy,vz=0.,dbdx,dbdy;
334  IssmDouble gmb,fmb,phi,basalmeltingvalue;
335 
336  if(!element->IsOnBase()) return NULL;
337 
338  /*Fetch number of nodes for this finite element*/
339  int numnodes = element->GetNumberOfNodes();
340 
341  /*Initialize Element vector*/
342  ElementVector* pe = element->NewElementVector();
343  IssmDouble* basis = xNew<IssmDouble>(numnodes);
344 
345  /*Retrieve all inputs and parameters*/
346  element->GetVerticesCoordinates(&xyz_list);
347  element->GetVerticesCoordinatesBase(&xyz_list_base);
348  element->GetInputValue(&approximation,ApproximationEnum);
349  Input2* base_input=element->GetInput2(BaseEnum); _assert_(base_input);
350  Input2* groundedice_input=element->GetInput2(MaskOceanLevelsetEnum); _assert_(groundedice_input);
351  Input2* groundedice_melting_input=element->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(groundedice_melting_input);
352  Input2* floatingice_melting_input=element->GetInput2(BasalforcingsFloatingiceMeltingRateEnum); _assert_(floatingice_melting_input);
353  Input2* vx_input=element->GetInput2(VxEnum); _assert_(vx_input);
354  Input2* vy_input=element->GetInput2(VyEnum); _assert_(vy_input);
355  Input2* vzFS_input=NULL;
356  if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
357  vzFS_input=element->GetInput2(VzFSEnum); _assert_(vzFS_input);
358  }
359 
360  /* Start looping on the number of gaussian points: */
361  Gauss* gauss=element->NewGaussBase(2);
362  for(int ig=gauss->begin();ig<gauss->end();ig++){
363  gauss->GaussPoint(ig);
364 
365  groundedice_melting_input->GetInputValue(&gmb,gauss);
366  floatingice_melting_input->GetInputValue(&fmb,gauss);
367  groundedice_input->GetInputValue(&phi,gauss);
368  base_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
369  vx_input->GetInputValue(&vx,gauss);
370  vy_input->GetInputValue(&vy,gauss);
371  if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
372  vzFS_input->GetInputValue(&vz,gauss);
373  }
374  dbdx=slope[0];
375  dbdy=slope[1];
376  if(phi>0.) basalmeltingvalue=gmb;
377  else basalmeltingvalue=fmb;
378 
379  element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
380  element->NodalFunctions(basis,gauss);
381 
382  for(int i=0;i<numnodes;i++) pe->values[i]+=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-vz-basalmeltingvalue)*basis[i];
383  }
384 
385  /*Clean up and return*/
386  delete gauss;
387  xDelete<IssmDouble>(basis);
388  xDelete<IssmDouble>(xyz_list);
389  xDelete<IssmDouble>(xyz_list_base);
390  return pe;
391 }/*}}}*/
392 ElementVector* StressbalanceVerticalAnalysis::CreatePVectorSurface(Element* element){/*{{{*/
393 
394  /*Intermediaries */
395  int approximation;
396  IssmDouble *xyz_list = NULL;
397  IssmDouble *xyz_list_surface= NULL;
398  IssmDouble Jdet,slope[3];
399  IssmDouble vx,vy,vz=0.,dsdx,dsdy;
400  IssmDouble smb,smbvalue;
401 
402  if(!element->IsOnSurface()) return NULL;
403 
404  /*Fetch number of nodes for this finite element*/
405  int numnodes = element->GetNumberOfNodes();
406 
407  /*Initialize Element vector*/
408  ElementVector* pe = element->NewElementVector();
409  IssmDouble* basis = xNew<IssmDouble>(numnodes);
410 
411  /*Retrieve all inputs and parameters*/
412  element->GetVerticesCoordinates(&xyz_list);
413  element->GetVerticesCoordinatesTop(&xyz_list_surface);
414  element->GetInputValue(&approximation,ApproximationEnum);
415  Input2* surface_input =element->GetInput2(SurfaceEnum); _assert_(surface_input);
416  Input2* smb_input=element->GetInput2(SmbMassBalanceEnum); _assert_(smb_input);
417  Input2* vx_input=element->GetInput2(VxEnum); _assert_(vx_input);
418  Input2* vy_input=element->GetInput2(VyEnum); _assert_(vy_input);
419  Input2* vzFS_input=NULL;
420  if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
421  vzFS_input=element->GetInput2(VzFSEnum); _assert_(vzFS_input);
422  }
423 
424  /* Start looping on the number of gaussian points: */
425  Gauss* gauss=element->NewGaussTop(2);
426  for(int ig=gauss->begin();ig<gauss->end();ig++){
427  gauss->GaussPoint(ig);
428 
429  smb_input->GetInputValue(&smb,gauss);
430  surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
431  vx_input->GetInputValue(&vx,gauss);
432  vy_input->GetInputValue(&vy,gauss);
433  if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
434  vzFS_input->GetInputValue(&vz,gauss);
435  }
436  dsdx=slope[0];
437  dsdy=slope[1];
438 
439  element->JacobianDeterminantTop(&Jdet,xyz_list_surface,gauss);
440  element->NodalFunctions(basis,gauss);
441 
442  for(int i=0;i<numnodes;i++) pe->values[i]+=-Jdet*gauss->weight*(vx*dsdx+vy*dsdy-vz+smb)*basis[i];
443  }
444 
445  /*Clean up and return*/
446  delete gauss;
447  xDelete<IssmDouble>(basis);
448  xDelete<IssmDouble>(xyz_list);
449  xDelete<IssmDouble>(xyz_list_surface);
450  return pe;
451 }/*}}}*/
452 ElementVector* StressbalanceVerticalAnalysis::CreatePVectorVolume(Element* element){/*{{{*/
453 
454  /*Intermediaries*/
455  int approximation;
456  IssmDouble Jdet,dudx,dvdy,dwdz;
457  IssmDouble du[3],dv[3],dw[3];
458  IssmDouble* xyz_list = NULL;
459 
460  /*Fetch number of nodes for this finite element*/
461  int numnodes = element->GetNumberOfNodes();
462 
463  /*Initialize Element vector and basis functions*/
464  ElementVector* pe = element->NewElementVector();
465  IssmDouble* basis = xNew<IssmDouble>(numnodes);
466 
467  /*Retrieve all inputs and parameters*/
468  element->GetVerticesCoordinates(&xyz_list);
469  element->GetInputValue(&approximation,ApproximationEnum);
470  Input2* vx_input=element->GetInput2(VxEnum); _assert_(vx_input);
471  Input2* vy_input=element->GetInput2(VyEnum); _assert_(vy_input);
472  Input2* vzFS_input=NULL;
473  if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
474  vzFS_input=element->GetInput2(VzFSEnum); _assert_(vzFS_input);
475  }
476 
477  /* Start looping on the number of gaussian points: */
478  Gauss* gauss=element->NewGauss(2);
479  for(int ig=gauss->begin();ig<gauss->end();ig++){
480  gauss->GaussPoint(ig);
481 
482  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
483  element->NodalFunctions(basis,gauss);
484 
485  vx_input->GetInputDerivativeValue(&du[0],xyz_list,gauss);
486  vy_input->GetInputDerivativeValue(&dv[0],xyz_list,gauss);
487  if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
488  vzFS_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
489  dwdz=dw[2];
490  }
491  else dwdz=0;
492  dudx=du[0];
493  dvdy=dv[1];
494 
495  for(int i=0;i<numnodes;i++) pe->values[i] += (dudx+dvdy+dwdz)*Jdet*gauss->weight*basis[i];
496  }
497 
498  /*Clean up and return*/
499  delete gauss;
500  xDelete<IssmDouble>(basis);
501  xDelete<IssmDouble>(xyz_list);
502  return pe;
503 }/*}}}*/
504 void StressbalanceVerticalAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
505  element->GetSolutionFromInputsOneDof(solution,VzEnum);
506 }/*}}}*/
507 void StressbalanceVerticalAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
508  _error_("Not implemented yet");
509 }/*}}}*/
510 void StressbalanceVerticalAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
511 
512  int numnodes = element->GetNumberOfNodes();
513  int numdof=numnodes*1;
514 
515  int i;
516  int approximation;
517  int* doflist = NULL;
518  IssmDouble* xyz_list = NULL;
519  IssmDouble rho_ice,g;
520 
521  /*Get the approximation and do nothing if the element in FS or None*/
522  element->GetInputValue(&approximation,ApproximationEnum);
523  if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
524  return;
525  }
526 
527  /*Get dof list and vertices coordinates: */
528  element->GetVerticesCoordinates(&xyz_list);
529  element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
530  IssmDouble* values = xNew<IssmDouble>(numdof);
531  IssmDouble* vx = xNew<IssmDouble>(numnodes);
532  IssmDouble* vy = xNew<IssmDouble>(numnodes);
533  IssmDouble* vz = xNew<IssmDouble>(numnodes);
534  IssmDouble* vzSSA = xNew<IssmDouble>(numnodes);
535  IssmDouble* vzHO = xNew<IssmDouble>(numnodes);
536  IssmDouble* vzFS = xNew<IssmDouble>(numnodes);
537  IssmDouble* vel = xNew<IssmDouble>(numnodes);
538  IssmDouble* pressure = xNew<IssmDouble>(numnodes);
539  IssmDouble* surface = xNew<IssmDouble>(numnodes);
540 
541  /*Use the dof list to index into the solution vector vz: */
542  for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
543  for(i=0;i<numdof;i++){
544  vz[i]=values[i*1+0];
545 
546  /*Check solution*/
547  if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
548  if(xIsInf<IssmDouble>(vz[i])) _error_("Inf found in solution vector");
549  }
550 
551  /*Get Vx and Vy*/
552  element->GetInputListOnNodes(&vx[0],VxEnum,0.0); //default is 0
553  element->GetInputListOnNodes(&vy[0],VyEnum,0.0); //default is 0
554 
555  /*Do some modifications if we actually have a HOFS or SSAFS element*/
556  if(approximation==HOFSApproximationEnum){
557  Input2* vzFS_input=element->GetInput2(VzFSEnum);
558  if (vzFS_input){
559  if (vzFS_input->ObjectEnum()!=PentaInput2Enum) _error_("Cannot compute Vel as VzFS is of type " << EnumToStringx(vzFS_input->ObjectEnum()));
560  element->GetInputListOnNodes(&vzFS[0],VzFSEnum,0.);
561  }
562  else _error_("Cannot compute Vz as VzFS in not present in HOFS element");
563  for(i=0;i<numnodes;i++){
564  vzHO[i]=vz[i];
565  vz[i]=vzHO[i]+vzFS[i];
566  }
567  }
568  else if(approximation==SSAFSApproximationEnum){
569  Input2* vzFS_input=element->GetInput2(VzFSEnum);
570  if (vzFS_input){
571  if (vzFS_input->ObjectEnum()!=PentaInput2Enum) _error_("Cannot compute Vel as VzFS is of type " << EnumToStringx(vzFS_input->ObjectEnum()));
572  element->GetInputListOnNodes(&vzFS[0],VzFSEnum,0.);
573  }
574  else _error_("Cannot compute Vz as VzFS in not present in SSAFS element");
575  for(i=0;i<numnodes;i++){
576  vzSSA[i]=vz[i];
577  vz[i]=vzSSA[i]+vzFS[i];
578  }
579  }
580 
581  /*Now Compute vel*/
582  for(i=0;i<numnodes;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
583 
584  /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D,
585  *so the pressure is just the pressure at the z elevation: except it this is a HOFS element */
586  if(approximation!=HOFSApproximationEnum && approximation!=SSAFSApproximationEnum){
587  rho_ice = element->FindParam(MaterialsRhoIceEnum);
588  g = element->FindParam(ConstantsGEnum);
589  element->GetInputListOnNodes(&surface[0],SurfaceEnum,0.);
590  for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
591  }
592  if(approximation!=HOFSApproximationEnum && approximation!=SSAFSApproximationEnum){
593  element->AddInput2(PressureEnum,pressure,element->GetElementType());
594  }
595  else if(approximation==HOFSApproximationEnum){
596  element->AddInput2(VzHOEnum,vzHO,P1Enum);
597  }
598  else if(approximation==SSAFSApproximationEnum){
599  element->AddInput2(VzSSAEnum,vzSSA,P1Enum);
600  }
601  element->AddInput2(VzEnum,vz,P1Enum);
602  element->AddInput2(VelEnum,vel,P1Enum);
603 
604  /*Free ressources:*/
605  xDelete<IssmDouble>(surface);
606  xDelete<IssmDouble>(pressure);
607  xDelete<IssmDouble>(vel);
608  xDelete<IssmDouble>(vz);
609  xDelete<IssmDouble>(vzSSA);
610  xDelete<IssmDouble>(vzHO);
611  xDelete<IssmDouble>(vzFS);
612  xDelete<IssmDouble>(vy);
613  xDelete<IssmDouble>(vx);
614  xDelete<IssmDouble>(values);
615  xDelete<IssmDouble>(xyz_list);
616  xDelete<int>(doflist);
617 }/*}}}*/
618 void StressbalanceVerticalAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
619  /*Default, do nothing*/
620  SetActiveNodesLSMx(femmodel);
621  return;
622 }/*}}}*/
BasalforcingsPicoBasinIdEnum
@ BasalforcingsPicoBasinIdEnum
Definition: EnumDefinitions.h:486
BasalforcingsOceanTempEnum
@ BasalforcingsOceanTempEnum
Definition: EnumDefinitions.h:485
BaseEnum
@ BaseEnum
Definition: EnumDefinitions.h:495
Element::GetElementType
virtual int GetElementType(void)=0
SmbMassBalanceEnum
@ SmbMassBalanceEnum
Definition: EnumDefinitions.h:748
Element::GetInputListOnNodes
void GetInputListOnNodes(IssmDouble *pvalue, int enumtype)
Definition: Element.cpp:1106
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
Element::IsOnBase
bool IsOnBase()
Definition: Element.cpp:1984
Nodes
Declaration of Nodes class.
Definition: Nodes.h:19
StressbalanceVerticalAnalysis::CreateJacobianMatrix
ElementMatrix * CreateJacobianMatrix(Element *element)
Definition: StressbalanceVerticalAnalysis.cpp:167
Element::GetNumberOfNodes
virtual int GetNumberOfNodes(void)=0
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
StressbalanceVerticalAnalysis::Core
void Core(FemModel *femmodel)
Definition: StressbalanceVerticalAnalysis.cpp:157
Element::NewGaussBase
virtual Gauss * NewGaussBase(int order)=0
DataSet::AddObject
int AddObject(Object *object)
Definition: DataSet.cpp:252
MaskOceanLevelsetEnum
@ MaskOceanLevelsetEnum
Definition: EnumDefinitions.h:640
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
MaskIceLevelsetEnum
@ MaskIceLevelsetEnum
Definition: EnumDefinitions.h:641
StressbalanceVerticalAnalysis::CreateLoads
void CreateLoads(Loads *loads, IoModel *iomodel)
Definition: StressbalanceVerticalAnalysis.cpp:71
Constraints
Declaration of Constraints class.
Definition: Constraints.h:13
MeshVertexonbaseEnum
@ MeshVertexonbaseEnum
Definition: EnumDefinitions.h:653
Elements
Declaration of Elements class.
Definition: Elements.h:17
StressbalanceVerticalAnalysis::CreateConstraints
void CreateConstraints(Constraints *constraints, IoModel *iomodel)
Definition: StressbalanceVerticalAnalysis.cpp:9
PressureEnum
@ PressureEnum
Definition: EnumDefinitions.h:664
MaterialsRhoIceEnum
@ MaterialsRhoIceEnum
Definition: EnumDefinitions.h:264
MismipFloatingMeltRateEnum
@ MismipFloatingMeltRateEnum
Definition: EnumDefinitions.h:1190
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
VzFSEnum
@ VzFSEnum
Definition: EnumDefinitions.h:854
BasalforcingsPicoEnum
@ BasalforcingsPicoEnum
Definition: EnumDefinitions.h:990
StressbalanceVerticalAnalysis::InputUpdateFromSolution
void InputUpdateFromSolution(IssmDouble *solution, Element *element)
Definition: StressbalanceVerticalAnalysis.cpp:510
Element::JacobianDeterminantBase
virtual void JacobianDeterminantBase(IssmDouble *Jdet, IssmDouble *xyz_list_base, Gauss *gauss)=0
Element::JacobianDeterminantTop
virtual void JacobianDeterminantTop(IssmDouble *Jdet, IssmDouble *xyz_list_base, Gauss *gauss)=0
Element::NormalBase
virtual void NormalBase(IssmDouble *normal, IssmDouble *xyz_list)=0
BasalforcingsUpperwaterElevationEnum
@ BasalforcingsUpperwaterElevationEnum
Definition: EnumDefinitions.h:94
StressbalanceVerticalAnalysis::CreateDVector
ElementVector * CreateDVector(Element *element)
Definition: StressbalanceVerticalAnalysis.cpp:163
FSApproximationEnum
@ FSApproximationEnum
Definition: EnumDefinitions.h:1060
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
BasalforcingsDeepwaterMeltingRateEnum
@ BasalforcingsDeepwaterMeltingRateEnum
Definition: EnumDefinitions.h:62
solutionsequence_linear
void solutionsequence_linear(FemModel *femmodel)
Definition: solutionsequence_linear.cpp:10
VzSSAEnum
@ VzSSAEnum
Definition: EnumDefinitions.h:857
SpcStatic
Definition: SpcStatic.h:13
IoModel::my_vertices
bool * my_vertices
Definition: IoModel.h:72
Input2::GetInputDerivativeValue
virtual void GetInputDerivativeValue(IssmDouble *derivativevalues, IssmDouble *xyz_list, Gauss *gauss)
Definition: Input2.h:37
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
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
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
Element::NewElementVector
ElementVector * NewElementVector(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2505
IoModel::DeleteData
void DeleteData(int num,...)
Definition: IoModel.cpp:500
VzHOEnum
@ VzHOEnum
Definition: EnumDefinitions.h:855
IoModel::numberofelements
int numberofelements
Definition: IoModel.h:96
ConstantsGEnum
@ ConstantsGEnum
Definition: EnumDefinitions.h:102
BasalforcingsGroundediceMeltingRateEnum
@ BasalforcingsGroundediceMeltingRateEnum
Definition: EnumDefinitions.h:478
NoneApproximationEnum
@ NoneApproximationEnum
Definition: EnumDefinitions.h:1201
Element::AddInput2
virtual void AddInput2(int input_enum, IssmDouble *values, int interpolation_enum)
Definition: Element.h:216
Element::NewGauss
virtual Gauss * NewGauss(void)=0
ApproximationEnum
@ ApproximationEnum
Definition: EnumDefinitions.h:470
StressbalanceVerticalAnalysis::CreatePVectorSurface
ElementVector * CreatePVectorSurface(Element *element)
Definition: StressbalanceVerticalAnalysis.cpp:392
FloatingMeltRateEnum
@ FloatingMeltRateEnum
Definition: EnumDefinitions.h:1069
BasalforcingsOceanSalinityEnum
@ BasalforcingsOceanSalinityEnum
Definition: EnumDefinitions.h:484
VzEnum
@ VzEnum
Definition: EnumDefinitions.h:853
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
LinearFloatingMeltRateEnum
@ LinearFloatingMeltRateEnum
Definition: EnumDefinitions.h:1143
BeckmannGoosseFloatingMeltRateEnum
@ BeckmannGoosseFloatingMeltRateEnum
Definition: EnumDefinitions.h:991
MantlePlumeGeothermalFluxEnum
@ MantlePlumeGeothermalFluxEnum
Definition: EnumDefinitions.h:1158
StressbalanceVerticalAnalysis::GradientJ
void GradientJ(Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
Definition: StressbalanceVerticalAnalysis.cpp:507
GsetEnum
@ GsetEnum
Definition: EnumDefinitions.h:1093
BasalforcingsFloatingiceMeltingRateEnum
@ BasalforcingsFloatingiceMeltingRateEnum
Definition: EnumDefinitions.h:476
IoModel::FindConstant
void FindConstant(bool *pvalue, const char *constant_name)
Definition: IoModel.cpp:2362
Element::NormalTop
virtual void NormalTop(IssmDouble *normal, IssmDouble *xyz_list)=0
StressbalanceVerticalAnalysis::CreateKMatrixBase
ElementMatrix * CreateKMatrixBase(Element *element)
Definition: StressbalanceVerticalAnalysis.cpp:190
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
SpatialLinearFloatingMeltRateEnum
@ SpatialLinearFloatingMeltRateEnum
Definition: EnumDefinitions.h:1278
SurfaceEnum
@ SurfaceEnum
Definition: EnumDefinitions.h:823
IoModel::FetchData
void FetchData(bool *pboolean, const char *data_name)
Definition: IoModel.cpp:933
Inputs2
Declaration of Inputs class.
Definition: Inputs2.h:23
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
BasalforcingsIsmip6Enum
@ BasalforcingsIsmip6Enum
Definition: EnumDefinitions.h:989
StressbalanceVerticalAnalysis::DofsPerNode
int DofsPerNode(int **doflist, int domaintype, int approximation)
Definition: StressbalanceVerticalAnalysis.cpp:85
SSAFSApproximationEnum
@ SSAFSApproximationEnum
Definition: EnumDefinitions.h:1256
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
SealevelEnum
@ SealevelEnum
Definition: EnumDefinitions.h:675
StressbalanceVerticalAnalysis.h
: header file for generic external result object
BasalforcingsIsmip6BasinIdEnum
@ BasalforcingsIsmip6BasinIdEnum
Definition: EnumDefinitions.h:480
Element::IsIceInElement
bool IsIceInElement()
Definition: Element.cpp:2021
StressbalanceVerticalAnalysis::CreateKMatrixSurface
ElementMatrix * CreateKMatrixSurface(Element *element)
Definition: StressbalanceVerticalAnalysis.cpp:227
IoModel::Data
IssmDouble * Data(const char *data_name)
Definition: IoModel.cpp:437
PentaInput2Enum
@ PentaInput2Enum
Definition: EnumDefinitions.h:1125
Loads
Declaration of Loads class.
Definition: Loads.h:16
Element::GetSolutionFromInputsOneDof
void GetSolutionFromInputsOneDof(Vector< IssmDouble > *solution, int solutionenum)
Definition: Element.cpp:1281
Element::GetVerticesCoordinatesBase
virtual void GetVerticesCoordinatesBase(IssmDouble **xyz_list)=0
StressbalanceVerticalAnalysisEnum
@ StressbalanceVerticalAnalysisEnum
Definition: EnumDefinitions.h:1289
BasalforcingsDeepwaterElevationEnum
@ BasalforcingsDeepwaterElevationEnum
Definition: EnumDefinitions.h:61
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
SetActiveNodesLSMx
void SetActiveNodesLSMx(FemModel *femmodel)
Definition: SetActiveNodesLSMx.cpp:12
StressbalanceVerticalAnalysis::UpdateConstraints
void UpdateConstraints(FemModel *femmodel)
Definition: StressbalanceVerticalAnalysis.cpp:618
StressbalanceVerticalAnalysis::CreatePVectorBase
ElementVector * CreatePVectorBase(Element *element)
Definition: StressbalanceVerticalAnalysis.cpp:326
StressbalanceVerticalAnalysis::CreatePVector
ElementVector * CreatePVector(Element *element)
Definition: StressbalanceVerticalAnalysis.cpp:307
Gauss::begin
virtual int begin(void)=0
VerboseSolution
bool VerboseSolution(void)
Definition: Verbosity.cpp:24
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
MeshVertexonsurfaceEnum
@ MeshVertexonsurfaceEnum
Definition: EnumDefinitions.h:655
StressbalanceVerticalAnalysis::CreatePVectorVolume
ElementVector * CreatePVectorVolume(Element *element)
Definition: StressbalanceVerticalAnalysis.cpp:452
StressbalanceVerticalAnalysis::CreateKMatrixVolume
ElementMatrix * CreateKMatrixVolume(Element *element)
Definition: StressbalanceVerticalAnalysis.cpp:264
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
FemModel::SetCurrentConfiguration
void SetCurrentConfiguration(int configuration_type)
Definition: FemModel.cpp:634
Element::GetInputValue
void GetInputValue(bool *pvalue, int enum_type)
Definition: Element.cpp:1177
HOFSApproximationEnum
@ HOFSApproximationEnum
Definition: EnumDefinitions.h:1096
ThicknessEnum
@ ThicknessEnum
Definition: EnumDefinitions.h:840
StressbalanceVerticalAnalysis::UpdateParameters
void UpdateParameters(Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
Definition: StressbalanceVerticalAnalysis.cpp:150
IoModel::FetchDataToInput
void FetchDataToInput(Inputs2 *inputs2, Elements *elements, const char *vector_name, int input_enum)
Definition: IoModel.cpp:1651
VelEnum
@ VelEnum
Definition: EnumDefinitions.h:844
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
Input2::GetInputValue
virtual void GetInputValue(IssmDouble *pvalue, Gauss *gauss)
Definition: Input2.h:38
Gauss::weight
IssmDouble weight
Definition: Gauss.h:11
Element::GetVerticesCoordinatesTop
virtual void GetVerticesCoordinatesTop(IssmDouble **xyz_list)=0
Element::IsOnSurface
bool IsOnSurface()
Definition: Element.cpp:1981
IoModel
Definition: IoModel.h:48
Element::NewGaussTop
virtual Gauss * NewGaussTop(int order)=0
StressbalanceVerticalAnalysis::GetSolutionFromInputs
void GetSolutionFromInputs(Vector< IssmDouble > *solution, Element *element)
Definition: StressbalanceVerticalAnalysis.cpp:504
IoModel::domaintype
int domaintype
Definition: IoModel.h:78
StressbalanceVerticalAnalysis::CreateKMatrix
ElementMatrix * CreateKMatrix(Element *element)
Definition: StressbalanceVerticalAnalysis.cpp:170
ElementMatrix
Definition: ElementMatrix.h:19
Vector< IssmDouble >
Gauss
Definition: Gauss.h:8
Element::NodalFunctionsDerivatives
virtual void NodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
StressbalanceVerticalAnalysis::UpdateElements
void UpdateElements(Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
Definition: StressbalanceVerticalAnalysis.cpp:88
StressbalanceVerticalAnalysis::CreateNodes
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
Definition: StressbalanceVerticalAnalysis.cpp:76
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
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