Ice Sheet System Model  4.18
Code documentation
Balancethickness2Analysis.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 /*Model processing*/
9 
10  int finiteelement = P1Enum;
11  IoModelToConstraintsx(constraints,iomodel,"md.balancethickness.spcthickness",Balancethickness2AnalysisEnum,finiteelement);
12 
13 }/*}}}*/
15 
16 }/*}}}*/
17 void Balancethickness2Analysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
18 
19  int finiteelement = P1Enum;
20  ::CreateNodes(nodes,iomodel,Balancethickness2AnalysisEnum,finiteelement);
21 }/*}}}*/
22 int Balancethickness2Analysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
23  return 1;
24 }/*}}}*/
25 void Balancethickness2Analysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
26 
27  /*Finite element type*/
28  int finiteelement = P1Enum;
29 
30  /*Load variables in element*/
31  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.thickness",ThicknessEnum);
32  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.surface",SurfaceEnum);
33  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.base",BaseEnum);
34  iomodel->FetchDataToInput(inputs2,elements,"md.solidearth.sealevel",SealevelEnum,0);
35  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
36  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vx",VxEnum);
37  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vy",VyEnum);
38  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
39  iomodel->FetchDataToInput(inputs2,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
40  iomodel->FetchDataToInput(inputs2,elements,"md.balancethickness.thickening_rate",BalancethicknessThickeningRateEnum);
41 
42  /*Update elements: */
43  int counter=0;
44  for(int i=0;i<iomodel->numberofelements;i++){
45  if(iomodel->my_elements[i]){
46  Element* element=(Element*)elements->GetObjectByOffset(counter);
47  element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,finiteelement);
48 
49  counter++;
50  }
51  }
52 
53 }/*}}}*/
54 void Balancethickness2Analysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
55 }/*}}}*/
56 
57 /*Finite Element Analysis*/
59  _error_("not implemented yet");
60 }/*}}}*/
62  /*Default, return NULL*/
63  return NULL;
64 }/*}}}*/
66 _error_("Not implemented");
67 }/*}}}*/
69 
70  /*Intermediaries */
71  IssmDouble yts = 365*24*3600.;
72  IssmDouble Jdet,vx,vy,vel;
73  IssmDouble* xyz_list = NULL;
74 
75  /*Fetch number of nodes and dof for this finite element*/
76  int numnodes = element->GetNumberOfNodes();
77 
78  /*Initialize Element vector and other vectors*/
79  ElementMatrix* Ke = element->NewElementMatrix();
80  IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);
81 
82  /*Retrieve all inputs and parameters*/
83  element->GetVerticesCoordinates(&xyz_list);
84  Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
85  Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
86 
87  /*Get element characteristic length*/
88  IssmDouble h = element->CharacteristicLength();
89 
90  /* Start looping on the number of gaussian points: */
91  Gauss* gauss=element->NewGauss(2);
92  for(int ig=gauss->begin();ig<gauss->end();ig++){
93  gauss->GaussPoint(ig);
94  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
95  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
96  vx_input->GetInputValue(&vx,gauss);
97  vy_input->GetInputValue(&vy,gauss);
98 
99  /*make sure are diffusivisty is large enough*/
100  vel = sqrt(vx*vx+vy*vy);
101  if(sqrt(vx*vx+vy*vy)==0.){
102  vx = 0.1/yts;
103  vy = 0.1/yts;
104  vel = sqrt(vx*vx+vy*vy);
105  }
106  else if(vel<30./yts){
107  vx = 0.;//vx/vel*0.01;
108  vy = 0.;//vy/vel*0.01;
109  vel = sqrt(vx*vx+vy*vy);
110  vel = 30./yts*500000.;
111  }
112 
113  for(int i=0;i<numnodes;i++){
114  for(int j=0;j<numnodes;j++){
115  Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
116  (vx*dbasis[0*numnodes+i] + vy*dbasis[1*numnodes+i])*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j])
117  + vel/500000.*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]));
118  }
119  }
120  }
121 
122  /*Clean up and return*/
123  xDelete<IssmDouble>(xyz_list);
124  xDelete<IssmDouble>(dbasis);
125  delete gauss;
126  return Ke;
127 }/*}}}*/
129 
130  return NULL;
131  /*Intermediaries */
132  IssmDouble dhdt[2],mb[2],ms[2],Jdet;
133  IssmDouble* xyz_list = NULL;
134 
135  /*Fetch number of nodes and dof for this finite element*/
136  int numnodes = element->GetNumberOfNodes();
137 
138  /*Initialize Element vector and other vectors*/
139  ElementVector* pe = element->NewElementVector();
140  IssmDouble* basis = xNew<IssmDouble>(numnodes);
141 
142  /*Retrieve all inputs and parameters*/
143  element->GetVerticesCoordinates(&xyz_list);
144  Input2* ms_input = element->GetInput2(SmbMassBalanceEnum); _assert_(ms_input);
145  Input2* mb_input = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(mb_input);
146  Input2* dhdt_input = element->GetInput2(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
147 
148  /* Start looping on the number of gaussian points: */
149  Gauss* gauss=element->NewGauss(2);
150  for(int ig=gauss->begin();ig<gauss->end();ig++){
151  gauss->GaussPoint(ig);
152 
153  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
154  element->NodalFunctions(basis,gauss);
155 
156  ms_input->GetInputDerivativeValue(&ms[0],xyz_list,gauss);
157 
158  ms_input->GetInputDerivativeValue(&ms[0],xyz_list,gauss);
159  mb_input->GetInputDerivativeValue(&mb[0],xyz_list,gauss);
160  dhdt_input->GetInputDerivativeValue(&dhdt[0],xyz_list,gauss);
161 
162  for(int i=0;i<numnodes;i++) pe->values[i]+=0*Jdet*gauss->weight*(
163  (ms[0]+ms[1]-mb[0]-mb[1]-dhdt[0]-dhdt[1])*basis[i]
164  );
165  }
166 
167  /*Clean up and return*/
168  xDelete<IssmDouble>(xyz_list);
169  xDelete<IssmDouble>(basis);
170  delete gauss;
171  return pe;
172 }/*}}}*/
174  element->GetSolutionFromInputsOneDof(solution,ThicknessEnum);
175 }/*}}}*/
176 void Balancethickness2Analysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
177  _error_("Not implemented yet");
178 }/*}}}*/
180 
181  element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
182 
183 }/*}}}*/
185  /*Default, do nothing*/
186  return;
187 }/*}}}*/
BaseEnum
@ BaseEnum
Definition: EnumDefinitions.h:495
SmbMassBalanceEnum
@ SmbMassBalanceEnum
Definition: EnumDefinitions.h:748
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
Nodes
Declaration of Nodes class.
Definition: Nodes.h:19
Element::GetNumberOfNodes
virtual int GetNumberOfNodes(void)=0
Balancethickness2Analysis::InputUpdateFromSolution
void InputUpdateFromSolution(IssmDouble *solution, Element *element)
Definition: Balancethickness2Analysis.cpp:179
BalancethicknessThickeningRateEnum
@ BalancethicknessThickeningRateEnum
Definition: EnumDefinitions.h:474
Balancethickness2Analysis::CreateNodes
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
Definition: Balancethickness2Analysis.cpp:17
Balancethickness2AnalysisEnum
@ Balancethickness2AnalysisEnum
Definition: EnumDefinitions.h:979
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
MaskIceLevelsetEnum
@ MaskIceLevelsetEnum
Definition: EnumDefinitions.h:641
Balancethickness2Analysis::CreateConstraints
void CreateConstraints(Constraints *constraints, IoModel *iomodel)
Definition: Balancethickness2Analysis.cpp:8
Balancethickness2Analysis::CreateKMatrix
ElementMatrix * CreateKMatrix(Element *element)
Definition: Balancethickness2Analysis.cpp:68
Constraints
Declaration of Constraints class.
Definition: Constraints.h:13
Balancethickness2Analysis::Core
void Core(FemModel *femmodel)
Definition: Balancethickness2Analysis.cpp:58
Elements
Declaration of Elements class.
Definition: Elements.h:17
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
Balancethickness2Analysis::CreateDVector
ElementVector * CreateDVector(Element *element)
Definition: Balancethickness2Analysis.cpp:61
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
Balancethickness2Analysis::DofsPerNode
int DofsPerNode(int **doflist, int domaintype, int approximation)
Definition: Balancethickness2Analysis.cpp:22
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
Input2::GetInputDerivativeValue
virtual void GetInputDerivativeValue(IssmDouble *derivativevalues, IssmDouble *xyz_list, Gauss *gauss)
Definition: Input2.h:37
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
Element
Definition: Element.h:41
Element::NodalFunctions
virtual void NodalFunctions(IssmDouble *basis, Gauss *gauss)=0
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
Balancethickness2Analysis::UpdateConstraints
void UpdateConstraints(FemModel *femmodel)
Definition: Balancethickness2Analysis.cpp:184
Element::NewElementVector
ElementVector * NewElementVector(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2505
IoModel::numberofelements
int numberofelements
Definition: IoModel.h:96
Balancethickness2Analysis.h
: header file for generic external result object
BasalforcingsGroundediceMeltingRateEnum
@ BasalforcingsGroundediceMeltingRateEnum
Definition: EnumDefinitions.h:478
Element::NewGauss
virtual Gauss * NewGauss(void)=0
Element::InputUpdateFromSolutionOneDof
virtual void InputUpdateFromSolutionOneDof(IssmDouble *solution, int inputenum)=0
Balancethickness2Analysis::GetSolutionFromInputs
void GetSolutionFromInputs(Vector< IssmDouble > *solution, Element *element)
Definition: Balancethickness2Analysis.cpp:173
Balancethickness2Analysis::UpdateElements
void UpdateElements(Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
Definition: Balancethickness2Analysis.cpp:25
Balancethickness2Analysis::UpdateParameters
void UpdateParameters(Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
Definition: Balancethickness2Analysis.cpp:54
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
SurfaceEnum
@ SurfaceEnum
Definition: EnumDefinitions.h:823
Inputs2
Declaration of Inputs class.
Definition: Inputs2.h:23
Balancethickness2Analysis::GradientJ
void GradientJ(Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
Definition: Balancethickness2Analysis.cpp:176
Input2
Definition: Input2.h:18
FemModel
Definition: FemModel.h:31
SealevelEnum
@ SealevelEnum
Definition: EnumDefinitions.h:675
Loads
Declaration of Loads class.
Definition: Loads.h:16
Element::GetSolutionFromInputsOneDof
void GetSolutionFromInputsOneDof(Vector< IssmDouble > *solution, int solutionenum)
Definition: Element.cpp:1281
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Balancethickness2Analysis::CreateJacobianMatrix
ElementMatrix * CreateJacobianMatrix(Element *element)
Definition: Balancethickness2Analysis.cpp:65
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
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
Balancethickness2Analysis::CreatePVector
ElementVector * CreatePVector(Element *element)
Definition: Balancethickness2Analysis.cpp:128
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
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
ElementMatrix
Definition: ElementMatrix.h:19
Vector< IssmDouble >
Gauss
Definition: Gauss.h:8
Balancethickness2Analysis::CreateLoads
void CreateLoads(Loads *loads, IoModel *iomodel)
Definition: Balancethickness2Analysis.cpp:14
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