source: issm/trunk/src/c/analyses/SamplingAnalysis.cpp

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

merged trunk-jpl and trunk for revision 27230
File size: 14.4 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.sampling.tau",SamplingTauEnum);
79 iomodel->FetchDataToInput(inputs,elements,"md.initialization.sample",SampleEnum,0.);
80 iomodel->FetchDataToInput(inputs,elements,"md.initialization.sample",SampleNoiseEnum,0.);
81 if(iomodel->solution_enum==TransientSolutionEnum) iomodel->FetchDataToInput(inputs,elements,"md.sampling.phi",SamplingPhiEnum,0.);
82
83}/*}}}*/
84void SamplingAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
85
86 int numoutputs;
87 char** requestedoutputs = NULL;
88
89 parameters->AddObject(iomodel->CopyConstantObject("md.sampling.alpha",SamplingAlphaEnum));
90 parameters->AddObject(iomodel->CopyConstantObject("md.sampling.robin",SamplingRobinEnum));
91 parameters->AddObject(iomodel->CopyConstantObject("md.sampling.seed",SamplingSeedEnum));
92
93 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.sampling.requested_outputs");
94 parameters->AddObject(new IntParam(SamplingNumRequestedOutputsEnum,numoutputs));
95 if(numoutputs)parameters->AddObject(new StringArrayParam(SamplingRequestedOutputsEnum,requestedoutputs,numoutputs));
96 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.sampling.requested_outputs");
97
98}/*}}}*/
99
100/*Finite Element Analysis*/
101void SamplingAnalysis::Core(FemModel* femmodel){/*{{{*/
102 _error_("not implemented");
103}/*}}}*/
104void SamplingAnalysis::PreCore(FemModel* femmodel){/*{{{*/
105 _error_("not implemented");
106}/*}}}*/
107ElementVector* SamplingAnalysis::CreateDVector(Element* element){/*{{{*/
108 /*Default, return NULL*/
109 return NULL;
110}/*}}}*/
111ElementMatrix* SamplingAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
112_error_("not implemented");
113}/*}}}*/
114ElementMatrix* SamplingAnalysis::CreateKMatrix(Element* element){/*{{{*/
115
116 /* Check if ice in element */
117 //if(!element->IsIceInElement()) return NULL; // Ice in element not required for sampling capability
118
119 /*Intermediaries*/
120 bool robin;
121 int domaintype;
122
123 /*compute all stiffness matrices for this element*/
124
125 ElementMatrix* Ke = NULL;
126 element->FindParam(&robin,SamplingRobinEnum);
127
128 if(!robin) Ke=CreateKMatrixModifiedHelmholtz(element);
129 else{
130 ElementMatrix* Ke1=CreateKMatrixModifiedHelmholtz(element);
131 ElementMatrix* Ke2=CreateKMatrixRobinBC(element);
132 Ke =new ElementMatrix(Ke1,Ke2);
133
134 delete Ke1;
135 delete Ke2;
136 }
137
138 /*clean-up and return*/
139 return Ke;
140
141}/*}}}*/
142ElementMatrix* SamplingAnalysis::CreateKMatrixModifiedHelmholtz(Element* element){/*{{{*/
143
144 /* Check if ice in element */
145 //if(!element->IsIceInElement()) return NULL;
146
147 /*Intermediaries*/
148 int dim;
149 IssmDouble D,Jdet;
150 IssmDouble* xyz_list = NULL;
151 IssmDouble kappa;
152
153 /*Get problem dimension*/
154 dim=2; // Only 2D horizontal domain is implemented so far
155
156 /*Fetch number of nodes and dof for this finite element*/
157 int numnodes = element->GetNumberOfNodes();
158
159 /*Initialize Element vector and other vectors*/
160 ElementMatrix* Ke = element->NewElementMatrix();
161 IssmDouble* basis = xNew<IssmDouble>(numnodes);
162 IssmDouble* dbasis = xNew<IssmDouble>(dim*numnodes);
163
164 /*Retrieve all inputs and parameters*/
165 element->GetVerticesCoordinates(&xyz_list);
166 Input* kappa_input=element->GetInput(SamplingKappaEnum); _assert_(kappa_input);
167
168 /* Start looping on the number of gaussian points: */
169 Gauss* gauss=element->NewGauss(2);
170 while(gauss->next()){
171
172 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
173 element->NodalFunctions(basis,gauss);
174 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
175
176 kappa_input->GetInputValue(&kappa,gauss); _assert_(kappa>0);
177
178 D=gauss->weight*Jdet;
179
180 /* Laplacian operator */
181 for(int i=0;i<numnodes;i++){
182 for(int j=0;j<numnodes;j++){
183 Ke->values[i*numnodes+j] += D*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
184 }
185 }
186
187 /* Identity identity */
188 for(int i=0;i<numnodes;i++){
189 for(int j=0;j<numnodes;j++){
190 Ke->values[i*numnodes+j] += D*kappa*kappa*(basis[j]*basis[i]);
191 }
192 }
193
194 }
195
196 /*Clean up and return*/
197 xDelete<IssmDouble>(xyz_list);
198 xDelete<IssmDouble>(basis);
199 xDelete<IssmDouble>(dbasis);
200 delete gauss;
201 return Ke;
202
203}/*}}}*/
204ElementMatrix* SamplingAnalysis::CreateKMatrixRobinBC(Element* element){/*{{{*/
205
206 /* Check if ice in element */
207 //if(!element->IsIceInElement()) return NULL;
208
209 /*If no boundary, return NULL*/
210 if(!IsFaceOnBoundary(element)) return NULL;
211
212 /*Intermediaries*/
213 IssmDouble D,Jdet;
214 IssmDouble* xyz_list = NULL;
215 IssmDouble* xyz_list_boundary = NULL;
216 IssmDouble beta;
217
218 /*Fetch number of nodes and dof for this finite element*/
219 int numnodes = element->GetNumberOfNodes();
220
221 /*Initialize Element vector and other vectors*/
222 ElementMatrix* Ke = element->NewElementMatrix();
223 IssmDouble* basis = xNew<IssmDouble>(numnodes);
224
225 /*Retrieve all inputs and parameters*/
226 element->GetVerticesCoordinates(&xyz_list);
227 element->GetLevelCoordinates(&xyz_list_boundary,xyz_list,MeshVertexonboundaryEnum,1.);
228 Input* beta_input=element->GetInput(SamplingBetaEnum); _assert_(beta_input);
229
230 /* Start looping on the number of gaussian points: */
231 Gauss* gauss=element->NewGauss(xyz_list,xyz_list_boundary,3);
232 while(gauss->next()){
233
234 element->JacobianDeterminantSurface(&Jdet,xyz_list_boundary,gauss);
235 element->NodalFunctions(basis,gauss);
236
237 beta_input->GetInputValue(&beta,gauss); _assert_(beta>=0);
238
239 D=gauss->weight*Jdet;
240
241 for(int i=0;i<numnodes;i++){
242 for(int j=0;j<numnodes;j++){
243 Ke->values[i*numnodes+j] += D*beta*(basis[j]*basis[i]);
244 }
245 }
246
247 }
248
249 /*Clean up and return*/
250 xDelete<IssmDouble>(xyz_list);
251 xDelete<IssmDouble>(xyz_list_boundary);
252 xDelete<IssmDouble>(basis);
253 delete gauss;
254 return Ke;
255
256}/*}}}*/
257ElementVector* SamplingAnalysis::CreatePVector(Element* element){/*{{{*/
258 //_error_("not supported");
259 return NULL;
260}/*}}}*/
261void SamplingAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
262 element->GetSolutionFromInputsOneDof(solution,SampleEnum);
263}/*}}}*/
264void SamplingAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
265 _error_("not supported");
266}/*}}}*/
267void SamplingAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
268
269 /*Fetch number of nodes and dof for this finite element*/
270 int numnodes = element->GetNumberOfNodes();
271
272 /*Fetch dof list and allocate solution vector*/
273 int* doflist=NULL;
274 element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
275 IssmDouble* newsample = xNew<IssmDouble>(numnodes);
276 IssmDouble* tau = xNew<IssmDouble>(numnodes);
277 element->GetInputListOnNodes(&tau[0],SamplingTauEnum);
278
279 /*Use the dof list to index into the solution vector: */
280 for(int i=0;i<numnodes;i++){
281 newsample[i]=solution[doflist[i]] / tau[i]; // new
282
283 /*Check solution*/
284 if(xIsNan<IssmDouble>(newsample[i])) _error_("NaN found in solution vector");
285 if(xIsInf<IssmDouble>(newsample[i])) _error_("Inf found in solution vector");
286 }
287
288 /*Add sample inputs to the tria element: */
289 element->AddInput(SampleEnum,newsample,element->GetElementType());
290
291 /*Free resources:*/
292 xDelete<IssmDouble>(newsample);
293 xDelete<IssmDouble>(tau);
294 xDelete<int>(doflist);
295 }/*}}}*/
296void SamplingAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
297 /*Default, do nothing*/
298 return;
299}/*}}}*/
300
301ElementMatrix* SamplingAnalysis::CreateMassMatrix(Element* element){/*{{{*/
302
303 /* Check if ice in element */
304 //if(!element->IsIceInElement()) return NULL;
305
306 /*Intermediaries*/
307 IssmDouble D,Jdet;
308 IssmDouble* xyz_list = NULL;
309
310 /*Fetch number of nodes and dof for this finite element*/
311 int numnodes = element->GetNumberOfNodes();
312
313 /*Initialize Element vector and other vectors*/
314 ElementMatrix* Me = element->NewElementMatrix();
315 IssmDouble* basis = xNew<IssmDouble>(numnodes);
316
317 /*Retrieve all inputs and parameters*/
318 element->GetVerticesCoordinates(&xyz_list);
319
320 /* Start looping on the number of gaussian points: */
321 Gauss* gauss=element->NewGauss(2);
322 while(gauss->next()){
323
324 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
325 element->NodalFunctions(basis,gauss);
326
327 D=gauss->weight*Jdet;
328 TripleMultiply(basis,1,numnodes,1,
329 &D,1,1,0,
330 basis,1,numnodes,0,
331 &Me->values[0],1);
332 }
333
334 /*Clean up and return*/
335 xDelete<IssmDouble>(xyz_list);
336 xDelete<IssmDouble>(basis);
337 delete gauss;
338 return Me;
339}/*}}}*/
340void SamplingAnalysis::LumpedKMatrix(Vector<IssmDouble>** pKlff,FemModel* femmodel){/*{{{*/
341
342 /*Initialize Lumped mass matrix (actually we just save its diagonal)*/
343 int fsize = femmodel->nodes->NumberOfDofs(FsetEnum);
344 int flocalsize = femmodel->nodes->NumberOfDofsLocal(FsetEnum);
345 Vector<IssmDouble>* Klff = new Vector<IssmDouble>(flocalsize,fsize);
346
347 /*Create and assemble matrix*/
348 for(int i=0;i<femmodel->elements->Size();i++){
349 Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
350 ElementMatrix* KLe = this->CreateKMatrix(element);
351 if(KLe){
352 KLe->Lump();
353 KLe->AddDiagonalToGlobal(Klff);
354 }
355 delete KLe;
356 }
357 Klff->Assemble();
358
359 /*Assign output pointer*/
360 *pKlff=Klff;
361}/*}}}*/
362void SamplingAnalysis::KMatrix(Matrix<IssmDouble>** pKff,FemModel* femmodel){/*{{{*/
363
364 /*Initialize K matrix*/
365 Matrix<IssmDouble> *Kff = NULL;
366 AllocateSystemMatricesx(&Kff,NULL,NULL,NULL,femmodel);
367
368 /*Create and assemble matrix*/
369 for(int i=0;i<femmodel->elements->Size();i++){
370 Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
371 ElementMatrix* Ke = this->CreateKMatrix(element);
372 if(Ke){
373 Ke->AddToGlobal(Kff);
374 }
375 delete Ke;
376 }
377 Kff->Assemble();
378
379 /*Assign output pointer*/
380 *pKff=Kff;
381
382}/*}}}*/
383void SamplingAnalysis::LumpedMassMatrix(Vector<IssmDouble>** pMlff,FemModel* femmodel){/*{{{*/
384
385 /*Initialize Lumped mass matrix (actually we just save its diagonal)*/
386 int fsize = femmodel->nodes->NumberOfDofs(FsetEnum);
387 int flocalsize = femmodel->nodes->NumberOfDofsLocal(FsetEnum);
388 Vector<IssmDouble>* Mlff = new Vector<IssmDouble>(flocalsize,fsize);
389
390 /*Create and assemble matrix*/
391 for(int i=0;i<femmodel->elements->Size();i++){
392 Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
393 ElementMatrix* MLe = this->CreateMassMatrix(element);
394 if(MLe){
395 MLe->Lump();
396 MLe->AddDiagonalToGlobal(Mlff);
397 }
398 delete MLe;
399 }
400 Mlff->Assemble();
401
402 /*Assign output pointer*/
403 *pMlff=Mlff;
404}/*}}}*/
405void SamplingAnalysis::MassMatrix(Matrix<IssmDouble>** pMff,FemModel* femmodel){/*{{{*/
406
407 /*Initialize Mass matrix*/
408 Matrix<IssmDouble> *Mff = NULL;
409 AllocateSystemMatricesx(&Mff,NULL,NULL,NULL,femmodel);
410
411 /*Create and assemble matrix*/
412 for(int i=0;i<femmodel->elements->Size();i++){
413 Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
414 ElementMatrix* MLe = this->CreateMassMatrix(element);
415 if(MLe){
416 MLe->AddToGlobal(Mff);
417 }
418 delete MLe;
419 }
420 Mff->Assemble();
421
422 /*Assign output pointer*/
423 *pMff=Mff;
424}/*}}}*/
425void SamplingAnalysis::UpdateTransientSample(FemModel * femmodel){
426
427 for(int j=0;j<femmodel->elements->Size();j++){
428 Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
429 UpdateTransientSample(element);
430 }
431
432}
433void SamplingAnalysis::UpdateTransientSample(Element * element){
434
435 /*Intermediaries */
436 IssmDouble phi, sample, noise;
437
438 /*Fetch number vertices for this element*/
439 int numvertices = element->GetNumberOfVertices();
440
441 /*Initialize new sample*/
442 IssmDouble* sample_new = xNew<IssmDouble>(numvertices);
443
444 /*Retrieve all inputs and parameters*/
445 Input* sample_input=element->GetInput(SampleOldEnum); _assert_(sample_input);
446 Input* noise_input=element->GetInput(SampleNoiseEnum); _assert_(noise_input);
447 Input* phi_input=element->GetInput(SamplingPhiEnum); _assert_(phi_input);
448
449 /* Start looping on the number of gaussian points: */
450 Gauss* gauss=element->NewGauss();
451 for(int iv=0;iv<numvertices;iv++){
452 gauss->GaussVertex(iv);
453
454 /*Get input values at gauss points*/
455 sample_input->GetInputValue(&sample,gauss);
456 noise_input->GetInputValue(&noise,gauss);
457 phi_input->GetInputValue(&phi,gauss);
458
459 /*Get new sample*/
460 sample_new[iv] = phi*sample + noise;
461
462 }
463
464 element->AddInput(SampleEnum,sample_new,element->GetElementType());
465
466 /*Clean up and return*/
467 xDelete<IssmDouble>(sample_new);
468 delete gauss;
469
470}
Note: See TracBrowser for help on using the repository browser.