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

#include <UzawaPressureAnalysis.h>

Inheritance diagram for UzawaPressureAnalysis:
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)
 
ElementVectorCreateDVector (Element *element)
 
ElementMatrixCreateJacobianMatrix (Element *element)
 
ElementMatrixCreateKMatrix (Element *element)
 
ElementVectorCreatePVector (Element *element)
 
void GetM (IssmDouble *M, Element *element, Gauss *gauss)
 
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 UzawaPressureAnalysis.h.

Member Function Documentation

◆ CreateConstraints()

void UzawaPressureAnalysis::CreateConstraints ( Constraints constraints,
IoModel iomodel 
)
virtual

Implements Analysis.

Definition at line 8 of file UzawaPressureAnalysis.cpp.

8  {/*{{{*/
9  return;
10 }/*}}}*/

◆ CreateLoads()

void UzawaPressureAnalysis::CreateLoads ( Loads loads,
IoModel iomodel 
)
virtual

Implements Analysis.

Definition at line 11 of file UzawaPressureAnalysis.cpp.

11  {/*{{{*/
12  return;
13 }/*}}}*/

◆ CreateNodes()

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

Implements Analysis.

Definition at line 14 of file UzawaPressureAnalysis.cpp.

14  {/*{{{*/
15 
16  int finiteelement;
17  int fe_FS;
18 
19  iomodel->FindConstant(&fe_FS,"md.flowequation.fe_FS");
20  if(fe_FS==LATaylorHoodEnum) finiteelement = P1Enum;
21  else if(fe_FS==LACrouzeixRaviartEnum) finiteelement = P1DGEnum;
22  else _error_("solution not supported yet");
23 
24  ::CreateNodes(nodes,iomodel,UzawaPressureAnalysisEnum,finiteelement);
25 }/*}}}*/

◆ DofsPerNode()

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

Implements Analysis.

Definition at line 26 of file UzawaPressureAnalysis.cpp.

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

◆ UpdateElements()

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

Implements Analysis.

Definition at line 29 of file UzawaPressureAnalysis.cpp.

29  {/*{{{*/
30 
31  /*Update elements: */
32  int finiteelement;
33  int counter=0;
34  int fe_FS;
35 
36  iomodel->FindConstant(&fe_FS,"md.flowequation.fe_FS");
37  if(fe_FS==LATaylorHoodEnum) finiteelement = P1Enum;
38  else if(fe_FS==LACrouzeixRaviartEnum) finiteelement = P1DGEnum;
39  else _error_("solution not supported yet");
40 
41  for(int i=0;i<iomodel->numberofelements;i++){
42  if(iomodel->my_elements[i]){
43  Element* element=(Element*)elements->GetObjectByOffset(counter);
44  element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,finiteelement);
45  counter++;
46  }
47  }
48 
49  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vx",VxEnum,0.);
50  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vy",VyEnum,0.);
51  if(iomodel->domaintype==Domain3DEnum) iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vz",VzEnum,0.);
52  iomodel->FetchDataToInput(inputs2,elements,"md.initialization.pressure",PressureEnum,0.);
53  InputUpdateFromConstantx(inputs2,elements,0.,SigmaNNEnum);
54 }/*}}}*/

◆ UpdateParameters()

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

Implements Analysis.

Definition at line 55 of file UzawaPressureAnalysis.cpp.

55  {/*{{{*/
56 
57  parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.augmented_lagrangian_rhop",AugmentedLagrangianRhopEnum));
58  parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.augmented_lagrangian_rholambda",AugmentedLagrangianRholambdaEnum));
59 }/*}}}*/

◆ Core()

void UzawaPressureAnalysis::Core ( FemModel femmodel)
virtual

Implements Analysis.

Definition at line 62 of file UzawaPressureAnalysis.cpp.

62  {/*{{{*/
63  _error_("not implemented");
64 }/*}}}*/

◆ CreateDVector()

ElementVector * UzawaPressureAnalysis::CreateDVector ( Element element)
virtual

Implements Analysis.

