Ice Sheet System Model  4.18
Code documentation
StressbalanceSIAAnalysis.cpp
Go to the documentation of this file.
1 #include "./StressbalanceSIAAnalysis.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 StressbalanceSIAAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
10 
11  /*Intermediaries*/
12  bool isSIA,isSSA,isL1L2,isHO,isFS,iscoupling;
13 
14  /*Fetch parameters: */
15  iomodel->FindConstant(&isSIA,"md.flowequation.isSIA");
16  iomodel->FindConstant(&isSSA,"md.flowequation.isSSA");
17  iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2");
18  iomodel->FindConstant(&isHO,"md.flowequation.isHO");
19  iomodel->FindConstant(&isFS,"md.flowequation.isFS");
20 
21  /*Now, is the flag isSIA on? otherwise, do nothing: */
22  if (!isSIA) return;
23 
24  /*Do we have coupling*/
25  if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
26  iscoupling = true;
27  else
28  iscoupling = false;
29 
30  /*If no coupling, call Regular IoModelToConstraintsx, else, OLD stuff, keep for now*/
31  if(!iscoupling){
32  IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceSIAAnalysisEnum,P1Enum,0);
33  IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceSIAAnalysisEnum,P1Enum,1);
34  }
35  else{
36  /*Fetch data: */
37  iomodel->FetchData(3,"md.stressbalance.spcvx","md.stressbalance.spcvy","md.flowequation.vertex_equation");
38 
39  /*Initialize conunter*/
40  int count = 0;
41 
42  /*vx and vy are spc'd if we are not on nodeonSIA: */
43  for(int i=0;i<iomodel->numberofvertices;i++){
44  /*keep only this partition's nodes:*/
45  if((iomodel->my_vertices[i])){
46  if (IoCodeToEnumVertexEquation(reCast<int>(iomodel->Data("md.flowequation.vertex_equation")[i]))!=SIAApproximationEnum){
47 
48  constraints->AddObject(new SpcStatic(count+1,i+1,0,0,StressbalanceSIAAnalysisEnum));
49  count++;
50 
51  constraints->AddObject(new SpcStatic(count+1,i+1,1,0,StressbalanceSIAAnalysisEnum));
52  count++;
53  }
54  else{
55  if (!xIsNan<IssmDouble>(iomodel->Data("md.stressbalance.spcvx")[i])){
56  constraints->AddObject(new SpcStatic(count+1,i+1,0,iomodel->Data("md.stressbalance.spcvx")[i],StressbalanceSIAAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
57  count++;
58  }
59 
60  if (!xIsNan<IssmDouble>(iomodel->Data("md.stressbalance.spcvy")[i])){
61  constraints->AddObject(new SpcStatic(count+1,i+1,1,iomodel->Data("md.stressbalance.spcvy")[i],StressbalanceSIAAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
62  count++;
63  }
64  }
65  }
66  }
67  }
68 
69  /*Free data: */
70  iomodel->DeleteData(3,"md.stressbalance.spcvx","md.stressbalance.spcvy","md.flowequation.vertex_equation");
71 
72 }/*}}}*/
73 void StressbalanceSIAAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
74 
75  /*No loads*/
76 
77 }/*}}}*/
78 void StressbalanceSIAAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
79 
80  /*Intermediaries*/
81  bool isSIA;
82 
83  /*Fetch parameters: */
84  iomodel->FindConstant(&isSIA,"md.flowequation.isSIA");
85 
86  /*Now, is the flag isSIA on? otherwise, do nothing: */
87  if(!isSIA) return;
88 
89  /*First create nodes*/
90  int lid=0;
91  iomodel->FetchData(4,"md.flowequation.borderSSA","md.flowequation.borderFS","md.flowequation.vertex_equation","md.stressbalance.referential");
92  if(iomodel->domaintype!=Domain2DhorizontalEnum){
93  iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
94  }
95 
96  ::CreateNodes(nodes,iomodel,StressbalanceSIAAnalysisEnum,P1Enum,SIAApproximationEnum);
97  for(int i=0;i<nodes->Size();i++){
98  Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
99  int sid = node->Sid();
100  int approximation=IoCodeToEnumVertexEquation(reCast<int>(iomodel->Data("md.flowequation.vertex_equation")[sid]));
101  if(approximation!=SIAApproximationEnum) node->Deactivate();
102  }
103  iomodel->DeleteData(6,"md.mesh.vertexonbase","md.mesh.vertexonsurface","md.flowequation.borderSSA","md.flowequation.borderFS","md.flowequation.vertex_equation","md.stressbalance.referential");
104 
105 }/*}}}*/
106 int StressbalanceSIAAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
107  return 2;
108 }/*}}}*/
109 void StressbalanceSIAAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
110 
111  /*Fetch data needed: */
112  bool isSIA;
113  bool ismovingfront;
114  int frictionlaw;
115  iomodel->FindConstant(&isSIA,"md.flowequation.isSIA");
116  iomodel->FindConstant(&ismovingfront,"md.transient.ismovingfront");
117  iomodel->FindConstant(&frictionlaw,"md.friction.law");
118 
119  /*Now, is the flag SIA on? otherwise, do nothing: */
120  if (!isSIA)return;
121 
122  iomodel->FetchData(1,"md.flowequation.element_equation");
123 
124  /*Update elements: */
125  int counter=0;
126  for(int i=0;i<iomodel->numberofelements;i++){
127  if(iomodel->my_elements[i]){
128  Element* element=(Element*)elements->GetObjectByOffset(counter);
129  element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,P1Enum);
130  /*Need to know the type of approximation for this element*/
131  if(iomodel->Data("md.flowequation.element_equation")){
132  inputs2->SetInput(ApproximationEnum,counter,IoCodeToEnumElementEquation(reCast<int>(iomodel->Data("md.flowequation.element_equation")[i])));
133  }
134  counter++;
135  }
136  }
137 
138  /*Free data: */
139  iomodel->DeleteData(1,"md.flowequation.element_equation");
140 
141  /*Friction law variables*/
142  switch(frictionlaw){
143  case 1:
144  iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficient",FrictionCoefficientEnum);
145  iomodel->FetchDataToInput(inputs2,elements,"md.friction.p",FrictionPEnum);
146  iomodel->FetchDataToInput(inputs2,elements,"md.friction.q",FrictionQEnum);
147  break;
148  case 2:
149  iomodel->FetchDataToInput(inputs2,elements,"md.friction.C",FrictionCEnum);
150  iomodel->FetchDataToInput(inputs2,elements,"md.friction.m",FrictionMEnum);
151  break;
152  case 6:
153  iomodel->FetchDataToInput(inputs2,elements,"md.friction.C",FrictionCEnum);
154  iomodel->FetchDataToInput(inputs2,elements,"md.friction.m",FrictionMEnum);
155  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.pressure",PressureEnum);
156  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.temperature",TemperatureEnum);
157  break;
158  default:
159  _error_("not supported");
160  }
161 
162  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.thickness",ThicknessEnum);
163  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
164  if(ismovingfront){
165  if(iomodel->domaintype!=Domain2DhorizontalEnum)
166  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum); // required for updating active nodes
167  }
168 
169 }/*}}}*/
170 void StressbalanceSIAAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
171 
172  /*No specific parameters*/
173 
174 }/*}}}*/
175 
176 /*Finite Element Analysis*/
177 void StressbalanceSIAAnalysis::Core(FemModel* femmodel){/*{{{*/
178 
179  if(VerboseSolution()) _printf0_(" computing SIA velocities\n");
180  femmodel->SetCurrentConfiguration(StressbalanceSIAAnalysisEnum);
181  solutionsequence_linear(femmodel);
182 }/*}}}*/
183 ElementVector* StressbalanceSIAAnalysis::CreateDVector(Element* element){/*{{{*/
184  /*Default, return NULL*/
185  return NULL;
186 }/*}}}*/
187 ElementMatrix* StressbalanceSIAAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
188 _error_("Not implemented");
189 }/*}}}*/
190 ElementMatrix* StressbalanceSIAAnalysis::CreateKMatrix(Element* element){/*{{{*/
191 
192  /* Check if ice in element */
193  if(!element->IsIceInElement()) return NULL;
194 
195  int domaintype;
196  element->FindParam(&domaintype,DomainTypeEnum);
197  switch(domaintype){
198  case Domain2DhorizontalEnum:
199  return CreateKMatrix2D(element);
200  case Domain3DEnum:
201  return CreateKMatrix3D(element);
202  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
203  }
204 }/*}}}*/
205 ElementMatrix* StressbalanceSIAAnalysis::CreateKMatrix2D(Element* element){/*{{{*/
206 
207  /* Check if ice in element */
208  if(!element->IsIceInElement()) return NULL;
209 
210  /*Intermediaries */
211  IssmDouble connectivity;
212 
213  /*Fetch number vertices for this element*/
214  int numvertices = element->GetNumberOfVertices();
215 
216  /*Initialize Element vector*/
217  ElementMatrix* Ke=element->NewElementMatrix();
218  for(int iv=0;iv<numvertices;iv++){
219  connectivity=(IssmDouble)element->VertexConnectivity(iv);
220  Ke->values[(2*iv+0)*2*numvertices+(2*iv+0)]=1./connectivity;
221  Ke->values[(2*iv+1)*2*numvertices+(2*iv+1)]=1./connectivity;
222  }
223 
224  /*Clean up and return*/
225  return Ke;
226 }/*}}}*/
227 ElementMatrix* StressbalanceSIAAnalysis::CreateKMatrix3D(Element* element){/*{{{*/
228 
229  /* Check if ice in element */
230  if(!element->IsIceInElement()) return NULL;
231 
232  /*Intermediaries */
233  int i0,i1,j0,j1,nodeup,nodedown,numsegments;
234  IssmDouble slope[2],connectivity[2],one0,one1;
235  int *pairindices = NULL;
236 
237  /*Fetch number vertices for this element*/
238  int numvertices = element->GetNumberOfVertices();
239  int numdof = 2*numvertices;
240 
241  /*Initialize Element vector*/
242  ElementMatrix* Ke=element->NewElementMatrix();
243 
244  element->VerticalSegmentIndices(&pairindices,&numsegments);
245  for(int is=0;is<numsegments;is++){
246  nodedown = pairindices[is*2+0];
247  nodeup = pairindices[is*2+1];
248  connectivity[0]=(IssmDouble)element->VertexConnectivity(nodedown);
249  connectivity[1]=(IssmDouble)element->VertexConnectivity(nodeup);
250  one0=1./connectivity[0];
251  one1=1./connectivity[1];
252 
253  /*2 dofs of first node*/
254  i0=2*nodedown; i1=2*nodedown+1;
255  /*2 dofs of second node*/
256  j0=2*nodeup; j1=2*nodeup+1;
257 
258  /*Create matrix for these two nodes*/
259  if(element->IsOnBase() && element->IsOnSurface()){
260  Ke->values[i0*numdof+i0] = +one0;
261  Ke->values[i1*numdof+i1] = +one0;
262  Ke->values[j0*numdof+i0] = -one1;
263  Ke->values[j0*numdof+j0] = +one1;
264  Ke->values[j1*numdof+i1] = -one1;
265  Ke->values[j1*numdof+j1] = +one1;
266  }
267  else if(element->IsOnBase()){
268  Ke->values[i0*numdof+i0] = one0;
269  Ke->values[i1*numdof+i1] = one0;
270  Ke->values[j0*numdof+i0] = -2.*one1;
271  Ke->values[j0*numdof+j0] = +2.*one1;
272  Ke->values[j1*numdof+i1] = -2.*one1;
273  Ke->values[j1*numdof+j1] = +2.*one1;
274  }
275  else if(element->IsOnSurface()){
276  Ke->values[j0*numdof+i0] = -one1;
277  Ke->values[j0*numdof+j0] = +one1;
278  Ke->values[j1*numdof+i1] = -one1;
279  Ke->values[j1*numdof+j1] = +one1;
280  }
281  else{ //node is on two horizontal layers and beams include the values only once, so the have to use half of the connectivity
282  Ke->values[j0*numdof+i0] = -2.*one1;
283  Ke->values[j0*numdof+j0] = +2.*one1;
284  Ke->values[j1*numdof+i1] = -2.*one1;
285  Ke->values[j1*numdof+j1] = +2.*one1;
286  }
287  }
288 
289  /*Clean up and return*/
290  xDelete<int>(pairindices);
291  return Ke;
292 }/*}}}*/
293 ElementVector* StressbalanceSIAAnalysis::CreatePVector(Element* element){/*{{{*/
294 
295  /* Check if ice in element */
296  if(!element->IsIceInElement()) return NULL;
297 
298  int domaintype;
299  element->FindParam(&domaintype,DomainTypeEnum);
300  switch(domaintype){
301  case Domain2DhorizontalEnum:
302  return CreatePVector2D(element);
303  case Domain3DEnum:
304  return CreatePVector3D(element);
305  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
306  }
307 }/*}}}*/
308 ElementVector* StressbalanceSIAAnalysis::CreatePVector2D(Element* element){/*{{{*/
309 
310  /* Check if ice in element */
311  if(!element->IsIceInElement()) return NULL;
312 
313  /*Intermediaries */
314  int frictionlaw = 1;
315  IssmDouble ub,vb,slope2,drag,thickness,surface,connectivity;
316  IssmDouble slope[2];
317 
318  /*Fetch number vertices for this element*/
319  int numvertices = element->GetNumberOfVertices();
320 
321  /*Initialize Element vector*/
322  ElementVector* pe=element->NewElementVector();
323 
324  /*Retrieve all inputs and parameters*/
325  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
326  IssmDouble gravity = element->FindParam(ConstantsGEnum);
327  IssmDouble B,n;
328  Input2* B_input = element->GetInput2(MaterialsRheologyBbarEnum);_assert_(B_input);
329  Input2* n_input = element->GetInput2(MaterialsRheologyNEnum); _assert_(n_input);
330  Input2* slopex_input = element->GetInput2(SurfaceSlopeXEnum); _assert_(slopex_input);
331  Input2* slopey_input = element->GetInput2(SurfaceSlopeYEnum); _assert_(slopey_input);
332  Input2* thickness_input = element->GetInput2(ThicknessEnum); _assert_(thickness_input);
333  Input2* surface_input = element->GetInput2(SurfaceEnum); _assert_(surface_input);
334  Input2* drag_input = NULL;
335  if(frictionlaw!=5 && frictionlaw!=1){
336  drag_input = element->GetInput2(FrictionCoefficientEnum); _assert_(drag_input);
337  }
338 
339  Gauss* gauss=element->NewGauss();
340  for(int iv=0;iv<numvertices;iv++){
341  gauss->GaussVertex(iv);
342 
343  connectivity=(IssmDouble)element->VertexConnectivity(iv);
344 
345  B_input->GetInputValue(&B,gauss);
346  n_input->GetInputValue(&n,gauss);
347  thickness_input->GetInputValue(&thickness,gauss);
348  surface_input->GetInputValue(&surface,gauss);
349  slopex_input->GetInputValue(&slope[0],gauss);
350  slopey_input->GetInputValue(&slope[1],gauss);
351  slope2=slope[0]*slope[0]+slope[1]*slope[1];
352 
353  switch(frictionlaw){
354  case 1:
355  /*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10 m/s/Pa*/
356  ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
357  vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
358  break;
359  case 2:
360  /*Ritz et al. 1996*/
361  drag_input->GetInputValue(&drag,gauss);
362  ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
363  vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
364  break;
365  case 3:
366  /*Rutt et al. 2009*/
367  drag_input->GetInputValue(&drag,gauss);
368  ub=-drag*rho_ice*gravity*thickness*slope[0];
369  vb=-drag*rho_ice*gravity*thickness*slope[1];
370  break;
371  case 4:
372  /*Henning Akesson*/
373  drag = -4e-15 * surface + 8.6e-12;
374  ub=-drag*rho_ice*gravity*thickness*slope[0];
375  vb=-drag*rho_ice*gravity*thickness*slope[1];
376  break;
377  case 6:
378  /*No sliding*/
379  ub=0.;
380  vb=0.;
381  break;
382  default:
383  _error_("Not supported yet");
384  }
385 
386  pe->values[2*iv+0]=(ub-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n+1.)/(pow(B,n)*(n+2))*slope[0])/connectivity;
387  pe->values[2*iv+1]=(vb-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n+1.)/(pow(B,n)*(n+2))*slope[1])/connectivity;
388  }
389 
390  /*Clean up and return*/
391  delete gauss;
392  return pe;
393 }/*}}}*/
394 ElementVector* StressbalanceSIAAnalysis::CreatePVector3D(Element* element){/*{{{*/
395 
396  /* Check if ice in element */
397  if(!element->IsIceInElement()) return NULL;
398 
399  /*Intermediaries */
400  int frictionlaw = 1;
401  int nodeup,nodedown,numsegments;
402  IssmDouble ub,vb,slope2,drag,surface,thickness,constant_part,z,Jdet;
403  IssmDouble slope[2],connectivity[2],xyz_list_line[2][3];
404  IssmDouble *xyz_list = NULL;
405  int *pairindices = NULL;
406 
407  /*Fetch number vertices for this element*/
408  int numvertices = element->GetNumberOfVertices();
409 
410  /*Initialize Element vector*/
411  ElementVector* pe=element->NewElementVector();
412 
413  /*Retrieve all inputs and parameters*/
414  element->GetVerticesCoordinates(&xyz_list);
415  IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
416  IssmDouble gravity = element->FindParam(ConstantsGEnum);
417  IssmDouble B,n;
418  Input2* B_input = element->GetInput2(MaterialsRheologyBEnum); _assert_(B_input);
419  Input2* n_input = element->GetInput2(MaterialsRheologyNEnum); _assert_(n_input);
420  Input2* surface_input = element->GetInput2(SurfaceEnum); _assert_(surface_input);
421  Input2* slopex_input = element->GetInput2(SurfaceSlopeXEnum); _assert_(slopex_input);
422  Input2* slopey_input = element->GetInput2(SurfaceSlopeYEnum); _assert_(slopey_input);
423  Input2* thickness_input = element->GetInput2(ThicknessEnum); _assert_(thickness_input);
424  Input2* drag_input = NULL;
425  Friction* friction = NULL;
426  if(frictionlaw!=5 && frictionlaw!=1){
427  drag_input = element->GetInput2(FrictionCoefficientEnum); _assert_(drag_input);
428  }
429  else if(frictionlaw==5){
430  friction=new Friction(element,3);
431  }
432 
433  /*Get Vertical segment indices*/
434  element->VerticalSegmentIndices(&pairindices,&numsegments);
435  for(int is=0;is<numsegments;is++){
436  nodedown = pairindices[is*2+0];
437  nodeup = pairindices[is*2+1];
438  connectivity[0]=(IssmDouble)element->VertexConnectivity(nodedown);
439  connectivity[1]=(IssmDouble)element->VertexConnectivity(nodeup);
440  for(int i=0;i<3;i++){
441  xyz_list_line[0][i]=xyz_list[nodedown*3+i];
442  xyz_list_line[1][i]=xyz_list[nodeup*3+i];
443  }
444 
445  Gauss* gauss=element->NewGaussLine(nodedown,nodeup,3);
446  for(int ig=gauss->begin();ig<gauss->end();ig++){
447  gauss->GaussPoint(ig);
448 
449  B_input->GetInputValue(&B,gauss);
450  n_input->GetInputValue(&n,gauss);
451  slopex_input->GetInputValue(&slope[0],gauss);
452  slopey_input->GetInputValue(&slope[1],gauss);
453  surface_input->GetInputValue(&surface,gauss);
454  thickness_input->GetInputValue(&thickness,gauss);
455 
456  slope2=slope[0]*slope[0]+slope[1]*slope[1];
457  constant_part=-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.));
458 
459  z = element->GetZcoord(xyz_list,gauss);
460  element->JacobianDeterminantLine(&Jdet,&xyz_list_line[0][0],gauss);
461 
462  if(element->IsOnSurface()){
463  pe->values[2*nodeup+0]+=constant_part*pow((surface-z)/B,n)*slope[0]*Jdet*gauss->weight/connectivity[1];
464  pe->values[2*nodeup+1]+=constant_part*pow((surface-z)/B,n)*slope[1]*Jdet*gauss->weight/connectivity[1];
465  }
466  else{/*connectivity is too large, should take only half on it*/
467  pe->values[2*nodeup+0]+=constant_part*pow((surface-z)/B,n)*slope[0]*Jdet*gauss->weight*2./connectivity[1];
468  pe->values[2*nodeup+1]+=constant_part*pow((surface-z)/B,n)*slope[1]*Jdet*gauss->weight*2./connectivity[1];
469  }
470  }
471 
472  /*Deal with basal velocities*/
473  if(element->IsOnBase()){
474  /*Make sure to now push the gauss point to the base*/
475  delete gauss; gauss=element->NewGauss();
476  gauss->GaussVertex(nodedown);
477  switch(frictionlaw){
478  case 1:
479  /*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10 m/s/Pa*/
480  ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
481  vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
482  break;
483  case 2:
484  /*Ritz et al. 1996*/
485  drag_input->GetInputValue(&drag,gauss);
486  ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
487  vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
488  break;
489  case 3:
490  /*Rutt et al. 2009*/
491  drag_input->GetInputValue(&drag,gauss);
492  ub=-drag*rho_ice*gravity*thickness*slope[0];
493  vb=-drag*rho_ice*gravity*thickness*slope[1];
494  break;
495  case 4:
496  /*Henning Akesson*/
497  drag_input->GetInputValue(&drag,gauss);
498  drag = -4e-15 * surface + 8.6e-12;
499  ub=-drag*rho_ice*gravity*thickness*slope[0];
500  vb=-drag*rho_ice*gravity*thickness*slope[1];
501  break;
502  case 5: /*Weertman temp for Kevin*/{
503  friction->GetAlpha2WeertmanTemp(&drag,gauss);
504  ub = -1./drag * rho_ice*gravity*thickness*slope[0];
505  vb = -1./drag * rho_ice*gravity*thickness*slope[1];
506  }
507  break;
508  default:
509  _error_("Not supported yet");
510  }
511 
512  pe->values[2*nodedown+0]+=ub/connectivity[0];
513  pe->values[2*nodedown+1]+=vb/connectivity[0];
514  }
515  delete gauss;
516  }
517 
518  /*Clean up and return*/
519  xDelete<int>(pairindices);
520  xDelete<IssmDouble>(xyz_list);
521  if(frictionlaw==5) delete friction;
522  return pe;
523 
524 }/*}}}*/
525 void StressbalanceSIAAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
526 
527  IssmDouble vx,vy;
528  int *doflist = NULL;
529 
530  /*Fetch number of nodes and initialize values*/
531  int numnodes = element->GetNumberOfNodes();
532  int numdof = numnodes*2;
533  IssmDouble* values = xNew<IssmDouble>(numdof);
534 
535  /*Get dof list and inputs */
536  element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
537  Input2* vx_input=element->GetInput2(VxEnum); _assert_(vx_input);
538  Input2* vy_input=element->GetInput2(VyEnum); _assert_(vy_input);
539 
540  /*Ok, we have the velocities in inputs, fill in solution */
541  Gauss* gauss=element->NewGauss();
542  for(int i=0;i<numnodes;i++){
543  gauss->GaussVertex(i);
544  vx_input->GetInputValue(&vx,gauss);
545  vy_input->GetInputValue(&vy,gauss);
546  values[i*2+0]=vx;
547  values[i*2+1]=vy;
548  }
549 
550  /*Add value to global vector*/
551  solution->SetValues(numdof,doflist,values,INS_VAL);
552 
553  /*Free ressources:*/
554  delete gauss;
555  xDelete<int>(doflist);
556  xDelete<IssmDouble>(values);
557 }/*}}}*/
558 void StressbalanceSIAAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
559  _error_("Not implemented yet");
560 }/*}}}*/
561 void StressbalanceSIAAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
562 
563  int i,domaintype;
564  IssmDouble rho_ice,g;
565  int* doflist=NULL;
566  IssmDouble* xyz_list=NULL;
567 
568  /*Fetch number of nodes and dof for this finite element*/
569  int numnodes = element->GetNumberOfNodes();
570  int numdof = numnodes*2;
571 
572  /*Fetch dof list and allocate solution vectors*/
573  element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
574  IssmDouble* values = xNew<IssmDouble>(numdof);
575  IssmDouble* vx = xNew<IssmDouble>(numdof);
576  IssmDouble* vy = xNew<IssmDouble>(numdof);
577  IssmDouble* vz = xNew<IssmDouble>(numdof);
578  IssmDouble* vel = xNew<IssmDouble>(numdof);
579  IssmDouble* pressure = xNew<IssmDouble>(numdof);
580  IssmDouble* thickness = xNew<IssmDouble>(numdof);
581  IssmDouble* surface = xNew<IssmDouble>(numnodes);
582 
583  /*Use the dof list to index into the solution vector: */
584  for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
585 
586  /*Transform solution in Cartesian Space*/
587  element->TransformSolutionCoord(&values[0],XYEnum);
588 
589  /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
590  for(i=0;i<numnodes;i++){
591  vx[i]=values[i*2+0];
592  vy[i]=values[i*2+1];
593 
594  /*Check solution*/
595  if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
596  if(xIsInf<IssmDouble>(vx[i])) _error_("Inf found in solution vector");
597  if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
598  if(xIsInf<IssmDouble>(vy[i])) _error_("Inf found in solution vector");
599  }
600 
601  /*Get Vz and compute vel*/
602  element->GetInputListOnNodes(&vz[0],VzEnum,0.);
603  for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
604 
605  /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
606  *so the pressure is just the pressure at the bedrock: */
607  rho_ice = element->FindParam(MaterialsRhoIceEnum);
608  g = element->FindParam(ConstantsGEnum);
609  element->FindParam(&domaintype,DomainTypeEnum);
610  switch(domaintype){
611  case Domain2DhorizontalEnum:
612  element->GetInputListOnNodes(&thickness[0],ThicknessEnum);
613  for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*thickness[i];
614  break;
615  case Domain3DEnum:
616  element->GetVerticesCoordinates(&xyz_list);
617  element->GetInputListOnNodes(&surface[0],SurfaceEnum);
618  for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
619  break;
620  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
621  }
622 
623  /*Add vx and vy as inputs to the tria element: */
624  element->AddInput2(VxEnum,vx,P1Enum);
625  element->AddInput2(VyEnum,vy,P1Enum);
626  element->AddInput2(VelEnum,vel,P1Enum);
627  element->AddInput2(PressureEnum,pressure,P1Enum);
628 
629  /*Free ressources:*/
630  xDelete<IssmDouble>(thickness);
631  xDelete<IssmDouble>(surface);
632  xDelete<IssmDouble>(pressure);
633  xDelete<IssmDouble>(vel);
634  xDelete<IssmDouble>(vz);
635  xDelete<IssmDouble>(vy);
636  xDelete<IssmDouble>(vx);
637  xDelete<IssmDouble>(values);
638  xDelete<IssmDouble>(xyz_list);
639  xDelete<int>(doflist);
640 }/*}}}*/
641 void StressbalanceSIAAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
642  SetActiveNodesLSMx(femmodel);
643 }/*}}}*/
DataSet::Size
int Size()
Definition: DataSet.cpp:399
FrictionCoefficientEnum
@ FrictionCoefficientEnum
Definition: EnumDefinitions.h:574
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::TransformSolutionCoord
void TransformSolutionCoord(IssmDouble *solution, int cs_enum)
Definition: Element.cpp:4397
Element::IsOnBase
bool IsOnBase()
Definition: Element.cpp:1984
Nodes
Declaration of Nodes class.
Definition: Nodes.h:19
Element::GetDofList
void GetDofList(int **pdoflist, int approximation_enum, int setenum)
Definition: Element.cpp:961
Element::GetNumberOfNodes
virtual int GetNumberOfNodes(void)=0
FrictionPEnum
@ FrictionPEnum
Definition: EnumDefinitions.h:578
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
SurfaceSlopeXEnum
@ SurfaceSlopeXEnum
Definition: EnumDefinitions.h:829
StressbalanceSIAAnalysis::UpdateParameters
void UpdateParameters(Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
Definition: StressbalanceSIAAnalysis.cpp:170
StressbalanceSIAAnalysis::CreateKMatrix3D
ElementMatrix * CreateKMatrix3D(Element *element)
Definition: StressbalanceSIAAnalysis.cpp:227
DataSet::AddObject
int AddObject(Object *object)
Definition: DataSet.cpp:252
MaskOceanLevelsetEnum
@ MaskOceanLevelsetEnum
Definition: EnumDefinitions.h:640
Element::VerticalSegmentIndices
virtual void VerticalSegmentIndices(int **pindices, int *pnumseg)=0
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
StressbalanceSIAAnalysis::CreateNodes
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
Definition: StressbalanceSIAAnalysis.cpp:78
Gauss::end
virtual int end(void)=0
Constraints
Declaration of Constraints class.
Definition: Constraints.h:13
IoCodeToEnumElementEquation
int IoCodeToEnumElementEquation(int enum_in)
Definition: IoCodeConversions.cpp:311
MeshVertexonbaseEnum
@ MeshVertexonbaseEnum
Definition: EnumDefinitions.h:653
Elements
Declaration of Elements class.
Definition: Elements.h:17
PressureEnum
@ PressureEnum
Definition: EnumDefinitions.h:664
MaterialsRhoIceEnum
@ MaterialsRhoIceEnum
Definition: EnumDefinitions.h:264
StressbalanceSIAAnalysis::CreateKMatrix
ElementMatrix * CreateKMatrix(Element *element)
Definition: StressbalanceSIAAnalysis.cpp:190
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
Element::GetZcoord
IssmDouble GetZcoord(IssmDouble *xyz_list, Gauss *gauss)
Definition: Element.cpp:1494
FrictionQEnum
@ FrictionQEnum
Definition: EnumDefinitions.h:580
StressbalanceSIAAnalysis::DofsPerNode
int DofsPerNode(int **doflist, int domaintype, int approximation)
Definition: StressbalanceSIAAnalysis.cpp:106
MaterialsRheologyNEnum
@ MaterialsRheologyNEnum
Definition: EnumDefinitions.h:651
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
solutionsequence_linear
void solutionsequence_linear(FemModel *femmodel)
Definition: solutionsequence_linear.cpp:10
StressbalanceSIAAnalysis::UpdateElements
void UpdateElements(Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
Definition: StressbalanceSIAAnalysis.cpp:109
SpcStatic
Definition: SpcStatic.h:13
IoModel::my_vertices
bool * my_vertices
Definition: IoModel.h:72
MaterialsRheologyBbarEnum
@ MaterialsRheologyBbarEnum
Definition: EnumDefinitions.h:644
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
Element
Definition: Element.h:41
IoModel::numberofvertices
int numberofvertices
Definition: IoModel.h:99
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
StressbalanceSIAAnalysis::CreateDVector
ElementVector * CreateDVector(Element *element)
Definition: StressbalanceSIAAnalysis.cpp:183
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
XYEnum
@ XYEnum
Definition: EnumDefinitions.h:1330
ConstantsGEnum
@ ConstantsGEnum
Definition: EnumDefinitions.h:102
NoneApproximationEnum
@ NoneApproximationEnum
Definition: EnumDefinitions.h:1201
Vector::SetValues
void SetValues(int ssize, int *list, doubletype *values, InsMode mode)
Definition: Vector.h:153
DomainTypeEnum
@ DomainTypeEnum
Definition: EnumDefinitions.h:124
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
VzEnum
@ VzEnum
Definition: EnumDefinitions.h:853
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
StressbalanceSIAAnalysis::CreateConstraints
void CreateConstraints(Constraints *constraints, IoModel *iomodel)
Definition: StressbalanceSIAAnalysis.cpp:9
GsetEnum
@ GsetEnum
Definition: EnumDefinitions.h:1093
IoModel::FindConstant
void FindConstant(bool *pvalue, const char *constant_name)
Definition: IoModel.cpp:2362
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
StressbalanceSIAAnalysis::GetSolutionFromInputs
void GetSolutionFromInputs(Vector< IssmDouble > *solution, Element *element)
Definition: StressbalanceSIAAnalysis.cpp:525
StressbalanceSIAAnalysis::UpdateConstraints
void UpdateConstraints(FemModel *femmodel)
Definition: StressbalanceSIAAnalysis.cpp:641
INS_VAL
@ INS_VAL
Definition: toolkitsenums.h:14
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
Friction
Definition: Friction.h:13
StressbalanceSIAAnalysis::GradientJ
void GradientJ(Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
Definition: StressbalanceSIAAnalysis.cpp:558
Element::NewGaussLine
virtual Gauss * NewGaussLine(int vertex1, int vertex2, int order)=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
StressbalanceSIAAnalysis::CreatePVector3D
ElementVector * CreatePVector3D(Element *element)
Definition: StressbalanceSIAAnalysis.cpp:394
TemperatureEnum
@ TemperatureEnum
Definition: EnumDefinitions.h:831
Gauss::GaussVertex
virtual void GaussVertex(int iv)=0
StressbalanceSIAAnalysisEnum
@ StressbalanceSIAAnalysisEnum
Definition: EnumDefinitions.h:1287
Element::IsIceInElement
bool IsIceInElement()
Definition: Element.cpp:2021
IoModel::Data
IssmDouble * Data(const char *data_name)
Definition: IoModel.cpp:437
Loads
Declaration of Loads class.
Definition: Loads.h:16
StressbalanceSIAAnalysis.h
: header file for generic external result object
Inputs2::SetInput
void SetInput(int enum_in, int index, bool value)
Definition: Inputs2.cpp:572
Node
Definition: Node.h:23
Node::Sid
int Sid(void)
Definition: Node.cpp:622
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Node::Deactivate
void Deactivate(void)
Definition: Node.cpp:681
SetActiveNodesLSMx
void SetActiveNodesLSMx(FemModel *femmodel)
Definition: SetActiveNodesLSMx.cpp:12
MaterialsRheologyBEnum
@ MaterialsRheologyBEnum
Definition: EnumDefinitions.h:643
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
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
StressbalanceSIAAnalysis::CreateJacobianMatrix
ElementMatrix * CreateJacobianMatrix(Element *element)
Definition: StressbalanceSIAAnalysis.cpp:187
FemModel::SetCurrentConfiguration
void SetCurrentConfiguration(int configuration_type)
Definition: FemModel.cpp:634
ThicknessEnum
@ ThicknessEnum
Definition: EnumDefinitions.h:840
IoModel::FetchDataToInput
void FetchDataToInput(Inputs2 *inputs2, Elements *elements, const char *vector_name, int input_enum)
Definition: IoModel.cpp:1651
Element::VertexConnectivity
virtual int VertexConnectivity(int vertexindex)=0
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
StressbalanceSIAAnalysis::CreatePVector
ElementVector * CreatePVector(Element *element)
Definition: StressbalanceSIAAnalysis.cpp:293
Gauss::weight
IssmDouble weight
Definition: Gauss.h:11
FrictionMEnum
@ FrictionMEnum
Definition: EnumDefinitions.h:577
Friction::GetAlpha2WeertmanTemp
void GetAlpha2WeertmanTemp(IssmDouble *palpha2, Gauss *gauss)
Definition: Friction.cpp:517
Element::IsOnSurface
bool IsOnSurface()
Definition: Element.cpp:1981
IoModel
Definition: IoModel.h:48
IoCodeToEnumVertexEquation
int IoCodeToEnumVertexEquation(int enum_in)
Definition: IoCodeConversions.cpp:297
StressbalanceSIAAnalysis::CreatePVector2D
ElementVector * CreatePVector2D(Element *element)
Definition: StressbalanceSIAAnalysis.cpp:308
StressbalanceSIAAnalysis::Core
void Core(FemModel *femmodel)
Definition: StressbalanceSIAAnalysis.cpp:177
Element::GetNumberOfVertices
virtual int GetNumberOfVertices(void)=0
Element::JacobianDeterminantLine
virtual void JacobianDeterminantLine(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)=0
IoModel::domaintype
int domaintype
Definition: IoModel.h:78
ElementMatrix
Definition: ElementMatrix.h:19
FrictionCEnum
@ FrictionCEnum
Definition: EnumDefinitions.h:572
Vector< IssmDouble >
SIAApproximationEnum
@ SIAApproximationEnum
Definition: EnumDefinitions.h:1241
StressbalanceSIAAnalysis::CreateLoads
void CreateLoads(Loads *loads, IoModel *iomodel)
Definition: StressbalanceSIAAnalysis.cpp:73
Gauss
Definition: Gauss.h:8
StressbalanceSIAAnalysis::CreateKMatrix2D
ElementMatrix * CreateKMatrix2D(Element *element)
Definition: StressbalanceSIAAnalysis.cpp:205
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
SurfaceSlopeYEnum
@ SurfaceSlopeYEnum
Definition: EnumDefinitions.h:830
StressbalanceSIAAnalysis::InputUpdateFromSolution
void InputUpdateFromSolution(IssmDouble *solution, Element *element)
Definition: StressbalanceSIAAnalysis.cpp:561
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:15 for Ice Sheet System Model by   doxygen 1.8.19