Ice Sheet System Model  4.18
Code documentation
LevelsetAnalysis.cpp
Go to the documentation of this file.
1 #ifdef HAVE_CONFIG_H
2 #include <config.h>
3 #else
4 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
5 #endif
6 #include "./LevelsetAnalysis.h"
7 #include "../toolkits/toolkits.h"
8 #include "../classes/classes.h"
9 #include "../shared/shared.h"
10 #include "../modules/modules.h"
11 #include "../solutionsequences/solutionsequences.h"
12 
13 void LevelsetAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
14  int finiteelement;
15  iomodel->FindConstant(&finiteelement,"md.levelset.fe");
16  IoModelToConstraintsx(constraints,iomodel,"md.levelset.spclevelset",LevelsetAnalysisEnum,finiteelement);
17 }
18 /*}}}*/
19 void LevelsetAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
20  return;
21 }/*}}}*/
22 void LevelsetAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
23  int finiteelement;
24  iomodel->FindConstant(&finiteelement,"md.levelset.fe");
25  if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
26  ::CreateNodes(nodes,iomodel,LevelsetAnalysisEnum,finiteelement);
27  iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
28 }
29 /*}}}*/
30 int LevelsetAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
31  return 1;
32 }
33 /*}}}*/
34 void LevelsetAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
35 
36  /*Finite element type*/
37  int finiteelement;
38  iomodel->FindConstant(&finiteelement,"md.levelset.fe");
39 
40  /*Update elements: */
41  int counter=0;
42  for(int i=0;i<iomodel->numberofelements;i++){
43  if(iomodel->my_elements[i]){
44  Element* element=(Element*)elements->GetObjectByOffset(counter);
45  element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,finiteelement);
46  counter++;
47  }
48  }
49 
50  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
51  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vx",VxEnum);
52  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vy",VyEnum);
53 
54  /*Get moving front parameters*/
55  int calvinglaw;
56  iomodel->FindConstant(&calvinglaw,"md.calving.law");
57  switch(calvinglaw){
58  case DefaultCalvingEnum:
59  iomodel->FetchDataToInput(inputs2,elements,"md.calving.calvingrate",CalvingCalvingrateEnum);
60  break;
61  case CalvingLevermannEnum:
62  iomodel->FetchDataToInput(inputs2,elements,"md.calving.coeff",CalvinglevermannCoeffEnum);
63  break;
64  case CalvingVonmisesEnum:
65  iomodel->FetchDataToInput(inputs2,elements,"md.calving.stress_threshold_groundedice",CalvingStressThresholdGroundediceEnum);
66  iomodel->FetchDataToInput(inputs2,elements,"md.calving.stress_threshold_floatingice",CalvingStressThresholdFloatingiceEnum);
67  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.bed",BedEnum);
68  break;
69  case CalvingMinthicknessEnum:
70  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.bed",BedEnum);
71  break;
72  case CalvingHabEnum:
73  iomodel->FetchDataToInput(inputs2,elements,"md.calving.flotation_fraction",CalvingHabFractionEnum);
74  break;
75  case CalvingCrevasseDepthEnum:
76  iomodel->FetchDataToInput(inputs2,elements,"md.calving.water_height",WaterheightEnum);
77  break;
78  case CalvingDev2Enum:
79  iomodel->FetchDataToInput(inputs2,elements,"md.calving.stress_threshold_groundedice",CalvingStressThresholdGroundediceEnum);
80  iomodel->FetchDataToInput(inputs2,elements,"md.calving.stress_threshold_floatingice",CalvingStressThresholdFloatingiceEnum);
81  break;
82  default:
83  _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
84  }
85 
86  /*Get frontal melt parameters*/
87  int melt_parameterization;
88  iomodel->FindConstant(&melt_parameterization,"md.frontalforcings.parameterization");
89  switch(melt_parameterization){
90  case FrontalForcingsDefaultEnum:
91  iomodel->FetchDataToInput(inputs2,elements,"md.frontalforcings.meltingrate",CalvingMeltingrateEnum);
92  break;
93  case FrontalForcingsRignotEnum:
94  iomodel->FetchDataToInput(inputs2,elements,"md.frontalforcings.basin",FrontalForcingsBasinIdEnum);
95  iomodel->FetchDataToInput(inputs2,elements,"md.frontalforcings.subglacial_discharge",FrontalForcingsSubglacialDischargeEnum);
96  iomodel->FetchDataToInput(inputs2,elements,"md.frontalforcings.thermalforcing",FrontalForcingsThermalForcingEnum);
97  break;
98  default:
99  _error_("Frontal forcings"<<EnumToStringx(melt_parameterization)<<" not supported yet");
100  }
101 }
102 /*}}}*/
103 void LevelsetAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
104 
105  parameters->AddObject(iomodel->CopyConstantObject("md.levelset.stabilization",LevelsetStabilizationEnum));
106  parameters->AddObject(iomodel->CopyConstantObject("md.levelset.reinit_frequency",LevelsetReinitFrequencyEnum));
107  parameters->AddObject(iomodel->CopyConstantObject("md.levelset.kill_icebergs",LevelsetKillIcebergsEnum));
108  parameters->AddObject(iomodel->CopyConstantObject("md.levelset.calving_max",CalvingMaxEnum));
109 
110  int calvinglaw;
111  iomodel->FindConstant(&calvinglaw,"md.calving.law");
112  switch(calvinglaw){
113  case DefaultCalvingEnum:
114  case CalvingLevermannEnum:
115  break;
116  case CalvingVonmisesEnum:
117  parameters->AddObject(iomodel->CopyConstantObject("md.calving.min_thickness",CalvingMinthicknessEnum));
118  break;
119  case CalvingMinthicknessEnum:
120  parameters->AddObject(iomodel->CopyConstantObject("md.calving.min_thickness",CalvingMinthicknessEnum));
121  break;
122  case CalvingHabEnum:
123  break;
124  case CalvingCrevasseDepthEnum:
125  parameters->AddObject(iomodel->CopyConstantObject("md.calving.crevasse_opening_stress",CalvingCrevasseDepthEnum));
126  break;
127  case CalvingDev2Enum:
128  parameters->AddObject(iomodel->CopyConstantObject("md.calving.height_above_floatation",CalvingHeightAboveFloatationEnum));
129  break;
130  default:
131  _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
132  }
133 
134  /*Get frontal melt parameters*/
135  int melt_parameterization;
136  iomodel->FindConstant(&melt_parameterization,"md.frontalforcings.parameterization");
137  switch(melt_parameterization){
138  case FrontalForcingsDefaultEnum:
139  break;
140  case FrontalForcingsRignotEnum:
141  parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.numberofbasins",FrontalForcingsNumberofBasinsEnum));
142  break;
143  default:
144  _error_("Frontal forcings "<<EnumToStringx(melt_parameterization)<<" not supported yet");
145  }
146 }
147 /*}}}*/
148 
149 /*Finite element Analysis*/
150 void LevelsetAnalysis::Core(FemModel* femmodel){/*{{{*/
151 
152  /*parameters: */
153  int stabilization;
154  bool save_results;
155  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
156  femmodel->parameters->FindParam(&stabilization,LevelsetStabilizationEnum);
157 
158  /*activate formulation: */
159  femmodel->SetCurrentConfiguration(LevelsetAnalysisEnum);
160 
161  if(VerboseSolution()) _printf0_(" call computational core:\n");
162  if(stabilization==4){
163  solutionsequence_fct(femmodel);
164  }
165  else{
166  solutionsequence_linear(femmodel);
167  }
168 
169  if(save_results){
170  if(VerboseSolution()) _printf0_(" saving levelset results\n");
171  int outputs[1] = {MaskIceLevelsetEnum};
172  femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
173  }
174 }/*}}}*/
175 ElementVector* LevelsetAnalysis::CreateDVector(Element* element){/*{{{*/
176  /*Default, return NULL*/
177  return NULL;
178 }/*}}}*/
179 ElementMatrix* LevelsetAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
180  /* Jacobian required for the Newton solver */
181  _error_("not implemented yet");
182 }/*}}}*/
183 ElementMatrix* LevelsetAnalysis::CreateKMatrix(Element* element){/*{{{*/
184 
185  if(!element->IsOnBase()) return NULL;
186  Element* basalelement = element->SpawnBasalElement();
187 
188  /*Intermediaries */
189  int stabilization,dim, domaintype, calvinglaw;
190  int i,j,k,row, col;
191  IssmDouble kappa;
192  IssmDouble Jdet, dt, D_scalar;
193  IssmDouble h,hx,hy,hz;
194  IssmDouble vel,v[3],w[3],c[3],m[3],dlsf[3];
195  IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate, groundedice;
196  IssmDouble calvingmax, calvinghaf, heaviside, haf_eps;
197  IssmDouble* xyz_list = NULL;
198 
199  /*Get problem dimension and whether there is moving front or not*/
200  basalelement->FindParam(&domaintype,DomainTypeEnum);
201  basalelement->FindParam(&calvinglaw,CalvingLawEnum);
202  basalelement->FindParam(&stabilization,LevelsetStabilizationEnum);
203  switch(domaintype){
204  case Domain2DverticalEnum: dim = 1; break;
205  case Domain2DhorizontalEnum: dim = 2; break;
206  case Domain3DEnum: dim = 2; break;
207  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
208  }
209 
210  /*Calving threshold*/
211 
212  /*Fetch number of nodes and dof for this finite element*/
213  int numnodes = basalelement->GetNumberOfNodes();
214 
215  /*Initialize Element vector and other vectors*/
216  ElementMatrix* Ke = basalelement->NewElementMatrix();
217  IssmDouble* basis = xNew<IssmDouble>(numnodes);
218  IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);
219  IssmDouble* Bprime = NULL;
220  if(stabilization==2){
221  Bprime = xNew<IssmDouble>(dim*numnodes);
222  }
223 
224  /*Retrieve all inputs and parameters*/
225  basalelement->GetVerticesCoordinates(&xyz_list);
226  basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
227  basalelement->FindParam(&calvingmax,CalvingMaxEnum);
228  Input2* vx_input = NULL;
229  Input2* vy_input = NULL;
230  Input2* calvingratex_input = NULL;
231  Input2* calvingratey_input = NULL;
232  Input2* lsf_slopex_input = NULL;
233  Input2* lsf_slopey_input = NULL;
234  Input2* calvingrate_input = NULL;
235  Input2* meltingrate_input = NULL;
236  Input2* gr_input = NULL;
237 
238  /*Load velocities*/
239  switch(domaintype){
240  case Domain2DverticalEnum:
241  vx_input=basalelement->GetInput2(VxEnum); _assert_(vx_input);
242  break;
243  case Domain2DhorizontalEnum:
244  vx_input=basalelement->GetInput2(VxEnum); _assert_(vx_input);
245  vy_input=basalelement->GetInput2(VyEnum); _assert_(vy_input);
246  gr_input=basalelement->GetInput2(MaskOceanLevelsetEnum); _assert_(gr_input);
247  break;
248  case Domain3DEnum:
249  vx_input=basalelement->GetInput2(VxAverageEnum); _assert_(vx_input);
250  vy_input=basalelement->GetInput2(VyAverageEnum); _assert_(vy_input);
251  gr_input=basalelement->GetInput2(MaskOceanLevelsetEnum); _assert_(gr_input);
252  break;
253  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
254  }
255 
256  /*Load calving inputs*/
257  switch(calvinglaw){
258  case DefaultCalvingEnum:
259  case CalvingVonmisesEnum:
260  lsf_slopex_input = basalelement->GetInput2(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
261  if(dim==2) lsf_slopey_input = basalelement->GetInput2(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
262  calvingrate_input = basalelement->GetInput2(CalvingCalvingrateEnum); _assert_(calvingrate_input);
263  meltingrate_input = basalelement->GetInput2(CalvingMeltingrateEnum); _assert_(meltingrate_input);
264  break;
265  case CalvingLevermannEnum:
266  switch(domaintype){
267  case Domain2DverticalEnum:
268  calvingratex_input=basalelement->GetInput2(CalvingratexEnum); _assert_(calvingratex_input);
269  break;
270  case Domain2DhorizontalEnum:
271  calvingratex_input=basalelement->GetInput2(CalvingratexEnum); _assert_(calvingratex_input);
272  calvingratey_input=basalelement->GetInput2(CalvingrateyEnum); _assert_(calvingratey_input);
273  break;
274  case Domain3DEnum:
275  calvingratex_input=basalelement->GetInput2(CalvingratexAverageEnum); _assert_(calvingratex_input);
276  calvingratey_input=basalelement->GetInput2(CalvingrateyAverageEnum); _assert_(calvingratey_input);
277  break;
278  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
279  }
280  meltingrate_input = basalelement->GetInput2(CalvingMeltingrateEnum); _assert_(meltingrate_input);
281  break;
282  case CalvingMinthicknessEnum:
283  lsf_slopex_input = basalelement->GetInput2(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
284  if(dim==2) lsf_slopey_input = basalelement->GetInput2(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
285  meltingrate_input = basalelement->GetInput2(CalvingMeltingrateEnum); _assert_(meltingrate_input);
286  break;
287  case CalvingHabEnum:
288  lsf_slopex_input = basalelement->GetInput2(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
289  if(dim==2) lsf_slopey_input = basalelement->GetInput2(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
290  meltingrate_input = basalelement->GetInput2(CalvingMeltingrateEnum); _assert_(meltingrate_input);
291  break;
292  case CalvingCrevasseDepthEnum:
293  lsf_slopex_input = basalelement->GetInput2(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
294  if(dim==2) lsf_slopey_input = basalelement->GetInput2(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
295  meltingrate_input = basalelement->GetInput2(CalvingMeltingrateEnum); _assert_(meltingrate_input);
296  break;
297  case CalvingDev2Enum:
298  basalelement->FindParam(&calvinghaf,CalvingHeightAboveFloatationEnum);
299  lsf_slopex_input = basalelement->GetInput2(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
300  if(dim==2) lsf_slopey_input = basalelement->GetInput2(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
301  calvingrate_input = basalelement->GetInput2(CalvingCalvingrateEnum); _assert_(calvingrate_input);
302  meltingrate_input = basalelement->GetInput2(CalvingMeltingrateEnum); _assert_(meltingrate_input);
303  break;
304  default:
305  _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
306  }
307 
308  /* Start looping on the number of gaussian points: */
309  Gauss* gauss=basalelement->NewGauss(2);
310  for(int ig=gauss->begin();ig<gauss->end();ig++){
311  gauss->GaussPoint(ig);
312 
313  basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
314  basalelement->NodalFunctions(basis,gauss);
315  basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
316  D_scalar=gauss->weight*Jdet;
317 
318  /* Transient */
319  if(dt!=0.){
320  for(i=0;i<numnodes;i++){
321  for(j=0;j<numnodes;j++){
322  Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
323  }
324  }
325  D_scalar=D_scalar*dt;
326  }
327 
328  /* Advection */
329  vx_input->GetInputValue(&v[0],gauss);
330  vy_input->GetInputValue(&v[1],gauss);
331  gr_input->GetInputValue(&groundedice,gauss);
332 
333  /*Get calving speed*/
334  switch(calvinglaw){
335  case DefaultCalvingEnum:
336  case CalvingVonmisesEnum:
337  lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
338  if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
339  calvingrate_input->GetInputValue(&calvingrate,gauss);
340  meltingrate_input->GetInputValue(&meltingrate,gauss);
341 
342  /*Limit calving rate to c <= v + 3 km/yr */
343  vel=sqrt(v[0]*v[0] + v[1]*v[1]);
344  if(calvingrate>calvingmax+vel) calvingrate = vel+calvingmax;
345  if(groundedice<0) meltingrate = 0.;
346 
347  norm_dlsf=0.;
348  for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
349  norm_dlsf=sqrt(norm_dlsf);
350 
351  if(norm_dlsf>1.e-10)
352  for(i=0;i<dim;i++){
353  c[i]=calvingrate*dlsf[i]/norm_dlsf; m[i]=meltingrate*dlsf[i]/norm_dlsf;
354  }
355  else
356  for(i=0;i<dim;i++){
357  c[i]=0.; m[i]=0.;
358  }
359  break;
360 
361  case CalvingLevermannEnum:
362  calvingratex_input->GetInputValue(&c[0],gauss);
363  if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
364  meltingrate_input->GetInputValue(&meltingrate,gauss);
365  norm_calving=0.;
366  for(i=0;i<dim;i++) norm_calving+=pow(c[i],2);
367  norm_calving=sqrt(norm_calving)+1.e-14;
368  for(i=0;i<dim;i++) m[i]=meltingrate*c[i]/norm_calving;
369  break;
370 
371  case CalvingMinthicknessEnum:
372  lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
373  if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
374  meltingrate_input->GetInputValue(&meltingrate,gauss);
375 
376  norm_dlsf=0.;
377  for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
378  norm_dlsf=sqrt(norm_dlsf);
379 
380  if(norm_dlsf>1.e-10)
381  for(i=0;i<dim;i++){
382  c[i]=0.;
383  m[i]=meltingrate*dlsf[i]/norm_dlsf;
384  }
385  else
386  for(i=0;i<dim;i++){
387  c[i]=0.;
388  m[i]=0.;
389  }
390  break;
391 
392  case CalvingHabEnum:
393  lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
394  if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
395  meltingrate_input->GetInputValue(&meltingrate,gauss);
396 
397  norm_dlsf=0.;
398  for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
399  norm_dlsf=sqrt(norm_dlsf);
400 
401  if(norm_dlsf>1.e-10)
402  for(i=0;i<dim;i++){
403  c[i]=0.;
404  m[i]=meltingrate*dlsf[i]/norm_dlsf;
405  }
406  else
407  for(i=0;i<dim;i++){
408  c[i]=0.;
409  m[i]=0.;
410  }
411  break;
412 
413  case CalvingCrevasseDepthEnum:
414  lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
415  if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
416  meltingrate_input->GetInputValue(&meltingrate,gauss);
417 
418  if(groundedice<0) meltingrate = 0.;
419 
420  norm_dlsf=0.;
421  for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
422  norm_dlsf=sqrt(norm_dlsf);
423 
424  if(norm_dlsf>1.e-10)
425  for(i=0;i<dim;i++){
426  c[i]=0.;
427  m[i]=meltingrate*dlsf[i]/norm_dlsf;
428  }
429  else
430  for(i=0;i<dim;i++){
431  c[i]=0.;
432  m[i]=0.;
433  }
434  break;
435 
436  case CalvingDev2Enum:
437  {
438  lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
439  if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
440  calvingrate_input->GetInputValue(&calvingrate,gauss);
441  meltingrate_input->GetInputValue(&meltingrate,gauss);
442  gr_input->GetInputValue(&groundedice,gauss);
443 
444  //idea: no retreat on ice above critical calving height "calvinghaf" . Limit using regularized Heaviside function.
445  vel=sqrt(v[0]*v[0] + v[1]*v[1]);
446  haf_eps=10.;
447  if(groundedice-calvinghaf<=-haf_eps){
448  // ice floats freely below calvinghaf: calve freely
449  // undercutting has no effect:
450  meltingrate=0.;
451  }
452  else if(groundedice-calvinghaf>=haf_eps){
453  // ice is well above calvinghaf -> no calving back, i.e. limit calving rate to ice velocity
454  calvingrate=min(calvingrate,vel);
455  // ice is almost grounded: frontal undercutting has maximum effect (do nothing).
456  }
457  else{ // ice is close to calvinghaf: smooth transition between limitation and free calving.
458  //heaviside: 0 for floating, 1 for grounded
459  heaviside=(groundedice-calvinghaf+haf_eps)/(2.*haf_eps) + sin(PI*(groundedice-calvinghaf)/haf_eps)/(2.*PI);
460  calvingrate=heaviside*(min(calvingrate,vel)-calvingrate)+calvingrate;
461  meltingrate=heaviside*meltingrate+0.;
462  }
463 
464  norm_dlsf=0.;
465  for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
466  norm_dlsf=sqrt(norm_dlsf);
467 
468  if(norm_dlsf>1.e-10)
469  for(i=0;i<dim;i++){
470  c[i]=calvingrate*dlsf[i]/norm_dlsf;
471  m[i]=meltingrate*dlsf[i]/norm_dlsf;
472  }
473  else
474  for(i=0;i<dim;i++){
475  c[i]=0.;
476  m[i]=0.;
477  }
478  break;
479  }
480 
481  default:
482  _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
483  }
484 
485  /*Levelset speed is ice velocity - calving rate*/
486  for(i=0;i<dim;i++) w[i]=v[i]-c[i]-m[i];
487 
488  /*Compute D*/
489  for(i=0;i<numnodes;i++){
490  for(j=0;j<numnodes;j++){
491  for(k=0;k<dim;k++){
492  Ke->values[i*numnodes+j] += D_scalar*w[k]*dbasis[k*numnodes+j]*basis[i];
493  }
494  }
495  }
496 
497  /* Stabilization */
498  vel=0.;
499  for(i=0;i<dim;i++) vel+=w[i]*w[i];
500  vel=sqrt(vel)+1.e-14;
501  switch(stabilization){
502  case 0:
503  // no stabilization, do nothing
504  break;
505  case 1:
506  /* Artificial Diffusion */
507  basalelement->ElementSizes(&hx,&hy,&hz);
508  h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
509  kappa=h*vel/2.;
510  for(i=0;i<numnodes;i++){
511  for(j=0;j<numnodes;j++){
512  for(k=0;k<dim;k++){
513  Ke->values[i*numnodes+j] += D_scalar*kappa*dbasis[k*numnodes+j]*dbasis[k*numnodes+i];
514  }
515  }
516  }
517  break;
518  case 2:
519  {
520  /* Streamline Upwinding */
521  basalelement->ElementSizes(&hx,&hy,&hz);
522  h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
523  for(int i=0;i<numnodes;i++){
524  for(int j=0;j<numnodes;j++){
525  Ke->values[i*numnodes+j] += D_scalar*h/(2.*vel)*(
526  dbasis[0*numnodes+i] *(w[0]*w[0]*dbasis[0*numnodes+j] + w[0]*w[1]*dbasis[1*numnodes+j]) +
527  dbasis[1*numnodes+i] *(w[1]*w[0]*dbasis[0*numnodes+j] + w[1]*w[1]*dbasis[1*numnodes+j])
528  );
529  }
530  }
531  }
532  break;
533  default:
534  _error_("unknown type of stabilization in LevelsetAnalysis.cpp");
535  }
536  }
537 
538  /*Clean up and return*/
539  xDelete<IssmDouble>(xyz_list);
540  xDelete<IssmDouble>(basis);
541  xDelete<IssmDouble>(dbasis);
542  xDelete<IssmDouble>(Bprime);
543  delete gauss;
544  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
545  return Ke;
546 }/*}}}*/
547 ElementVector* LevelsetAnalysis::CreatePVector(Element* element){/*{{{*/
548 
549  if(!element->IsOnBase()) return NULL;
550  Element* basalelement = element->SpawnBasalElement();
551 
552  /*Intermediaries */
553  int i, ig, domaintype;
554  IssmDouble Jdet,dt;
555  IssmDouble lsf;
556  IssmDouble* xyz_list = NULL;
557 
558  /*Fetch number of nodes and dof for this finite element*/
559  int numnodes = basalelement->GetNumberOfNodes();
560 
561  /*Initialize Element vector*/
562  ElementVector* pe = basalelement->NewElementVector();
563  basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
564 
565  if(dt!=0.){
566  /*Initialize basis vector*/
567  IssmDouble* basis = xNew<IssmDouble>(numnodes);
568 
569  /*Retrieve all inputs and parameters*/
570  basalelement->GetVerticesCoordinates(&xyz_list);
571  Input2* levelset_input = basalelement->GetInput2(MaskIceLevelsetEnum); _assert_(levelset_input);
572 
573  /* Start looping on the number of gaussian points: */
574  Gauss* gauss=basalelement->NewGauss(2);
575  for(ig=gauss->begin();ig<gauss->end();ig++){
576  gauss->GaussPoint(ig);
577 
578  basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
579  basalelement->NodalFunctions(basis,gauss);
580 
581  /* old function value */
582  levelset_input->GetInputValue(&lsf,gauss);
583  for(i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*lsf*basis[i];
584  }
585 
586  /*Clean up and return*/
587  xDelete<IssmDouble>(xyz_list);
588  xDelete<IssmDouble>(basis);
589  basalelement->FindParam(&domaintype,DomainTypeEnum);
590  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
591  delete gauss;
592  }
593 
594  return pe;
595 }/*}}}*/
596 IssmDouble LevelsetAnalysis::GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1){/*{{{*/
597  // returns distance d of point q to straight going through points s0, s1
598  // d=|a x b|/|b|
599  // with a=q-s0, b=s1-s0
600 
601  /* Intermediaries */
602  const int dim=2;
603  int i;
604  IssmDouble a[dim], b[dim];
605  IssmDouble norm_b;
606 
607  for(i=0;i<dim;i++){
608  a[i]=q[i]-s0[i];
609  b[i]=s1[i]-s0[i];
610  }
611 
612  norm_b=0.;
613  for(i=0;i<dim;i++)
614  norm_b+=b[i]*b[i];
615  norm_b=sqrt(norm_b);
616  _assert_(norm_b>0.);
617 
618  return fabs(a[0]*b[1]-a[1]*b[0])/norm_b;
619 }/*}}}*/
620 void LevelsetAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
621  _error_("not implemented yet");
622 }/*}}}*/
623 void LevelsetAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
624  _error_("Not implemented yet");
625 }/*}}}*/
626 void LevelsetAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
627 
628  int domaintype;
629  element->FindParam(&domaintype,DomainTypeEnum);
630  switch(domaintype){
631  case Domain2DhorizontalEnum:
632  element->InputUpdateFromSolutionOneDof(solution,MaskIceLevelsetEnum);
633  break;
634  case Domain3DEnum:
635  element->InputUpdateFromSolutionOneDofCollapsed(solution,MaskIceLevelsetEnum);
636  break;
637  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
638  }
639 }/*}}}*/
640 void LevelsetAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
641 
642  /*Intermediaries*/
643  int calvinglaw;
644  IssmDouble min_thickness,thickness,hab_fraction;
645  IssmDouble crevassedepth,surface_crevasse,surface,critical_fraction;
646  IssmDouble rho_ice,rho_water;
647  IssmDouble bed,water_depth;
648  IssmDouble levelset,sealevel;
649 
650  femmodel->parameters->FindParam(&calvinglaw,CalvingLawEnum);
651 
652  if(calvinglaw==CalvingMinthicknessEnum || calvinglaw==CalvingVonmisesEnum){
653 
654  /*Get minimum thickness threshold*/
655  femmodel->parameters->FindParam(&min_thickness,CalvingMinthicknessEnum);
656 
657  /*Loop over all elements of this partition*/
658  for(int i=0;i<femmodel->elements->Size();i++){
659  Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
660 
661  int numnodes = element->GetNumberOfNodes();
662  Gauss* gauss = element->NewGauss();
663  Input2* H_input = element->GetInput2(ThicknessEnum); _assert_(H_input);
664  Input2* b_input = element->GetInput2(BedEnum); _assert_(b_input);
665  Input2* sl_input = element->GetInput2(SealevelEnum); _assert_(sl_input);
666 
667  /*Potentially constrain nodes of this element*/
668  for(int in=0;in<numnodes;in++){
669  gauss->GaussNode(element->GetElementType(),in);
670  Node* node=element->GetNode(in);
671  if(!node->IsActive()) continue;
672 
673  H_input->GetInputValue(&thickness,gauss);
674  b_input->GetInputValue(&bed,gauss);
675  sl_input->GetInputValue(&sealevel,gauss);
676  if(thickness<min_thickness && bed<sealevel){
677  node->ApplyConstraint(0,+1.);
678  }
679  else {
680  /* no ice, set no spc */
681  node->DofInFSet(0);
682  }
683  }
684  delete gauss;
685  }
686  }
687 
688  if(calvinglaw==CalvingHabEnum){
689 
690  /*Get the fraction of the flotation thickness at the terminus*/
691  InputDuplicatex(femmodel,MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);
692  femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);
693 
694  /*Loop over all elements of this partition*/
695  for(int i=0;i<femmodel->elements->Size();i++){
696  Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
697 
698  rho_ice = element->FindParam(MaterialsRhoIceEnum);
699  rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
700 
701  int numnodes = element->GetNumberOfNodes();
702  Gauss* gauss = element->NewGauss();
703  Input2* H_input = element->GetInput2(ThicknessEnum); _assert_(H_input);
704  Input2* bed_input = element->GetInput2(BedEnum); _assert_(bed_input);
705  Input2* hab_fraction_input = element->GetInput2(CalvingHabFractionEnum); _assert_(hab_fraction_input);
706  Input2* ls_input = element->GetInput2(DistanceToCalvingfrontEnum); _assert_(ls_input);
707 
708  /*Potentially constrain nodes of this element*/
709  for(int in=0;in<numnodes;in++){
710  gauss->GaussNode(element->GetElementType(),in);
711  Node* node=element->GetNode(in);
712  if(!node->IsActive()) continue;
713 
714  H_input->GetInputValue(&thickness,gauss);
715  bed_input->GetInputValue(&water_depth,gauss);
716  ls_input->GetInputValue(&levelset,gauss);
717  hab_fraction_input->GetInputValue(&hab_fraction,gauss);
718 
719  if(thickness<((rho_water/rho_ice)*(1+hab_fraction)*-water_depth) && levelset>-300. && levelset<0.){
720  node->ApplyConstraint(0,+1.);
721  }
722  else {
723  /* no ice, set no spc */
724  node->DofInFSet(0);
725  }
726  }
727  delete gauss;
728  }
729  }
730 
731  if(calvinglaw==CalvingCrevasseDepthEnum){
732 
733  int nflipped,local_nflipped;
734  Vector<IssmDouble>* vec_constraint_nodes = NULL;
735  IssmDouble* constraint_nodes = NULL;
736 
737  /*Get the DistanceToCalvingfront*/
738  InputDuplicatex(femmodel,MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);
739  femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);
740 
741  /*Vector of size number of nodes*/
742  vec_constraint_nodes=new Vector<IssmDouble>(femmodel->nodes->NumberOfNodes());
743 
744  for(int i=0;i<femmodel->elements->Size();i++){
745  Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
746  int numnodes = element->GetNumberOfNodes();
747  Gauss* gauss = element->NewGauss();
748  Input2* crevassedepth_input = element->GetInput2(CrevasseDepthEnum); _assert_(crevassedepth_input);
749  Input2* bed_input = element->GetInput2(BedEnum); _assert_(bed_input);
750  Input2* surface_crevasse_input = element->GetInput2(SurfaceCrevasseEnum); _assert_(surface_crevasse_input);
751  Input2* thickness_input = element->GetInput2(ThicknessEnum); _assert_(thickness_input);
752  Input2* surface_input = element->GetInput2(SurfaceEnum); _assert_(surface_input);
753 
754  /*First, look at ice front and figure out if any of the nodes will be calved*/
755  if(element->IsIcefront()){
756  for(int in=0;in<numnodes;in++){
757  gauss->GaussNode(element->GetElementType(),in);
758  Node* node=element->GetNode(in);
759  if(!node->IsActive()) continue;
760 
761  crevassedepth_input->GetInputValue(&crevassedepth,gauss);
762  bed_input->GetInputValue(&bed,gauss);
763  surface_crevasse_input->GetInputValue(&surface_crevasse,gauss);
764  thickness_input->GetInputValue(&thickness,gauss);
765  surface_input->GetInputValue(&surface,gauss);
766 
767  if((surface_crevasse-surface>0. || crevassedepth-thickness>0.) && bed<0.){
768  vec_constraint_nodes->SetValue(node->Sid(),1.0,INS_VAL);
769  }
770  }
771  }
772  delete gauss;
773  }
774 
775  /*Assemble vector and serialize: */
776  vec_constraint_nodes->Assemble();
777  constraint_nodes=vec_constraint_nodes->ToMPISerial();
778 
779  nflipped=1;
780  while(nflipped){
781  local_nflipped=0;
782  for(int i=0;i<femmodel->elements->Size();i++){
783  Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
784  int numnodes = element->GetNumberOfNodes();
785  Gauss* gauss = element->NewGauss();
786  Input2* levelset_input = element->GetInput2(DistanceToCalvingfrontEnum); _assert_(levelset_input);
787  Input2* crevassedepth_input = element->GetInput2(CrevasseDepthEnum); _assert_(crevassedepth_input);
788  Input2* bed_input = element->GetInput2(BedEnum); _assert_(bed_input);
789  Input2* surface_crevasse_input = element->GetInput2(SurfaceCrevasseEnum); _assert_(surface_crevasse_input);
790  Input2* thickness_input = element->GetInput2(ThicknessEnum); _assert_(thickness_input);
791  Input2* surface_input = element->GetInput2(SurfaceEnum); _assert_(surface_input);
792 
793  /*Is this element connected to a node that should be calved*/
794  bool isconnected = false;
795  for(int in=0;in<numnodes;in++){
796  Node* node=element->GetNode(in);
797  if(constraint_nodes[node->Sid()]==1.){
798  isconnected = true;
799  break;
800  }
801  }
802 
803  /*Check status if connected*/
804  if(isconnected){
805  for(int in=0;in<numnodes;in++){
806  gauss->GaussNode(element->GetElementType(),in);
807  Node* node=element->GetNode(in);
808  levelset_input->GetInputValue(&levelset,gauss);
809  crevassedepth_input->GetInputValue(&crevassedepth,gauss);
810  bed_input->GetInputValue(&bed,gauss);
811  surface_crevasse_input->GetInputValue(&surface_crevasse,gauss);
812  thickness_input->GetInputValue(&thickness,gauss);
813  surface_input->GetInputValue(&surface,gauss);
814 
815  if((surface_crevasse-surface>0. || crevassedepth-thickness>0.) && bed<0. && levelset>-300. && levelset<0. && constraint_nodes[node->Sid()]==0.){
816  local_nflipped++;
817  vec_constraint_nodes->SetValue(node->Sid(),1.0,INS_VAL);
818  }
819  }
820  }
821  }
822 
823  /*Count how many new nodes were found*/
824  ISSM_MPI_Allreduce(&local_nflipped,&nflipped,1,ISSM_MPI_INT,ISSM_MPI_SUM,IssmComm::GetComm());
825  //_printf0_("Found "<<nflipped<<" to flip\n");
826 
827  /*Assemble and serialize flag vector*/
828  vec_constraint_nodes->Assemble();
829  xDelete<IssmDouble>(constraint_nodes);
830  constraint_nodes=vec_constraint_nodes->ToMPISerial();
831  }
832  /*Free ressources:*/
833  delete vec_constraint_nodes;
834 
835  /*Contrain the nodes that will be calved*/
836  for(int i=0;i<femmodel->elements->Size();i++){
837  Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
838  int numnodes = element->GetNumberOfNodes();
839  Gauss* gauss = element->NewGauss();
840  /*Potentially constrain nodes of this element*/
841  for(int in=0;in<numnodes;in++){
842  gauss->GaussNode(element->GetElementType(),in);
843  Node* node=element->GetNode(in);
844  if(!node->IsActive()) continue;
845 
846  if(constraint_nodes[node->Sid()]>0.){
847  node->ApplyConstraint(0,+1.);
848  }
849  else {
850  /* no ice, set no spc */
851  node->DofInFSet(0);
852  }
853  }
854  delete gauss;
855  }
856  xDelete<IssmDouble>(constraint_nodes);
857  }
858 
859  /*Default, do nothing*/
860  return;
861 }/*}}}*/
DataSet::Size
int Size()
Definition: DataSet.cpp:399
Node::IsActive
bool IsActive(void)
Definition: Node.cpp:795
SaveResultsEnum
@ SaveResultsEnum
Definition: EnumDefinitions.h:302
Element::GetElementType
virtual int GetElementType(void)=0
CalvingStressThresholdGroundediceEnum
@ CalvingStressThresholdGroundediceEnum
Definition: EnumDefinitions.h:506
CalvingVonmisesEnum
@ CalvingVonmisesEnum
Definition: EnumDefinitions.h:1004
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
CalvingHabEnum
@ CalvingHabEnum
Definition: EnumDefinitions.h:1002
Element::IsOnBase
bool IsOnBase()
Definition: Element.cpp:1984
Nodes
Declaration of Nodes class.
Definition: Nodes.h:19
LevelsetAnalysis::GetDistanceToStraight
IssmDouble GetDistanceToStraight(IssmDouble *q, IssmDouble *s0, IssmDouble *s1)
Definition: LevelsetAnalysis.cpp:596
LevelsetAnalysis::CreateJacobianMatrix
ElementMatrix * CreateJacobianMatrix(Element *element)
Definition: LevelsetAnalysis.cpp:179
Element::GetNumberOfNodes
virtual int GetNumberOfNodes(void)=0
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
ISSM_MPI_Allreduce
int ISSM_MPI_Allreduce(void *sendbuf, void *recvbuf, int count, ISSM_MPI_Datatype datatype, ISSM_MPI_Op op, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:116
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
LevelsetReinitFrequencyEnum
@ LevelsetReinitFrequencyEnum
Definition: EnumDefinitions.h:228
MaskOceanLevelsetEnum
@ MaskOceanLevelsetEnum
Definition: EnumDefinitions.h:640
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
LevelsetfunctionSlopeXEnum
@ LevelsetfunctionSlopeXEnum
Definition: EnumDefinitions.h:635
MaskIceLevelsetEnum
@ MaskIceLevelsetEnum
Definition: EnumDefinitions.h:641
ISSM_MPI_SUM
#define ISSM_MPI_SUM
Definition: issmmpi.h:134
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
IssmComm::GetComm
static ISSM_MPI_Comm GetComm(void)
Definition: IssmComm.cpp:30
CalvinglevermannCoeffEnum
@ CalvinglevermannCoeffEnum
Definition: EnumDefinitions.h:507
TimesteppingTimeStepEnum
@ TimesteppingTimeStepEnum
Definition: EnumDefinitions.h:433
LevelsetfunctionSlopeYEnum
@ LevelsetfunctionSlopeYEnum
Definition: EnumDefinitions.h:636
Constraints
Declaration of Constraints class.
Definition: Constraints.h:13
BedEnum
@ BedEnum
Definition: EnumDefinitions.h:499
Element::SpawnBasalElement
virtual Element * SpawnBasalElement(void)=0
Elements
Declaration of Elements class.
Definition: Elements.h:17
FemModel::results
Results * results
Definition: FemModel.h:48
CalvingCrevasseDepthEnum
@ CalvingCrevasseDepthEnum
Definition: EnumDefinitions.h:96
MaterialsRhoIceEnum
@ MaterialsRhoIceEnum
Definition: EnumDefinitions.h:264
LevelsetAnalysis::CreateLoads
void CreateLoads(Loads *loads, IoModel *iomodel)
Definition: LevelsetAnalysis.cpp:19
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
CalvingMaxEnum
@ CalvingMaxEnum
Definition: EnumDefinitions.h:99
FrontalForcingsThermalForcingEnum
@ FrontalForcingsThermalForcingEnum
Definition: EnumDefinitions.h:587
Gauss::GaussNode
virtual void GaussNode(int finitelement, int iv)=0
LevelsetAnalysis::CreateDVector
ElementVector * CreateDVector(Element *element)
Definition: LevelsetAnalysis.cpp:175
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
solutionsequence_linear
void solutionsequence_linear(FemModel *femmodel)
Definition: solutionsequence_linear.cpp:10
Parameters::AddObject
void AddObject(Param *newparam)
Definition: Parameters.cpp:67
CalvingLawEnum
@ CalvingLawEnum
Definition: EnumDefinitions.h:98
CalvingMeltingrateEnum
@ CalvingMeltingrateEnum
Definition: EnumDefinitions.h:504
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
Element::DeleteMaterials
void DeleteMaterials(void)
Definition: Element.cpp:429
Element
Definition: Element.h:41
Element::NodalFunctions
virtual void NodalFunctions(IssmDouble *basis, Gauss *gauss)=0
LevelsetAnalysis.h
: header file for generic external result object
LevelsetStabilizationEnum
@ LevelsetStabilizationEnum
Definition: EnumDefinitions.h:229
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
Element::ElementSizes
virtual void ElementSizes(IssmDouble *phx, IssmDouble *phy, IssmDouble *phz)=0
CalvingrateyEnum
@ CalvingrateyEnum
Definition: EnumDefinitions.h:511
Nodes::NumberOfNodes
int NumberOfNodes(void)
Definition: Nodes.cpp:354
Element::NewElementVector
ElementVector * NewElementVector(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2505
ISSM_MPI_INT
#define ISSM_MPI_INT
Definition: issmmpi.h:127
IoModel::DeleteData
void DeleteData(int num,...)
Definition: IoModel.cpp:500
IoModel::numberofelements
int numberofelements
Definition: IoModel.h:96
IoModel::CopyConstantObject
Param * CopyConstantObject(const char *constant_name, int param_enum)
Definition: IoModel.cpp:418
CalvingratexAverageEnum
@ CalvingratexAverageEnum
Definition: EnumDefinitions.h:508
DomainTypeEnum
@ DomainTypeEnum
Definition: EnumDefinitions.h:124
VxAverageEnum
@ VxAverageEnum
Definition: EnumDefinitions.h:845
FemModel::nodes
Nodes * nodes
Definition: FemModel.h:56
FrontalForcingsDefaultEnum
@ FrontalForcingsDefaultEnum
Definition: EnumDefinitions.h:1073
Element::NewGauss
virtual Gauss * NewGauss(void)=0
FrontalForcingsNumberofBasinsEnum
@ FrontalForcingsNumberofBasinsEnum
Definition: EnumDefinitions.h:153
FrontalForcingsRignotEnum
@ FrontalForcingsRignotEnum
Definition: EnumDefinitions.h:1074
SurfaceCrevasseEnum
@ SurfaceCrevasseEnum
Definition: EnumDefinitions.h:822
Element::InputUpdateFromSolutionOneDof
virtual void InputUpdateFromSolutionOneDof(IssmDouble *solution, int inputenum)=0
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
LevelsetAnalysisEnum
@ LevelsetAnalysisEnum
Definition: EnumDefinitions.h:1141
LevelsetAnalysis::UpdateParameters
void UpdateParameters(Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
Definition: LevelsetAnalysis.cpp:103
DistanceToCalvingfrontEnum
@ DistanceToCalvingfrontEnum
Definition: EnumDefinitions.h:532
IoModel::FindConstant
void FindConstant(bool *pvalue, const char *constant_name)
Definition: IoModel.cpp:2362
CalvingHabFractionEnum
@ CalvingHabFractionEnum
Definition: EnumDefinitions.h:503
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
Element::GetNode
Node * GetNode(int nodeindex)
Definition: Element.cpp:1207
LevelsetAnalysis::UpdateConstraints
void UpdateConstraints(FemModel *femmodel)
Definition: LevelsetAnalysis.cpp:640
INS_VAL
@ INS_VAL
Definition: toolkitsenums.h:14
InputDuplicatex
void InputDuplicatex(FemModel *femmodel, int original_enum, int new_enum)
Definition: InputDuplicatex.cpp:10
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
LevelsetAnalysis::GetSolutionFromInputs
void GetSolutionFromInputs(Vector< IssmDouble > *solution, Element *element)
Definition: LevelsetAnalysis.cpp:620
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
MaterialsRhoSeawaterEnum
@ MaterialsRhoSeawaterEnum
Definition: EnumDefinitions.h:265
PI
const double PI
Definition: constants.h:11
LevelsetAnalysis::Core
void Core(FemModel *femmodel)
Definition: LevelsetAnalysis.cpp:150
FemModel::elements
Elements * elements
Definition: FemModel.h:44
CalvingHeightAboveFloatationEnum
@ CalvingHeightAboveFloatationEnum
Definition: EnumDefinitions.h:97
Vector::Assemble
void Assemble(void)
Definition: Vector.h:142
Input2
Definition: Input2.h:18
FemModel
Definition: FemModel.h:31
CalvingDev2Enum
@ CalvingDev2Enum
Definition: EnumDefinitions.h:1001
LevelsetKillIcebergsEnum
@ LevelsetKillIcebergsEnum
Definition: EnumDefinitions.h:227
SealevelEnum
@ SealevelEnum
Definition: EnumDefinitions.h:675
CalvingLevermannEnum
@ CalvingLevermannEnum
Definition: EnumDefinitions.h:1003
Loads
Declaration of Loads class.
Definition: Loads.h:16
Node
Definition: Node.h:23
Node::Sid
int Sid(void)
Definition: Node.cpp:622
LevelsetAnalysis::DofsPerNode
int DofsPerNode(int **doflist, int domaintype, int approximation)
Definition: LevelsetAnalysis.cpp:30
CalvingStressThresholdFloatingiceEnum
@ CalvingStressThresholdFloatingiceEnum
Definition: EnumDefinitions.h:505
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Node::DofInFSet
void DofInFSet(int dof)
Definition: Node.cpp:694
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
LevelsetAnalysis::CreateConstraints
void CreateConstraints(Constraints *constraints, IoModel *iomodel)
Definition: LevelsetAnalysis.cpp:13
LevelsetAnalysis::InputUpdateFromSolution
void InputUpdateFromSolution(IssmDouble *solution, Element *element)
Definition: LevelsetAnalysis.cpp:626
min
IssmDouble min(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:14
LevelsetAnalysis::CreateKMatrix
ElementMatrix * CreateKMatrix(Element *element)
Definition: LevelsetAnalysis.cpp:183
CalvingCalvingrateEnum
@ CalvingCalvingrateEnum
Definition: EnumDefinitions.h:502
Element::JacobianDeterminant
virtual void JacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)=0
Gauss::GaussPoint
virtual void GaussPoint(int ig)=0
LevelsetAnalysis::GradientJ
void GradientJ(Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
Definition: LevelsetAnalysis.cpp:623
Element::Update
virtual void Update(Inputs2 *inputs2, int index, IoModel *iomodel, int analysis_counter, int analysis_type, int finite_element)=0
FemModel::RequestedOutputsx
void RequestedOutputsx(Results **presults, char **requested_outputs, int numoutputs, bool save_results=true)
Definition: FemModel.cpp:2267
FemModel::SetCurrentConfiguration
void SetCurrentConfiguration(int configuration_type)
Definition: FemModel.cpp:634
VyAverageEnum
@ VyAverageEnum
Definition: EnumDefinitions.h:849
LevelsetAnalysis::UpdateElements
void UpdateElements(Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
Definition: LevelsetAnalysis.cpp:34
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
ThicknessEnum
@ ThicknessEnum
Definition: EnumDefinitions.h:840
WaterheightEnum
@ WaterheightEnum
Definition: EnumDefinitions.h:863
LevelsetAnalysis::CreateNodes
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
Definition: LevelsetAnalysis.cpp:22
IoModel::FetchDataToInput
void FetchDataToInput(Inputs2 *inputs2, Elements *elements, const char *vector_name, int input_enum)
Definition: IoModel.cpp:1651
FrontalForcingsBasinIdEnum
@ FrontalForcingsBasinIdEnum
Definition: EnumDefinitions.h:585
IoModelToConstraintsx
void IoModelToConstraintsx(Constraints *constraints, IoModel *iomodel, const char *spc_name, int analysis_type, int finite_element, int dof)
Definition: IoModelToConstraintsx.cpp:10
Node::ApplyConstraint
void ApplyConstraint(int dof, IssmDouble value)
Definition: Node.cpp:646
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
IoModel
Definition: IoModel.h:48
Vector::ToMPISerial
doubletype * ToMPISerial(void)
Definition: Vector.h:277
solutionsequence_fct
void solutionsequence_fct(FemModel *femmodel)
Definition: solutionsequence_fct.cpp:388
LevelsetAnalysis::CreatePVector
ElementVector * CreatePVector(Element *element)
Definition: LevelsetAnalysis.cpp:547
IoModel::domaintype
int domaintype
Definition: IoModel.h:78
FemModel::DistanceToFieldValue
void DistanceToFieldValue(int fieldenum, IssmDouble fieldvalue, int distanceenum)
Definition: FemModel.cpp:1103
CalvingMinthicknessEnum
@ CalvingMinthicknessEnum
Definition: EnumDefinitions.h:100
Domain2DverticalEnum
@ Domain2DverticalEnum
Definition: EnumDefinitions.h:535
ElementMatrix
Definition: ElementMatrix.h:19
FrontalForcingsSubglacialDischargeEnum
@ FrontalForcingsSubglacialDischargeEnum
Definition: EnumDefinitions.h:586
Vector< IssmDouble >
CalvingrateyAverageEnum
@ CalvingrateyAverageEnum
Definition: EnumDefinitions.h:510
CrevasseDepthEnum
@ CrevasseDepthEnum
Definition: EnumDefinitions.h:515
Element::InputUpdateFromSolutionOneDofCollapsed
virtual void InputUpdateFromSolutionOneDofCollapsed(IssmDouble *solution, int inputenum)=0
Element::IsIcefront
virtual bool IsIcefront(void)=0
Gauss
Definition: Gauss.h:8
CalvingratexEnum
@ CalvingratexEnum
Definition: EnumDefinitions.h:509
Vector::SetValue
void SetValue(int dof, doubletype value, InsMode mode)
Definition: Vector.h:163
Element::NodalFunctionsDerivatives
virtual void NodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
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
DefaultCalvingEnum
@ DefaultCalvingEnum
Definition: EnumDefinitions.h:1033
Generated on Thu Jul 2 2020 08:09:15 for Ice Sheet System Model by   doxygen 1.8.19