Definition at line 65 of file UzawaPressureAnalysis.cpp.

65  {/*{{{*/
66  /*Default, return NULL*/
67  return NULL;
68 }/*}}}*/

◆ CreateJacobianMatrix()

ElementMatrix * UzawaPressureAnalysis::CreateJacobianMatrix ( Element element)
virtual

Implements Analysis.

Definition at line 69 of file UzawaPressureAnalysis.cpp.

69  {/*{{{*/
70 _error_("Not implemented");
71 }/*}}}*/

◆ CreateKMatrix()

ElementMatrix * UzawaPressureAnalysis::CreateKMatrix ( Element element)
virtual

Implements Analysis.

Definition at line 72 of file UzawaPressureAnalysis.cpp.

72  {/*{{{*/
73 
74  /*Intermediaries*/
75  IssmDouble D_scalar,Jdet;
76  IssmDouble *xyz_list = NULL;
77  int numnodes = element->GetNumberOfNodes();
78 
79  /*Initialize Element matrix and vectors*/
80  ElementMatrix* Ke = element->NewElementMatrix();
81  IssmDouble* M = xNew<IssmDouble>(numnodes);
82 
83  IssmDouble connectivity;
84  /*Retrieve all inputs and parameters*/
85  element->GetVerticesCoordinates(&xyz_list);
86  int numvertices = element->GetNumberOfVertices();
87 
88  //Gauss* gauss = element->NewGauss(5);
89  //for(int ig=gauss->begin();ig<gauss->end();ig++){
90  // gauss->GaussPoint(ig);
91 
92  // element->JacobianDeterminant(&Jdet,xyz_list,gauss);
93  // this->GetM(M,element,gauss);
94  // D_scalar=gauss->weight*Jdet;
95  // //TripleMultiply(M,1,numnodes,1,
96  // // &D_scalar,1,1,0,
97  // // M,1,numnodes,0,
98  // // &Ke->values[0],1);
99 
100  //}
101  for(int iv=0;iv<numvertices;iv++){
102  connectivity=(IssmDouble)element->VertexConnectivity(iv);
103  Ke->values[iv*numvertices+iv]=1./connectivity;
104  }
105 
106  /*Clean up and return*/
107  //delete gauss;
108  xDelete<IssmDouble>(xyz_list);
109  xDelete<IssmDouble>(M);
110  return Ke;
111 }/*}}}*/

◆ CreatePVector()

ElementVector * UzawaPressureAnalysis::CreatePVector ( Element element)
virtual

Implements Analysis.

Definition at line 112 of file UzawaPressureAnalysis.cpp.

112  {/*{{{*/
113 
114  /*Intermediaries*/
115  int dim;
116  IssmDouble Jdet,rhop,divu;
117  IssmDouble *xyz_list = NULL;
118  int numnodes = element->GetNumberOfNodes();
119 
120  /*Retrieve all inputs and parameters*/
121  element->FindParam(&dim,DomainDimensionEnum);
122  element->FindParam(&rhop,AugmentedLagrangianRhopEnum);
123  element->GetVerticesCoordinates(&xyz_list);
124 
125  /*Initialize Element matrix and vectors*/
126  ElementVector *pe = element->NewElementVector();
127  IssmDouble *basis = xNew<IssmDouble>(numnodes);
128  IssmDouble dvx[3];
129  IssmDouble dvy[3];
130  IssmDouble dvz[3];
131 
132  Input2* vx_input=element->GetInput2(VxEnum); _assert_(vx_input);
133  Input2* vy_input=element->GetInput2(VyEnum); _assert_(vy_input);
134  Input2* vz_input = NULL;
135  if(dim==3){vz_input=element->GetInput2(VzEnum); _assert_(vz_input);}
136 
137  Gauss* gauss = element->NewGauss(5);
138  for(int ig=gauss->begin();ig<gauss->end();ig++){
139  gauss->GaussPoint(ig);
140 
141  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
142  element->NodalFunctions(basis, gauss);
143  vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
144  vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
145  if(dim==3){
146  vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
147  }
148 
149  divu=dvx[0]+dvy[1];
150  if (dim==3) divu=divu+dvz[2];
151 
152  for(int i=0;i<numnodes;i++){
153  pe->values[i] += - rhop * divu * Jdet * gauss->weight * basis[i];
154  }
155  }
156 
157  /*Clean up and return*/
158  delete gauss;
159  xDelete<IssmDouble>(xyz_list);
160  xDelete<IssmDouble>(basis);
161  return pe;
162 }/*}}}*/

