source: issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/L2ProjectionEPLAnalysis.cpp

Last change on this file was 20690, checked in by Mathieu Morlighem, 10 years ago

NEW: marhsall strings instead of enums
File size: 8.4 KB
Line 
1#include "./L2ProjectionEPLAnalysis.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 L2ProjectionEPLAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
9
10 /*No constraints*/
11}/*}}}*/
12void L2ProjectionEPLAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
13
14 /*No loads*/
15}/*}}}*/
16void L2ProjectionEPLAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
17 /*Now, do we really want DC?*/
18 int hydrology_model;
19 iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
20 if(hydrology_model!=HydrologydcEnum) return;
21
22 /*Do we want an efficient layer*/
23 bool isefficientlayer;
24 iomodel->FindConstant(&isefficientlayer,"md.hydrology.isefficientlayer");
25 if(!isefficientlayer) return;
26
27 if(iomodel->domaintype==Domain3DEnum){
28 iomodel->FetchData(1,"md.mesh.vertexonbase");
29 }
30 else if(iomodel->domaintype==Domain2DverticalEnum){
31 iomodel->FetchData(1,"md.mesh.vertexonbase");
32 }
33 ::CreateNodes(nodes,iomodel,L2ProjectionEPLAnalysisEnum,P1Enum);
34 iomodel->DeleteData(1,"md.mesh.vertexonbase");
35}/*}}}*/
36int L2ProjectionEPLAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
37 return 1;
38}/*}}}*/
39void L2ProjectionEPLAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
40
41 bool isefficientlayer;
42 int hydrology_model;
43
44 /*Now, do we really want DC?*/
45 iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
46 if(hydrology_model!=HydrologydcEnum) return;
47
48 /*Do we want an efficient layer*/
49 iomodel->FindConstant(&isefficientlayer,"md.hydrology.isefficientlayer");
50 if(!isefficientlayer) return;
51
52 /*Update elements: */
53 int counter=0;
54 for(int i=0;i<iomodel->numberofelements;i++){
55 if(iomodel->my_elements[i]){
56 Element* element=(Element*)elements->GetObjectByOffset(counter);
57 element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
58 counter++;
59 }
60 }
61
62 iomodel->FetchDataToInput(elements,"md.initialization.epl_head",EplHeadEnum);
63 iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
64 if(iomodel->domaintype!=Domain2DhorizontalEnum){
65 iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
66 iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
67 }
68}/*}}}*/
69void L2ProjectionEPLAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
70}/*}}}*/
71
72/*Finite Element Analysis*/
73void L2ProjectionEPLAnalysis::Core(FemModel* femmodel){/*{{{*/
74 _error_("not implemented");
75}/*}}}*/
76ElementVector* L2ProjectionEPLAnalysis::CreateDVector(Element* element){/*{{{*/
77 /*Default, return NULL*/
78 return NULL;
79}/*}}}*/
80ElementMatrix* L2ProjectionEPLAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
81_error_("Not implemented");
82}/*}}}*/
83ElementMatrix* L2ProjectionEPLAnalysis::CreateKMatrix(Element* element){/*{{{*/
84
85 /*Intermediaries*/
86 int domaintype;
87 bool active_element;
88 Element* basalelement;
89
90 /*Get basal element*/
91 element->FindParam(&domaintype,DomainTypeEnum);
92 switch(domaintype){
93 case Domain2DhorizontalEnum:
94 basalelement = element;
95 break;
96 case Domain2DverticalEnum:
97 if(!element->IsOnBase()) return NULL;
98 basalelement = element->SpawnBasalElement();
99 break;
100 case Domain3DEnum:
101 if(!element->IsOnBase()) return NULL;
102 basalelement = element->SpawnBasalElement();
103 break;
104 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
105 }
106
107 Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
108 active_element_input->GetInputValue(&active_element);
109
110 /* Check that all nodes are active, else return empty matrix */
111 if(!active_element){
112 if(domaintype!=Domain2DhorizontalEnum){
113 basalelement->DeleteMaterials();
114 delete basalelement;
115 }
116 return NULL;
117 }
118
119 /*Intermediaries */
120 IssmDouble D,Jdet;
121 IssmDouble *xyz_list = NULL;
122
123 /*Fetch number of nodes and dof for this finite element*/
124 int numnodes = basalelement->GetNumberOfNodes();
125
126 /*Initialize Element vector*/
127 ElementMatrix* Ke = basalelement->NewElementMatrix(NoneApproximationEnum);
128 IssmDouble* basis = xNew<IssmDouble>(numnodes);
129
130 /*Retrieve all inputs and parameters*/
131 basalelement->GetVerticesCoordinates(&xyz_list);
132
133 /* Start looping on the number of gaussian points: */
134 Gauss* gauss=basalelement->NewGauss(2);
135 for(int ig=gauss->begin();ig<gauss->end();ig++){
136 gauss->GaussPoint(ig);
137
138 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
139 basalelement->NodalFunctions(basis,gauss);
140 D=gauss->weight*Jdet;
141
142 TripleMultiply(basis,1,numnodes,1,
143 &D,1,1,0,
144 basis,1,numnodes,0,
145 &Ke->values[0],1);
146 }
147
148 /*Clean up and return*/
149 xDelete<IssmDouble>(xyz_list);
150 xDelete<IssmDouble>(basis);
151 delete gauss;
152 if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
153 return Ke;
154}/*}}}*/
155ElementVector* L2ProjectionEPLAnalysis::CreatePVector(Element* element){/*{{{*/
156
157 /*Intermediaries*/
158 int domaintype;
159 bool active_element;
160 Element* basalelement;
161
162 /*Get basal element*/
163 element->FindParam(&domaintype,DomainTypeEnum);
164 switch(domaintype){
165 case Domain2DhorizontalEnum:
166 basalelement = element;
167 break;
168 case Domain3DEnum:
169 if(!element->IsOnBase()) return NULL;
170 basalelement = element->SpawnBasalElement();
171 break;
172 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
173 }
174
175 Input* active_element_input = basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
176 active_element_input->GetInputValue(&active_element);
177
178 /*Check that all nodes are active, else return empty matrix*/
179 if(!active_element) {
180 if(domaintype!=Domain2DhorizontalEnum){
181 basalelement->DeleteMaterials();
182 delete basalelement;
183 }
184 return NULL;
185 }
186
187 /*Intermediaries */
188 int input_enum,index;
189 IssmDouble Jdet,slopes[2];
190 Input *input = NULL;
191 IssmDouble *xyz_list = NULL;
192
193 /*Fetch number of nodes and dof for this finite element*/
194 int numnodes = basalelement->GetNumberOfNodes();
195
196 /*Initialize Element vector*/
197 ElementVector* pe = basalelement->NewElementVector();
198 IssmDouble* basis = xNew<IssmDouble>(numnodes);
199
200 /*Retrieve all inputs and parameters*/
201 basalelement->GetVerticesCoordinates(&xyz_list);
202 basalelement->FindParam(&input_enum,InputToL2ProjectEnum);
203 switch(input_enum){
204 case EplHeadSlopeXEnum: input = basalelement->GetInput(EplHeadEnum); index = 0; _assert_(input); break;
205 case EplHeadSlopeYEnum: input = basalelement->GetInput(EplHeadEnum); index = 1; _assert_(input); break;
206 default: _error_("not implemented");
207 }
208
209 /* Start looping on the number of gaussian points: */
210 Gauss* gauss=basalelement->NewGauss(2);
211 for(int ig=gauss->begin();ig<gauss->end();ig++){
212 gauss->GaussPoint(ig);
213
214 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
215 basalelement->NodalFunctions(basis,gauss);
216
217 input->GetInputDerivativeValue(&slopes[0],xyz_list,gauss);
218 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*slopes[index]*basis[i];
219 }
220
221 /*Clean up and return*/
222 xDelete<IssmDouble>(xyz_list);
223 xDelete<IssmDouble>(basis);
224 delete gauss;
225 if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
226 return pe;
227}/*}}}*/
228void L2ProjectionEPLAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
229 _error_("not implemented yet");
230}/*}}}*/
231void L2ProjectionEPLAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
232 _error_("Not implemented yet");
233}/*}}}*/
234void L2ProjectionEPLAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
235 int inputenum,domaintype;
236
237 element->FindParam(&inputenum,InputToL2ProjectEnum);
238 element->FindParam(&domaintype,DomainTypeEnum);
239 switch(domaintype){
240 case Domain2DhorizontalEnum:
241 element->InputUpdateFromSolutionOneDof(solution,inputenum);
242 break;
243 case Domain2DverticalEnum:
244 element->InputUpdateFromSolutionOneDof(solution,inputenum);
245 break;
246 case Domain3DEnum:
247 element->InputUpdateFromSolutionOneDofCollapsed(solution,inputenum);
248 break;
249 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
250 }
251}/*}}}*/
252void L2ProjectionEPLAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
253 /*Default, do nothing*/
254 return;
255}/*}}}*/
256
Note: See TracBrowser for help on using the repository browser.