Ice Sheet System Model  4.18
Code documentation
SmoothAnalysis.cpp
Go to the documentation of this file.
1 #include "./SmoothAnalysis.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 SmoothAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
9 }/*}}}*/
10 void SmoothAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
11 }/*}}}*/
12 void SmoothAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
13 
15 
16 }/*}}}*/
17 int SmoothAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
18  return 1;
19 }/*}}}*/
20 void SmoothAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
21 
22  /*Update elements: */
23  int counter=0;
24  for(int i=0;i<iomodel->numberofelements;i++){
25  if(iomodel->my_elements[i]){
26  Element* element=(Element*)elements->GetObjectByOffset(counter);
27  element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,P1Enum);
28  counter++;
29  }
30  }
31 }/*}}}*/
32 void SmoothAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
33 }/*}}}*/
34 
35 /*Finite Element Analysis*/
37  _error_("not implemented");
38 }/*}}}*/
40  /*Default, return NULL*/
41  return NULL;
42 }/*}}}*/
44 _error_("Not implemented");
45 }/*}}}*/
47 
48  /* Intermediaries */
49  int domaintype;
50  IssmDouble Jdet,thickness,l;
51  IssmDouble *xyz_list = NULL;
52 
53  /*Check dimension*/
54  element->FindParam(&domaintype,DomainTypeEnum);
55  switch(domaintype){
57  break;
58  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
59  }
60 
61  /*Fetch number of nodes and dof for this finite element*/
62  int numnodes = element->GetNumberOfNodes();
63 
64  /*Initialize Element matrix and vectors*/
65  ElementMatrix* Ke = element->NewElementMatrix();
66  IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);
67  IssmDouble* basis = xNew<IssmDouble>(numnodes);
68 
69  /*Retrieve all inputs and parameters*/
71  element->GetVerticesCoordinates(&xyz_list);
72  Input2* thickness_input = element->GetInput2(ThicknessEnum); _assert_(thickness_input);
73 
74  /* Start looping on the number of gaussian points: */
75  Gauss* gauss=element->NewGauss(2);
76  for(int ig=gauss->begin();ig<gauss->end();ig++){
77  gauss->GaussPoint(ig);
78 
79  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
80  thickness_input->GetInputValue(&thickness,gauss);
81  if(thickness<50.) thickness=50.;
82 
83  element->NodalFunctions(basis,gauss);
84  element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
85 
86  for(int i=0;i<numnodes;i++){
87  for(int j=0;j<numnodes;j++){
88  Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
89  basis[i]*basis[j]
90  +(l*thickness)*(l*thickness)*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])
91  );
92  }
93  }
94  }
95 
96  /*Clean up and return*/
97  delete gauss;
98  xDelete<IssmDouble>(dbasis);
99  xDelete<IssmDouble>(basis);
100  xDelete<IssmDouble>(xyz_list);
101  return Ke;
102 }/*}}}*/
104 
105  /*Get basal element*/
106  int domaintype;
107  element->FindParam(&domaintype,DomainTypeEnum);
108  switch(domaintype){
110  break;
111  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
112  }
113 
114  /*Intermediaries */
115  int input_enum;
116  IssmDouble Jdet,value;
117  IssmDouble *xyz_list = NULL;
118  Input2 *input = NULL;
119 
120  /*SPECIFICS: Driving stress for balance velocities*/
121  Input2* H_input = NULL, *surface_input = NULL, *vx_input = NULL, *vy_input = NULL;
122  IssmDouble taud_x,norms,normv,vx,vy;
123  IssmDouble rho_ice,gravity,slope[2],thickness;
124 
125  /*Fetch number of nodes and dof for this finite element*/
126  int numnodes = element->GetNumberOfNodes();
127 
128  /*Initialize Element vector*/
129  ElementVector* pe = element->NewElementVector();
130  IssmDouble* basis = xNew<IssmDouble>(numnodes);
131 
132  /*Retrieve all inputs and parameters*/
133  element->GetVerticesCoordinates(&xyz_list);
134  element->FindParam(&input_enum,InputToSmoothEnum);
135 
136  switch(input_enum){
137  case DrivingStressXEnum:
138  case DrivingStressYEnum:{
139  rho_ice = element->FindParam(MaterialsRhoIceEnum);
140  gravity = element->FindParam(ConstantsGEnum);
141  H_input = element->GetInput2(ThicknessEnum); _assert_(H_input);
142  surface_input = element->GetInput2(SurfaceEnum); _assert_(surface_input);
143  vx_input = element->GetInput2(VxEnum);
144  vy_input = element->GetInput2(VyEnum);
145  }
146  break;
147  case SurfaceSlopeXEnum:
148  case SurfaceSlopeYEnum:{
149  surface_input = element->GetInput2(SurfaceEnum); _assert_(surface_input);
150  }
151  break;
152  default: input = element->GetInput2(input_enum);
153  }
154 
155  /* Start looping on the number of gaussian points: */
156  Gauss* gauss=element->NewGauss(2);
157  for(int ig=gauss->begin();ig<gauss->end();ig++){
158  gauss->GaussPoint(ig);
159 
160  element->JacobianDeterminant(&Jdet,xyz_list,gauss);
161  element->NodalFunctions(basis,gauss);
162 
163  switch(input_enum){
164  case DrivingStressXEnum:
165  case DrivingStressYEnum:{
166  H_input->GetInputValue(&thickness,gauss);
167  surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
168  if(vx_input && vy_input){
169  vx_input->GetInputValue(&vx,gauss);
170  vy_input->GetInputValue(&vy,gauss);
171  norms = sqrt(slope[0]*slope[0]+slope[1]*slope[1]+1.e-10);
172  normv = sqrt(vx*vx + vy*vy);
173  if(normv>15./(365.*24.*3600.)){
174  slope[0] = -vx/normv*norms;
175  slope[1] = -vy/normv*norms;
176  }
177  }
178  if(input_enum==DrivingStressXEnum)
179  value = rho_ice*gravity*thickness*slope[0];
180  else
181  value = rho_ice*gravity*thickness*slope[1];
182  }
183  break;
184  case SurfaceSlopeXEnum:
185  surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
186  value = slope[0];
187  break;
188  case SurfaceSlopeYEnum:
189  surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
190  value = slope[1];
191  break;
192  default:
193  input->GetInputValue(&value,gauss);
194  }
195 
196  for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*value*basis[i];
197  }
198 
199  /*Clean up and return*/
200  xDelete<IssmDouble>(xyz_list);
201  xDelete<IssmDouble>(basis);
202  delete gauss;
203  return pe;
204 }/*}}}*/
206  _error_("not implemented yet");
207 }/*}}}*/
208 void SmoothAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
209  _error_("Not implemented yet");
210 }/*}}}*/
212  int inputenum,domaintype,elementtype;
213 
214  element->FindParam(&inputenum,InputToSmoothEnum);
215  element->FindParam(&domaintype,DomainTypeEnum);
216  switch(domaintype){
218  element->InputUpdateFromSolutionOneDof(solution,inputenum);
219  break;
220  default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
221  }
222 }/*}}}*/
224  /*Default, do nothing*/
225  return;
226 }/*}}}*/
SmoothAnalysis.h
: header file for generic external result object
_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
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
SurfaceSlopeXEnum
@ SurfaceSlopeXEnum
Definition: EnumDefinitions.h:829
DrivingStressYEnum
@ DrivingStressYEnum
Definition: EnumDefinitions.h:539
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
Constraints
Declaration of Constraints class.
Definition: Constraints.h:13
Elements
Declaration of Elements class.
Definition: Elements.h:17
DrivingStressXEnum
@ DrivingStressXEnum
Definition: EnumDefinitions.h:538
MaterialsRhoIceEnum
@ MaterialsRhoIceEnum
Definition: EnumDefinitions.h:264
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
SmoothAnalysis::GradientJ
void GradientJ(Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
Definition: SmoothAnalysis.cpp:208
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
SmoothAnalysisEnum
@ SmoothAnalysisEnum
Definition: EnumDefinitions.h:1276
SmoothAnalysis::GetSolutionFromInputs
void GetSolutionFromInputs(Vector< IssmDouble > *solution, Element *element)
Definition: SmoothAnalysis.cpp:205
SmoothAnalysis::DofsPerNode
int DofsPerNode(int **doflist, int domaintype, int approximation)
Definition: SmoothAnalysis.cpp:17
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
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
SmoothAnalysis::InputUpdateFromSolution
void InputUpdateFromSolution(IssmDouble *solution, Element *element)
Definition: SmoothAnalysis.cpp:211
Element::NewElementVector
ElementVector * NewElementVector(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2505
IoModel::numberofelements
int numberofelements
Definition: IoModel.h:96
SmoothAnalysis::Core
void Core(FemModel *femmodel)
Definition: SmoothAnalysis.cpp:36
ConstantsGEnum
@ ConstantsGEnum
Definition: EnumDefinitions.h:102
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
InputToSmoothEnum
@ InputToSmoothEnum
Definition: EnumDefinitions.h:207
SmoothThicknessMultiplierEnum
@ SmoothThicknessMultiplierEnum
Definition: EnumDefinitions.h:397
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
SurfaceEnum
@ SurfaceEnum
Definition: EnumDefinitions.h:823
SmoothAnalysis::CreateDVector
ElementVector * CreateDVector(Element *element)
Definition: SmoothAnalysis.cpp:39
Inputs2
Declaration of Inputs class.
Definition: Inputs2.h:23
Input2
Definition: Input2.h:18
FemModel
Definition: FemModel.h:31
SmoothAnalysis::CreateNodes
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
Definition: SmoothAnalysis.cpp:12
Loads
Declaration of Loads class.
Definition: Loads.h:16
_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::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
SmoothAnalysis::CreateConstraints
void CreateConstraints(Constraints *constraints, IoModel *iomodel)
Definition: SmoothAnalysis.cpp:8
ThicknessEnum
@ ThicknessEnum
Definition: EnumDefinitions.h:840
SmoothAnalysis::CreateJacobianMatrix
ElementMatrix * CreateJacobianMatrix(Element *element)
Definition: SmoothAnalysis.cpp:43
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
SmoothAnalysis::UpdateElements
void UpdateElements(Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
Definition: SmoothAnalysis.cpp:20
IoModel
Definition: IoModel.h:48
SmoothAnalysis::UpdateConstraints
void UpdateConstraints(FemModel *femmodel)
Definition: SmoothAnalysis.cpp:223
SmoothAnalysis::CreateLoads
void CreateLoads(Loads *loads, IoModel *iomodel)
Definition: SmoothAnalysis.cpp:10
ElementMatrix
Definition: ElementMatrix.h:19
Vector< IssmDouble >
SmoothAnalysis::CreateKMatrix
ElementMatrix * CreateKMatrix(Element *element)
Definition: SmoothAnalysis.cpp:46
SmoothAnalysis::UpdateParameters
void UpdateParameters(Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
Definition: SmoothAnalysis.cpp:32
Gauss
Definition: Gauss.h:8
Element::NodalFunctionsDerivatives
virtual void NodalFunctionsDerivatives(IssmDouble *dbasis, IssmDouble *xyz_list, Gauss *gauss)=0
SmoothAnalysis::CreatePVector
ElementVector * CreatePVector(Element *element)
Definition: SmoothAnalysis.cpp:103
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
SurfaceSlopeYEnum
@ SurfaceSlopeYEnum
Definition: EnumDefinitions.h:830
Element::NewElementMatrix
ElementMatrix * NewElementMatrix(int approximation_enum=NoneApproximationEnum)
Definition: Element.cpp:2497
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16