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

Last change on this file since 26744 was 26744, checked in by Mathieu Morlighem, 3 years ago

merged trunk-jpl and trunk for revision 26742
File size: 6.9 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.initialization.sealevel",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}/*}}}*/
61void Balancethickness2Analysis::PreCore(FemModel* femmodel){/*{{{*/
62 _error_("not implemented yet");
63}/*}}}*/
64ElementVector* Balancethickness2Analysis::CreateDVector(Element* element){/*{{{*/
65 /*Default, return NULL*/
66 return NULL;
67}/*}}}*/
68ElementMatrix* Balancethickness2Analysis::CreateJacobianMatrix(Element* element){/*{{{*/
69_error_("Not implemented");
70}/*}}}*/
71ElementMatrix* Balancethickness2Analysis::CreateKMatrix(Element* element){/*{{{*/
72
73 /*Intermediaries */
74 IssmDouble yts = 365*24*3600.;
75 IssmDouble Jdet,vx,vy,vel;
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 vector and other vectors*/
82 ElementMatrix* Ke = element->NewElementMatrix();
83 IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);
84
85 /*Retrieve all inputs and parameters*/
86 element->GetVerticesCoordinates(&xyz_list);
87 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
88 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
89
90 /*Get element characteristic length*/
91 IssmDouble h = element->CharacteristicLength();
92
93 /* Start looping on the number of gaussian points: */
94 Gauss* gauss=element->NewGauss(2);
95 while(gauss->next()){
96 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
97 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
98 vx_input->GetInputValue(&vx,gauss);
99 vy_input->GetInputValue(&vy,gauss);
100
101 /*make sure are diffusivisty is large enough*/
102 vel = sqrt(vx*vx+vy*vy);
103 if(sqrt(vx*vx+vy*vy)==0.){
104 vx = 0.1/yts;
105 vy = 0.1/yts;
106 vel = sqrt(vx*vx+vy*vy);
107 }
108 else if(vel<30./yts){
109 vx = 0.;//vx/vel*0.01;
110 vy = 0.;//vy/vel*0.01;
111 vel = sqrt(vx*vx+vy*vy);
112 vel = 30./yts*500000.;
113 }
114
115 for(int i=0;i<numnodes;i++){
116 for(int j=0;j<numnodes;j++){
117 Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
118 (vx*dbasis[0*numnodes+i] + vy*dbasis[1*numnodes+i])*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j])
119 + vel/500000.*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]));
120 }
121 }
122 }
123
124 /*Clean up and return*/
125 xDelete<IssmDouble>(xyz_list);
126 xDelete<IssmDouble>(dbasis);
127 delete gauss;
128 return Ke;
129}/*}}}*/
130ElementVector* Balancethickness2Analysis::CreatePVector(Element* element){/*{{{*/
131
132 return NULL;
133 /*Intermediaries */
134 IssmDouble dhdt[2],mb[2],ms[2],Jdet;
135 IssmDouble* xyz_list = NULL;
136
137 /*Fetch number of nodes and dof for this finite element*/
138 int numnodes = element->GetNumberOfNodes();
139
140 /*Initialize Element vector and other vectors*/
141 ElementVector* pe = element->NewElementVector();
142 IssmDouble* basis = xNew<IssmDouble>(numnodes);
143
144 /*Retrieve all inputs and parameters*/
145 element->GetVerticesCoordinates(&xyz_list);
146 Input* ms_input = element->GetInput(SmbMassBalanceEnum); _assert_(ms_input);
147 Input* mb_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(mb_input);
148 Input* dhdt_input = element->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
149
150 /* Start looping on the number of gaussian points: */
151 Gauss* gauss=element->NewGauss(2);
152 while(gauss->next()){
153
154 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
155 element->NodalFunctions(basis,gauss);
156
157 ms_input->GetInputDerivativeValue(&ms[0],xyz_list,gauss);
158
159 ms_input->GetInputDerivativeValue(&ms[0],xyz_list,gauss);
160 mb_input->GetInputDerivativeValue(&mb[0],xyz_list,gauss);
161 dhdt_input->GetInputDerivativeValue(&dhdt[0],xyz_list,gauss);
162
163 for(int i=0;i<numnodes;i++) pe->values[i]+=0*Jdet*gauss->weight*(
164 (ms[0]+ms[1]-mb[0]-mb[1]-dhdt[0]-dhdt[1])*basis[i]
165 );
166 }
167
168 /*Clean up and return*/
169 xDelete<IssmDouble>(xyz_list);
170 xDelete<IssmDouble>(basis);
171 delete gauss;
172 return pe;
173}/*}}}*/
174void Balancethickness2Analysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
175 element->GetSolutionFromInputsOneDof(solution,ThicknessEnum);
176}/*}}}*/
177void Balancethickness2Analysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
178 _error_("Not implemented yet");
179}/*}}}*/
180void Balancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
181
182 element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
183
184}/*}}}*/
185void Balancethickness2Analysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
186 /*Default, do nothing*/
187 return;
188}/*}}}*/
Note: See TracBrowser for help on using the repository browser.