◆ GetM()

void UzawaPressureAnalysis::GetM ( IssmDouble M,
Element element,
Gauss gauss 
)

Definition at line 163 of file UzawaPressureAnalysis.cpp.

163  {/*{{{*/
164  /*Compute B matrix. M=[M1 M2 M3] */
165 
166  /*Fetch number of nodes for this finite element*/
167  int numnodes = element->GetNumberOfNodes();
168 
169  /*Get nodal functions*/
170  IssmDouble* basis=xNew<IssmDouble>(numnodes);
171  element->NodalFunctions(basis,gauss);
172 
173  /*Build B: */
174  for(int i=0;i<numnodes;i++){
175  M[i] = basis[i];
176  }
177 
178  /*Clean-up*/
179  xDelete<IssmDouble>(basis);
180 }/*}}}*/

◆ GetSolutionFromInputs()

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

Implements Analysis.

Definition at line 181 of file UzawaPressureAnalysis.cpp.

181  {/*{{{*/
182  _error_("not implemented yet");
183 }/*}}}*/

◆ GradientJ()

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

Implements Analysis.

Definition at line 184 of file UzawaPressureAnalysis.cpp.

184  {/*{{{*/
185  _error_("Not implemented yet");
186 }/*}}}*/

◆ InputUpdateFromSolution()

void UzawaPressureAnalysis::InputUpdateFromSolution ( IssmDouble solution,
Element element 
)
virtual

Implements Analysis.

Definition at line 187 of file UzawaPressureAnalysis.cpp.

