source: issm/trunk/src/c/analyses/Balancethickness2Analysis.cpp@ 25836

Last change on this file since 25836 was 25836, checked in by Mathieu Morlighem, 4 years ago

merged trunk-jpl and trunk for revision 25834
File size: 6.8 KB
Line 
1#include "./Balancethickness2Analysis.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*/
8void Balancethickness2Analysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
9
10 int finiteelement = P1Enum;
11 IoModelToConstraintsx(constraints,iomodel,"md.balancethickness.spcthickness",Balancethickness2AnalysisEnum,finiteelement);
12
13}/*}}}*/
14void Balancethickness2Analysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
15
16}/*}}}*/
17void Balancethickness2Analysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
18
19 int finiteelement = P1Enum;
20 ::CreateNodes(nodes,iomodel,Balancethickness2AnalysisEnum,finiteelement);
21}/*}}}*/
22int Balancethickness2Analysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
23 return 1;
24}/*}}}*/
25void Balancethickness2Analysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
26
27 /*Finite element type*/
28 int finiteelement = P1Enum;
29
30 /*Load variables in element*/
31 iomodel->FetchDataToInput(inputs,elements,"md.geometry.thickness",ThicknessEnum);
32 iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
33 iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
34 iomodel->FetchDataToInput(inputs,elements,"md.solidearth.initialsealevel",SealevelEnum,0);
35 iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
36 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum);
37 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum);
38 iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
39 iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
40 iomodel->FetchDataToInput(inputs,elements,"md.balancethickness.thickening_rate",BalancethicknessThickeningRateEnum);
41
42 /*Update elements: */
43 int counter=0;
44 for(int i=0;i<iomodel->numberofelements;i++){
45 if(iomodel->my_elements[i]){
46 Element* element=(Element*)elements->GetObjectByOffset(counter);
47 element->Update(inputs,i,iomodel,analysis_counter,analysis_type,finiteelement);
48
49 counter++;
50 }
51 }
52
53}/*}}}*/
54void Balancethickness2Analysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
55}/*}}}*/
56
57/*Finite Element Analysis*/
58void Balancethickness2Analysis::Core(FemModel* femmodel){/*{{{*/
59 _error_("not implemented yet");
60}/*}}}*/
61ElementVector* Balancethickness2Analysis::CreateDVector(Element* element){/*{{{*/
62 /*Default, return NULL*/
63 return NULL;
64}/*}}}*/
65ElementMatrix* Balancethickness2Analysis::CreateJacobianMatrix(Element* element){/*{{{*/
66_error_("Not implemented");
67}/*}}}*/
68ElementMatrix* Balancethickness2Analysis::CreateKMatrix(Element* element){/*{{{*/
69
70 /*Intermediaries */
71 IssmDouble yts = 365*24*3600.;
72 IssmDouble Jdet,vx,vy,vel;
73 IssmDouble* xyz_list = NULL;
74
75 /*Fetch number of nodes and dof for this finite element*/
76 int numnodes = element->GetNumberOfNodes();
77
78 /*Initialize Element vector and other vectors*/
79 ElementMatrix* Ke = element->NewElementMatrix();
80 IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);
81
82 /*Retrieve all inputs and parameters*/
83 element->GetVerticesCoordinates(&xyz_list);
84 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
85 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
86
87 /*Get element characteristic length*/
88 IssmDouble h = element->CharacteristicLength();
89
90 /* Start looping on the number of gaussian points: */
91 Gauss* gauss=element->NewGauss(2);
92 while(gauss->next()){
93 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
94 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
95 vx_input->GetInputValue(&vx,gauss);
96 vy_input->GetInputValue(&vy,gauss);
97
98 /*make sure are diffusivisty is large enough*/
99 vel = sqrt(vx*vx+vy*vy);
100 if(sqrt(vx*vx+vy*vy)==0.){
101 vx = 0.1/yts;
102 vy = 0.1/yts;
103 vel = sqrt(vx*vx+vy*vy);
104 }
105 else if(vel<30./yts){
106 vx = 0.;//vx/vel*0.01;
107 vy = 0.;//vy/vel*0.01;
108 vel = sqrt(vx*vx+vy*vy);
109 vel = 30./yts*500000.;
110 }
111
112 for(int i=0;i<numnodes;i++){
113 for(int j=0;j<numnodes;j++){
114 Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
115 (vx*dbasis[0*numnodes+i] + vy*dbasis[1*numnodes+i])*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j])
116 + vel/500000.*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]));
117 }
118 }
119 }
120
121 /*Clean up and return*/
122 xDelete<IssmDouble>(xyz_list);
123 xDelete<IssmDouble>(dbasis);
124 delete gauss;
125 return Ke;
126}/*}}}*/
127ElementVector* Balancethickness2Analysis::CreatePVector(Element* element){/*{{{*/
128
129 return NULL;
130 /*Intermediaries */
131 IssmDouble dhdt[2],mb[2],ms[2],Jdet;
132 IssmDouble* xyz_list = NULL;
133
134 /*Fetch number of nodes and dof for this finite element*/
135 int numnodes = element->GetNumberOfNodes();
136
137 /*Initialize Element vector and other vectors*/
138 ElementVector* pe = element->NewElementVector();
139 IssmDouble* basis = xNew<IssmDouble>(numnodes);
140
141 /*Retrieve all inputs and parameters*/
142 element->GetVerticesCoordinates(&xyz_list);
143 Input* ms_input = element->GetInput(SmbMassBalanceEnum); _assert_(ms_input);
144 Input* mb_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(mb_input);
145 Input* dhdt_input = element->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
146
147 /* Start looping on the number of gaussian points: */
148 Gauss* gauss=element->NewGauss(2);
149 while(gauss->next()){
150
151 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
152 element->NodalFunctions(basis,gauss);
153
154 ms_input->GetInputDerivativeValue(&ms[0],xyz_list,gauss);
155
156 ms_input->GetInputDerivativeValue(&ms[0],xyz_list,gauss);
157 mb_input->GetInputDerivativeValue(&mb[0],xyz_list,gauss);
158 dhdt_input->GetInputDerivativeValue(&dhdt[0],xyz_list,gauss);
159
160 for(int i=0;i<numnodes;i++) pe->values[i]+=0*Jdet*gauss->weight*(
161 (ms[0]+ms[1]-mb[0]-mb[1]-dhdt[0]-dhdt[1])*basis[i]
162 );
163 }
164
165 /*Clean up and return*/
166 xDelete<IssmDouble>(xyz_list);
167 xDelete<IssmDouble>(basis);
168 delete gauss;
169 return pe;
170}/*}}}*/
171void Balancethickness2Analysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
172 element->GetSolutionFromInputsOneDof(solution,ThicknessEnum);
173}/*}}}*/
174void Balancethickness2Analysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
175 _error_("Not implemented yet");
176}/*}}}*/
177void Balancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
178
179 element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
180
181}/*}}}*/
182void Balancethickness2Analysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
183 /*Default, do nothing*/
184 return;
185}/*}}}*/
Note: See TracBrowser for help on using the repository browser.