Ice Sheet System Model  4.18
Code documentation
HydrologyDCInefficientAnalysis.cpp
Go to the documentation of this file.
2 #include "../toolkits/toolkits.h"
3 #include "../classes/classes.h"
4 #include "../classes/Node.h"
5 #include "../shared/shared.h"
6 #include "../modules/modules.h"
7 #include "../classes/Inputs2/TransientInput2.h"
8 
9 /*Model processing*/
10 int HydrologyDCInefficientAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
11  return 1;
12 }/*}}}*/
13 void HydrologyDCInefficientAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
14 
15  int hydrology_model;
16  int sedimentlimit_flag;
17  int transfer_flag;
18  int unconfined_flag;
19  int penalty_lock;
20  int hydro_maxiter;
21  int hydroslices;
22  int averaging_method;
23  int numoutputs;
24  bool isefficientlayer;
25  IssmDouble penalty_factor;
26  IssmDouble rel_tol;
27  IssmDouble leakagefactor;
28  IssmDouble sedimentlimit;
29  char** requestedoutputs = NULL;
30 
31  /*retrieve some parameters: */
32  bool issmb;
33  iomodel->FindConstant(&issmb,"md.transient.issmb");
34  iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
35 
36  /*Now, do we really want DC?*/
37  if(hydrology_model!=HydrologydcEnum) return;
38 
39  iomodel->FetchData(&sedimentlimit_flag, "md.hydrology.sedimentlimit_flag" );
40  iomodel->FetchData(&transfer_flag, "md.hydrology.transfer_flag" );
41  iomodel->FetchData(&unconfined_flag, "md.hydrology.unconfined_flag" );
42  iomodel->FetchData(&penalty_lock, "md.hydrology.penalty_lock" );
43  iomodel->FetchData(&hydro_maxiter, "md.hydrology.max_iter" );
44  iomodel->FetchData(&hydroslices, "md.hydrology.steps_per_step");
45  iomodel->FetchData(&averaging_method, "md.hydrology.averaging");
46  iomodel->FetchData(&isefficientlayer, "md.hydrology.isefficientlayer");
47  iomodel->FetchData(&penalty_factor, "md.hydrology.penalty_factor" );
48  iomodel->FetchData(&rel_tol, "md.hydrology.rel_tol" );
49 
50  parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
51  parameters->AddObject(new IntParam(HydrologydcSedimentlimitFlagEnum,sedimentlimit_flag));
52  parameters->AddObject(new IntParam(HydrologydcTransferFlagEnum,transfer_flag));
53  parameters->AddObject(new IntParam(HydrologydcUnconfinedFlagEnum,unconfined_flag));
54  parameters->AddObject(new IntParam(HydrologydcPenaltyLockEnum,penalty_lock));
55  parameters->AddObject(new IntParam(HydrologydcMaxIterEnum,hydro_maxiter));
56  parameters->AddObject(new IntParam(HydrologyStepsPerStepEnum,hydroslices));
57  parameters->AddObject(new IntParam(HydrologyAveragingEnum,averaging_method));
58 
59  parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer));
60  parameters->AddObject(new DoubleParam(HydrologydcPenaltyFactorEnum,penalty_factor));
61  parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol));
62  if(transfer_flag==1){
63  iomodel->FetchData(&leakagefactor,"md.hydrology.leakage_factor");
64  parameters->AddObject(new DoubleParam(HydrologydcLeakageFactorEnum,leakagefactor));
65  }
66  if(sedimentlimit_flag==1){
67  iomodel->FetchData(&sedimentlimit,"md.hydrology.sedimentlimit");
68  parameters->AddObject(new DoubleParam(HydrologydcSedimentlimitEnum,sedimentlimit));
69  }
70  if(!issmb){
71  parameters->AddObject(iomodel->CopyConstantObject("md.smb.model",SmbEnum));
72  parameters->AddObject(iomodel->CopyConstantObject("md.smb.averaging",SmbAveragingEnum));
73  }
74 
75  /*Requested outputs*/
76  iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.hydrology.requested_outputs");
77  parameters->AddObject(new IntParam(HydrologyNumRequestedOutputsEnum,numoutputs));
78  if(numoutputs)parameters->AddObject(new StringArrayParam(HydrologyRequestedOutputsEnum,requestedoutputs,numoutputs));
79  iomodel->DeleteData(&requestedoutputs,numoutputs,"md.hydrology.requested_outputs");
80 }/*}}}*/
81 void HydrologyDCInefficientAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
82 
83  bool isefficientlayer;
84  int hydrology_model;
85 
86  /*Fetch data needed: */
87  iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
88 
89  /*Now, do we really want DC?*/
90  if(hydrology_model!=HydrologydcEnum) return;
91 
92  /*Fetch data needed: */
93  iomodel->FindConstant(&isefficientlayer,"md.hydrology.isefficientlayer");
94 
95  /*Update elements: */
96  int counter=0;
97  for(int i=0;i<iomodel->numberofelements;i++){
98  if(iomodel->my_elements[i]){
99  Element* element=(Element*)elements->GetObjectByOffset(counter);
100  element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,P1Enum);
101  counter++;
102  }
103  }
104  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.thickness",ThicknessEnum);
105  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.base",BaseEnum);
106  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
107  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
108  iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.basal_moulin_input",HydrologydcBasalMoulinInputEnum);
109  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.sediment_head",SedimentHeadSubstepEnum);
110  iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.sediment_transmitivity",HydrologydcSedimentTransmitivityEnum);
111  iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.mask_thawed_node",HydrologydcMaskThawedNodeEnum);
112  if(iomodel->domaintype!=Domain2DhorizontalEnum){
113  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
114  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
115  }
116  if(isefficientlayer){
117  iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.mask_eplactive_node",HydrologydcMaskEplactiveNodeEnum);
118  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.epl_head",EplHeadSubstepEnum);
119  }
120 
121 }/*}}}*/
122 void HydrologyDCInefficientAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
123 
124  /*Fetch parameters: */
125  int hydrology_model;
126  iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
127 
128  /*Now, do we really want DC?*/
129  if(hydrology_model!=HydrologydcEnum) return;
130 
131  if(iomodel->domaintype!=Domain2DhorizontalEnum){
132  iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
133  }
135  iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
136 }/*}}}*/
138 
139  /*retrieve some parameters: */
140  int hydrology_model;
141  iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
142  if(hydrology_model!=HydrologydcEnum) return;
143 
144  IoModelToConstraintsx(constraints,iomodel,"md.hydrology.spcsediment_head",HydrologyDCInefficientAnalysisEnum,P1Enum);
145 }/*}}}*/
147 
148  /*Fetch parameters: */
149  int hydrology_model;
150  iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
151  if(hydrology_model!=HydrologydcEnum) return;
152 
153  if(iomodel->domaintype==Domain3DEnum){
154  iomodel->FetchData(1,"md.mesh.vertexonbase");
155  }
156  //create penalties for nodes: no node can have water above the max
158  for(int i=0;i<iomodel->numberofvertices;i++){
159  if (iomodel->domaintype!=Domain3DEnum){
160  /*keep only this partition's nodes:*/
161  if(iomodel->my_vertices[i]){
162  loads->AddObject(new Pengrid(i+1,i,iomodel));
163  loads->AddObject(new Moulin(i+1,i,iomodel));
164  }
165  }
166  else if(reCast<int>(iomodel->Data("md.mesh.vertexonbase")[i])){
167  if(iomodel->my_vertices[i]){
168  loads->AddObject(new Pengrid(i+1,i,iomodel));
169  loads->AddObject(new Moulin(i+1,i,iomodel));
170  }
171  }
172  }
173  iomodel->DeleteData(1,"md.mesh.vertexonbase");
174 }/*}}}*/
175 
176 /*Finite Element Analysis*/
178  _error_("not implemented");
179 }/*}}}*/
181  /*Default, return NULL*/
182  return NULL;
183 }/*}}}*/
185 _error_("Not implemented");
186 }/*}}}*/
188 
189  /*Intermediaries*/
190  bool thawed_element;
191  int domaintype;
192  Element* basalelement;
193 
194  /*Get basal element*/
195  element->FindParam(&domaintype,DomainTypeEnum);
196  switch(domaintype){
198  basalelement = element;
199  break;
200  case Domain3DEnum:
201  if(!element->IsOnBase()) return NULL;
202  basalelement = element->SpawnBasalElement();
203  break;
204  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
205  }
206 
207  basalelement->GetInput2Value(&thawed_element,HydrologydcMaskThawedEltEnum);
208 
209  /*Check that all nodes are active, else return empty matrix*/
210  if(!thawed_element) {
211  if(domaintype!=Domain2DhorizontalEnum){
212  basalelement->DeleteMaterials();
213  delete basalelement;
214  }
215  return NULL;
216  }
217 
218  /*Intermediaries */
219  bool active_element,isefficientlayer;
220  IssmDouble D_scalar,Jdet,dt;
221  IssmDouble sediment_transmitivity;
222  IssmDouble transfer,sediment_storing;
223  IssmDouble *xyz_list = NULL;
224 
225  /*Fetch number of nodes and dof for this finite element*/
226  int numnodes = basalelement->GetNumberOfNodes();
227 
228  /*Initialize Element vector*/
229  ElementMatrix* Ke = basalelement->NewElementMatrix();
230  IssmDouble* B = xNew<IssmDouble>(2*numnodes);
231  IssmDouble* basis = xNew<IssmDouble>(numnodes);
232  IssmDouble D[2][2]= {0.};
233 
234  /*Retrieve all inputs and parameters*/
235  basalelement ->GetVerticesCoordinates(&xyz_list);
236  basalelement ->FindParam(&dt,TimesteppingTimeStepEnum);
237  basalelement ->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
238  Input2* SedTrans_input = basalelement->GetInput2(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input);
239  Input2* sed_head_input = basalelement->GetInput2(SedimentHeadSubstepEnum);
240  Input2* base_input = basalelement->GetInput2(BaseEnum);
241  Input2* old_wh_input = basalelement->GetInput2(SedimentHeadOldEnum); _assert_(old_wh_input);
242 
243  /*Transfer related Inputs*/
244  if(isefficientlayer){
245  basalelement->GetInput2Value(&active_element,HydrologydcMaskEplactiveEltEnum);
246  }
247 
248  /* Start looping on the number of gaussian points: */
249  Gauss* gauss=basalelement->NewGauss(2);
250 
251  for(int ig=gauss -> begin();ig<gauss->end();ig++){
252  gauss -> GaussPoint(ig);
253  basalelement -> JacobianDeterminant(&Jdet,xyz_list,gauss);
254 
255  sediment_transmitivity = SedimentTransmitivity(basalelement,gauss,sed_head_input,base_input,SedTrans_input);
256  sediment_storing = SedimentStoring(basalelement,gauss,sed_head_input,base_input);
257 
258  /*Diffusivity*/
259  D_scalar=sediment_transmitivity*gauss->weight*Jdet;
260  //D_scalar=gauss->weight*Jdet;
261  if(dt!=0.) D_scalar=D_scalar*dt;
262  D[0][0]=D_scalar;
263  D[1][1]=D_scalar;
264  GetB(B,basalelement,xyz_list,gauss);
265  TripleMultiply(B,2,numnodes,1,
266  &D[0][0],2,2,0,
267  B,2,numnodes,0,
268  &Ke->values[0],1);
269 
270  /*Transient*/
271  if(dt!=0.){
272  basalelement->NodalFunctions(&basis[0],gauss);
273  D_scalar=sediment_storing*gauss->weight*Jdet;
274  //D_scalar=(sediment_storing/sediment_transmitivity)*gauss->weight*Jdet;
275  TripleMultiply(basis,numnodes,1,0,
276  &D_scalar,1,1,0,
277  basis,1,numnodes,0,
278  &Ke->values[0],1);
279 
280  /*Transfer EPL part*/
281  if(isefficientlayer){
282  if(active_element){
283  transfer=GetHydrologyKMatrixTransfer(basalelement);
284  basalelement->NodalFunctions(&basis[0],gauss);
285  D_scalar=dt*transfer*gauss->weight*Jdet;
286  //D_scalar=dt*(transfer/sediment_transmitivity)*gauss->weight*Jdet;
287  TripleMultiply(basis,numnodes,1,0,
288  &D_scalar,1,1,0,
289  basis,1,numnodes,0,
290  &Ke->values[0],1);
291  }
292  }
293  }
294  }
295  /*Clean up and return*/
296  xDelete<IssmDouble>(xyz_list);
297  xDelete<IssmDouble>(B);
298  xDelete<IssmDouble>(basis);
299  delete gauss;
300  if(domaintype!=Domain2DhorizontalEnum){
301  basalelement->DeleteMaterials();
302  delete basalelement;
303  }
304  return Ke;
305 }/*}}}*/
307 
308  /*Intermediaries*/
309  bool thawed_element;
310  int domaintype;
311  Element* basalelement;
312 
313  /*Get basal element*/
314  element->FindParam(&domaintype,DomainTypeEnum);
315  switch(domaintype){
317  basalelement = element;
318  break;
319  case Domain3DEnum:
320  if(!element->IsOnBase()) return NULL;
321  basalelement = element->SpawnBasalElement();
322  break;
323  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
324  }
325 
326  basalelement->GetInput2Value(&thawed_element,HydrologydcMaskThawedEltEnum);
327 
328  /*Check that all nodes are active, else return empty matrix*/
329  if(!thawed_element) {
330  if(domaintype!=Domain2DhorizontalEnum){
331  basalelement->DeleteMaterials();
332  delete basalelement;
333  }
334  return NULL;
335  }
336 
337  /*Intermediaries */
338  bool active_element,isefficientlayer;
339  int smb_model;
340  int smbsubstepping;
341  int hydrologysubstepping;
342  int smb_averaging;
343  IssmDouble dt,scalar,sediment_storing;
344  IssmDouble water_head,sediment_transmitivity;
345  IssmDouble water_load,runoff_value,transfer;
346  IssmDouble Jdet;
347 
348  IssmDouble *xyz_list = NULL;
349  Input2 *active_element_input = NULL;
350  Input2 *old_wh_input = NULL;
351  Input2 *dummy_input = NULL;
352  Input2 *surface_runoff_input = NULL;
353 
354  /*Fetch number of nodes and dof for this finite element*/
355  int numnodes = basalelement->GetNumberOfNodes();
356 
357  /*Initialize Element vector*/
358  ElementVector* pe = basalelement->NewElementVector();
359  IssmDouble* basis = xNew<IssmDouble>(numnodes);
360 
361  /*Retrieve all inputs and parameters*/
362  basalelement->GetVerticesCoordinates(&xyz_list);
363  basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
364  basalelement->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
365  basalelement->FindParam(&smb_model,SmbEnum);
366  basalelement->FindParam(&smb_averaging,SmbAveragingEnum);
367 
368  Input2* sed_head_input = basalelement->GetInput2(SedimentHeadSubstepEnum);
369  Input2* epl_head_input = basalelement->GetInput2(EplHeadSubstepEnum);
370  Input2* base_input = basalelement->GetInput2(BaseEnum);
371  Input2* basal_melt_input = basalelement->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(basal_melt_input);
372  Input2* SedTrans_input = basalelement->GetInput2(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input);
373 
374  IssmDouble time;
375  basalelement->FindParam(&time,TimeEnum);
376 
377  if(dt!= 0.){
378  old_wh_input = basalelement->GetInput2(SedimentHeadOldEnum); _assert_(old_wh_input);
379  }
380  if(smb_model==SMBgradientscomponentsEnum){
381  basalelement->FindParam(&smbsubstepping,SmbStepsPerStepEnum);
382  basalelement->FindParam(&hydrologysubstepping,HydrologyStepsPerStepEnum);
383 
384  if(smbsubstepping==1){
385  //no substeping for the smb we take the result from there
386  dummy_input = basalelement->GetInput2(SmbRunoffEnum); _assert_(dummy_input);
387  }
388  else if(smbsubstepping>1 && smbsubstepping<=hydrologysubstepping){
389  //finer hydro stepping, we take the value at the needed time
390  dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum, time); _assert_(dummy_input);
391  }
392  else{
393  //finer stepping in smb, we average the runoff from transient input
394  dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum,time-dt,time,smb_averaging); _assert_(dummy_input);
395  }
396  surface_runoff_input=xDynamicCast<Input2*>(dummy_input); _assert_(surface_runoff_input);
397  }
398 
399  /*Transfer related Inputs*/
400  if(isefficientlayer){
401  basalelement->GetInput2Value(&active_element,HydrologydcMaskEplactiveEltEnum);
402  }
403 
404  /* Start looping on the number of gaussian points: */
405  Gauss* gauss=basalelement->NewGauss(2);
406 
407  IssmDouble yts;
408  basalelement->FindParam(&yts,ConstantsYtsEnum);
409 
410  for(int ig=gauss->begin();ig<gauss->end();ig++){
411  gauss->GaussPoint(ig);
412  basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
413  basalelement->NodalFunctions(basis,gauss);
414  sediment_transmitivity = SedimentTransmitivity(basalelement,gauss,sed_head_input,base_input,SedTrans_input);
415 
416  /*Loading term*/
417  if(!isefficientlayer){
418  basal_melt_input->GetInputValue(&water_load,gauss);
419  if(surface_runoff_input) surface_runoff_input->GetInputValue(&runoff_value,gauss);
420  else runoff_value = 0.;
421  scalar = Jdet*gauss->weight*(water_load+runoff_value);
422  //scalar = Jdet*gauss->weight*(water_load)/sediment_transmitivity;
423  if(dt!=0.) scalar = scalar*dt;
424  for(int i=0;i<numnodes;i++){
425  pe->values[i]+=scalar*basis[i];
426  }
427  }
428  else{
429  /*if EPL is present and active input is there not here*/
430  if(!active_element){
431  basal_melt_input->GetInputValue(&water_load,gauss);
432  if(surface_runoff_input)surface_runoff_input->GetInputValue(&runoff_value,gauss);
433  else runoff_value = 0.;
434  scalar = Jdet*gauss->weight*(water_load+runoff_value);
435  //scalar = Jdet*gauss->weight*(water_load)/sediment_transmitivity;
436  if(dt!=0.) scalar = scalar*dt;
437  for(int i=0;i<numnodes;i++){
438  pe->values[i]+=scalar*basis[i];
439  }
440  }
441  }
442 
443  /*Transient and transfer terms*/
444  if(dt!=0.){
445  old_wh_input->GetInputValue(&water_head,gauss);
446  sediment_storing = SedimentStoring(basalelement,gauss,sed_head_input,base_input);
447  if(isefficientlayer){
448  /*Dealing with the sediment part of the transfer term*/
449  if(active_element){
450  transfer=GetHydrologyPVectorTransfer(basalelement,gauss,epl_head_input);
451  }
452  else{
453  transfer=0.0;
454  }
455  scalar = Jdet*gauss->weight*((water_head*sediment_storing)+(dt*transfer));
456  //scalar = Jdet*gauss->weight*((water_head*sediment_storing)+(dt*transfer))/sediment_transmitivity;
457  for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
458  }
459  else{
460  scalar = Jdet*gauss->weight*(water_head*sediment_storing);
461  //scalar = Jdet*gauss->weight*(water_head*sediment_storing)/sediment_transmitivity;
462  for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
463  }
464  }
465  }
466  /*Clean up and return*/
467  xDelete<IssmDouble>(xyz_list);
468  xDelete<IssmDouble>(basis);
469  delete gauss;
470  if(domaintype!=Domain2DhorizontalEnum){
471  basalelement->DeleteMaterials();
472  delete basalelement;
473  }
474  return pe;
475 }/*}}}*/
476 void HydrologyDCInefficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
477  /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*2.
478  * For node i, Bi can be expressed in the actual coordinate system
479  * by:
480  * Bi=[ dN/dx ]
481  * [ dN/dy ]
482  * where N is the finiteelement function for node i.
483  *
484  * We assume B has been allocated already, of size: 3x(2*numnodes)
485  */
486 
487  /*Fetch number of nodes for this finite element*/
488  int numnodes = element->GetNumberOfNodes();
489 
490  /*Get nodal functions derivatives*/
491  IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
492  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
493 
494  /*Build B: */
495  for(int i=0;i<numnodes;i++){
496  B[numnodes*0+i] = dbasis[0*numnodes+i];
497  B[numnodes*1+i] = dbasis[1*numnodes+i];
498  }
499 
500  /*Clean-up*/
501  xDelete<IssmDouble>(dbasis);
502 }/*}}}*/
505 }/*}}}*/
506 void HydrologyDCInefficientAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
507  _error_("Not implemented yet");
508 }/*}}}*/
510 
511  /*Intermediaries*/
512  int domaintype;
513  Element* basalelement=NULL;
514  bool converged;
515  int* doflist = NULL;
516 
517  /*Get basal element*/
518  element->FindParam(&domaintype,DomainTypeEnum);
519  switch(domaintype){
521  basalelement = element;
522  break;
523  case Domain3DEnum:
524  if(!element->IsOnBase()) return;
525  basalelement = element->SpawnBasalElement();
526  break;
527  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
528  }
529  /*Fetch number of nodes for this finite element*/
530  int numnodes = basalelement->GetNumberOfNodes();
531 
532  /*Fetch dof list and allocate solution vector*/
533  basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
534  IssmDouble* values = xNew<IssmDouble>(numnodes);
535  IssmDouble* pressure = xNew<IssmDouble>(numnodes);
536  IssmDouble* residual = xNew<IssmDouble>(numnodes);
537 
538  for(int i=0;i<numnodes;i++){
539  values[i] =solution[doflist[i]];
540  if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
541  if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
542 
543  }
544 
545  /*If converged keep the residual in mind, also compute effective pressure*/
546  basalelement->GetInputValue(&converged,ConvergedEnum);
547 
548  /*Get inputs*/
549  if(converged){
550  IssmDouble penalty_factor,kmax,kappa,h_max;
551  IssmDouble* thickness = xNew<IssmDouble>(numnodes);
552  IssmDouble* base = xNew<IssmDouble>(numnodes);
553  IssmDouble* transmitivity = xNew<IssmDouble>(numnodes);
554 
555  basalelement->FindParam(&kmax,HydrologySedimentKmaxEnum);
556  basalelement->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
557  IssmDouble rho_freshwater = basalelement->FindParam(MaterialsRhoFreshwaterEnum);
558  IssmDouble rho_ice = basalelement->FindParam(MaterialsRhoIceEnum);
559  IssmDouble g = basalelement->FindParam(ConstantsGEnum);
560 
561  basalelement->GetInputListOnVertices(&thickness[0],ThicknessEnum);
562  basalelement->GetInputListOnVertices(&base[0],BaseEnum);
563  basalelement->GetInputListOnVertices(&transmitivity[0],HydrologydcSedimentTransmitivityEnum);
564 
565  kappa=kmax*pow(10.,penalty_factor);
566 
567  for(int i=0;i<numnodes;i++){
568  GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->GetNode(i));
569  if(values[i]>h_max) {
570  residual[i] = kappa*(values[i]-h_max);
571  //residual[i] = kappa*(values[i]-h_max)*transmitivity[i];
572  }
573  else{
574  residual[i] = 0.;
575  }
576  pressure[i]=(rho_ice*g*thickness[i])-(rho_freshwater*g*(values[i]-base[i]));
577  }
578  xDelete<IssmDouble>(thickness);
579  xDelete<IssmDouble>(base);
580  xDelete<IssmDouble>(transmitivity);
581  }
582 
583  /*Add input to the element: */
587 
588  /*Free ressources:*/
589  xDelete<IssmDouble>(values);
590  xDelete<IssmDouble>(residual);
591  xDelete<IssmDouble>(pressure);
592  xDelete<int>(doflist);
593  if(domaintype!=Domain2DhorizontalEnum){
594  basalelement->DeleteMaterials();
595  delete basalelement;
596  }
597 }/*}}}*/
599  /*Default, do nothing*/
600  return;
601 }/*}}}*/
602 /*Intermediaries*/
603 IssmDouble HydrologyDCInefficientAnalysis::SedimentStoring(Element* element,Gauss* gauss,Input2* sed_head_input, Input2* base_input){/*{{{*/
604  int unconf_scheme;
605  IssmDouble expfac;
606  IssmDouble sediment_storing;
607  IssmDouble storing,yield;
608  IssmDouble base_elev,prestep_head,water_sheet;
609  IssmDouble rho_freshwater = element->FindParam(MaterialsRhoFreshwaterEnum);
610  IssmDouble g = element->FindParam(ConstantsGEnum);
611  IssmDouble sediment_porosity = element->FindParam(HydrologydcSedimentPorosityEnum);
612  IssmDouble sediment_thickness = element->FindParam(HydrologydcSedimentThicknessEnum);
613  IssmDouble sediment_compressibility = element->FindParam(HydrologydcSedimentCompressibilityEnum);
614  IssmDouble water_compressibility = element->FindParam(HydrologydcWaterCompressibilityEnum);
615  element->FindParam(&unconf_scheme,HydrologydcUnconfinedFlagEnum);
616  switch(unconf_scheme){
617  case 0:
618  sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
619  break;
620  case 1:
621  base_input->GetInputValue(&base_elev,gauss);
622  sed_head_input->GetInputValue(&prestep_head,gauss);
623  water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
624 
625  /* if (water_sheet<sediment_thickness){ */
626  /* sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); */
627  /* } */
628  /* else{ */
629  /* sediment_storing=sediment_porosity; */
630  /* } */
631  storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
632  //using logistic function for heavyside approximation
633  expfac=10.;
634  yield=sediment_porosity;
635  sediment_storing=yield+(storing-yield)/(1+exp(-2*expfac*(water_sheet-0.99*sediment_thickness)));
636  break;
637  default:
638  _error_("UnconfinedFlag is 0 or 1");
639  }
640  return sediment_storing;
641 }/*}}}*/
642 IssmDouble HydrologyDCInefficientAnalysis::SedimentTransmitivity(Element* element,Gauss* gauss,Input2* sed_head_input, Input2* base_input,Input2* SedTrans_input){/*{{{*/
643  int unconf_scheme;
644  IssmDouble ratio,expfac;
645  IssmDouble sediment_transmitivity;
646  IssmDouble FullLayer_transmitivity;
647  IssmDouble meltingrate;
648  IssmDouble groundedice;
649  IssmDouble base_elev,prestep_head,water_sheet;
650  IssmDouble sediment_thickness = element->FindParam(HydrologydcSedimentThicknessEnum);
651 
652  element->FindParam(&unconf_scheme,HydrologydcUnconfinedFlagEnum);
653  SedTrans_input->GetInputValue(&FullLayer_transmitivity,gauss);
654 
655  switch(unconf_scheme){
656  case 0:
657  sediment_transmitivity=FullLayer_transmitivity;
658  break;
659  case 1:
660  base_input->GetInputValue(&base_elev,gauss);
661  sed_head_input->GetInputValue(&prestep_head,gauss);
662  water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
663 
664  //min definition of the if test
665  sediment_transmitivity=FullLayer_transmitivity*min(water_sheet/sediment_thickness,1.);
666  if (sediment_transmitivity==0){
667  sediment_transmitivity=1.0e-20;
668  }
669 
670  break;
671  default:
672  _error_("UnconfinedFlag is 0 or 1");
673  }
674  return sediment_transmitivity;
675 }/*}}}*/
677 
678  int hmax_flag;
679  IssmDouble h_max;
680  IssmDouble rho_ice,rho_water;
681  IssmDouble thickness,bed;
682  /*Get the flag to the limitation method*/
683  element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);
684 
685  /*Switch between the different cases*/
686  switch(hmax_flag){
687  case 0:
688  h_max=1.0e+10;
689  break;
690  case 1:
691  element->FindParam(&h_max,HydrologydcSedimentlimitEnum);
692  break;
693  case 2:
694  /*Compute max*/
695  rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
696  rho_ice = element->FindParam(MaterialsRhoIceEnum);
697  element->GetInputValue(&thickness,innode,ThicknessEnum);
698  element->GetInputValue(&bed,innode,BaseEnum);
699  h_max=((rho_ice*thickness)/rho_water)+bed;
700  break;
701  case 3:
702  _error_("Using normal stress not supported yet");
703  break;
704  default:
705  _error_("no case higher than 3 for SedimentlimitFlag");
706  }
707  /*Assign output pointer*/
708  *ph_max=h_max;
709 }
710 /*}}}*/
712 
713  int transfermethod;
714  IssmDouble leakage,transfer;
715  element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
716 
717  /*Switch between the different transfer methods cases*/
718  switch(transfermethod){
719  case 0:
720  /*Just keepping the transfer to zero*/
721  transfer=0.0;
722  break;
723  case 1:
724  element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
725  transfer=+leakage;
726  break;
727  default:
728  _error_("no case higher than 1 for the Transfer method");
729  }
730  return transfer;
731 }/*}}}*/
733 
734  int transfermethod;
735  IssmDouble epl_head;
736  IssmDouble leakage,transfer;
737  element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
738 
739  /*Switch between the different transfer methods cases*/
740  switch(transfermethod){
741  case 0:
742  /*Just keepping the transfer to zero*/
743  transfer=0.0;
744  break;
745  case 1:
746  _assert_(epl_head_input);
747  epl_head_input->GetInputValue(&epl_head,gauss);
748  element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
749  transfer=+epl_head*leakage;
750  break;
751  default:
752  _error_("no case higher than 1 for the Transfer method");
753  }
754  return transfer;
755 }/*}}}*/
757 
758  bool element_active;
759  Element* element=NULL;
760 
761  for(int i=0;i<femmodel->elements->Size();i++){
762  element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
763 
765  if(input->GetInputMax()>0.){
766  element_active = true;
767  }
768  else{
769  element_active = false;
770  }
771  element->SetBoolInput(element->inputs2,HydrologydcMaskEplactiveEltEnum,element_active);
772  }
773 }/*}}}*/
775  bool active_element;
776  int domaintype;
777  Element* basalelement=NULL;
778 
779  /*Get basal element*/
780  element->FindParam(&domaintype,DomainTypeEnum);
781  switch(domaintype){
783  basalelement = element;
784  break;
785  case Domain3DEnum:
786  if(!element->IsOnBase()) return;
787  basalelement = element->SpawnBasalElement();
788  break;
789  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
790  }
791  /*Intermediaries*/
792  int numnodes = basalelement->GetNumberOfNodes();
793  IssmDouble* meltingrate = xNew<IssmDouble>(numnodes);
794  IssmDouble* groundedice = xNew<IssmDouble>(numnodes);
795 
797  basalelement->GetInputListOnVertices(&groundedice[0],MaskOceanLevelsetEnum);
798 
799  /*if melting rate is not positive and node is not floating, deactivate*/
800  for(int i=0;i<numnodes;i++){
801  if ((meltingrate[i]<=0.0) && (groundedice[i]>0)){
802  vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
803  }
804  else{
805  vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
806  }
807  }
808 
809  if(domaintype!=Domain2DhorizontalEnum){
810  basalelement->DeleteMaterials();
811  delete basalelement;
812  }
813  xDelete<IssmDouble>(meltingrate);
814  xDelete<IssmDouble>(groundedice);
815 }/*}}}*/
817 
818  bool element_active;
819  Element* element=NULL;
820 
821  for(int i=0;i<femmodel->elements->Size();i++){
822  element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
823 
824  Input2* input=element->GetInput2(HydrologydcMaskThawedNodeEnum); _assert_(input);
825  if(input->GetInputMax()>0.){
826  element_active = true;
827  }
828  else{
829  element_active = false;
830  }
831  element->SetBoolInput(element->inputs2,HydrologydcMaskThawedEltEnum,element_active);
832  }
833 }/*}}}*/
835  /*Constants*/
836  int domaintype;
837  Element* basalelement=NULL;
838 
839  /*Get basal element*/
840  element->FindParam(&domaintype,DomainTypeEnum);
841  switch(domaintype){
843  basalelement = element;
844  break;
845  case Domain3DEnum:
846  if(!element->IsOnBase()) return;
847  basalelement = element->SpawnBasalElement();
848  break;
849  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
850  }
851 
852  const int numnodes = basalelement->GetNumberOfNodes();
853  IssmDouble flag = 0.;
854  IssmDouble* active = xNew<IssmDouble>(numnodes);
855  bool active_element;
856 
857  /*Pass the activity mask from elements to nodes*/
858  basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskThawedNodeEnum);
859  basalelement->GetInput2Value(&active_element,HydrologydcMaskThawedEltEnum);
860 
861  for(int i=0;i<numnodes;i++) flag+=active[i];
862 
863  /*If any node is active all the node in the element are active*/
864  if(flag>0.){
865  for(int i=0;i<numnodes;i++){
866  active_vec->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
867  }
868  }
869  /*If the element is active all its nodes are active*/
870  else if(active_element){
871  for(int i=0;i<numnodes;i++){
872  active_vec->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
873  }
874  }
875  else{
876  /*Do not do anything: at least one node is active for this element but this element is not solved for*/
877  }
878  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
879  xDelete<IssmDouble>(active);
880 }
881 /*}}}*/
HydrologyDCInefficientAnalysis::UpdateConstraints
void UpdateConstraints(FemModel *femmodel)
Definition: HydrologyDCInefficientAnalysis.cpp:598
DataSet::Size
int Size()
Definition: DataSet.cpp:399
HydrologyDCInefficientAnalysis::CreateKMatrix
ElementMatrix * CreateKMatrix(Element *element)
Definition: HydrologyDCInefficientAnalysis.cpp:187
SedimentHeadOldEnum
@ SedimentHeadOldEnum
Definition: EnumDefinitions.h:699
EffectivePressureSubstepEnum
@ EffectivePressureSubstepEnum
Definition: EnumDefinitions.h:549
BaseEnum
@ BaseEnum
Definition: EnumDefinitions.h:495
SmbEnum
@ SmbEnum
Definition: EnumDefinitions.h:358
_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
HydrologyDCInefficientAnalysis::GetSolutionFromInputs
void GetSolutionFromInputs(Vector< IssmDouble > *solution, Element *element)
Definition: HydrologyDCInefficientAnalysis.cpp:503
HydrologyDCInefficientAnalysis::GradientJ
void GradientJ(Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
Definition: HydrologyDCInefficientAnalysis.cpp:506
Element::GetNumberOfNodes
virtual int GetNumberOfNodes(void)=0
HydrologyDCInefficientAnalysis::UpdateElements
void UpdateElements(Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
Definition: HydrologyDCInefficientAnalysis.cpp:81
SmbAveragingEnum
@ SmbAveragingEnum
Definition: EnumDefinitions.h:349
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
Element::SetBoolInput
void SetBoolInput(Inputs2 *inputs2, int enum_in, bool value)
Definition: Element.cpp:3355
HydrologydcEnum
@ HydrologydcEnum
Definition: EnumDefinitions.h:1106
HydrologyDCInefficientAnalysis::CreateJacobianMatrix
ElementMatrix * CreateJacobianMatrix(Element *element)
Definition: HydrologyDCInefficientAnalysis.cpp:184
DataSet::AddObject
int AddObject(Object *object)
Definition: DataSet.cpp:252
HydrologyRequestedOutputsEnum
@ HydrologyRequestedOutputsEnum
Definition: EnumDefinitions.h:173
HydrologyDCInefficientAnalysis::SedimentStoring
IssmDouble SedimentStoring(Element *element, Gauss *gauss, Input2 *sed_head_input, Input2 *base_input)
Definition: HydrologyDCInefficientAnalysis.cpp:603
HydrologydcMaskEplactiveNodeEnum
@ HydrologydcMaskEplactiveNodeEnum
Definition: EnumDefinitions.h:608
Moulin
Definition: Moulin.h:21
MaskOceanLevelsetEnum
@ MaskOceanLevelsetEnum
Definition: EnumDefinitions.h:640
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
MaskIceLevelsetEnum
@ MaskIceLevelsetEnum
Definition: EnumDefinitions.h:641
TimeEnum
@ TimeEnum
Definition: EnumDefinitions.h:427
MaterialsRhoFreshwaterEnum
@ MaterialsRhoFreshwaterEnum
Definition: EnumDefinitions.h:263
TimesteppingTimeStepEnum
@ TimesteppingTimeStepEnum
Definition: EnumDefinitions.h:433
Gauss::end
virtual int end(void)=0
Constraints
Declaration of Constraints class.
Definition: Constraints.h:13
ConvergedEnum
@ ConvergedEnum
Definition: EnumDefinitions.h:514
Element::SpawnBasalElement
virtual Element * SpawnBasalElement(void)=0
HydrologydcMaskThawedNodeEnum
@ HydrologydcMaskThawedNodeEnum
Definition: EnumDefinitions.h:610
HydrologyDCInefficientAnalysis::CreateLoads
void CreateLoads(Loads *loads, IoModel *iomodel)
Definition: HydrologyDCInefficientAnalysis.cpp:146
MeshVertexonbaseEnum
@ MeshVertexonbaseEnum
Definition: EnumDefinitions.h:653
HydrologydcIsefficientlayerEnum
@ HydrologydcIsefficientlayerEnum
Definition: EnumDefinitions.h:185
Elements
Declaration of Elements class.
Definition: Elements.h:17
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
HydrologydcPenaltyFactorEnum
@ HydrologydcPenaltyFactorEnum
Definition: EnumDefinitions.h:188
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
HydrologydcSedimentPorosityEnum
@ HydrologydcSedimentPorosityEnum
Definition: EnumDefinitions.h:194
HydrologyDCInefficientAnalysis::GetB
void GetB(IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss)
Definition: HydrologyDCInefficientAnalysis.cpp:476
ConstantsYtsEnum
@ ConstantsYtsEnum
Definition: EnumDefinitions.h:104
HydrologydcRelTolEnum
@ HydrologydcRelTolEnum
Definition: EnumDefinitions.h:190
HydrologyDCInefficientAnalysis::GetHydrologyPVectorTransfer
IssmDouble GetHydrologyPVectorTransfer(Element *element, Gauss *gauss, Input2 *epl_head_input)
Definition: HydrologyDCInefficientAnalysis.cpp:732
Parameters::AddObject
void AddObject(Param *newparam)
Definition: Parameters.cpp:67
IoModel::my_vertices
bool * my_vertices
Definition: IoModel.h:72
Input2::GetInputMax
virtual IssmDouble GetInputMax(void)
Definition: Input2.h:34
HydrologyDCInefficientAnalysisEnum
@ HydrologyDCInefficientAnalysisEnum
Definition: EnumDefinitions.h:1099
SmbStepsPerStepEnum
@ SmbStepsPerStepEnum
Definition: EnumDefinitions.h:389
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
Element::DeleteMaterials
void DeleteMaterials(void)
Definition: Element.cpp:429
Element::nodes
Node ** nodes
Definition: Element.h:48
Element
Definition: Element.h:41
Element::NodalFunctions
virtual void NodalFunctions(IssmDouble *basis, Gauss *gauss)=0
SedimentHeadResidualEnum
@ SedimentHeadResidualEnum
Definition: EnumDefinitions.h:702
IoModel::numberofvertices
int numberofvertices
Definition: IoModel.h:99
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
HydrologydcMaskThawedEltEnum
@ HydrologydcMaskThawedEltEnum
Definition: EnumDefinitions.h:609
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
HydrologydcSedimentCompressibilityEnum
@ HydrologydcSedimentCompressibilityEnum
Definition: EnumDefinitions.h:191
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
DoubleParam
Definition: DoubleParam.h:20
IoModel::CopyConstantObject
Param * CopyConstantObject(const char *constant_name, int param_enum)
Definition: IoModel.cpp:418
ConstantsGEnum
@ ConstantsGEnum
Definition: EnumDefinitions.h:102
HydrologySedimentKmaxEnum
@ HydrologySedimentKmaxEnum
Definition: EnumDefinitions.h:174
BasalforcingsGroundediceMeltingRateEnum
@ BasalforcingsGroundediceMeltingRateEnum
Definition: EnumDefinitions.h:478
Element::inputs2
Inputs2 * inputs2
Definition: Element.h:47
NoneApproximationEnum
@ NoneApproximationEnum
Definition: EnumDefinitions.h:1201
DomainTypeEnum
@ DomainTypeEnum
Definition: EnumDefinitions.h:124
HydrologyDCInefficientAnalysis::ElementizeIdsMask
void ElementizeIdsMask(FemModel *femmodel)
Definition: HydrologyDCInefficientAnalysis.cpp:816
Element::NewGauss
virtual Gauss * NewGauss(void)=0
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
Pengrid
Definition: Pengrid.h:21
SMBgradientscomponentsEnum
@ SMBgradientscomponentsEnum
Definition: EnumDefinitions.h:1248
GsetEnum
@ GsetEnum
Definition: EnumDefinitions.h:1093
IoModel::FindConstant
void FindConstant(bool *pvalue, const char *constant_name)
Definition: IoModel.cpp:2362
HydrologydcUnconfinedFlagEnum
@ HydrologydcUnconfinedFlagEnum
Definition: EnumDefinitions.h:197
StringArrayParam
Definition: StringArrayParam.h:20
Element::GetInput2Value
void GetInput2Value(bool *pvalue, int enum_type)
Definition: Element.cpp:1185
HydrologydcMaxIterEnum
@ HydrologydcMaxIterEnum
Definition: EnumDefinitions.h:187
HydrologydcLeakageFactorEnum
@ HydrologydcLeakageFactorEnum
Definition: EnumDefinitions.h:186
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
Element::GetNode
Node * GetNode(int nodeindex)
Definition: Element.cpp:1207
HydrologyDCInefficientAnalysis::SedimentTransmitivity
IssmDouble SedimentTransmitivity(Element *element, Gauss *gauss, Input2 *sed_head_input, Input2 *base_input, Input2 *SedTrans_input)
Definition: HydrologyDCInefficientAnalysis.cpp:642
HydrologyDCInefficientAnalysis::CreatePVector
ElementVector * CreatePVector(Element *element)
Definition: HydrologyDCInefficientAnalysis.cpp:306
IntParam
Definition: IntParam.h:20
INS_VAL
@ INS_VAL
Definition: toolkitsenums.h:14
HydrologyDCInefficientAnalysis::GetHydrologyDCInefficientHmax
void GetHydrologyDCInefficientHmax(IssmDouble *ph_max, Element *element, Node *innode)
Definition: HydrologyDCInefficientAnalysis.cpp:676
HydrologyDCInefficientAnalysis::GetHydrologyKMatrixTransfer
IssmDouble GetHydrologyKMatrixTransfer(Element *element)
Definition: HydrologyDCInefficientAnalysis.cpp:711
HydrologyAveragingEnum
@ HydrologyAveragingEnum
Definition: EnumDefinitions.h:162
HydrologydcSedimentlimitEnum
@ HydrologydcSedimentlimitEnum
Definition: EnumDefinitions.h:192
IoModel::FetchData
void FetchData(bool *pboolean, const char *data_name)
Definition: IoModel.cpp:933
Inputs2
Declaration of Inputs class.
Definition: Inputs2.h:23
HydrologyNumRequestedOutputsEnum
@ HydrologyNumRequestedOutputsEnum
Definition: EnumDefinitions.h:170
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
HydrologyDCInefficientAnalysis::HydrologyIdsGetActive
void HydrologyIdsGetActive(Vector< IssmDouble > *active_vec, Element *element)
Definition: HydrologyDCInefficientAnalysis.cpp:834
CreateSingleNodeToElementConnectivity
void CreateSingleNodeToElementConnectivity(IoModel *iomodel)
Definition: CreateSingleNodeToElementConnectivity.cpp:10
FemModel::elements
Elements * elements
Definition: FemModel.h:44
HydrologydcWaterCompressibilityEnum
@ HydrologydcWaterCompressibilityEnum
Definition: EnumDefinitions.h:198
HydrologydcSedimentTransmitivityEnum
@ HydrologydcSedimentTransmitivityEnum
Definition: EnumDefinitions.h:611
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
HydrologydcPenaltyLockEnum
@ HydrologydcPenaltyLockEnum
Definition: EnumDefinitions.h:189
HydrologydcSedimentlimitFlagEnum
@ HydrologydcSedimentlimitFlagEnum
Definition: EnumDefinitions.h:193
Element::AddBasalInput2
virtual void AddBasalInput2(int input_enum, IssmDouble *values, int interpolation_enum)
Definition: Element.h:215
SedimentHeadSubstepEnum
@ SedimentHeadSubstepEnum
Definition: EnumDefinitions.h:700
SmbRunoffEnum
@ SmbRunoffEnum
Definition: EnumDefinitions.h:772
HydrologyDCInefficientAnalysis::DofsPerNode
int DofsPerNode(int **doflist, int domaintype, int approximation)
Definition: HydrologyDCInefficientAnalysis.cpp:10
IoModel::Data
IssmDouble * Data(const char *data_name)
Definition: IoModel.cpp:437
HydrologyDCInefficientAnalysis.h
: header file for generic external result object
Loads
Declaration of Loads class.
Definition: Loads.h:16
Node
Definition: Node.h:23
HydrologyModelEnum
@ HydrologyModelEnum
Definition: EnumDefinitions.h:169
Element::GetSolutionFromInputsOneDof
void GetSolutionFromInputsOneDof(Vector< IssmDouble > *solution, int solutionenum)
Definition: Element.cpp:1281
Node::Sid
int Sid(void)
Definition: Node.cpp:622
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Gauss::begin
virtual int begin(void)=0
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
MeshVertexonsurfaceEnum
@ MeshVertexonsurfaceEnum
Definition: EnumDefinitions.h:655
HydrologydcSedimentThicknessEnum
@ HydrologydcSedimentThicknessEnum
Definition: EnumDefinitions.h:195
min
IssmDouble min(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:14
HydrologydcBasalMoulinInputEnum
@ HydrologydcBasalMoulinInputEnum
Definition: EnumDefinitions.h:602
HydrologyDCInefficientAnalysis::CreateNodes
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
Definition: HydrologyDCInefficientAnalysis.cpp:122
SmbRunoffTransientEnum
@ SmbRunoffTransientEnum
Definition: EnumDefinitions.h:774
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
Element::GetInputValue
void GetInputValue(bool *pvalue, int enum_type)
Definition: Element.cpp:1177
HydrologyStepsPerStepEnum
@ HydrologyStepsPerStepEnum
Definition: EnumDefinitions.h:175
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
HydrologydcTransferFlagEnum
@ HydrologydcTransferFlagEnum
Definition: EnumDefinitions.h:196
HydrologyDCInefficientAnalysis::Core
void Core(FemModel *femmodel)
Definition: HydrologyDCInefficientAnalysis.cpp:177
IoModelToConstraintsx
void IoModelToConstraintsx(Constraints *constraints, IoModel *iomodel, const char *spc_name, int analysis_type, int finite_element, int dof)
Definition: IoModelToConstraintsx.cpp:10
ElementVector
Definition: ElementVector.h:20
Input2::GetInputValue
virtual void GetInputValue(IssmDouble *pvalue, Gauss *gauss)
Definition: Input2.h:38
Gauss::weight
IssmDouble weight
Definition: Gauss.h:11
IoModel
Definition: IoModel.h:48
HydrologyDCInefficientAnalysis::ElementizeEplMask
void ElementizeEplMask(FemModel *femmodel)
Definition: HydrologyDCInefficientAnalysis.cpp:756
EplHeadSubstepEnum
@ EplHeadSubstepEnum
Definition: EnumDefinitions.h:557
BoolParam
Definition: BoolParam.h:20
max
IssmDouble max(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:24
HydrologyDCInefficientAnalysis::CreateConstraints
void CreateConstraints(Constraints *constraints, IoModel *iomodel)
Definition: HydrologyDCInefficientAnalysis.cpp:137
Element::GetInputListOnVertices
void GetInputListOnVertices(IssmDouble *pvalue, int enumtype)
Definition: Element.cpp:1131
IoModel::domaintype
int domaintype
Definition: IoModel.h:78
ElementMatrix
Definition: ElementMatrix.h:19
Vector< IssmDouble >
HydrologyDCInefficientAnalysis::UpdateParameters
void UpdateParameters(Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
Definition: HydrologyDCInefficientAnalysis.cpp:13
HydrologydcMaskEplactiveEltEnum
@ HydrologydcMaskEplactiveEltEnum
Definition: EnumDefinitions.h:607
Gauss
Definition: Gauss.h:8
Vector::SetValue
void SetValue(int dof, doubletype value, InsMode mode)
Definition: Vector.h:163
HydrologyDCInefficientAnalysis::InputUpdateFromSolution
void InputUpdateFromSolution(IssmDouble *solution, Element *element)
Definition: HydrologyDCInefficientAnalysis.cpp:509
HydrologyDCInefficientAnalysis::CreateDVector
ElementVector * CreateDVector(Element *element)
Definition: HydrologyDCInefficientAnalysis.cpp:180
HydrologyDCInefficientAnalysis::HydrologyIDSGetMask
void HydrologyIDSGetMask(Vector< IssmDouble > *vec_mask, Element *element)
Definition: HydrologyDCInefficientAnalysis.cpp:774
Element::NodalFunctionsDerivatives
virtual void NodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
Element::NewElementMatrix
ElementMatrix * NewElementMatrix(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2497
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16