187  {/*{{{*/
188 
189  int dim;
190  int *doflist = NULL;
191  IssmDouble rholambda,un,vx,vy,vz,sigmann;
192  IssmDouble *xyz_list_base = NULL;
193 
194  /*Fetch number of nodes and dof for this finite element*/
195  int numnodes = element->GetNumberOfNodes();
196  int numnodessigma;
197  if(element->element_type==P1Enum) numnodessigma=element->GetNumberOfNodes(P2Enum);
198  else if(element->element_type==P1DGEnum) numnodessigma=element->GetNumberOfNodes(P2Enum);
199  else _error_("finite element not supported yet");
200 
201  element->FindParam(&dim,DomainDimensionEnum);
202 
203  /*Fetch dof list and allocate solution vector*/
205  IssmDouble* values = xNew<IssmDouble>(numnodes);
206  IssmDouble* valueslambda = xNewZeroInit<IssmDouble>(numnodessigma);
207  IssmDouble* pressure = xNew<IssmDouble>(numnodes);
208  Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
209  Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input);
210  Input2* vz_input = NULL;
211  if(dim==3){vz_input = element->GetInput2(VzEnum); _assert_(vz_input);}
212  element->GetInputListOnNodes(&pressure[0],PressureEnum);
213 
214  /*Update pressure enum first*/
215  for(int i=0;i<numnodes;i++){
216  values[i] = pressure[i] + solution[doflist[i]];
217  }
218  element->AddInput2(PressureEnum,values,element->GetElementType());
219 
220  /*Now compute sigmann if on base*/
221  if(element->IsOnBase() && 0){
222  Input2* sigmann_input = element->GetInput2(SigmaNNEnum); _assert_(sigmann_input);
223  if(dim==3) _error_("not implemented yet");
224 
225  int baselist[3];
226  int onbase=0;
227  IssmDouble Jdet;
228  IssmDouble bed_normal[3];
229  IssmDouble Jlambda[3][3] = {0.0};
230  IssmDouble Cuk[3] = {0.0};
231  IssmDouble deltalambda[3] = {0.0};
232  IssmDouble* vertexonbase = xNew<IssmDouble>(numnodessigma);
233  Input2* vertexonbase_input = element->GetInput2(MeshVertexonbaseEnum); _assert_(vertexonbase_input);
234  Gauss* gauss = element->NewGauss();
235 
236  IssmDouble* basis = xNewZeroInit<IssmDouble>(numnodessigma);
237  element->GetVerticesCoordinatesBase(&xyz_list_base);
238  element->NormalBase(&bed_normal[0],xyz_list_base);
239  element->FindParam(&rholambda,AugmentedLagrangianRholambdaEnum);
240 
241  for(int i=0;i<numnodessigma;i++){
242  gauss->GaussNode(P2Enum,i);
243  vertexonbase_input->GetInputValue(&vertexonbase[i], gauss);
244  if(vertexonbase[i]==1){
245  baselist[onbase]=i;
246  onbase += 1;
247  }
248  }
249  if(onbase!=3) _error_("basal nodes of element not found");
250 
251  delete gauss;
252  gauss = element->NewGaussBase(3);
253  for(int ig=gauss->begin();ig<gauss->end();ig++){
254  gauss->GaussPoint(ig);
255 
256  /*Compute Jlambda*/
257  element->NodalFunctionsP2(basis,gauss);
258  element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
259  for(int i=0;i<3;i++){
260  for(int j=0;j<3;j++){
261  Jlambda[i][j] += Jdet*gauss->weight*basis[baselist[i]]*basis[baselist[j]];
262  }
263  }
264 
265  /*Compute rho_lambd C u^k*/
266  vx_input->GetInputValue(&vx, gauss);
267  vy_input->GetInputValue(&vy, gauss);
268  un=bed_normal[0]*vx + bed_normal[1]*vy;
269  for(int i=0;i<3;i++) Cuk[i] += - un*rholambda*Jdet*gauss->weight*basis[baselist[i]];
270  }
271 
272  /*Now update sigmann*/
273  Matrix3x3Solve(&deltalambda[0],&Jlambda[0][0],&Cuk[0]);
274  delete gauss;
275  gauss = element->NewGauss();
276  for(int i=0;i<3;i++){
277  gauss->GaussNode(P2Enum,baselist[i]);
278  sigmann_input->GetInputValue(&sigmann, gauss);
279  valueslambda[baselist[i]] = sigmann + deltalambda[i];
280  }
281 
282  delete gauss;
283  xDelete<IssmDouble>(vertexonbase);
284  xDelete<IssmDouble>(xyz_list_base);
285  xDelete<IssmDouble>(basis);
286 
287  element->AddInput2(SigmaNNEnum,valueslambda,P2Enum);
288  }
289 
290  /*Free ressources:*/
291  xDelete<IssmDouble>(values);
292  xDelete<IssmDouble>(valueslambda);
293  xDelete<IssmDouble>(pressure);
294  xDelete<int>(doflist);
295 }/*}}}*/

◆ UpdateConstraints()

void UzawaPressureAnalysis::UpdateConstraints ( FemModel femmodel)
virtual

Implements Analysis.

Definition at line 296 of file UzawaPressureAnalysis.cpp.

296  {/*{{{*/
297  /*Default, do nothing*/
298  return;
299 }/*}}}*/

