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

#include <ExtrapolationAnalysis.h>

Inheritance diagram for ExtrapolationAnalysis:
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)
 
int GetExtrapolationCase (Element *element)
 
void SetConstraintsOnIce (Element *element)
 
void UpdateConstraints (FemModel *femmodel)
 
- Public Member Functions inherited from Analysis
virtual ~Analysis ()
 

Detailed Description

Definition at line 11 of file ExtrapolationAnalysis.h.

Member Function Documentation

◆ CreateConstraints()

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

Implements Analysis.

Definition at line 8 of file ExtrapolationAnalysis.cpp.

8  {/*{{{*/
9  // do nothing for now
10  return;
11 }

◆ CreateLoads()

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

Implements Analysis.

Definition at line 13 of file ExtrapolationAnalysis.cpp.

13  {/*{{{*/
14  // do nothing for now
15  return;
16 }/*}}}*/

◆ CreateNodes()

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

Implements Analysis.

Definition at line 17 of file ExtrapolationAnalysis.cpp.

17  {/*{{{*/
18  int finiteelement=P1Enum;
19  if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
20  ::CreateNodes(nodes,iomodel,ExtrapolationAnalysisEnum,finiteelement);
21  iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
22 }

◆ DofsPerNode()

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

Implements Analysis.

Definition at line 24 of file ExtrapolationAnalysis.cpp.

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

◆ UpdateElements()

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

Implements Analysis.

Definition at line 28 of file ExtrapolationAnalysis.cpp.

28  {/*{{{*/
29  int finiteelement;
30 
31  /*Finite element type*/
32  finiteelement = P1Enum;
33 
34  /*Update elements: */
35  int counter=0;
36  for(int i=0;i<iomodel->numberofelements;i++){
37  if(iomodel->my_elements[i]){
38  Element* element=(Element*)elements->GetObjectByOffset(counter);
39  element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,finiteelement);
40  counter++;
41  }
42  }
43  if(iomodel->domaintype!=Domain2DhorizontalEnum){
44  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
45  iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
46  }
47 }

◆ UpdateParameters()

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

Implements Analysis.

Definition at line 49 of file ExtrapolationAnalysis.cpp.

49  {/*{{{*/
50  //do nothing for now
51  return;
52 }

◆ Core()

void ExtrapolationAnalysis::Core ( FemModel *  femmodel)
virtual

Implements Analysis.

Definition at line 56 of file ExtrapolationAnalysis.cpp.

56  {/*{{{*/
57 
58  /* Intermediaries */
59  bool save_results;
60  int extvar_enum;
61  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
62  femmodel->parameters->FindParam(&extvar_enum, ExtrapolationVariableEnum);
63 
64  /*activate formulation: */
65  femmodel->SetCurrentConfiguration(ExtrapolationAnalysisEnum);
66 
67  if(VerboseSolution()) _printf0_(" extrapolation of " << EnumToStringx(extvar_enum) << ":\n");
68  solutionsequence_linear(femmodel);
69 
70  if(save_results){
71  if(VerboseSolution()) _printf0_(" saving results\n");
72  femmodel->RequestedOutputsx(&femmodel->results,&extvar_enum,1);
73  }
74 }/*}}}*/

◆ CreateDVector()

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

Implements Analysis.

Definition at line 75 of file ExtrapolationAnalysis.cpp.

75  {/*{{{*/
76  /*Default, return NULL*/
77  return NULL;
78 }/*}}}*/

◆ CreateJacobianMatrix()

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

Implements Analysis.

Definition at line 79 of file ExtrapolationAnalysis.cpp.

79  {/*{{{*/
80  /* Jacobian required for the Newton solver */
81  _error_("not implemented yet");
82 }/*}}}*/

◆ CreateKMatrix()

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

Implements Analysis.

Definition at line 83 of file ExtrapolationAnalysis.cpp.

