Ice Sheet System Model  4.18
Code documentation
BalancevelocityAnalysis.cpp
Go to the documentation of this file.
1 #include "./BalancevelocityAnalysis.h"
2 #include "../toolkits/toolkits.h"
3 #include "../classes/classes.h"
4 #include "../shared/shared.h"
5 #include "../modules/modules.h"
6 
7 /*Model processing*/
8 void BalancevelocityAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
9 
10  /*No constraints for now*/
11  //IoModelToConstraintsx(constraints,iomodel,"md.balancethickness.spcthickness",BalancevelocityAnalysisEnum,P1Enum);
12 }/*}}}*/
13 void BalancevelocityAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
14 
15  /*No loads*/
16 }/*}}}*/
17 void BalancevelocityAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
18 
19  /*Check in 3d*/
20  if(iomodel->domaintype==Domain3DEnum) _error_("DG 3d not implemented yet");
21 
22  /*First fetch data: */
23  ::CreateNodes(nodes,iomodel,BalancevelocityAnalysisEnum,P1Enum);
24 }/*}}}*/
25 int BalancevelocityAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
26  return 1;
27 }/*}}}*/
28 void BalancevelocityAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
29 
30  /*Update elements: */
31  int counter=0;
32  for(int i=0;i<iomodel->numberofelements;i++){
33  if(iomodel->my_elements[i]){
34  Element* element=(Element*)elements->GetObjectByOffset(counter);
35  element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,P1Enum);
36  counter++;
37  }
38  }
39 
40  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.thickness",ThicknessEnum);
41  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.surface",SurfaceEnum);
42  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.base",BaseEnum);
43  iomodel->FetchDataToInput(inputs2,elements,"md.solidearth.sealevel",SealevelEnum,0);
44  iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
45  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vx",VxEnum);
46  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vy",VyEnum);
47  iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
48  iomodel->FetchDataToInput(inputs2,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
49  iomodel->FetchDataToInput(inputs2,elements,"md.balancethickness.thickening_rate",BalancethicknessThickeningRateEnum);
50 
51  if(iomodel->domaintype!=Domain2DhorizontalEnum){
52  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
53  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
54  }
55 }/*}}}*/
56 void BalancevelocityAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
57 }/*}}}*/
58 
59 /*Finite Element Analysis*/
60 void BalancevelocityAnalysis::Core(FemModel* femmodel){/*{{{*/
61  _error_("not implemented");
62 }/*}}}*/
63 ElementVector* BalancevelocityAnalysis::CreateDVector(Element* element){/*{{{*/
64  /*Default, return NULL*/
65  return NULL;
66 }/*}}}*/
67 ElementMatrix* BalancevelocityAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
68 _error_("Not implemented");
69 }/*}}}*/
70 ElementMatrix* BalancevelocityAnalysis::CreateKMatrix(Element* element){/*{{{*/
71 
72  /*Intermediaries */
73  IssmDouble dhdt,mb,ms,Jdet;
74  IssmDouble h,gamma,thickness;
75  IssmDouble hnx,hny,dhnx[2],dhny[2];
76  IssmDouble *xyz_list = NULL;
77 
78  /*Fetch number of nodes and dof for this finite element*/
79  int numnodes = element->GetNumberOfNodes();
80 
81  /*Initialize Element matrix and vectors*/
82  ElementMatrix* Ke = element->NewElementMatrix();
83  IssmDouble* basis = xNew<IssmDouble>(numnodes);
84  IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);
85  IssmDouble* HNx = xNew<IssmDouble>(numnodes);
86  IssmDouble* HNy = xNew<IssmDouble>(numnodes);
87  IssmDouble* H = xNew<IssmDouble>(numnodes);
88  IssmDouble* Nx = xNew<IssmDouble>(numnodes);
89  IssmDouble* Ny = xNew<IssmDouble>(numnodes);
90 
91  /*Retrieve all Inputs and parameters: */
92  element->GetVerticesCoordinates(&xyz_list);
93  Input2* H_input = element->GetInput2(ThicknessEnum); _assert_(H_input);
94  h = element->CharacteristicLength();
95 
96  /*Get vector N for all nodes and build HNx and HNy*/
97  element->GetInputListOnNodes(Nx,DrivingStressXEnum);
98  element->GetInputListOnNodes(Ny,DrivingStressYEnum);
99  element->GetInputListOnNodes(H,ThicknessEnum);
100  for(int i=0;i<numnodes;i++){
101  IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]+1.e-10);
102  HNx[i] = -H[i]*Nx[i]/norm;
103  HNy[i] = -H[i]*Ny[i]/norm;
104  }
105 
106  /*Start looping on the number of gaussian points:*/
107  Gauss* gauss=element->NewGauss(2);
108  for(int ig=gauss->begin();ig<gauss->end();ig++){
109  gauss->GaussPoint(ig);
110 
111  H_input->GetInputValue(&thickness,gauss);
112  if(thickness<50.) thickness=50.;
113  element->ValueP1DerivativesOnGauss(&dhnx[0],HNx,xyz_list,gauss);
114  element->ValueP1DerivativesOnGauss(&dhny[0],HNy,xyz_list,gauss);
115  element->ValueP1OnGauss(&hnx,HNx,gauss);
116  element->ValueP1OnGauss(&hny,HNy,gauss);
117 
118  gamma=h/(2.*thickness+1.e-10);
119 
120  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
121  element->NodalFunctions(basis,gauss);
122  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
123 
124  for(int i=0;i<numnodes;i++){
125  for(int j=0;j<numnodes;j++){
126  Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
127  (basis[i]+gamma*(basis[i]*(dhnx[0]+dhny[1]) + dbasis[0*numnodes+i]*hnx + dbasis[1*numnodes+i]*hny))*
128  (basis[j]*(dhnx[0]+dhny[1]) + dbasis[0*numnodes+j]*hnx + dbasis[1*numnodes+j]*hny)
129  );
130  }
131  }
132  }
133 
134  /*Clean up and return*/
135  xDelete<IssmDouble>(xyz_list);
136  xDelete<IssmDouble>(basis);
137  xDelete<IssmDouble>(dbasis);
138  xDelete<IssmDouble>(H);
139  xDelete<IssmDouble>(Nx);
140  xDelete<IssmDouble>(Ny);
141  xDelete<IssmDouble>(HNx);
142  xDelete<IssmDouble>(HNy);
143  delete gauss;
144  return Ke;
145 }/*}}}*/
146 ElementVector* BalancevelocityAnalysis::CreatePVector(Element* element){/*{{{*/
147 
148  /*Intermediaries*/
149  int domaintype;
150  Element* basalelement;
151 
152  /*Get basal element*/
153  element->FindParam(&domaintype,DomainTypeEnum);
154  switch(domaintype){
155  case Domain2DhorizontalEnum:
156  basalelement = element;
157  break;
158  case Domain3DEnum:
159  if(!element->IsOnBase()) return NULL;
160  basalelement = element->SpawnBasalElement();
161  break;
162  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
163  }
164 
165  /*Intermediaries */
166  IssmDouble dhdt,mb,ms,Jdet;
167  IssmDouble gamma,thickness;
168  IssmDouble hnx,hny,dhnx[2],dhny[2];
169  IssmDouble *xyz_list = NULL;
170 
171  /*Fetch number of nodes and dof for this finite element*/
172  int numnodes = basalelement->GetNumberOfNodes();
173 
174  /*Initialize Element vector*/
175  ElementVector* pe = basalelement->NewElementVector();
176  IssmDouble* basis = xNew<IssmDouble>(numnodes);
177  IssmDouble* dbasis = xNew<IssmDouble>(numnodes*2);
178  IssmDouble* H = xNew<IssmDouble>(numnodes);
179  IssmDouble* Nx = xNew<IssmDouble>(numnodes);
180  IssmDouble* Ny = xNew<IssmDouble>(numnodes);
181 
182  /*Retrieve all inputs and parameters*/
183  basalelement->GetVerticesCoordinates(&xyz_list);
184  Input2* ms_input = basalelement->GetInput2(SmbMassBalanceEnum); _assert_(ms_input);
185  Input2* mb_input = basalelement->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(mb_input);
186  Input2* dhdt_input = basalelement->GetInput2(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
187  Input2* H_input = basalelement->GetInput2(ThicknessEnum); _assert_(H_input);
188  IssmDouble h = basalelement->CharacteristicLength();
189 
190  /*Get vector N for all nodes*/
191  basalelement->GetInputListOnNodes(Nx,DrivingStressXEnum);
192  basalelement->GetInputListOnNodes(Ny,DrivingStressYEnum);
193  basalelement->GetInputListOnNodes(H,ThicknessEnum);
194  for(int i=0;i<numnodes;i++){
195  IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]+1.e-10);
196  Nx[i] = -H[i]*Nx[i]/norm;
197  Ny[i] = -H[i]*Ny[i]/norm;
198  }
199 
200  /* Start looping on the number of gaussian points: */
201  Gauss* gauss=basalelement->NewGauss(2);
202  for(int ig=gauss->begin();ig<gauss->end();ig++){
203  gauss->GaussPoint(ig);
204 
205  basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
206  basalelement->NodalFunctions(basis,gauss);
207  basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
208 
209  element->ValueP1DerivativesOnGauss(&dhnx[0],Nx,xyz_list,gauss);
210  element->ValueP1DerivativesOnGauss(&dhny[0],Ny,xyz_list,gauss);
211  element->ValueP1OnGauss(&hnx,Nx,gauss);
212  element->ValueP1OnGauss(&hny,Ny,gauss);
213 
214  ms_input->GetInputValue(&ms,gauss);
215  mb_input->GetInputValue(&mb,gauss);
216  dhdt_input->GetInputValue(&dhdt,gauss);
217  H_input->GetInputValue(&thickness,gauss);
218  if(thickness<50.) thickness=50.;
219 
220  gamma=h/(2.*thickness+1.e-10);
221 
222  for(int i=0;i<numnodes;i++){
223  pe->values[i]+=Jdet*gauss->weight*(ms-mb-dhdt)*( basis[i] + gamma*(basis[i]*(dhnx[0]+dhny[1])+hnx*dbasis[0*numnodes+i] + hny*dbasis[1*numnodes+i]));
224  }
225  }
226 
227  /*Clean up and return*/
228  xDelete<IssmDouble>(xyz_list);
229  xDelete<IssmDouble>(basis);
230  xDelete<IssmDouble>(dbasis);
231  xDelete<IssmDouble>(H);
232  xDelete<IssmDouble>(Nx);
233  xDelete<IssmDouble>(Ny);
234  delete gauss;
235  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
236  return pe;
237 }/*}}}*/
238 void BalancevelocityAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
239  _error_("not implemented yet");
240 }/*}}}*/
241 void BalancevelocityAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
242  _error_("Not implemented yet");
243 }/*}}}*/
244 void BalancevelocityAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
245 
246  int domaintype;
247  element->FindParam(&domaintype,DomainTypeEnum);
248  switch(domaintype){
249  case Domain2DhorizontalEnum:
250  element->InputUpdateFromSolutionOneDof(solution,VelEnum);
251  break;
252  case Domain3DEnum:
253  element->InputUpdateFromSolutionOneDofCollapsed(solution,VelEnum);
254  break;
255  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
256  }
257 }/*}}}*/
258 void BalancevelocityAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
259  /*Default, do nothing*/
260  return;
261 }/*}}}*/
BalancevelocityAnalysis::DofsPerNode
int DofsPerNode(int **doflist, int domaintype, int approximation)
Definition: BalancevelocityAnalysis.cpp:25
BaseEnum
@ BaseEnum
Definition: EnumDefinitions.h:495
SmbMassBalanceEnum
@ SmbMassBalanceEnum
Definition: EnumDefinitions.h:748
Element::GetInputListOnNodes
void GetInputListOnNodes(IssmDouble *pvalue, int enumtype)
Definition: Element.cpp:1106
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
Element::IsOnBase
bool IsOnBase()
Definition: Element.cpp:1984
Nodes
Declaration of Nodes class.
Definition: Nodes.h:19
BalancevelocityAnalysis::CreateKMatrix
ElementMatrix * CreateKMatrix(Element *element)
Definition: BalancevelocityAnalysis.cpp:70
Element::GetNumberOfNodes
virtual int GetNumberOfNodes(void)=0
BalancevelocityAnalysis::GetSolutionFromInputs
void GetSolutionFromInputs(Vector< IssmDouble > *solution, Element *element)
Definition: BalancevelocityAnalysis.cpp:238
BalancevelocityAnalysisEnum
@ BalancevelocityAnalysisEnum
Definition: EnumDefinitions.h:987
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
BalancethicknessThickeningRateEnum
@ BalancethicknessThickeningRateEnum
Definition: EnumDefinitions.h:474
BalancevelocityAnalysis::UpdateParameters
void UpdateParameters(Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
Definition: BalancevelocityAnalysis.cpp:56
DrivingStressYEnum
@ DrivingStressYEnum
Definition: EnumDefinitions.h:539
BalancevelocityAnalysis.h
: header file for generic external result object
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
BalancevelocityAnalysis::Core
void Core(FemModel *femmodel)
Definition: BalancevelocityAnalysis.cpp:60
MaskIceLevelsetEnum
@ MaskIceLevelsetEnum
Definition: EnumDefinitions.h:641
Constraints
Declaration of Constraints class.
Definition: Constraints.h:13
Element::SpawnBasalElement
virtual Element * SpawnBasalElement(void)=0
MeshVertexonbaseEnum
@ MeshVertexonbaseEnum
Definition: EnumDefinitions.h:653
Elements
Declaration of Elements class.
Definition: Elements.h:17
DrivingStressXEnum
@ DrivingStressXEnum
Definition: EnumDefinitions.h:538
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
BalancevelocityAnalysis::CreateLoads
void CreateLoads(Loads *loads, IoModel *iomodel)
Definition: BalancevelocityAnalysis.cpp:13
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
BalancevelocityAnalysis::UpdateElements
void UpdateElements(Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
Definition: BalancevelocityAnalysis.cpp:28
Element::DeleteMaterials
void DeleteMaterials(void)
Definition: Element.cpp:429
Element
Definition: Element.h:41
Element::NodalFunctions
virtual void NodalFunctions(IssmDouble *basis, Gauss *gauss)=0
Element::ValueP1OnGauss
virtual void ValueP1OnGauss(IssmDouble *pvalue, IssmDouble *values, Gauss *gauss)=0
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
Element::NewElementVector
ElementVector * NewElementVector(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2505
IoModel::numberofelements
int numberofelements
Definition: IoModel.h:96
BasalforcingsGroundediceMeltingRateEnum
@ BasalforcingsGroundediceMeltingRateEnum
Definition: EnumDefinitions.h:478
DomainTypeEnum
@ DomainTypeEnum
Definition: EnumDefinitions.h:124
BalancevelocityAnalysis::CreateJacobianMatrix
ElementMatrix * CreateJacobianMatrix(Element *element)
Definition: BalancevelocityAnalysis.cpp:67
Element::NewGauss
virtual Gauss * NewGauss(void)=0
Element::InputUpdateFromSolutionOneDof
virtual void InputUpdateFromSolutionOneDof(IssmDouble *solution, int inputenum)=0
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
BalancevelocityAnalysis::CreatePVector
ElementVector * CreatePVector(Element *element)
Definition: BalancevelocityAnalysis.cpp:146
BalancevelocityAnalysis::UpdateConstraints
void UpdateConstraints(FemModel *femmodel)
Definition: BalancevelocityAnalysis.cpp:258
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
BalancevelocityAnalysis::CreateConstraints
void CreateConstraints(Constraints *constraints, IoModel *iomodel)
Definition: BalancevelocityAnalysis.cpp:8
SurfaceEnum
@ SurfaceEnum
Definition: EnumDefinitions.h:823
BalancevelocityAnalysis::CreateDVector
ElementVector * CreateDVector(Element *element)
Definition: BalancevelocityAnalysis.cpp:63
Inputs2
Declaration of Inputs class.
Definition: Inputs2.h:23
BalancevelocityAnalysis::CreateNodes
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
Definition: BalancevelocityAnalysis.cpp:17
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
Input2
Definition: Input2.h:18
FemModel
Definition: FemModel.h:31
SealevelEnum
@ SealevelEnum
Definition: EnumDefinitions.h:675
Element::ValueP1DerivativesOnGauss
virtual void ValueP1DerivativesOnGauss(IssmDouble *dvalue, IssmDouble *values, IssmDouble *xyz_list, Gauss *gauss)=0
Loads
Declaration of Loads class.
Definition: Loads.h:16
BalancevelocityAnalysis::GradientJ
void GradientJ(Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
Definition: BalancevelocityAnalysis.cpp:241
_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
Element::CharacteristicLength
virtual IssmDouble CharacteristicLength(void)=0
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
MeshVertexonsurfaceEnum
@ MeshVertexonsurfaceEnum
Definition: EnumDefinitions.h:655
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
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
VelEnum
@ VelEnum
Definition: EnumDefinitions.h:844
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
IoModel::domaintype
int domaintype
Definition: IoModel.h:78
ElementMatrix
Definition: ElementMatrix.h:19
Vector< IssmDouble >
Element::InputUpdateFromSolutionOneDofCollapsed
virtual void InputUpdateFromSolutionOneDofCollapsed(IssmDouble *solution, int inputenum)=0
Gauss
Definition: Gauss.h:8
Element::NodalFunctionsDerivatives
virtual void NodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
BalancevelocityAnalysis::InputUpdateFromSolution
void InputUpdateFromSolution(IssmDouble *solution, Element *element)
Definition: BalancevelocityAnalysis.cpp:244
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
Generated on Thu Jul 2 2020 08:09:14 for Ice Sheet System Model by   doxygen 1.8.19