source: issm/trunk/src/c/analyses/SamplingAnalysis.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: 12.8 KB
Line 
1#include "./SamplingAnalysis.h" //
2#include "../toolkits/toolkits.h"
3#include "../classes/classes.h"
4#include "../shared/shared.h"
5#include "../modules/modules.h"
6#include <random>
7
8#define FINITEELEMENT P1Enum
9#define NUMVERTICES 3
10
11/* Reimplementation of IsFaceOnBoundary with two edges on boundary supported */
12bool IsFaceOnBoundary(Element* element)
13{
14
15 IssmDouble values[NUMVERTICES];
16 IssmDouble sum;
17
18 /*Retrieve all inputs and parameters*/
19 element->GetInputListOnVertices(&values[0],MeshVertexonboundaryEnum);
20 sum = values[0]+values[1]+values[2];
21
22 _assert_(sum==0. || sum==1. || sum==2. || sum==3.);
23
24 if(sum>1.){
25 return true;
26 }
27 else{
28 return false;
29 }
30}
31
32/*Model processing*/
33void SamplingAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
34 /*Default, do nothing*/
35 return;
36}/*}}}*/
37void SamplingAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
38 /*Default, do nothing*/
39 return;
40}/*}}}*/
41void SamplingAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
42
43 int finiteelement;
44 finiteelement = FINITEELEMENT;
45
46 /*Check in 2d*/
47 if(iomodel->domaintype!=Domain2DhorizontalEnum) _error_("Only 2D horizontal domain is implemented.");
48
49 ::CreateNodes(nodes,iomodel,SamplingAnalysisEnum,finiteelement);
50
51}/*}}}*/
52int SamplingAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
53 return 1;
54}/*}}}*/
55void SamplingAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
56
57 int finiteelement;
58
59 /*Finite element type*/
60 finiteelement = FINITEELEMENT;
61
62 /*Update elements: */
63 int counter=0;
64 for(int i=0;i<iomodel->numberofelements;i++){
65 if(iomodel->my_elements[i]){
66 Element* element=(Element*)elements->GetObjectByOffset(counter);
67 element->Update(inputs,i,iomodel,analysis_counter,analysis_type,finiteelement);
68 counter++;
69 }
70 }
71
72 /*Create inputs: */
73 iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
74 iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonboundary",MeshVertexonboundaryEnum);
75
76 iomodel->FetchDataToInput(inputs,elements,"md.sampling.kappa",SamplingKappaEnum);
77 iomodel->FetchDataToInput(inputs,elements,"md.sampling.beta",SamplingBetaEnum,0.);
78 iomodel->FetchDataToInput(inputs,elements,"md.initialization.sample",SampleEnum,0.);
79
80}/*}}}*/
81void SamplingAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
82
83 int numoutputs;
84 char** requestedoutputs = NULL;
85
86 parameters->AddObject(iomodel->CopyConstantObject("md.sampling.tau",SamplingTauEnum));
87 parameters->AddObject(iomodel->CopyConstantObject("md.sampling.alpha",SamplingAlphaEnum));
88 parameters->AddObject(iomodel->CopyConstantObject("md.sampling.robin",SamplingRobinEnum));
89 if(solution_enum==TransientSolutionEnum) parameters->AddObject(iomodel->CopyConstantObject("md.sampling.phi",SamplingPhiEnum));
90 parameters->AddObject(iomodel->CopyConstantObject("md.sampling.seed",SamplingSeedEnum));
91
92 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.sampling.requested_outputs");
93 parameters->AddObject(new IntParam(SamplingNumRequestedOutputsEnum,numoutputs));
94 if(numoutputs)parameters->AddObject(new StringArrayParam(SamplingRequestedOutputsEnum,requestedoutputs,numoutputs));
95 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.sampling.requested_outputs");
96
97}/*}}}*/
98
99/*Finite Element Analysis*/
100void SamplingAnalysis::Core(FemModel* femmodel){/*{{{*/
101 _error_("not implemented");
102}/*}}}*/
103void SamplingAnalysis::PreCore(FemModel* femmodel){/*{{{*/
104 _error_("not implemented");
105}/*}}}*/
106ElementVector* SamplingAnalysis::CreateDVector(Element* element){/*{{{*/
107 /*Default, return NULL*/
108 return NULL;
109}/*}}}*/
110ElementMatrix* SamplingAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
111_error_("not implemented");
112}/*}}}*/
113ElementMatrix* SamplingAnalysis::CreateKMatrix(Element* element){/*{{{*/
114
115 /* Check if ice in element */
116 //if(!element->IsIceInElement()) return NULL; // Ice in element not required for sampling capability
117
118 /*Intermediaries*/
119 bool robin;
120 int domaintype;
121
122 /*compute all stiffness matrices for this element*/
123
124 ElementMatrix* Ke = NULL;
125 element->FindParam(&robin,SamplingRobinEnum);
126
127 if(!robin) Ke=CreateKMatrixModifiedHelmholtz(element);
128 else{
129 ElementMatrix* Ke1=CreateKMatrixModifiedHelmholtz(element);
130 ElementMatrix* Ke2=CreateKMatrixRobinBC(element);
131 Ke =new ElementMatrix(Ke1,Ke2);
132
133 delete Ke1;
134 delete Ke2;
135 }
136
137 /*clean-up and return*/
138 return Ke;
139
140}/*}}}*/
141ElementMatrix* SamplingAnalysis::CreateKMatrixModifiedHelmholtz(Element* element){/*{{{*/
142
143 /* Check if ice in element */
144 //if(!element->IsIceInElement()) return NULL;
145
146 /*Intermediaries*/
147 int dim;
148 IssmDouble D,Jdet;
149 IssmDouble* xyz_list = NULL;
150 IssmDouble kappa;
151
152 /*Get problem dimension*/
153 dim=2; // Only 2D horizontal domain is implemented so far
154
155 /*Fetch number of nodes and dof for this finite element*/
156 int numnodes = element->GetNumberOfNodes();
157
158 /*Initialize Element vector and other vectors*/
159 ElementMatrix* Ke = element->NewElementMatrix();
160 IssmDouble* basis = xNew<IssmDouble>(numnodes);
161 IssmDouble* dbasis = xNew<IssmDouble>(dim*numnodes);
162
163 /*Retrieve all inputs and parameters*/
164 element->GetVerticesCoordinates(&xyz_list);
165 Input* kappa_input=element->GetInput(SamplingKappaEnum); _assert_(kappa_input);
166
167 /* Start looping on the number of gaussian points: */
168 Gauss* gauss=element->NewGauss(2);
169 while(gauss->next()){
170
171 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
172 element->NodalFunctions(basis,gauss);
173 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
174
175 kappa_input->GetInputValue(&kappa,gauss); _assert_(kappa>0);
176
177 D=gauss->weight*Jdet;
178
179 /* Laplacian operator */
180 for(int i=0;i<numnodes;i++){
181 for(int j=0;j<numnodes;j++){
182 Ke->values[i*numnodes+j] += D*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
183 }
184 }
185
186 /* Identity identity */
187 for(int i=0;i<numnodes;i++){
188 for(int j=0;j<numnodes;j++){
189 Ke->values[i*numnodes+j] += D*kappa*kappa*(basis[j]*basis[i]);
190 }
191 }
192
193 }
194
195 /*Clean up and return*/
196 xDelete<IssmDouble>(xyz_list);
197 xDelete<IssmDouble>(basis);
198 xDelete<IssmDouble>(dbasis);
199 delete gauss;
200 return Ke;
201
202}/*}}}*/
203ElementMatrix* SamplingAnalysis::CreateKMatrixRobinBC(Element* element){/*{{{*/
204
205 /* Check if ice in element */
206 //if(!element->IsIceInElement()) return NULL;
207
208 /*If no boundary, return NULL*/
209 if(!IsFaceOnBoundary(element)) return NULL;
210
211 /*Intermediaries*/
212 IssmDouble D,Jdet;
213 IssmDouble* xyz_list = NULL;
214 IssmDouble* xyz_list_boundary = NULL;
215 IssmDouble beta;
216
217 /*Fetch number of nodes and dof for this finite element*/
218 int numnodes = element->GetNumberOfNodes();
219
220 /*Initialize Element vector and other vectors*/
221 ElementMatrix* Ke = element->NewElementMatrix();
222 IssmDouble* basis = xNew<IssmDouble>(numnodes);
223
224 /*Retrieve all inputs and parameters*/
225 element->GetVerticesCoordinates(&xyz_list);
226 element->GetLevelCoordinates(&xyz_list_boundary,xyz_list,MeshVertexonboundaryEnum,1.);
227 Input* beta_input=element->GetInput(SamplingBetaEnum); _assert_(beta_input);
228
229 /* Start looping on the number of gaussian points: */
230 Gauss* gauss=element->NewGauss(xyz_list,xyz_list_boundary,3);
231 while(gauss->next()){
232
233 element->JacobianDeterminantSurface(&Jdet,xyz_list_boundary,gauss);
234 element->NodalFunctions(basis,gauss);
235
236 beta_input->GetInputValue(&beta,gauss); _assert_(beta>=0);
237
238 D=gauss->weight*Jdet;
239
240 for(int i=0;i<numnodes;i++){
241 for(int j=0;j<numnodes;j++){
242 Ke->values[i*numnodes+j] += D*beta*(basis[j]*basis[i]);
243 }
244 }
245
246 }
247
248 /*Clean up and return*/
249 xDelete<IssmDouble>(xyz_list);
250 xDelete<IssmDouble>(xyz_list_boundary);
251 xDelete<IssmDouble>(basis);
252 delete gauss;
253 return Ke;
254
255}/*}}}*/
256ElementVector* SamplingAnalysis::CreatePVector(Element* element){/*{{{*/
257 //_error_("not supported");
258 return NULL;
259}/*}}}*/
260void SamplingAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
261 element->GetSolutionFromInputsOneDof(solution,SampleEnum);
262}/*}}}*/
263void SamplingAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
264 _error_("not supported");
265}/*}}}*/
266void SamplingAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
267
268 /*Fetch number of nodes and dof for this finite element*/
269 int numnodes = element->GetNumberOfNodes();
270
271 /*Fetch dof list and allocate solution vector*/
272 int* doflist=NULL;
273 element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
274 IssmDouble* newsample = xNew<IssmDouble>(numnodes);
275
276 /*Use the dof list to index into the solution vector: */
277 for(int i=0;i<numnodes;i++){
278 newsample[i]=solution[doflist[i]];
279 /*Check solution*/
280 if(xIsNan<IssmDouble>(newsample[i])) _error_("NaN found in solution vector");
281 if(xIsInf<IssmDouble>(newsample[i])) _error_("Inf found in solution vector");
282 }
283
284 /*Add sample inputs to the tria element: */
285 element->AddInput(SampleEnum,newsample,element->GetElementType());
286
287 /*Free ressources:*/
288 xDelete<IssmDouble>(newsample);
289 xDelete<int>(doflist);
290 }/*}}}*/
291void SamplingAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
292 /*Default, do nothing*/
293 return;
294}/*}}}*/
295
296ElementMatrix* SamplingAnalysis::CreateMassMatrix(Element* element){/*{{{*/
297
298 /* Check if ice in element */
299 //if(!element->IsIceInElement()) return NULL;
300
301 /*Intermediaries*/
302 IssmDouble D,Jdet;
303 IssmDouble* xyz_list = NULL;
304
305 /*Fetch number of nodes and dof for this finite element*/
306 int numnodes = element->GetNumberOfNodes();
307
308 /*Initialize Element vector and other vectors*/
309 ElementMatrix* Me = element->NewElementMatrix();
310 IssmDouble* basis = xNew<IssmDouble>(numnodes);
311
312 /*Retrieve all inputs and parameters*/
313 element->GetVerticesCoordinates(&xyz_list);
314
315 /* Start looping on the number of gaussian points: */
316 Gauss* gauss=element->NewGauss(2);
317 while(gauss->next()){
318
319 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
320 element->NodalFunctions(basis,gauss);
321
322 D=gauss->weight*Jdet;
323 TripleMultiply(basis,1,numnodes,1,
324 &D,1,1,0,
325 basis,1,numnodes,0,
326 &Me->values[0],1);
327 }
328
329 /*Clean up and return*/
330 xDelete<IssmDouble>(xyz_list);
331 xDelete<IssmDouble>(basis);
332 delete gauss;
333 return Me;
334}/*}}}*/
335void SamplingAnalysis::LumpedKMatrix(Vector<IssmDouble>** pKlff,FemModel* femmodel){/*{{{*/
336
337 /*Initialize Lumped mass matrix (actually we just save its diagonal)*/
338 int fsize = femmodel->nodes->NumberOfDofs(FsetEnum);
339 int flocalsize = femmodel->nodes->NumberOfDofsLocal(FsetEnum);
340 Vector<IssmDouble>* Klff = new Vector<IssmDouble>(flocalsize,fsize);
341
342 /*Create and assemble matrix*/
343 for(int i=0;i<femmodel->elements->Size();i++){
344 Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
345 ElementMatrix* KLe = this->CreateKMatrix(element);
346 if(KLe){
347 KLe->Lump();
348 KLe->AddDiagonalToGlobal(Klff);
349 }
350 delete KLe;
351 }
352 Klff->Assemble();
353
354 /*Assign output pointer*/
355 *pKlff=Klff;
356}/*}}}*/
357void SamplingAnalysis::KMatrix(Matrix<IssmDouble>** pKff,FemModel* femmodel){/*{{{*/
358
359 /*Initialize K matrix*/
360 Matrix<IssmDouble> *Kff = NULL;
361 AllocateSystemMatricesx(&Kff,NULL,NULL,NULL,femmodel);
362
363 /*Create and assemble matrix*/
364 for(int i=0;i<femmodel->elements->Size();i++){
365 Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
366 ElementMatrix* Ke = this->CreateKMatrix(element);
367 if(Ke){
368 Ke->AddToGlobal(Kff);
369 }
370 delete Ke;
371 }
372 Kff->Assemble();
373
374 /*Assign output pointer*/
375 *pKff=Kff;
376
377}/*}}}*/
378void SamplingAnalysis::LumpedMassMatrix(Vector<IssmDouble>** pMlff,FemModel* femmodel){/*{{{*/
379
380 /*Initialize Lumped mass matrix (actually we just save its diagonal)*/
381 int fsize = femmodel->nodes->NumberOfDofs(FsetEnum);
382 int flocalsize = femmodel->nodes->NumberOfDofsLocal(FsetEnum);
383 Vector<IssmDouble>* Mlff = new Vector<IssmDouble>(flocalsize,fsize);
384
385 /*Create and assemble matrix*/
386 for(int i=0;i<femmodel->elements->Size();i++){
387 Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
388 ElementMatrix* MLe = this->CreateMassMatrix(element);
389 if(MLe){
390 MLe->Lump();
391 MLe->AddDiagonalToGlobal(Mlff);
392 }
393 delete MLe;
394 }
395 Mlff->Assemble();
396
397 /*Assign output pointer*/
398 *pMlff=Mlff;
399}/*}}}*/
400void SamplingAnalysis::MassMatrix(Matrix<IssmDouble>** pMff,FemModel* femmodel){/*{{{*/
401
402 /*Initialize Mass matrix*/
403 Matrix<IssmDouble> *Mff = NULL;
404 AllocateSystemMatricesx(&Mff,NULL,NULL,NULL,femmodel);
405
406 /*Create and assemble matrix*/
407 for(int i=0;i<femmodel->elements->Size();i++){
408 Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
409 ElementMatrix* MLe = this->CreateMassMatrix(element);
410 if(MLe){
411 MLe->AddToGlobal(Mff);
412 }
413 delete MLe;
414 }
415 Mff->Assemble();
416
417 /*Assign output pointer*/
418 *pMff=Mff;
419}/*}}}*/
Note: See TracBrowser for help on using the repository browser.