Ice Sheet System Model  4.18
Code documentation
Public Member Functions
BalancevelocityAnalysis Class Reference

#include <BalancevelocityAnalysis.h>

Inheritance diagram for BalancevelocityAnalysis:
Analysis

Public Member Functions

void CreateConstraints (Constraints *constraints, IoModel *iomodel)
 
void CreateLoads (Loads *loads, IoModel *iomodel)
 
void CreateNodes (Nodes *nodes, IoModel *iomodel, bool isamr=false)
 
int DofsPerNode (int **doflist, int domaintype, int approximation)
 
void UpdateElements (Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
 
void UpdateParameters (Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
 
void Core (FemModel *femmodel)
 
ElementVector * CreateDVector (Element *element)
 
ElementMatrix * CreateJacobianMatrix (Element *element)
 
ElementMatrix * CreateKMatrix (Element *element)
 
ElementVector * CreatePVector (Element *element)
 
void GetSolutionFromInputs (Vector< IssmDouble > *solution, Element *element)
 
void GradientJ (Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
 
void InputUpdateFromSolution (IssmDouble *solution, Element *element)
 
void UpdateConstraints (FemModel *femmodel)
 
- Public Member Functions inherited from Analysis
virtual ~Analysis ()
 

Detailed Description

Definition at line 11 of file BalancevelocityAnalysis.h.

Member Function Documentation

◆ CreateConstraints()

void BalancevelocityAnalysis::CreateConstraints ( Constraints *  constraints,
IoModel *  iomodel 
)
virtual

Implements Analysis.

Definition at line 8 of file BalancevelocityAnalysis.cpp.

8  {/*{{{*/
9 
10  /*No constraints for now*/
11  //IoModelToConstraintsx(constraints,iomodel,"md.balancethickness.spcthickness",BalancevelocityAnalysisEnum,P1Enum);
12 }/*}}}*/

◆ CreateLoads()

void BalancevelocityAnalysis::CreateLoads ( Loads *  loads,
IoModel *  iomodel 
)
virtual

Implements Analysis.

Definition at line 13 of file BalancevelocityAnalysis.cpp.

13  {/*{{{*/
14 
15  /*No loads*/
16 }/*}}}*/

◆ CreateNodes()

void BalancevelocityAnalysis::CreateNodes ( Nodes *  nodes,
IoModel *  iomodel,
bool  isamr = false 
)
virtual

Implements Analysis.

Definition at line 17 of file BalancevelocityAnalysis.cpp.

17  {/*{{{*/
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 }/*}}}*/

◆ DofsPerNode()

int BalancevelocityAnalysis::DofsPerNode ( int **  doflist,
int  domaintype,
int  approximation 
)
virtual

Implements Analysis.

Definition at line 25 of file BalancevelocityAnalysis.cpp.

25  {/*{{{*/
26  return 1;
27 }/*}}}*/

◆ UpdateElements()

void BalancevelocityAnalysis::UpdateElements ( Elements *  elements,
Inputs2 *  inputs2,
IoModel *  iomodel,
int  analysis_counter,
int  analysis_type 
)
virtual

Implements Analysis.

Definition at line 28 of file BalancevelocityAnalysis.cpp.

28  {/*{{{*/
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 }/*}}}*/

◆ UpdateParameters()

void BalancevelocityAnalysis::UpdateParameters ( Parameters *  parameters,
IoModel *  iomodel,
int  solution_enum,
int  analysis_enum 
)
virtual

Implements Analysis.

Definition at line 56 of file BalancevelocityAnalysis.cpp.

56  {/*{{{*/
57 }/*}}}*/

◆ Core()

void BalancevelocityAnalysis::Core ( FemModel *  femmodel)
virtual

Implements Analysis.

Definition at line 60 of file BalancevelocityAnalysis.cpp.

60  {/*{{{*/
61  _error_("not implemented");
62 }/*}}}*/

◆ CreateDVector()

ElementVector * BalancevelocityAnalysis::CreateDVector ( Element *  element)
virtual

Implements Analysis.

Definition at line 63 of file BalancevelocityAnalysis.cpp.

63  {/*{{{*/
64  /*Default, return NULL*/
65  return NULL;
66 }/*}}}*/

◆ CreateJacobianMatrix()

ElementMatrix * BalancevelocityAnalysis::CreateJacobianMatrix ( Element *  element)
virtual

Implements Analysis.

Definition at line 67 of file BalancevelocityAnalysis.cpp.

67  {/*{{{*/
68 _error_("Not implemented");
69 }/*}}}*/

◆ CreateKMatrix()

ElementMatrix * BalancevelocityAnalysis::CreateKMatrix ( Element *  element)
virtual

Implements Analysis.

Definition at line 70 of file BalancevelocityAnalysis.cpp.

70  {/*{{{*/
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 }/*}}}*/

◆ CreatePVector()

ElementVector * BalancevelocityAnalysis::CreatePVector ( Element *  element)
virtual

Implements Analysis.

Definition at line 146 of file BalancevelocityAnalysis.cpp.

146  {/*{{{*/
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 }/*}}}*/

◆ GetSolutionFromInputs()

void BalancevelocityAnalysis::GetSolutionFromInputs ( Vector< IssmDouble > *  solution,
Element *  element 
)
virtual

Implements Analysis.

Definition at line 238 of file BalancevelocityAnalysis.cpp.

238  {/*{{{*/
239  _error_("not implemented yet");
240 }/*}}}*/

◆ GradientJ()

void BalancevelocityAnalysis::GradientJ ( Vector< IssmDouble > *  gradient,
Element *  element,
int  control_type,
int  control_index 
)
virtual

Implements Analysis.

Definition at line 241 of file BalancevelocityAnalysis.cpp.

241  {/*{{{*/
242  _error_("Not implemented yet");
243 }/*}}}*/

◆ InputUpdateFromSolution()

void BalancevelocityAnalysis::InputUpdateFromSolution ( IssmDouble *  solution,
Element *  element 
)
virtual

Implements Analysis.

Definition at line 244 of file BalancevelocityAnalysis.cpp.

244  {/*{{{*/
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 }/*}}}*/

◆ UpdateConstraints()

void BalancevelocityAnalysis::UpdateConstraints ( FemModel *  femmodel)
virtual

Implements Analysis.

Definition at line 258 of file BalancevelocityAnalysis.cpp.

258  {/*{{{*/
259  /*Default, do nothing*/
260  return;
261 }/*}}}*/
The documentation for this class was generated from the following files:
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
Element::GetNumberOfNodes
virtual int GetNumberOfNodes(void)=0
BalancevelocityAnalysisEnum
@ BalancevelocityAnalysisEnum
Definition: EnumDefinitions.h:987
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
BalancethicknessThickeningRateEnum
@ BalancethicknessThickeningRateEnum
Definition: EnumDefinitions.h:474
DrivingStressYEnum
@ DrivingStressYEnum
Definition: EnumDefinitions.h:539
MaskIceLevelsetEnum
@ MaskIceLevelsetEnum
Definition: EnumDefinitions.h:641
Element::SpawnBasalElement
virtual Element * SpawnBasalElement(void)=0
MeshVertexonbaseEnum
@ MeshVertexonbaseEnum
Definition: EnumDefinitions.h:653
DrivingStressXEnum
@ DrivingStressXEnum
Definition: EnumDefinitions.h:538
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
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
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
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
SurfaceEnum
@ SurfaceEnum
Definition: EnumDefinitions.h:823
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
SealevelEnum
@ SealevelEnum
Definition: EnumDefinitions.h:675
Element::ValueP1DerivativesOnGauss
virtual void ValueP1DerivativesOnGauss(IssmDouble *dvalue, IssmDouble *values, IssmDouble *xyz_list, Gauss *gauss)=0
_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::domaintype
int domaintype
Definition: IoModel.h:78
ElementMatrix
Definition: ElementMatrix.h:19
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
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
Element::NewElementMatrix
ElementMatrix * NewElementMatrix(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2497
Generated on Thu Jul 2 2020 08:09:26 for Ice Sheet System Model by   doxygen 1.8.19