83  {/*{{{*/
84 
85  /*Intermediaries */
86  bool extrapolatebydiffusion = true;
87  int dim, domaintype, extrapolationcase;
88  int i,row,col,stabilization;
89  IssmDouble Jdet,D_scalar,h;
90  IssmDouble norm_dlsf;
91  IssmDouble hx,hy,hz,kappa;
92  IssmDouble *xyz_list = NULL;
93  Element *workelement = NULL;
94 
95  /*Get problem case*/
96  extrapolationcase=GetExtrapolationCase(element);
97  switch(extrapolationcase){
98  case 0:
99  if(!element->IsOnBase()) return NULL;
100  workelement = element->SpawnBasalElement();
101  break;
102  case 1: case 2: case 3: workelement=element; break;
103  }
104 
105  /* get extrapolation dimension */
106  workelement->FindParam(&domaintype,DomainTypeEnum);
107  switch(domaintype){
108  case Domain2DverticalEnum: dim=1; break;
109  case Domain2DhorizontalEnum: dim=2; break;
110  case Domain3DEnum: dim=3; break;
111  default: _error_("not supported yet");
112  }
113 
114  /*Fetch number of nodes and dof for this finite element*/
115  int numnodes = workelement->GetNumberOfNodes();
116 
117  /*Initialize Element vector and other vectors*/
118  ElementMatrix *Ke = workelement->NewElementMatrix();
119  IssmDouble *basis = xNew<IssmDouble>(numnodes);
120  IssmDouble *dbasis = xNew<IssmDouble>(dim*numnodes);
121  IssmDouble dlsf[3];
122  IssmDouble normal[3];
123 
124  /*Retrieve all inputs and parameters*/
125  Input2* lsf_slopex_input=workelement->GetInput2(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
126  Input2* lsf_slopey_input=workelement->GetInput2(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
127  workelement->GetVerticesCoordinates(&xyz_list);
128 
129  /* In 3d we are going to extrude using horizontal diffusion. Since the layers are 3d,
130  * we change the geometry of the element to make it flat. That way, the diffusion is
131  * handled along the layers
132  */
133  if(element->ObjectEnum()==PentaEnum){
134  for(i=0;i<3;i++) xyz_list[3*i+2] = 0.;
135  for(i=3;i<6;i++) xyz_list[3*i+2] = 1.;
136  }
137 
138  /* Start looping on the number of gaussian points: */
139  Gauss* gauss=workelement->NewGauss(2);
140  for(int ig=gauss->begin();ig<gauss->end();ig++){
141  gauss->GaussPoint(ig);
142 
143  workelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
144  workelement->NodalFunctions(basis,gauss);
145  workelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
146 
147  D_scalar=gauss->weight*Jdet;
148 
149  if(extrapolatebydiffusion){
150  for(int i=0;i<numnodes;i++){
151  for(int j=0;j<numnodes;j++){
152  Ke->values[i*numnodes+j] += D_scalar*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
153  }
154  }
155  }
156  else{
157  /* extrapolate values along normal */
158  /* Get normal on ice boundary */
159  lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
160  if(dim>1)
161  lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
162  if(dim>2)
163  dlsf[2]=0.;
164  norm_dlsf=0.;
165  for(i=0;i<dim;i++) norm_dlsf+=dlsf[i]*dlsf[i];
166  norm_dlsf=sqrt(norm_dlsf); _assert_(norm_dlsf>0.);
167 
168  if(norm_dlsf>0.)
169  for(i=0;i<dim;i++) normal[i]=dlsf[i]/norm_dlsf;
170  else
171  for(i=0;i<dim;i++) normal[i]=0.;
172 
173  for(int i=0;i<numnodes;i++){
174  for(int j=0;j<numnodes;j++){
175  Ke->values[i*numnodes+j] += D_scalar*(normal[0]*dbasis[0*numnodes+j]*basis[i] + normal[1]*dbasis[1*numnodes+j]*basis[i]);
176  }
177  }
178 
179  /* stabilization */
180  /* do not use streamline upwinding for extrapolation: it yields oscillating results due to diffusion along normal direction, but none across */
181  stabilization=1;
182  if (stabilization==0){/* no stabilization, do nothing*/}
183  else if(stabilization==1){
184  /* Artificial Diffusion */
185  workelement->ElementSizes(&hx,&hy,&hz);
186  h=sqrt(pow(hx*normal[0],2) + pow(hy*normal[1],2));
187  kappa=h/2.+1.e-14;
188 
189  for(int i=0;i<numnodes;i++){
190  for(int j=0;j<numnodes;j++){
191  Ke->values[i*numnodes+j] += D_scalar*kappa*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
192  }
193  }
194  }
195  }
196  }
197 
198  /*Clean up and return*/
199  xDelete<IssmDouble>(xyz_list);
200  xDelete<IssmDouble>(basis);
201  xDelete<IssmDouble>(dbasis);
202  delete gauss;
203  if(extrapolationcase==0){workelement->DeleteMaterials(); delete workelement;};
204  return Ke;
205 
206 }/*}}}*/

◆ CreatePVector()

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

Implements Analysis.

Definition at line 207 of file ExtrapolationAnalysis.cpp.

207  {/*{{{*/
208  return NULL;
209 }/*}}}*/

◆ GetSolutionFromInputs()

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

Implements Analysis.

Definition at line 210 of file ExtrapolationAnalysis.cpp.

210  {/*{{{*/
211  _error_("not implemented yet");
212 }/*}}}*/

◆ GradientJ()

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

Implements Analysis.

Definition at line 213 of file ExtrapolationAnalysis.cpp.

213  {/*{{{*/
214  _error_("Not implemented yet");
215 }/*}}}*/

◆ InputUpdateFromSolution()

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

Implements Analysis.

Definition at line 216 of file ExtrapolationAnalysis.cpp.

216  {/*{{{*/
217 
218  int extrapolationvariable, extrapolationcase;
219  extrapolationcase=GetExtrapolationCase(element);
220  element->FindParam(&extrapolationvariable, ExtrapolationVariableEnum);
221  switch(extrapolationcase){
222  case 0:
223  element->InputUpdateFromSolutionOneDof(solution,extrapolationvariable);
224  break;
225  case 1:
226  element->InputUpdateFromSolutionOneDof(solution,extrapolationvariable);
227  break;
228  case 2:
229  element->InputUpdateFromSolutionOneDofCollapsed(solution,extrapolationvariable);
230  break;
231  case 3:
232  element->InputUpdateFromSolutionOneDof(solution,extrapolationvariable);
233  break;
234  }
235 }/*}}}*/

◆ GetExtrapolationCase()

int ExtrapolationAnalysis::GetExtrapolationCase ( Element *  element)

Definition at line 236 of file ExtrapolationAnalysis.cpp.

236  {/*{{{*/
237 
238  /* Get case of extrapolation, depending on domain quality, and extrapolation variable */
239  int domaintype, extrapolationvariable;
240  int extrapolationcase;
241 
242  element->FindParam(&domaintype,DomainTypeEnum);
243  switch(domaintype){
244  case Domain2DverticalEnum: extrapolationcase=0; break;
245  case Domain2DhorizontalEnum: extrapolationcase=1; break;
246  case Domain3DEnum:
247  element->FindParam(&extrapolationvariable, ExtrapolationVariableEnum);
248  if(extrapolationvariable==ThicknessEnum){
249  extrapolationcase=2; // scalar fields that are constant along z-axis
250  }
251  else{
252  extrapolationcase=3; // scalar fields that vary along z-axis
253  }
254  break;
255  }
256  return extrapolationcase;
257 }/*}}}*/

◆ SetConstraintsOnIce()

void ExtrapolationAnalysis::SetConstraintsOnIce ( Element *  element)

Definition at line 258 of file ExtrapolationAnalysis.cpp.

258  {/*{{{*/
259 
260  int numnodes=element->GetNumberOfNodes();
261 
262  /* Intermediaries */
263  int extvar_enum;
264  IssmDouble active,value;
265  Node* node = NULL;
266 
267  /* Get parameters */
268  element->FindParam(&extvar_enum, ExtrapolationVariableEnum);
269 
270  Input2* active_input=element->GetInput2(IceMaskNodeActivationEnum); _assert_(active_input);
271  Input2* extvar_input=element->GetInput2(extvar_enum); _assert_(extvar_input);
272 
273  Gauss* gauss=element->NewGauss();
274  for(int in=0;in<numnodes;in++){
275  gauss->GaussNode(element->GetElementType(),in);
276  node=element->GetNode(in);
277  active_input->GetInputValue(&active,gauss);
278  if(node->IsActive()){
279  if(active>0.5){
280  /* if ice, set dirichlet BC */
281  extvar_input->GetInputValue(&value,gauss);
282  node->ApplyConstraint(0,value);
283  }
284  else {
285  /* no ice, set no spc */
286  node->DofInFSet(0);
287  }
288  }
289  }
290  delete gauss;
291 }/*}}}*/

◆ UpdateConstraints()

void ExtrapolationAnalysis::UpdateConstraints ( FemModel *  femmodel)
virtual

Implements Analysis.

Definition at line 292 of file ExtrapolationAnalysis.cpp.

292  {/*{{{*/
293 
294  for(int i=0;i<femmodel->elements->Size();i++){
295  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
296  this->SetConstraintsOnIce(element);
297  }
298 }/*}}}*/
The documentation for this class was generated from the following files:
DataSet::Size
int Size()
Definition: DataSet.cpp:399
Node::IsActive
bool IsActive(void)
Definition: Node.cpp:795
SaveResultsEnum
@ SaveResultsEnum
Definition: EnumDefinitions.h:302
Element::GetElementType
virtual int GetElementType(void)=0
_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
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
LevelsetfunctionSlopeXEnum
@ LevelsetfunctionSlopeXEnum
Definition: EnumDefinitions.h:635
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
LevelsetfunctionSlopeYEnum
@ LevelsetfunctionSlopeYEnum
Definition: EnumDefinitions.h:636
Element::SpawnBasalElement
virtual Element * SpawnBasalElement(void)=0
MeshVertexonbaseEnum
@ MeshVertexonbaseEnum
Definition: EnumDefinitions.h:653
FemModel::results
Results * results
Definition: FemModel.h:48
IceMaskNodeActivationEnum
@ IceMaskNodeActivationEnum
Definition: EnumDefinitions.h:627
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
Gauss::GaussNode
virtual void GaussNode(int finitelement, int iv)=0
solutionsequence_linear
void solutionsequence_linear(FemModel *femmodel)
Definition: solutionsequence_linear.cpp:10
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
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
Element::ElementSizes
virtual void ElementSizes(IssmDouble *phx, IssmDouble *phy, IssmDouble *phz)=0
ExtrapolationAnalysisEnum
@ ExtrapolationAnalysisEnum
Definition: EnumDefinitions.h:1057
IoModel::DeleteData
void DeleteData(int num,...)
Definition: IoModel.cpp:500
IoModel::numberofelements
int numberofelements
Definition: IoModel.h:96
ExtrapolationVariableEnum
@ ExtrapolationVariableEnum
Definition: EnumDefinitions.h:135
DomainTypeEnum
@ DomainTypeEnum
Definition: EnumDefinitions.h:124
Element::NewGauss
virtual Gauss * NewGauss(void)=0
ExtrapolationAnalysis::SetConstraintsOnIce
void SetConstraintsOnIce(Element *element)
Definition: ExtrapolationAnalysis.cpp:258
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
Element::GetNode
Node * GetNode(int nodeindex)
Definition: Element.cpp:1207
ExtrapolationAnalysis::GetExtrapolationCase
int GetExtrapolationCase(Element *element)
Definition: ExtrapolationAnalysis.cpp:236
IoModel::FetchData
void FetchData(bool *pboolean, const char *data_name)
Definition: IoModel.cpp:933
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
FemModel::elements
Elements * elements
Definition: FemModel.h:44
Object::ObjectEnum
virtual int ObjectEnum()=0
Input2
Definition: Input2.h:18
Node
Definition: Node.h:23
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Node::DofInFSet
void DofInFSet(int dof)
Definition: Node.cpp:694
Gauss::begin
virtual int begin(void)=0
VerboseSolution
bool VerboseSolution(void)
Definition: Verbosity.cpp:24
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
MeshVertexonsurfaceEnum
@ MeshVertexonsurfaceEnum
Definition: EnumDefinitions.h:655
PentaEnum
@ PentaEnum
Definition: EnumDefinitions.h:1231
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
FemModel::RequestedOutputsx
void RequestedOutputsx(Results **presults, char **requested_outputs, int numoutputs, bool save_results=true)
Definition: FemModel.cpp:2267
FemModel::SetCurrentConfiguration
void SetCurrentConfiguration(int configuration_type)
Definition: FemModel.cpp:634
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
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
Node::ApplyConstraint
void ApplyConstraint(int dof, IssmDouble value)
Definition: Node.cpp:646
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
ExtrapolationAnalysis::CreateNodes
void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)
Definition: ExtrapolationAnalysis.cpp:17
Domain2DverticalEnum
@ Domain2DverticalEnum
Definition: EnumDefinitions.h:535
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
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:27 for Ice Sheet System Model by   doxygen 1.8.19