Ice Sheet System Model  4.18
Code documentation
MasstransportAnalysis.cpp
Go to the documentation of this file.
2 #include "../toolkits/toolkits.h"
3 #include "../classes/classes.h"
4 #include "../shared/shared.h"
5 #include "../modules/modules.h"
6 
7 #define FINITEELEMENT P1Enum
8 
9 /*Model processing*/
11 
12  /*Fetch parameters: */
13  int stabilization;
14  iomodel->FindConstant(&stabilization,"md.masstransport.stabilization");
15 
16  /*Do not add constraints in DG, they are weakly imposed*/
17  if(stabilization!=3){
18  IoModelToConstraintsx(constraints,iomodel,"md.masstransport.spcthickness",MasstransportAnalysisEnum,FINITEELEMENT);
19  }
20 
21  /*FCT, constraints are imposed using penalties*/
22  if(stabilization==4){
24  }
25 }/*}}}*/
26 void MasstransportAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
27 
28  /*Intermediaries*/
29  int penpair_ids[2];
30  int count=0;
31  int stabilization;
32  int numvertex_pairing;
33 
34  /*Fetch parameters: */
35  iomodel->FindConstant(&stabilization,"md.masstransport.stabilization");
36 
37  /*Loads only in DG*/
38  if(stabilization==3){
39 
40  /*Get faces and elements*/
41  CreateFaces(iomodel);
42  iomodel->FetchData(1,"md.geometry.thickness");
43 
44  /*First load data:*/
45  for(int i=0;i<iomodel->numberoffaces;i++){
46 
47  /*Get left and right elements*/
48  int element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
49 
50  /*Now, if this element is not in the partition, pass: */
51  if(!iomodel->my_elements[element]) continue;
52  loads->AddObject(new Numericalflux(i+1,i,i,iomodel));
53  }
54 
55  /*Free data: */
56  iomodel->DeleteData(1,"md.geometry.thickness");
57  }
58 
59  /*Create Penpair for vertex_pairing: */
60  IssmDouble *vertex_pairing=NULL;
61  IssmDouble *nodeonbase=NULL;
62  iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,"md.masstransport.vertex_pairing");
63  if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(&nodeonbase,NULL,NULL,"md.mesh.vertexonbase");
64 
65  for(int i=0;i<numvertex_pairing;i++){
66 
67  if(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+0])-1]){
68 
69  /*In debugging mode, check that the second node is in the same cpu*/
70  _assert_(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]);
71 
72  /*Skip if one of the two is not on the bed*/
73  if(iomodel->domaintype!=Domain2DhorizontalEnum){
74  if(!(reCast<bool>(nodeonbase[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonbase[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
75  }
76 
77  /*Get node ids*/
78  penpair_ids[0]=reCast<int>(vertex_pairing[2*i+0]);
79  penpair_ids[1]=reCast<int>(vertex_pairing[2*i+1]);
80 
81  /*Create Load*/
82  loads->AddObject(new Penpair(count+1,&penpair_ids[0]));
83  count++;
84  }
85  }
86 
87  /*free ressources: */
88  iomodel->DeleteData(vertex_pairing,"md.masstransport.vertex_pairing");
89  iomodel->DeleteData(nodeonbase,"md.mesh.vertexonbase");
90 }/*}}}*/
91 void MasstransportAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
92 
93  /*Fetch parameters: */
94  int stabilization;
95  iomodel->FindConstant(&stabilization,"md.masstransport.stabilization");
96 
97  /*Check in 3d*/
98  if(stabilization==3 && iomodel->domaintype==Domain3DEnum) _error_("DG 3d not implemented yet");
99 
100  /*Create Nodes either DG or CG depending on stabilization*/
101  if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
102  if(stabilization!=3){
104  }
105  else{
106  ::CreateNodes(nodes,iomodel,MasstransportAnalysisEnum,P1DGEnum,isamr);
107  }
108  iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
109 }/*}}}*/
110 int MasstransportAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
111  return 1;
112 }/*}}}*/
113 void MasstransportAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
114 
115  int stabilization,finiteelement;
116  bool dakota_analysis;
117  bool isgroundingline;
118  bool ismovingfront;
119  bool isoceancoupling;
120  bool issmb;
121 
122  /*Fetch data needed: */
123  iomodel->FindConstant(&stabilization,"md.masstransport.stabilization");
124  iomodel->FindConstant(&dakota_analysis,"md.qmu.isdakota");
125  iomodel->FindConstant(&isgroundingline,"md.transient.isgroundingline");
126  iomodel->FindConstant(&ismovingfront,"md.transient.ismovingfront");
127  iomodel->FindConstant(&isoceancoupling,"md.transient.isoceancoupling");
128  iomodel->FindConstant(&issmb,"md.transient.issmb");
129 
130  /*Finite element type*/
131  finiteelement = FINITEELEMENT;
132  if(stabilization==3){
133  finiteelement = P1DGEnum;
134  //finiteelement = P0DGEnum;
135  }
136 
137  /*Update elements: */
138  int counter=0;
139  for(int i=0;i<iomodel->numberofelements;i++){
140  if(iomodel->my_elements[i]){
141  Element* element=(Element*)elements->GetObjectByOffset(counter);
142  element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,finiteelement);
143  counter++;
144  }
145  }
146 
147  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.thickness",ThicknessEnum);
148  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.surface",SurfaceEnum);
149  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.base",BaseEnum);
150  iomodel->FetchDataToInput(inputs2,elements,"md.solidearth.sealevel",SealevelEnum,0);
151  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
152  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
153  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
154  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vx",VxEnum);
155  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vy",VyEnum);
156  if(isgroundingline) iomodel->FetchDataToInput(inputs2,elements,"md.geometry.bed",BedEnum);
157  /*Initialize cumdeltalthickness input*/
159  /*Initialize ThicknessResidual input*/
160  InputUpdateFromConstantx(inputs2,elements,0.,ThicknessResidualEnum);
161 
162  /*Get what we need for ocean-induced basal melting*/
163  int basalforcing_model;
164  iomodel->FindConstant(&basalforcing_model,"md.basalforcings.model");
165  switch(basalforcing_model){
167  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.floatingice_melting_rate",BasalforcingsFloatingiceMeltingRateEnum);
168  break;
170  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.perturbation_melting_rate",BasalforcingsPerturbationMeltingRateEnum,0.);
171  break;
173  break;
175  break;
177  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum);
178  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum);
179  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum);
180  break;
182  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.basin_id",BasalforcingsPicoBasinIdEnum);
183  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.overturning_coeff",BasalforcingsPicoOverturningCoeffEnum);
184  break;
186  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum);
187  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.melt_anomaly",BasalforcingsIsmip6MeltAnomalyEnum,0.);
188  IssmDouble** array3d = NULL; int* Ms = NULL; int* Ns = NULL; int K;
189  iomodel->FetchData(&array3d,&Ms,&Ns,&K,"md.basalforcings.tf");
190  if(!array3d) _error_("md.basalforcings.tf not found in binary file");
191  for(int i=0;i<elements->Size();i++){
192  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
193  if(iomodel->domaintype!=Domain2DhorizontalEnum && !element->IsOnBase()) continue;
194  for(int kk=0;kk<K;kk++){
195  element->DatasetInputAdd(BasalforcingsIsmip6TfEnum,array3d[kk],inputs2,iomodel,Ms[kk],Ns[kk],1,BasalforcingsIsmip6TfEnum,7,kk);
196  }
197  }
198  xDelete<int>(Ms); xDelete<int>(Ns);
199  for(int i=0;i<K;i++) xDelete<IssmDouble>(array3d[i]);
200  xDelete<IssmDouble*>(array3d);
201  }
202  break;
204  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
205  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
206  break;
207  default:
208  _error_("Basal forcing model "<<EnumToStringx(basalforcing_model)<<" not supported yet");
209  }
210 
211  if(!issmb){
212  iomodel->FetchDataToInput(inputs2,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
213  }
214  if(stabilization==3){
215  iomodel->FetchDataToInput(inputs2,elements,"md.masstransport.spcthickness",MasstransportSpcthicknessEnum); //for DG, we need the spc in the element
216  }
217  if(stabilization==4){
218  iomodel->FetchDataToInput(inputs2,elements,"md.masstransport.spcthickness",MasstransportSpcthicknessEnum); //for FCT, we need the spc in the element (penlaties)
219  }
220 
221  if(iomodel->domaintype!=Domain2DhorizontalEnum){
222  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
223  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
224  }
225 
226 }/*}}}*/
227 void MasstransportAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
228 
229  int numoutputs;
230  char** requestedoutputs = NULL;
231 
232  parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isFS",FlowequationIsFSEnum));
233  parameters->AddObject(iomodel->CopyConstantObject("md.masstransport.isfreesurface",MasstransportIsfreesurfaceEnum));
234  parameters->AddObject(iomodel->CopyConstantObject("md.masstransport.hydrostatic_adjustment",MasstransportHydrostaticAdjustmentEnum));
235  parameters->AddObject(iomodel->CopyConstantObject("md.masstransport.stabilization",MasstransportStabilizationEnum));
236  parameters->AddObject(iomodel->CopyConstantObject("md.masstransport.min_thickness",MasstransportMinThicknessEnum));
237  parameters->AddObject(iomodel->CopyConstantObject("md.masstransport.penalty_factor",MasstransportPenaltyFactorEnum));
238 
239  iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.masstransport.requested_outputs");
240  parameters->AddObject(new IntParam(MasstransportNumRequestedOutputsEnum,numoutputs));
241  if(numoutputs)parameters->AddObject(new StringArrayParam(MasstransportRequestedOutputsEnum,requestedoutputs,numoutputs));
242  iomodel->DeleteData(&requestedoutputs,numoutputs,"md.masstransport.requested_outputs");
243 
244 }/*}}}*/
245 
246 /*Finite Element Analysis*/
248  _error_("not implemented");
249 }/*}}}*/
251  /*Default, return NULL*/
252  return NULL;
253 }/*}}}*/
255 _error_("Not implemented");
256 }/*}}}*/
258 
259  /* Check if ice in element */
260  if(!element->IsIceInElement()) return NULL;
261 
262  if(!element->IsOnBase()) return NULL;
263  Element* basalelement = element->SpawnBasalElement();
264 
265  ElementMatrix* Ke = NULL;
266  switch(element->FiniteElement()){
267  case P1Enum: case P2Enum:
268  Ke = CreateKMatrixCG(basalelement);
269  break;
270  case P0DGEnum:
271  case P1DGEnum:
272  Ke = CreateKMatrixDG(basalelement);
273  break;
274  default:
275  _error_("Element type " << EnumToStringx(element->FiniteElement()) << " not supported yet");
276  }
277 
278  int domaintype;
279  element->FindParam(&domaintype,DomainTypeEnum);
280  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
281  return Ke;
282 }/*}}}*/
284 
285  /* Check if ice in element */
286  if(!element->IsIceInElement()) return NULL;
287 
288  /*Intermediaries */
289  int stabilization;
290  int domaintype,dim;
291  IssmDouble Jdet,D_scalar,dt,h;
292  IssmDouble vel,vx,vy,dvxdx,dvydy;
293  IssmDouble xi,tau;
294  IssmDouble dvx[2],dvy[2];
295  IssmDouble D[4];
296  IssmDouble* xyz_list = NULL;
297 
298  /*Get problem dimension*/
299  element->FindParam(&domaintype,DomainTypeEnum);
300  switch(domaintype){
301  case Domain2DverticalEnum: dim = 1; break;
302  case Domain2DhorizontalEnum: dim = 2; break;
303  case Domain3DEnum: dim = 2; break;
304  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
305  }
306 
307  /*Fetch number of nodes and dof for this finite element*/
308  int numnodes = element->GetNumberOfNodes();
309 
310  /*Initialize Element vector and other vectors*/
311  ElementMatrix* Ke = element->NewElementMatrix();
312  IssmDouble* basis = xNew<IssmDouble>(numnodes);
313  IssmDouble* dbasis = xNew<IssmDouble>(dim*numnodes);
314 
315  /*Retrieve all inputs and parameters*/
316  element->GetVerticesCoordinates(&xyz_list);
317  element->FindParam(&dt,TimesteppingTimeStepEnum);
318  element->FindParam(&domaintype,DomainTypeEnum);
319  element->FindParam(&stabilization,MasstransportStabilizationEnum);
320  Input2* vxaverage_input=element->GetInput2(VxAverageEnum); _assert_(vxaverage_input);
321  Input2* vyaverage_input=NULL;
322  if(dim==2){
323  vyaverage_input=element->GetInput2(VyAverageEnum); _assert_(vyaverage_input);
324  }
325 
326  h = element->CharacteristicLength();
327 
328  /* Start looping on the number of gaussian points: */
329  Gauss* gauss=element->NewGauss(2);
330  for(int ig=gauss->begin();ig<gauss->end();ig++){
331  gauss->GaussPoint(ig);
332 
333  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
334  element->NodalFunctions(basis,gauss);
335  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
336 
337  /*Transient term*/
338  D_scalar=gauss->weight*Jdet;
339  for(int i=0;i<numnodes;i++){
340  for(int j=0;j<numnodes;j++){
341  Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j];
342  }
343  }
344 
345  /*Advection terms*/
346  vxaverage_input->GetInputValue(&vx,gauss);
347  vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
348  D_scalar=dt*gauss->weight*Jdet;
349  if(dim==2){
350  vyaverage_input->GetInputValue(&vy,gauss);
351  vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
352  dvxdx=dvx[0];
353  dvydy=dvy[1];
354  for(int i=0;i<numnodes;i++){
355  for(int j=0;j<numnodes;j++){
356  /*\phi_i \phi_j \nabla\cdot v*/
357  Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy);
358  /*\phi_i v\cdot\nabla\phi_j*/
359  Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
360  }
361  }
362  }
363  else{
364  dvxdx=dvx[0];
365  for(int i=0;i<numnodes;i++){
366  for(int j=0;j<numnodes;j++){
367  Ke->values[i*numnodes+j] += D_scalar*dvxdx*dbasis[0*numnodes+i]*dbasis[0*numnodes+j];
368  Ke->values[i*numnodes+j] += D_scalar*vx*dbasis[0*numnodes+j]*basis[i];
369  }
370  }
371  }
372 
373  for(int i=0;i<4;i++) D[i]=0.;
374  switch(stabilization){
375  case 0:
376  /*Nothing to be done*/
377  break;
378  case 1:
379  /*SSA*/
380  vxaverage_input->GetInputAverage(&vx);
381  if(dim==2) vyaverage_input->GetInputAverage(&vy);
382  D[0*dim+0]=h/2.0*fabs(vx);
383  if(dim==2) D[1*dim+1]=h/2.0*fabs(vy);
384  break;
385  case 2:
386  /*Streamline upwinding*/
387  vxaverage_input->GetInputAverage(&vx);
388  if(dim==1){
389  vel=fabs(vx)+1.e-8;
390  }
391  else{
392  vyaverage_input->GetInputAverage(&vy);
393  vel=sqrt(vx*vx+vy*vy)+1.e-8;
394  }
395  tau=h/(2*vel);
396  break;
397  case 5:
398  /*SUPG*/
399  if(dim!=2) _error_("Stabilization "<<stabilization<<" not supported yet for dim != 2");
400  vxaverage_input->GetInputAverage(&vx);
401  vyaverage_input->GetInputAverage(&vy);
402  vel=sqrt(vx*vx+vy*vy)+1.e-8;
403  //xi=0.3130;
404  xi=1;
405  tau=xi*h/(2*vel);
406  break;
407  default:
408  _error_("Stabilization "<<stabilization<<" not supported yet");
409  }
410  if(stabilization==1){
411  /*SSA*/
412  if(dim==1) D[0]=D_scalar*D[0];
413  else{
414  D[0*dim+0]=D_scalar*D[0*dim+0];
415  D[1*dim+0]=D_scalar*D[1*dim+0];
416  D[0*dim+1]=D_scalar*D[0*dim+1];
417  D[1*dim+1]=D_scalar*D[1*dim+1];
418  }
419 
420  if(dim==2){
421  for(int i=0;i<numnodes;i++){
422  for(int j=0;j<numnodes;j++){
423  Ke->values[i*numnodes+j] += (
424  dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
425  dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j])
426  );
427  }
428  }
429  }
430  else{
431  for(int i=0;i<numnodes;i++){
432  for(int j=0;j<numnodes;j++){
433  Ke->values[i*numnodes+j] += D_scalar*h/(2.*vel)*dbasis[0*numnodes+i] *D[0]*dbasis[0*numnodes+j];
434  }
435  }
436  }
437  }
438  if(stabilization==2){
439  /*Streamline upwind*/
440  _assert_(dim==2);
441  for(int i=0;i<numnodes;i++){
442  for(int j=0;j<numnodes;j++){
443  Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i])*(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j]);
444  }
445  }
446  }
447  if(stabilization==5){/*{{{*/
448  /*Mass matrix - part 2*/
449  for(int i=0;i<numnodes;i++){
450  for(int j=0;j<numnodes;j++){
451  Ke->values[i*numnodes+j]+=gauss->weight*Jdet*tau*basis[j]*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
452  }
453  }
454  /*Mass matrix - part 3*/
455  for(int i=0;i<numnodes;i++){
456  for(int j=0;j<numnodes;j++){
457  Ke->values[i*numnodes+j]+=gauss->weight*Jdet*tau*basis[j]*(basis[i]*dvxdx+basis[i]*dvydy);
458  }
459  }
460 
461  /*Advection matrix - part 2, A*/
462  for(int i=0;i<numnodes;i++){
463  for(int j=0;j<numnodes;j++){
464  Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j])*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
465  }
466  }
467  /*Advection matrix - part 3, A*/
468  for(int i=0;i<numnodes;i++){
469  for(int j=0;j<numnodes;j++){
470  Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j])*(basis[i]*dvxdx+basis[i]*dvydy);
471  }
472  }
473 
474  /*Advection matrix - part 2, B*/
475  for(int i=0;i<numnodes;i++){
476  for(int j=0;j<numnodes;j++){
477  Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx+basis[j]*dvydy)*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
478  }
479  }
480  /*Advection matrix - part 3, B*/
481  for(int i=0;i<numnodes;i++){
482  for(int j=0;j<numnodes;j++){
483  Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx+basis[j]*dvydy)*(basis[i]*dvxdx+basis[i]*dvydy);
484  }
485  }
486  }/*}}}*/
487  }
488 
489  /*Clean up and return*/
490  xDelete<IssmDouble>(xyz_list);
491  xDelete<IssmDouble>(basis);
492  xDelete<IssmDouble>(dbasis);
493  delete gauss;
494  return Ke;
495 }/*}}}*/
497 
498  /* Check if ice in element */
499  if(!element->IsIceInElement()) return NULL;
500 
501  /*Intermediaries */
502  int domaintype;
503  IssmDouble Jdet,D_scalar,dt,vx,vy;
504  IssmDouble* xyz_list = NULL;
505 
506  /*Fetch number of nodes and dof for this finite element*/
507  int numnodes = element->GetNumberOfNodes();
508 
509  /*Initialize Element vector and other vectors*/
510  ElementMatrix* Ke = element->NewElementMatrix();
511  IssmDouble* basis = xNew<IssmDouble>(numnodes);
512  IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
513 
514  /*Retrieve all inputs and parameters*/
515  element->GetVerticesCoordinates(&xyz_list);
516  element->FindParam(&dt,TimesteppingTimeStepEnum);
517  element->FindParam(&domaintype,DomainTypeEnum);
518  Input2* vxaverage_input=element->GetInput2(VxAverageEnum); _assert_(vxaverage_input);
519  Input2* vyaverage_input=element->GetInput2(VyAverageEnum); _assert_(vyaverage_input);
520 
521  /* Start looping on the number of gaussian points: */
522  Gauss* gauss=element->NewGauss(2);
523  for(int ig=gauss->begin();ig<gauss->end();ig++){
524  gauss->GaussPoint(ig);
525 
526  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
527  element->NodalFunctions(basis,gauss);
528  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
529 
530  vxaverage_input->GetInputValue(&vx,gauss);
531  vyaverage_input->GetInputValue(&vy,gauss);
532 
533  D_scalar=gauss->weight*Jdet;
534  for(int i=0;i<numnodes;i++){
535  for(int j=0;j<numnodes;j++){
536  Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j];
537  }
538  }
539 
540  /*WARNING: basis and dbasis are inverted compared to CG*/
541  D_scalar = - dt*gauss->weight*Jdet;
542  for(int i=0;i<numnodes;i++){
543  for(int j=0;j<numnodes;j++){
544  Ke->values[i*numnodes+j] += D_scalar*(vx*dbasis[0*numnodes+i]*basis[j] + vy*dbasis[1*numnodes+i]*basis[j]);
545  }
546  }
547  }
548 
549  /*Clean up and return*/
550  xDelete<IssmDouble>(xyz_list);
551  xDelete<IssmDouble>(basis);
552  xDelete<IssmDouble>(dbasis);
553  delete gauss;
554  return Ke;
555 }/*}}}*/
557 
558  /* Check if ice in element */
559  if(!element->IsIceInElement()) return NULL;
560 
561  if(!element->IsOnBase()) return NULL;
562  Element* basalelement = element->SpawnBasalElement();
563 
564  ElementVector* pe = NULL;
565  switch(element->FiniteElement()){
566  case P1Enum: case P2Enum:
567  pe = CreatePVectorCG(basalelement);
568  break;
569  case P0DGEnum:
570  case P1DGEnum:
571  pe = CreatePVectorDG(basalelement);
572  break;
573  default:
574  _error_("Element type " << EnumToStringx(element->FiniteElement()) << " not supported yet");
575  }
576 
577  int domaintype;
578  element->FindParam(&domaintype,DomainTypeEnum);
579  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
580  return pe;
581 }/*}}}*/
583 
584  /* Check if ice in element */
585  if(!element->IsIceInElement()) return NULL;
586 
587  /*Intermediaries */
588  int stabilization,dim,domaintype;
589  int melt_style,point1;
590  bool mainlyfloating;
591  IssmDouble fraction1,fraction2;
592  IssmDouble Jdet,dt;
593  IssmDouble ms,mb,gmb,fmb,thickness;
594  IssmDouble vx,vy,vel,dvxdx,dvydy,xi,h,tau;
595  IssmDouble dvx[2],dvy[2];
596  IssmDouble gllevelset,phi=1.;
597  IssmDouble* xyz_list = NULL;
598  Gauss* gauss = NULL;
599 
600  /*Get problem dimension*/
601  element->FindParam(&domaintype,DomainTypeEnum);
602  switch(domaintype){
603  case Domain2DverticalEnum: dim = 1; break;
604  case Domain2DhorizontalEnum: dim = 2; break;
605  case Domain3DEnum: dim = 2; break;
606  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
607  }
608 
609  /*Fetch number of nodes and dof for this finite element*/
610  int numnodes = element->GetNumberOfNodes();
611 
612  /*Initialize Element vector and other vectors*/
613  ElementVector* pe = element->NewElementVector();
614  IssmDouble* basis = xNew<IssmDouble>(numnodes);
615  IssmDouble* dbasis= xNew<IssmDouble>(dim*numnodes);
616 
617  /*Retrieve all inputs and parameters*/
618  element->GetVerticesCoordinates(&xyz_list);
619  element->FindParam(&melt_style,GroundinglineMeltInterpolationEnum);
620  element->FindParam(&dt,TimesteppingTimeStepEnum);
621  element->FindParam(&stabilization,MasstransportStabilizationEnum);
622  Input2* gmb_input = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input);
623  Input2* fmb_input = element->GetInput2(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input);
624  Input2* gllevelset_input = element->GetInput2(MaskOceanLevelsetEnum); _assert_(gllevelset_input);
625  Input2* ms_input = element->GetInput2(SmbMassBalanceEnum); _assert_(ms_input);
626  Input2* thickness_input = element->GetInput2(ThicknessEnum); _assert_(thickness_input);
627  Input2* vxaverage_input = element->GetInput2(VxAverageEnum); _assert_(vxaverage_input);
628  Input2* vyaverage_input = element->GetInput2(VyAverageEnum); _assert_(vyaverage_input);
629 
630 // if(element->Id()==9){
631 // gmb_input->Echo();
632 // _error_("S");
633 // }
634 
635  h=element->CharacteristicLength();
636 
637  /*Recover portion of element that is grounded*/
638  phi=element->GetGroundedPortion(xyz_list);
639  if(melt_style==SubelementMelt2Enum){
640  element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
641  gauss = element->NewGauss(point1,fraction1,fraction2,3);
642  }
643  else{
644  gauss = element->NewGauss(3);
645  }
646 
647  /* Start looping on the number of gaussian points: */
648  for(int ig=gauss->begin();ig<gauss->end();ig++){
649  gauss->GaussPoint(ig);
650 
651  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
652  element->NodalFunctions(basis,gauss);
653 
654  ms_input->GetInputValue(&ms,gauss);
655  gmb_input->GetInputValue(&gmb,gauss);
656  fmb_input->GetInputValue(&fmb,gauss);
657  gllevelset_input->GetInputValue(&gllevelset,gauss);
658  thickness_input->GetInputValue(&thickness,gauss);
659 
660  if(melt_style==SubelementMelt1Enum){
661  if (phi>0.999999999) mb=gmb;
662  else mb=(1-phi)*fmb+phi*gmb; // phi is the fraction of grounded ice so (1-phi) is floating
663  }
664  else if(melt_style==SubelementMelt2Enum){
665  if(gllevelset>0.) mb=gmb;
666  else mb=fmb;
667  }
668  else if(melt_style==NoMeltOnPartiallyFloatingEnum){
669  if (phi<0.00000001) mb=fmb;
670  else mb=gmb;
671  }
672  else if(melt_style==FullMeltOnPartiallyFloatingEnum){
673  if (phi<0.99999999) mb=fmb;
674  else mb=gmb;
675  }
676  else _error_("melt interpolation "<<EnumToStringx(melt_style)<<" not implemented yet");
677 
678  for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i];
679 
680  if(stabilization==5){ //SUPG
681  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
682  vxaverage_input->GetInputAverage(&vx);
683  vyaverage_input->GetInputAverage(&vy);
684  vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
685  vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
686  vel=sqrt(vx*vx+vy*vy)+1.e-8;
687  dvxdx=dvx[0];
688  dvydy=dvy[1];
689  //xi=0.3130;
690  xi=1;
691  tau=xi*h/(2*vel);
692 
693  /*Force vector - part 2*/
694  for(int i=0;i<numnodes;i++){
695  pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*(tau*vx*dbasis[0*numnodes+i]+tau*vy*dbasis[1*numnodes+i]);
696  }
697  /*Force vector - part 3*/
698  for(int i=0;i<numnodes;i++){
699  pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*(tau*basis[i]*dvxdx+tau*basis[i]*dvydy);
700  }
701  }
702 
703  }
704 
705  /*Clean up and return*/
706  xDelete<IssmDouble>(xyz_list);
707  xDelete<IssmDouble>(basis);
708  xDelete<IssmDouble>(dbasis);
709  delete gauss;
710  return pe;
711 }/*}}}*/
713 
714  /* Check if ice in element */
715  if(!element->IsIceInElement()) return NULL;
716 
717  /*Intermediaries */
718  int melt_style, point1;
719  bool mainlyfloating;
720  IssmDouble fraction1,fraction2,gllevelset;
721  IssmDouble Jdet,dt;
722  IssmDouble ms,mb,gmb,fmb,thickness,phi=1.;
723  IssmDouble* xyz_list = NULL;
724 
725  /*Fetch number of nodes and dof for this finite element*/
726  int numnodes = element->GetNumberOfNodes();
727 
728  /*Initialize Element vector and other vectors*/
729  ElementVector* pe = element->NewElementVector();
730  IssmDouble* basis = xNew<IssmDouble>(numnodes);
731 
732  /*Retrieve all inputs and parameters*/
733  element->GetVerticesCoordinates(&xyz_list);
734  element->FindParam(&dt,TimesteppingTimeStepEnum);
735  element->FindParam(&melt_style,GroundinglineMeltInterpolationEnum);
736  Input2* gmb_input = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input);
737  Input2* fmb_input = element->GetInput2(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input);
738  Input2* ms_input = element->GetInput2(SmbMassBalanceEnum); _assert_(ms_input);
739  Input2* gllevelset_input = element->GetInput2(MaskOceanLevelsetEnum); _assert_(gllevelset_input);
740  Input2* thickness_input = element->GetInput2(ThicknessEnum); _assert_(thickness_input);
741 
742  /*Recover portion of element that is grounded*/
743  Gauss* gauss=NULL;
744  phi=element->GetGroundedPortion(xyz_list);
745  if(melt_style==SubelementMelt2Enum){
746  element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
747  gauss = element->NewGauss(point1,fraction1,fraction2,3);
748  }
749  else{
750  gauss = element->NewGauss(3);
751  }
752 
753  /* Start looping on the number of gaussian points: */
754  for(int ig=gauss->begin();ig<gauss->end();ig++){
755  gauss->GaussPoint(ig);
756 
757  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
758  element->NodalFunctions(basis,gauss);
759 
760  ms_input->GetInputValue(&ms,gauss);
761  gmb_input->GetInputValue(&gmb,gauss);
762  fmb_input->GetInputValue(&fmb,gauss);
763  gllevelset_input->GetInputValue(&gllevelset,gauss);
764  thickness_input->GetInputValue(&thickness,gauss);
765 
766  if(melt_style==SubelementMelt1Enum){
767  if (phi>0.999999999) mb=gmb;
768  else mb=(1-phi)*fmb+phi*gmb; // phi is the fraction of grounded ice so (1-phi) is floating
769  }
770  else if(melt_style==SubelementMelt2Enum){
771  if(gllevelset>0.) mb=gmb;
772  else mb=fmb;
773  }
774  else if(melt_style==NoMeltOnPartiallyFloatingEnum){
775  if (phi<0.00000001) mb=fmb;
776  else mb=gmb;
777  }
778  else if(melt_style==FullMeltOnPartiallyFloatingEnum){
779  if (phi<0.99999999) mb=fmb;
780  else mb=gmb;
781  }
782  else _error_("melt interpolation "<<EnumToStringx(melt_style)<<" not implemented yet");
783 
784  for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i];
785  }
786 
787  /*Clean up and return*/
788  xDelete<IssmDouble>(xyz_list);
789  xDelete<IssmDouble>(basis);
790  delete gauss;
791  return pe;
792 }/*}}}*/
794  element->GetSolutionFromInputsOneDof(solution,ThicknessEnum);
795 }/*}}}*/
796 void MasstransportAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
797  _error_("Not implemented yet");
798 }/*}}}*/
800 
801  /*Only update if on base*/
802  if(!element->IsOnBase()) return;
803 
804  /*Fetch dof list and allocate solution vector*/
805  int *doflist = NULL;
807 
808  int numnodes = element->GetNumberOfNodes();
809  IssmDouble* newthickness = xNew<IssmDouble>(numnodes);
810  IssmDouble* thicknessresidual = xNew<IssmDouble>(numnodes);
811 
812  /*Use the dof list to index into the solution vector: */
813  IssmDouble minthickness = element->FindParam(MasstransportMinThicknessEnum);
814  for(int i=0;i<numnodes;i++){
815  newthickness[i]=solution[doflist[i]];
816  thicknessresidual[i]=0.;
817  /*Check solution*/
818  if(xIsNan<IssmDouble>(newthickness[i])) _error_("NaN found in solution vector");
819  if(xIsInf<IssmDouble>(newthickness[i])) _error_("Inf found in solution vector");
820  if(newthickness[i]<minthickness){
821  thicknessresidual[i]=minthickness-newthickness[i];
822  newthickness[i]=minthickness;
823  }
824  }
825  element->AddBasalInput2(ThicknessEnum,newthickness,element->GetElementType());
826  element->AddBasalInput2(ThicknessResidualEnum,thicknessresidual,element->GetElementType());
827 
828  xDelete<int>(doflist);
829  xDelete<IssmDouble>(newthickness);
830  xDelete<IssmDouble>(thicknessresidual);
831 
832  /*Update bed and surface accordingly*/
833 
834  /*Get basal element*/
835  int domaintype; element->FindParam(&domaintype,DomainTypeEnum);
836  Element* basalelement=element;
837  if(domaintype!=Domain2DhorizontalEnum) basalelement = element->SpawnBasalElement();
838 
839  /*Fetch number of nodes and dof for this finite element*/
840  int numvertices = basalelement->GetNumberOfVertices();
841 
842  /*Now, we need to do some "processing"*/
843  newthickness = xNew<IssmDouble>(numvertices);
844  IssmDouble* oldthickness = xNew<IssmDouble>(numvertices);
845  IssmDouble* cumdeltathickness = xNew<IssmDouble>(numvertices);
846  IssmDouble* deltathickness = xNew<IssmDouble>(numvertices);
847  IssmDouble* newbase = xNew<IssmDouble>(numvertices);
848  IssmDouble* bed = xNew<IssmDouble>(numvertices);
849  IssmDouble* newsurface = xNew<IssmDouble>(numvertices);
850  IssmDouble* oldbase = xNew<IssmDouble>(numvertices);
851  IssmDouble* oldsurface = xNew<IssmDouble>(numvertices);
852  IssmDouble* phi = xNew<IssmDouble>(numvertices);
853  IssmDouble* sealevel = xNew<IssmDouble>(numvertices);
854 
855  /*Get previous base, thickness, surfac and current sealevel and bed:*/
856  basalelement->GetInputListOnVertices(&newthickness[0],ThicknessEnum);
857  basalelement->GetInputListOnVertices(&oldthickness[0],ThicknessOldEnum);
858  basalelement->GetInputListOnVertices(&oldbase[0],BaseOldEnum);
859  basalelement->GetInputListOnVertices(&oldsurface[0],SurfaceOldEnum);
860  basalelement->GetInputListOnVertices(&phi[0],MaskOceanLevelsetEnum);
861  basalelement->GetInputListOnVertices(&sealevel[0],SealevelEnum);
862  basalelement->GetInputListOnVertices(&cumdeltathickness[0],SealevelriseCumDeltathicknessOldEnum);
863 
864  /*Do we do grounding line migration?*/
865  bool isgroundingline;
866  element->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
867  if(isgroundingline) basalelement->GetInputListOnVertices(&bed[0],BedEnum);
868 
869  /*What is the delta thickness forcing the sea-level rise core: cumulated over time, hence the +=:*/
870  for(int i=0;i<numvertices;i++){
871  cumdeltathickness[i] += newthickness[i]-oldthickness[i];
872  deltathickness[i] = newthickness[i]-oldthickness[i];
873  }
874 
875  /*Find MasstransportHydrostaticAdjustment to figure out how to update the geometry:*/
876  int hydroadjustment;
877  basalelement->FindParam(&hydroadjustment,MasstransportHydrostaticAdjustmentEnum);
878  IssmDouble rho_ice = basalelement->FindParam(MaterialsRhoIceEnum);
879  IssmDouble rho_water = basalelement->FindParam(MaterialsRhoSeawaterEnum);
880 
881  for(int i=0;i<numvertices;i++) {
882  if (phi[i]>0.){ //this is grounded ice: just add thickness to base.
883  if(isgroundingline){
884  newsurface[i] = bed[i]+newthickness[i]; //surface = bed + newthickness
885  newbase[i] = bed[i]; //new base at new bed
886  }
887  else{
888  newsurface[i] = oldbase[i]+newthickness[i]; //surface = oldbase + newthickness
889  newbase[i] = oldbase[i]; //same base: do nothing
890  }
891  }
892  else{ //this is an ice shelf: hydrostatic equilibrium*/
893  if(hydroadjustment==AbsoluteEnum){
894  newsurface[i] = newthickness[i]*(1.-rho_ice/rho_water)+sealevel[i];
895  newbase[i] = newthickness[i]*(-rho_ice/rho_water)+sealevel[i];
896  }
897  else if(hydroadjustment==IncrementalEnum){
898  newsurface[i] = oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i])+sealevel[i]; //surface = oldsurface + (1-di) * dH
899  newbase[i] = oldbase[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i])+sealevel[i]; //base = oldbed + di * dH
900  }
901  else _error_("Hydrostatic adjustment " << hydroadjustment << " (" << EnumToStringx(hydroadjustment) << ") not supported yet");
902  }
903  }
904 
905  /*Add input to the element: */
906  element->AddBasalInput2(SurfaceEnum,newsurface,P1Enum);
907  element->AddBasalInput2(BaseEnum,newbase,P1Enum);
908  element->AddBasalInput2(SealevelriseCumDeltathicknessEnum,cumdeltathickness,P1Enum);
910 
911  /*Free ressources:*/
912  xDelete<IssmDouble>(newthickness);
913  xDelete<IssmDouble>(cumdeltathickness);
914  xDelete<IssmDouble>(newbase);
915  xDelete<IssmDouble>(newsurface);
916  xDelete<IssmDouble>(oldthickness);
917  xDelete<IssmDouble>(deltathickness);
918  xDelete<IssmDouble>(oldbase);
919  xDelete<IssmDouble>(oldsurface);
920  xDelete<IssmDouble>(phi);
921  xDelete<IssmDouble>(sealevel);
922  xDelete<IssmDouble>(bed);
923  xDelete<int>(doflist);
924  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
925 }/*}}}*/
928 }/*}}}*/
929 
930 /*Flux Correction Transport*/
932 
933  /* Check if ice in element */
934  if(!element->IsIceInElement()) return NULL;
935 
936  /*Intermediaries */
937  IssmDouble Jdet,D_scalar;
938  IssmDouble vx,vy,dvxdx,dvydy;
939  IssmDouble* xyz_list = NULL;
940  IssmDouble dvx[2],dvy[2];
941 
942  /*Fetch number of nodes and dof for this finite element*/
943  int numnodes = element->GetNumberOfNodes();
944  int dim = 2;
945 
946  /*Initialize Element vector and other vectors*/
947  ElementMatrix* Ke = element->NewElementMatrix();
948  IssmDouble* basis = xNew<IssmDouble>(numnodes);
949  IssmDouble* dbasis = xNew<IssmDouble>(dim*numnodes);
950  IssmDouble* D = xNewZeroInit<IssmDouble>(dim*dim);
951 
952  /*Retrieve all inputs and parameters*/
953  element->GetVerticesCoordinates(&xyz_list);
954  Input2* vxaverage_input=element->GetInput2(VxEnum); _assert_(vxaverage_input);
955  Input2* vyaverage_input=element->GetInput2(VyEnum); _assert_(vyaverage_input);
956 
957  /* Start looping on the number of gaussian points: */
958  Gauss* gauss=element->NewGauss(2);
959  for(int ig=gauss->begin();ig<gauss->end();ig++){
960  gauss->GaussPoint(ig);
961 
962  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
963  element->NodalFunctions(basis,gauss);
964  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
965  vxaverage_input->GetInputValue(&vx,gauss);
966  vyaverage_input->GetInputValue(&vy,gauss);
967  vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
968  vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
969  dvxdx=dvx[0];
970  dvydy=dvy[1];
971 
972  D_scalar = gauss->weight*Jdet;
973  for(int i=0;i<numnodes;i++){
974  for(int j=0;j<numnodes;j++){
975  /*\phi_i \phi_j \nabla\cdot v*/
976  Ke->values[i*numnodes+j] += -D_scalar*basis[i]*basis[j]*(dvxdx+dvydy);
977  /*\phi_i v\cdot\nabla\phi_j*/
978  Ke->values[i*numnodes+j] += -D_scalar*(vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i]);
979  }
980  }
981  }
982 
983  /*Clean up and return*/
984  xDelete<IssmDouble>(xyz_list);
985  xDelete<IssmDouble>(D);
986  xDelete<IssmDouble>(basis);
987  xDelete<IssmDouble>(dbasis);
988  delete gauss;
989  return Ke;
990 }/*}}}*/
992 
993  /* Check if ice in element */
994  if(!element->IsIceInElement()) return NULL;
995 
996  /*Intermediaries*/
997  IssmDouble D,Jdet;
998  IssmDouble* xyz_list = NULL;
999 
1000  /*Fetch number of nodes and dof for this finite element*/
1001  int numnodes = element->GetNumberOfNodes();
1002 
1003  /*Initialize Element vector and other vectors*/
1004  ElementMatrix* Me = element->NewElementMatrix();
1005  IssmDouble* basis = xNew<IssmDouble>(numnodes);
1006 
1007  /*Retrieve all inputs and parameters*/
1008  element->GetVerticesCoordinates(&xyz_list);
1009 
1010  /* Start looping on the number of gaussian points: */
1011  Gauss* gauss=element->NewGauss(2);
1012  for(int ig=gauss->begin();ig<gauss->end();ig++){
1013  gauss->GaussPoint(ig);
1014 
1015  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
1016  element->NodalFunctions(basis,gauss);
1017 
1018  D=gauss->weight*Jdet;
1019  TripleMultiply(basis,1,numnodes,1,
1020  &D,1,1,0,
1021  basis,1,numnodes,0,
1022  &Me->values[0],1);
1023  }
1024 
1025  /*Clean up and return*/
1026  xDelete<IssmDouble>(xyz_list);
1027  xDelete<IssmDouble>(basis);
1028  delete gauss;
1029  return Me;
1030 }/*}}}*/
1032 
1033  /* Check if ice in element */
1034  if(!element->IsIceInElement()) return NULL;
1035 
1036  /*Intermediaries */
1037  int stabilization,dim,domaintype;
1038  int melt_style,point1;
1039  bool mainlyfloating;
1040  IssmDouble fraction1,fraction2;
1041  IssmDouble Jdet;
1042  IssmDouble ms,mb,gmb,fmb;
1043  IssmDouble vx,vy,dvxdx,dvydy;
1044  IssmDouble dvx[2],dvy[2];
1045  IssmDouble gllevelset,phi=1.;
1046  IssmDouble* xyz_list = NULL;
1047  Gauss* gauss = NULL;
1048 
1049  /*Get problem dimension*/
1050  element->FindParam(&domaintype,DomainTypeEnum);
1051  switch(domaintype){
1052  case Domain2DverticalEnum: dim = 1; break;
1053  case Domain2DhorizontalEnum: dim = 2; break;
1054  case Domain3DEnum: dim = 2; break;
1055  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
1056  }
1057 
1058  /*Fetch number of nodes and dof for this finite element*/
1059  int numnodes = element->GetNumberOfNodes();
1060 
1061  /*Initialize Element vector and other vectors*/
1062  ElementVector* pe = element->NewElementVector();
1063  IssmDouble* basis = xNew<IssmDouble>(numnodes);
1064  IssmDouble* dbasis= xNew<IssmDouble>(dim*numnodes);
1065 
1066  /*Retrieve all inputs and parameters*/
1067  element->GetVerticesCoordinates(&xyz_list);
1068  element->FindParam(&melt_style,GroundinglineMeltInterpolationEnum);
1069  element->FindParam(&stabilization,MasstransportStabilizationEnum);
1070  Input2* gmb_input = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input);
1071  Input2* fmb_input = element->GetInput2(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input);
1072  Input2* gllevelset_input = element->GetInput2(MaskOceanLevelsetEnum); _assert_(gllevelset_input);
1073  Input2* ms_input = element->GetInput2(SmbMassBalanceEnum); _assert_(ms_input);
1074  Input2* vxaverage_input = element->GetInput2(VxAverageEnum); _assert_(vxaverage_input);
1075  Input2* vyaverage_input = element->GetInput2(VyAverageEnum); _assert_(vyaverage_input);
1076 
1077  /*Recover portion of element that is grounded*/
1078  phi=element->GetGroundedPortion(xyz_list);
1079  if(melt_style==SubelementMelt2Enum){
1080  element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
1081  gauss = element->NewGauss(point1,fraction1,fraction2,3);
1082  }
1083  else{
1084  gauss = element->NewGauss(3);
1085  }
1086 
1087  /* Start looping on the number of gaussian points: */
1088  for(int ig=gauss->begin();ig<gauss->end();ig++){
1089  gauss->GaussPoint(ig);
1090 
1091  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
1092  element->NodalFunctions(basis,gauss);
1093 
1094  ms_input->GetInputValue(&ms,gauss);
1095  gmb_input->GetInputValue(&gmb,gauss);
1096  fmb_input->GetInputValue(&fmb,gauss);
1097  gllevelset_input->GetInputValue(&gllevelset,gauss);
1098 
1099  if(melt_style==SubelementMelt1Enum){
1100  if (phi>0.999999999) mb=gmb;
1101  else mb=(1-phi)*fmb+phi*gmb; // phi is the fraction of grounded ice so (1-phi) is floating
1102  }
1103  else if(melt_style==SubelementMelt2Enum){
1104  if(gllevelset>0.) mb=gmb;
1105  else mb=fmb;
1106  }
1107  else if(melt_style==NoMeltOnPartiallyFloatingEnum){
1108  if (phi<0.00000001) mb=fmb;
1109  else mb=gmb;
1110  }
1111  else if(melt_style==FullMeltOnPartiallyFloatingEnum){
1112  if (phi<0.99999999) mb=fmb;
1113  else mb=gmb;
1114  }
1115  else _error_("melt interpolation "<<EnumToStringx(melt_style)<<" not implemented yet");
1116 
1117  for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms-mb)*basis[i];
1118 
1119  }
1120 
1121  /*Clean up and return*/
1122  xDelete<IssmDouble>(xyz_list);
1123  xDelete<IssmDouble>(basis);
1124  xDelete<IssmDouble>(dbasis);
1125  delete gauss;
1126  return pe;
1127 }/*}}}*/
1129 
1130  /*Output*/
1131  Matrix<IssmDouble>* Kff = NULL;
1132  Matrix<IssmDouble>* Kfs = NULL;
1133 
1134  /*Initialize Jacobian Matrix*/
1135  AllocateSystemMatricesx(&Kff,&Kfs,NULL,NULL,femmodel);
1136 
1137  /*Create and assemble matrix*/
1138  for(int i=0;i<femmodel->elements->Size();i++){
1139  Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
1140  ElementMatrix* Ke = this->CreateFctKMatrix(element);
1141  if(Ke) Ke->AddToGlobal(Kff,Kfs);
1142  delete Ke;
1143  }
1144  Kff->Assemble();
1145  Kfs->Assemble();
1146 
1147  /*Assign output pointer*/
1148  *pKff=Kff;
1149  if(pKfs){
1150  *pKfs=Kfs;
1151  }
1152  else{
1153  delete Kfs;
1154  }
1155 }/*}}}*/
1157 
1158  /*Output*/
1159  Vector<IssmDouble>* pf = NULL;
1160 
1161  /*Initialize P vector*/
1162  AllocateSystemMatricesx(NULL,NULL,NULL,&pf,femmodel);
1163 
1164  /*Create and assemble matrix*/
1165  for(int i=0;i<femmodel->elements->Size();i++){
1166  Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
1167  ElementVector* pe = this->CreateFctPVector(element);
1168  if(pe) pe->AddToGlobal(pf);
1169  delete pe;
1170  }
1171  pf->Assemble();
1172 
1173  /*Assign output pointer*/
1174  *ppf=pf;
1175 }/*}}}*/
1177 
1178  /*Initialize Mass matrix*/
1179  Matrix<IssmDouble> *Mff = NULL;
1180  AllocateSystemMatricesx(&Mff,NULL,NULL,NULL,femmodel);
1181 
1182  /*Create and assemble matrix*/
1183  for(int i=0;i<femmodel->elements->Size();i++){
1184  Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
1185  ElementMatrix* MLe = this->CreateMassMatrix(element);
1186  if(MLe){
1187  MLe->AddToGlobal(Mff);
1188  }
1189  delete MLe;
1190  }
1191  Mff->Assemble();
1192 
1193  /*Assign output pointer*/
1194  *pMff=Mff;
1195 }/*}}}*/
DataSet::Size
int Size()
Definition: DataSet.cpp:399
Matrix< IssmDouble >
FINITEELEMENT
#define FINITEELEMENT
Definition: MasstransportAnalysis.cpp:7
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
MasstransportAnalysis::CreateFctKMatrix
ElementMatrix * CreateFctKMatrix(Element *element)
Definition: MasstransportAnalysis.cpp:931
_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
MasstransportAnalysis::GradientJ
void GradientJ(Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
Definition: MasstransportAnalysis.cpp:796
Element::GetNumberOfNodes
virtual int GetNumberOfNodes(void)=0
MasstransportAnalysisEnum
@ MasstransportAnalysisEnum
Definition: EnumDefinitions.h:1163
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
MasstransportAnalysis::UpdateConstraints
void UpdateConstraints(FemModel *femmodel)
Definition: MasstransportAnalysis.cpp:926
MasstransportAnalysis::GetSolutionFromInputs
void GetSolutionFromInputs(Vector< IssmDouble > *solution, Element *element)
Definition: MasstransportAnalysis.cpp:793
SurfaceOldEnum
@ SurfaceOldEnum
Definition: EnumDefinitions.h:824
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
TimesteppingTimeStepEnum
@ TimesteppingTimeStepEnum
Definition: EnumDefinitions.h:433
Constraints
Declaration of Constraints class.
Definition: Constraints.h:13
BedEnum
@ BedEnum
Definition: EnumDefinitions.h:499
Element::SpawnBasalElement
virtual Element * SpawnBasalElement(void)=0
MeshVertexonbaseEnum
@ MeshVertexonbaseEnum
Definition: EnumDefinitions.h:653
Elements
Declaration of Elements class.
Definition: Elements.h:17
Matrix::Assemble
void Assemble(void)
Definition: Matrix.h:185
TripleMultiply
int TripleMultiply(IssmDouble *a, int nrowa, int ncola, int itrna, IssmDouble *b, int nrowb, int ncolb, int itrnb, IssmDouble *c, int nrowc, int ncolc, int itrnc, IssmDouble *d, int iaddd)
Definition: MatrixUtils.cpp:20
MaterialsRhoIceEnum
@ MaterialsRhoIceEnum
Definition: EnumDefinitions.h:264
MismipFloatingMeltRateEnum
@ MismipFloatingMeltRateEnum
Definition: EnumDefinitions.h:1190
IoModel::numberoffaces
int numberoffaces
Definition: IoModel.h:97
P1DGEnum
@ P1DGEnum
Definition: EnumDefinitions.h:1215
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
AllocateSystemMatricesx
void AllocateSystemMatricesx(Matrix< IssmDouble > **pKff, Matrix< IssmDouble > **pKfs, Vector< IssmDouble > **pdf, Vector< IssmDouble > **ppf, FemModel *femmodel)
Definition: AllocateSystemMatricesx.cpp:9
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
BasalforcingsPicoEnum
@ BasalforcingsPicoEnum
Definition: EnumDefinitions.h:990
MasstransportStabilizationEnum
@ MasstransportStabilizationEnum
Definition: EnumDefinitions.h:249
MasstransportRequestedOutputsEnum
@ MasstransportRequestedOutputsEnum
Definition: EnumDefinitions.h:248
BasalforcingsUpperwaterElevationEnum
@ BasalforcingsUpperwaterElevationEnum
Definition: EnumDefinitions.h:94
TransientIsgroundinglineEnum
@ TransientIsgroundinglineEnum
Definition: EnumDefinitions.h:447
SubelementMelt1Enum
@ SubelementMelt1Enum
Definition: EnumDefinitions.h:1295
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
BasalforcingsDeepwaterMeltingRateEnum
@ BasalforcingsDeepwaterMeltingRateEnum
Definition: EnumDefinitions.h:62
Parameters::AddObject
void AddObject(Param *newparam)
Definition: Parameters.cpp:67
IoModel::my_vertices
bool * my_vertices
Definition: IoModel.h:72
MasstransportAnalysis.h
: header file for generic external result object
MasstransportAnalysis::CreateKMatrix
ElementMatrix * CreateKMatrix(Element *element)
Definition: MasstransportAnalysis.cpp:257
MasstransportAnalysis::CreateFctPVector
ElementVector * CreateFctPVector(Element *element)
Definition: MasstransportAnalysis.cpp:1031
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::DeleteMaterials
void DeleteMaterials(void)
Definition: Element.cpp:429
P0DGEnum
@ P0DGEnum
Definition: EnumDefinitions.h:1214
Element
Definition: Element.h:41
Element::NodalFunctions
virtual void NodalFunctions(IssmDouble *basis, Gauss *gauss)=0
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
SealevelriseCumDeltathicknessEnum
@ SealevelriseCumDeltathicknessEnum
Definition: EnumDefinitions.h:688
Element::NewElementVector
ElementVector * NewElementVector(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2505
IoModel::DeleteData
void DeleteData(int num,...)
Definition: IoModel.cpp:500
IoModel::numberofelements
int numberofelements
Definition: IoModel.h:96
Constraints::ActivatePenaltyMethod
void ActivatePenaltyMethod(int in_analysis)
Definition: Constraints.cpp:22
BasalforcingsPerturbationMeltingRateEnum
@ BasalforcingsPerturbationMeltingRateEnum
Definition: EnumDefinitions.h:479
IoModel::CopyConstantObject
Param * CopyConstantObject(const char *constant_name, int param_enum)
Definition: IoModel.cpp:418
MasstransportSpcthicknessEnum
@ MasstransportSpcthicknessEnum
Definition: EnumDefinitions.h:642
MasstransportAnalysis::MassMatrix
void MassMatrix(Matrix< IssmDouble > **pMff, FemModel *femmodel)
Definition: MasstransportAnalysis.cpp:1176
BasalforcingsGroundediceMeltingRateEnum
@ BasalforcingsGroundediceMeltingRateEnum
Definition: EnumDefinitions.h:478
NoneApproximationEnum
@ NoneApproximationEnum
Definition: EnumDefinitions.h:1201
Penpair
Definition: Penpair.h:16
Element::GetGroundedPart
virtual void GetGroundedPart(int *point1, IssmDouble *fraction1, IssmDouble *fraction2, bool *mainlyfloating)=0
DomainTypeEnum
@ DomainTypeEnum
Definition: EnumDefinitions.h:124
VxAverageEnum
@ VxAverageEnum
Definition: EnumDefinitions.h:845
NoMeltOnPartiallyFloatingEnum
@ NoMeltOnPartiallyFloatingEnum
Definition: EnumDefinitions.h:1197
FullMeltOnPartiallyFloatingEnum
@ FullMeltOnPartiallyFloatingEnum
Definition: EnumDefinitions.h:1076
MasstransportAnalysis::UpdateElements
void UpdateElements(Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
Definition: MasstransportAnalysis.cpp:113
Element::NewGauss
virtual Gauss * NewGauss(void)=0
FloatingMeltRateEnum
@ FloatingMeltRateEnum
Definition: EnumDefinitions.h:1069
BasalforcingsOceanSalinityEnum
@ BasalforcingsOceanSalinityEnum
Definition: EnumDefinitions.h:484
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
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
StringArrayParam
Definition: StringArrayParam.h:20
GroundinglineMeltInterpolationEnum
@ GroundinglineMeltInterpolationEnum
Definition: EnumDefinitions.h:160
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
BasalforcingsPicoOverturningCoeffEnum
@ BasalforcingsPicoOverturningCoeffEnum
Definition: EnumDefinitions.h:488
IntParam
Definition: IntParam.h:20
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
InputUpdateFromConstantx
void InputUpdateFromConstantx(FemModel *femmodel, bool constant, int name)
Definition: InputUpdateFromConstantx.cpp:10
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
MaterialsRhoSeawaterEnum
@ MaterialsRhoSeawaterEnum
Definition: EnumDefinitions.h:265
BasalforcingsIsmip6Enum
@ BasalforcingsIsmip6Enum
Definition: EnumDefinitions.h:989
MasstransportAnalysis::CreateConstraints
void CreateConstraints(Constraints *constraints, IoModel *iomodel)
Definition: MasstransportAnalysis.cpp:10
MasstransportAnalysis::CreateKMatrixCG
ElementMatrix * CreateKMatrixCG(Element *element)
Definition: MasstransportAnalysis.cpp:283
MasstransportMinThicknessEnum
@ MasstransportMinThicknessEnum
Definition: EnumDefinitions.h:245
FemModel::elements
Elements * elements
Definition: FemModel.h:44
MasstransportAnalysis::CreateDVector
ElementVector * CreateDVector(Element *element)
Definition: MasstransportAnalysis.cpp:250
Vector::Assemble
void Assemble(void)
Definition: Vector.h:142
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
MasstransportAnalysis::FctKMatrix
void FctKMatrix(Matrix< IssmDouble > **pKff, Matrix< IssmDouble > **pKfs, FemModel *femmodel)
Definition: MasstransportAnalysis.cpp:1128
BasalforcingsIsmip6TfEnum
@ BasalforcingsIsmip6TfEnum
Definition: EnumDefinitions.h:481
Element::AddBasalInput2
virtual void AddBasalInput2(int input_enum, IssmDouble *values, int interpolation_enum)
Definition: Element.h:215
SealevelEnum
@ SealevelEnum
Definition: EnumDefinitions.h:675
MasstransportAnalysis::CreateNodes
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
Definition: MasstransportAnalysis.cpp:91
ThicknessResidualEnum
@ ThicknessResidualEnum
Definition: EnumDefinitions.h:843
BasalforcingsIsmip6BasinIdEnum
@ BasalforcingsIsmip6BasinIdEnum
Definition: EnumDefinitions.h:480
Element::IsIceInElement
bool IsIceInElement()
Definition: Element.cpp:2021
MasstransportAnalysis::InputUpdateFromSolution
void InputUpdateFromSolution(IssmDouble *solution, Element *element)
Definition: MasstransportAnalysis.cpp:799
Loads
Declaration of Loads class.
Definition: Loads.h:16
MasstransportAnalysis::CreatePVectorCG
ElementVector * CreatePVectorCG(Element *element)
Definition: MasstransportAnalysis.cpp:582
MasstransportIsfreesurfaceEnum
@ MasstransportIsfreesurfaceEnum
Definition: EnumDefinitions.h:244
MasstransportAnalysis::DofsPerNode
int DofsPerNode(int **doflist, int domaintype, int approximation)
Definition: MasstransportAnalysis.cpp:110
Element::GetSolutionFromInputsOneDof
void GetSolutionFromInputsOneDof(Vector< IssmDouble > *solution, int solutionenum)
Definition: Element.cpp:1281
MasstransportAnalysis::CreatePVectorDG
ElementVector * CreatePVectorDG(Element *element)
Definition: MasstransportAnalysis.cpp:712
Numericalflux
Definition: Numericalflux.h:18
ThicknessOldEnum
@ ThicknessOldEnum
Definition: EnumDefinitions.h:841
BasalforcingsDeepwaterElevationEnum
@ BasalforcingsDeepwaterElevationEnum
Definition: EnumDefinitions.h:61
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
SetActiveNodesLSMx
void SetActiveNodesLSMx(FemModel *femmodel)
Definition: SetActiveNodesLSMx.cpp:12
FlowequationIsFSEnum
@ FlowequationIsFSEnum
Definition: EnumDefinitions.h:138
IoModel::faces
int * faces
Definition: IoModel.h:88
Gauss::begin
virtual int begin(void)=0
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
Element::CharacteristicLength
virtual IssmDouble CharacteristicLength(void)=0
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
MeshVertexonsurfaceEnum
@ MeshVertexonsurfaceEnum
Definition: EnumDefinitions.h:655
SurfaceloadIceThicknessChangeEnum
@ SurfaceloadIceThicknessChangeEnum
Definition: EnumDefinitions.h:691
Element::GetGroundedPortion
virtual IssmDouble GetGroundedPortion(IssmDouble *xyz_list)=0
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
MasstransportAnalysis::CreatePVector
ElementVector * CreatePVector(Element *element)
Definition: MasstransportAnalysis.cpp:556
MasstransportAnalysis::CreateKMatrixDG
ElementMatrix * CreateKMatrixDG(Element *element)
Definition: MasstransportAnalysis.cpp:496
BasalforcingsIsmip6MeltAnomalyEnum
@ BasalforcingsIsmip6MeltAnomalyEnum
Definition: EnumDefinitions.h:483
VyAverageEnum
@ VyAverageEnum
Definition: EnumDefinitions.h:849
Element::FiniteElement
virtual int FiniteElement(void)=0
ThicknessEnum
@ ThicknessEnum
Definition: EnumDefinitions.h:840
ElementVector::AddToGlobal
void AddToGlobal(Vector< IssmDouble > *pf)
Definition: ElementVector.cpp:155
IoModel::FetchDataToInput
void FetchDataToInput(Inputs2 *inputs2, Elements *elements, const char *vector_name, int input_enum)
Definition: IoModel.cpp:1651
IncrementalEnum
@ IncrementalEnum
Definition: EnumDefinitions.h:1118
P2Enum
@ P2Enum
Definition: EnumDefinitions.h:1223
IoModelToConstraintsx
void IoModelToConstraintsx(Constraints *constraints, IoModel *iomodel, const char *spc_name, int analysis_type, int finite_element, int dof)
Definition: IoModelToConstraintsx.cpp:10
ElementMatrix::AddToGlobal
void AddToGlobal(Matrix< IssmDouble > *Kff, Matrix< IssmDouble > *Kfs)
Definition: ElementMatrix.cpp:271
ElementVector
Definition: ElementVector.h:20
Input2::GetInputValue
virtual void GetInputValue(IssmDouble *pvalue, Gauss *gauss)
Definition: Input2.h:38
Gauss::weight
IssmDouble weight
Definition: Gauss.h:11
SubelementMelt2Enum
@ SubelementMelt2Enum
Definition: EnumDefinitions.h:1296
IoModel
Definition: IoModel.h:48
SealevelriseCumDeltathicknessOldEnum
@ SealevelriseCumDeltathicknessOldEnum
Definition: EnumDefinitions.h:689
CreateFaces
void CreateFaces(IoModel *iomodel)
Definition: CreateFaces.cpp:9
Element::DatasetInputAdd
void DatasetInputAdd(int enum_type, IssmDouble *vector, Inputs2 *inputs2, IoModel *iomodel, int M, int N, int vector_type, int vector_enum, int code, int input_enum)
Definition: Element.cpp:1831
MasstransportNumRequestedOutputsEnum
@ MasstransportNumRequestedOutputsEnum
Definition: EnumDefinitions.h:246
AbsoluteEnum
@ AbsoluteEnum
Definition: EnumDefinitions.h:968
MasstransportAnalysis::UpdateParameters
void UpdateParameters(Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
Definition: MasstransportAnalysis.cpp:227
Element::GetNumberOfVertices
virtual int GetNumberOfVertices(void)=0
MasstransportAnalysis::FctPVector
void FctPVector(Vector< IssmDouble > **ppf, FemModel *femmodel)
Definition: MasstransportAnalysis.cpp:1156
Element::GetInputListOnVertices
void GetInputListOnVertices(IssmDouble *pvalue, int enumtype)
Definition: Element.cpp:1131
IoModel::domaintype
int domaintype
Definition: IoModel.h:78
MasstransportPenaltyFactorEnum
@ MasstransportPenaltyFactorEnum
Definition: EnumDefinitions.h:247
Domain2DverticalEnum
@ Domain2DverticalEnum
Definition: EnumDefinitions.h:535
BaseOldEnum
@ BaseOldEnum
Definition: EnumDefinitions.h:496
ElementMatrix
Definition: ElementMatrix.h:19
Vector< IssmDouble >
Input2::GetInputAverage
virtual void GetInputAverage(IssmDouble *pvalue)
Definition: Input2.h:33
Gauss
Definition: Gauss.h:8
MasstransportAnalysis::CreateJacobianMatrix
ElementMatrix * CreateJacobianMatrix(Element *element)
Definition: MasstransportAnalysis.cpp:254
MasstransportAnalysis::Core
void Core(FemModel *femmodel)
Definition: MasstransportAnalysis.cpp:247
MasstransportAnalysis::CreateLoads
void CreateLoads(Loads *loads, IoModel *iomodel)
Definition: MasstransportAnalysis.cpp:26
Element::NodalFunctionsDerivatives
virtual void NodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
MasstransportAnalysis::CreateMassMatrix
ElementMatrix * CreateMassMatrix(Element *element)
Definition: MasstransportAnalysis.cpp:991
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
MasstransportHydrostaticAdjustmentEnum
@ MasstransportHydrostaticAdjustmentEnum
Definition: EnumDefinitions.h:243
Element::NewElementMatrix
ElementMatrix * NewElementMatrix(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2497
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16