The documentation for this class was generated from the following files:
Element::GetElementType
virtual int GetElementType(void)=0
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
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
Element::NewGaussBase
virtual Gauss * NewGaussBase(int order)=0
UzawaPressureAnalysis::CreateNodes
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
Definition: UzawaPressureAnalysis.cpp:14
MeshVertexonbaseEnum
@ MeshVertexonbaseEnum
Definition: EnumDefinitions.h:653
PressureEnum
@ PressureEnum
Definition: EnumDefinitions.h:664
P1DGEnum
@ P1DGEnum
Definition: EnumDefinitions.h:1215
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
Gauss::GaussNode
virtual void GaussNode(int finitelement, int iv)=0
Element::JacobianDeterminantBase
virtual void JacobianDeterminantBase(IssmDouble *Jdet, IssmDouble *xyz_list_base, Gauss *gauss)=0
Element::NormalBase
virtual void NormalBase(IssmDouble *normal, IssmDouble *xyz_list)=0
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
Parameters::AddObject
void AddObject(Param *newparam)
Definition: Parameters.cpp:67
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
UzawaPressureAnalysisEnum
@ UzawaPressureAnalysisEnum
Definition: EnumDefinitions.h:1320
Element::NewElementVector
ElementVector * NewElementVector(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2505
LATaylorHoodEnum
@ LATaylorHoodEnum
Definition: EnumDefinitions.h:1139
IoModel::numberofelements
int numberofelements
Definition: IoModel.h:96
IoModel::CopyConstantObject
Param * CopyConstantObject(const char *constant_name, int param_enum)
Definition: IoModel.cpp:418
NoneApproximationEnum
@ NoneApproximationEnum
Definition: EnumDefinitions.h:1201
Element::AddInput2
virtual void AddInput2(int input_enum, IssmDouble *values, int interpolation_enum)
Definition: Element.h:216
Element::NewGauss
virtual Gauss * NewGauss(void)=0
VzEnum
@ VzEnum
Definition: EnumDefinitions.h:853
DomainDimensionEnum
@ DomainDimensionEnum
Definition: EnumDefinitions.h:123
GsetEnum
@ GsetEnum
Definition: EnumDefinitions.h:1093
IoModel::FindConstant
void FindConstant(bool *pvalue, const char *constant_name)
Definition: IoModel.cpp:2362
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
InputUpdateFromConstantx
void InputUpdateFromConstantx(FemModel *femmodel, bool constant, int name)
Definition: InputUpdateFromConstantx.cpp:10
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
Input2
Definition: Input2.h:18
Element::GetDofListLocal
void GetDofListLocal(int **pdoflist, int approximation_enum, int setenum)
Definition: Element.cpp:984
LACrouzeixRaviartEnum
@ LACrouzeixRaviartEnum
Definition: EnumDefinitions.h:1138
Element::GetVerticesCoordinatesBase
virtual void GetVerticesCoordinatesBase(IssmDouble **xyz_list)=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
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
Element::NodalFunctionsP2
virtual void NodalFunctionsP2(IssmDouble *basis, Gauss *gauss)=0
Element::JacobianDeterminant
virtual void JacobianDeterminant(IssmDouble *Jdet, IssmDouble *xyz_list, Gauss *gauss)=0
Gauss::GaussPoint
virtual void GaussPoint(int ig)=0
Element::Update
virtual void Update(Inputs2 *inputs2, int index, IoModel *iomodel, int analysis_counter, int analysis_type, int finite_element)=0
IoModel::FetchDataToInput
void FetchDataToInput(Inputs2 *inputs2, Elements *elements, const char *vector_name, int input_enum)
Definition: IoModel.cpp:1651
Element::VertexConnectivity
virtual int VertexConnectivity(int vertexindex)=0
P2Enum
@ P2Enum
Definition: EnumDefinitions.h:1223
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
SigmaNNEnum
@ SigmaNNEnum
Definition: EnumDefinitions.h:704
Matrix3x3Solve
void Matrix3x3Solve(IssmDouble *X, IssmDouble *A, IssmDouble *B)
Definition: MatrixUtils.cpp:471
Element::GetNumberOfVertices
virtual int GetNumberOfVertices(void)=0
IoModel::domaintype
int domaintype
Definition: IoModel.h:78
ElementMatrix
Definition: ElementMatrix.h:19
AugmentedLagrangianRhopEnum
@ AugmentedLagrangianRhopEnum
Definition: EnumDefinitions.h:39
Element::element_type
int element_type
Definition: Element.h:56
Gauss
Definition: Gauss.h:8
AugmentedLagrangianRholambdaEnum
@ AugmentedLagrangianRholambdaEnum
Definition: EnumDefinitions.h:38
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
Element::NewElementMatrix
ElementMatrix * NewElementMatrix(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2497