source: issm/trunk-jpl/src/c/classes/Loads/Icefront.cpp@ 15423

Last change on this file since 15423 was 15423, checked in by Mathieu Morlighem, 12 years ago

CHG: dim can now be retrieved directly in iomodel
File size: 29.6 KB
Line 
1/*!\file Icefront.c
2 * \brief: implementation of the Icefront object
3 */
4
5/*Headers:*/
6/*{{{*/
7#ifdef HAVE_CONFIG_H
8 #include <config.h>
9#else
10#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
11#endif
12
13#include "../classes.h"
14#include "shared/shared.h"
15/*}}}*/
16
17/*Load macros*/
18#define NUMVERTICESSEG 2
19#define NUMVERTICESQUA 4
20
21/*Icefront constructors and destructor*/
22/*FUNCTION Icefront::Icefront() {{{*/
23Icefront::Icefront(){
24
25 this->inputs=NULL;
26 this->parameters=NULL;
27
28 this->hnodes=NULL;
29 this->nodes= NULL;
30 this->hvertices=NULL;
31 this->vertices= NULL;
32 this->helement=NULL;
33 this->element= NULL;
34 this->hmatpar=NULL;
35 this->matpar= NULL;
36}
37/*}}}*/
38/*FUNCTION Icefront::Icefront(int id, int i, IoModel* iomodel,int analysis_type) {{{*/
39Icefront::Icefront(int icefront_id,int i, IoModel* iomodel,int in_icefront_type, int in_analysis_type){
40
41 int segment_width;
42 int element;
43 int numnodes;
44 int numvertices;
45 int numberofelements;
46
47 /*icefront constructor data: */
48 int icefront_eid;
49 int icefront_mparid;
50 int icefront_node_ids[NUMVERTICESQUA]; //initialize with largest size
51 int icefront_vertex_ids[NUMVERTICESQUA]; //initialize with largest size
52 int icefront_fill;
53
54 /*find parameters: */
55 iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);
56
57 /*First, retrieve element index and element type: */
58 if(iomodel->dim==2){
59 segment_width=4;
60 }
61 else{
62 segment_width=6;
63 }
64 _assert_(iomodel->Data(DiagnosticIcefrontEnum));
65 element=reCast<int,IssmDouble>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+segment_width-2)-1); //element is in the penultimate column (node1 node2 ... elem fill)
66
67 /*Build ids for hook constructors: */
68 icefront_eid=reCast<int,IssmDouble>( *(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+segment_width-2)); //matlab indexing
69 icefront_mparid=numberofelements+1; //matlab indexing
70
71 if (in_icefront_type==MacAyeal2dIceFrontEnum || in_icefront_type==MacAyeal3dIceFrontEnum){
72 icefront_node_ids[0]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
73 icefront_node_ids[1]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
74 icefront_vertex_ids[0]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
75 icefront_vertex_ids[1]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
76 }
77 else if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){
78 icefront_node_ids[0]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
79 icefront_node_ids[1]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
80 icefront_node_ids[2]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2));
81 icefront_node_ids[3]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3));
82 icefront_vertex_ids[0]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
83 icefront_vertex_ids[1]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
84 icefront_vertex_ids[2]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2));
85 icefront_vertex_ids[3]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3));
86 }
87 else _error_("in_icefront_type " << EnumToStringx(in_icefront_type) << " not supported yet!");
88
89 if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){
90 numnodes=4;
91 numvertices=4;
92 }
93 else{
94 numnodes=2;
95 numvertices=2;
96 }
97
98 /*Fill*/
99 icefront_fill=reCast<int>(iomodel->Data(DiagnosticIcefrontEnum)[segment_width*i+segment_width-1]);
100
101 /*Ok, we have everything to build the object: */
102 this->id=icefront_id;
103 this->analysis_type=in_analysis_type;
104
105 /*Hooks: */
106 this->hnodes=new Hook(icefront_node_ids,numnodes);
107 this->hvertices=new Hook(icefront_vertex_ids,numvertices);
108 this->helement=new Hook(&icefront_eid,1);
109 this->hmatpar=new Hook(&icefront_mparid,1);
110
111 //intialize and add as many inputs per element as requested:
112 this->inputs=new Inputs();
113 this->inputs->AddInput(new IntInput(FillEnum,icefront_fill));
114 this->inputs->AddInput(new IntInput(IceFrontTypeEnum,in_icefront_type));
115
116 //parameters and hooked fields: we still can't point to them, they may not even exist. Configure will handle this.
117 this->parameters = NULL;
118 this->nodes = NULL;
119 this->vertices = NULL;
120 this->element = NULL;
121 this->matpar = NULL;
122}
123
124/*}}}*/
125/*FUNCTION Icefront::Icefront(const char* element_type_in,Inputs* inputs_in,Matpar* matpar_in, int icefront_type, int in_analysis_type) {{{*/
126Icefront::Icefront(const char* element_type_in,Inputs* inputs_in,Matpar* matpar_in,int in_icefront_type, int in_analysis_type){
127
128 int segment_width;
129 int element;
130 int numnodes;
131 int numvertices;
132 int numberofelements;
133
134 /*icefront constructor data: */
135 int icefront_eid;
136 int icefront_mparid;
137 int icefront_node_ids[NUMVERTICESQUA]; //initialize with largest size
138 int icefront_vertex_ids[NUMVERTICESQUA]; //initialize with largest size
139
140// /*find parameters: */
141// iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);
142//
143 /*First, retrieve element index and element type: */
144 if(strcmp(element_type_in,"2d")==0){
145 segment_width=4;
146 }
147 else{
148 segment_width=6;
149 }
150// element=reCast<int,IssmDouble>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+segment_width-2)-1); //element is in the penultimate column (node1 node2 ... elem fill)
151//
152// /*Build ids for hook constructors: */
153// icefront_eid=reCast<int,IssmDouble>( *(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+segment_width-2)); //matlab indexing
154// icefront_mparid=numberofelements+1; //matlab indexing
155//
156 if (in_icefront_type==MacAyeal2dIceFrontEnum || in_icefront_type==MacAyeal3dIceFrontEnum){
157// icefront_node_ids[0]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
158// icefront_node_ids[1]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
159// icefront_vertex_ids[0]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
160// icefront_vertex_ids[1]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
161 }
162 else if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){
163// icefront_node_ids[0]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
164// icefront_node_ids[1]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
165// icefront_node_ids[2]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2));
166// icefront_node_ids[3]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3));
167// icefront_vertex_ids[0]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
168// icefront_vertex_ids[1]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
169// icefront_vertex_ids[2]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2));
170// icefront_vertex_ids[3]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3));
171 }
172 else _error_("in_icefront_type " << EnumToStringx(in_icefront_type) << " not supported yet!");
173
174 if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){
175 numnodes=4;
176 numvertices=4;
177 }
178 else{
179 numnodes=2;
180 numvertices=2;
181 }
182
183 /*Ok, we have everything to build the object: */
184 this->id=1;
185 this->analysis_type=in_analysis_type;
186
187 /*Hooks: */
188 this->hnodes=new Hook(icefront_node_ids,numnodes);
189 this->hvertices=new Hook(icefront_vertex_ids,numvertices);
190 this->helement=new Hook(&icefront_eid,1);
191 this->hmatpar=new Hook(&icefront_mparid,1);
192
193 //intialize and add as many inputs per element as requested:
194 this->inputs=inputs_in;
195 this->inputs->AddInput(new IntInput(FillEnum,1)); //We always consider we have water, if above sea level, only air will be applied
196 this->inputs->AddInput(new IntInput(IceFrontTypeEnum,in_icefront_type));
197
198 //parameters and hooked fields: we still can't point to them, they may not even exist. Configure will handle this.
199 this->parameters = NULL;
200 this->nodes = NULL;
201 this->vertices = NULL;
202 this->element = NULL;
203 this->matpar = matpar_in;
204}
205
206/*}}}*/
207/*FUNCTION Icefront::~Icefront() {{{*/
208Icefront::~Icefront(){
209 delete inputs;
210 this->parameters=NULL;
211 delete hnodes;
212 delete hvertices;
213 delete helement;
214 delete hmatpar;
215}
216/*}}}*/
217
218/*Object virtual functions definitions:*/
219/*FUNCTION Icefront::Echo {{{*/
220void Icefront::Echo(void){
221 _printf_("Icefront:\n");
222 _printf_(" id: " << id << "\n");
223 _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n");
224 hnodes->Echo();
225 hvertices->Echo();
226 helement->Echo();
227 hmatpar->Echo();
228 _printf_(" parameters: " << parameters << "\n");
229 if(parameters)parameters->Echo();
230 _printf_(" inputs: " << inputs << "\n");
231 if(inputs)inputs->Echo();
232}
233/*}}}*/
234/*FUNCTION Icefront::DeepEcho{{{*/
235void Icefront::DeepEcho(void){
236
237 _printf_("Icefront:\n");
238 _printf_(" id: " << id << "\n");
239 _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n");
240 hnodes->DeepEcho();
241 hvertices->DeepEcho();
242 helement->DeepEcho();
243 hmatpar->DeepEcho();
244 _printf_(" parameters: " << parameters << "\n");
245 if(parameters)parameters->DeepEcho();
246 _printf_(" inputs: " << inputs << "\n");
247 if(inputs)inputs->DeepEcho();
248}
249/*}}}*/
250/*FUNCTION Icefront::Id {{{*/
251int Icefront::Id(void){ return id; }
252/*}}}*/
253/*FUNCTION Icefront::ObjectEnum{{{*/
254int Icefront::ObjectEnum(void){
255
256 return IcefrontEnum;
257
258}
259/*}}}*/
260/*FUNCTION Icefront::copy {{{*/
261Object* Icefront::copy() {
262
263 Icefront* icefront=NULL;
264
265 icefront=new Icefront();
266
267 /*copy fields: */
268 icefront->id=this->id;
269 icefront->analysis_type=this->analysis_type;
270 if(this->inputs){
271 icefront->inputs=(Inputs*)this->inputs->Copy();
272 }
273 else{
274 icefront->inputs=new Inputs();
275 }
276 /*point parameters: */
277 icefront->parameters=this->parameters;
278
279 /*now deal with hooks and objects: */
280 icefront->hnodes = (Hook*)this->hnodes->copy();
281 icefront->hvertices = (Hook*)this->hvertices->copy();
282 icefront->helement = (Hook*)this->helement->copy();
283 icefront->hmatpar = (Hook*)this->hmatpar->copy();
284
285 /*corresponding fields*/
286 icefront->nodes = (Node**)icefront->hnodes->deliverp();
287 icefront->vertices = (Vertex**)icefront->hvertices->deliverp();
288 icefront->element = (Element*)icefront->helement->delivers();
289 icefront->matpar = (Matpar*)icefront->hmatpar->delivers();
290
291 return icefront;
292
293}
294/*}}}*/
295
296/*Load virtual functions definitions:*/
297/*FUNCTION Icefront::Configure {{{*/
298void Icefront::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){
299
300 /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective
301 * datasets, using internal ids and offsets hidden in hooks: */
302 hnodes->configure((DataSet*)nodesin);
303 hvertices->configure((DataSet*)verticesin);
304 helement->configure((DataSet*)elementsin);
305 hmatpar->configure((DataSet*)materialsin);
306
307 /*Initialize hooked fields*/
308 this->nodes = (Node**)hnodes->deliverp();
309 this->vertices = (Vertex**)hvertices->deliverp();
310 this->element = (Element*)helement->delivers();
311 this->matpar = (Matpar*)hmatpar->delivers();
312
313 /*point parameters to real dataset: */
314 this->parameters=parametersin;
315}
316/*}}}*/
317/*FUNCTION Icefront::SetCurrentConfiguration {{{*/
318void Icefront::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){
319}
320/*}}}*/
321/*FUNCTION Icefront::CreateKMatrix {{{*/
322void Icefront::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){
323
324 /*No stiffness loads applied, do nothing: */
325 return;
326
327}
328/*}}}*/
329/*FUNCTION Icefront::CreatePVector {{{*/
330void Icefront::CreatePVector(Vector<IssmDouble>* pf){
331
332 /*Checks in debugging mode*/
333 _assert_(nodes);
334 _assert_(element);
335 _assert_(matpar);
336
337 /*Retrieve parameters: */
338 ElementVector* pe=NULL;
339 int analysis_type;
340 this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
341
342 /*Just branch to the correct element icefront vector generator, according to the type of analysis we are carrying out: */
343 switch(analysis_type){
344 #ifdef _HAVE_DIAGNOSTIC_
345 case DiagnosticHorizAnalysisEnum:
346 pe=CreatePVectorDiagnosticHoriz();
347 break;
348 #endif
349 #ifdef _HAVE_CONTROL_
350 case AdjointHorizAnalysisEnum:
351 pe=CreatePVectorAdjointHoriz();
352 break;
353 #endif
354 default:
355 _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
356 }
357
358 /*Add to global Vector*/
359 if(pe){
360 pe->AddToGlobal(pf);
361 delete pe;
362 }
363}
364/*}}}*/
365/*FUNCTION Icefront::CreateJacobianMatrix{{{*/
366void Icefront::CreateJacobianMatrix(Matrix<IssmDouble>* Jff){
367 this->CreateKMatrix(Jff,NULL);
368}
369/*}}}*/
370/*FUNCTION Icefront::GetNodesSidList{{{*/
371void Icefront::GetNodesSidList(int* sidlist){
372
373 int type;
374 inputs->GetInputValue(&type,IceFrontTypeEnum);
375 _assert_(sidlist);
376 _assert_(nodes);
377
378 switch(type){
379 case MacAyeal2dIceFrontEnum:
380 case MacAyeal3dIceFrontEnum:
381 for(int i=0;i<NUMVERTICESSEG;i++) sidlist[i]=nodes[i]->Sid();
382 return;
383#ifdef _HAVE_3D_
384 case PattynIceFrontEnum:
385 case StokesIceFrontEnum:
386 for(int i=0;i<NUMVERTICESQUA;i++) sidlist[i]=nodes[i]->Sid();
387 return;
388#endif
389 default:
390 _error_("Icefront type " << EnumToStringx(type) << " not supported yet");
391 }
392}
393/*}}}*/
394/*FUNCTION Icefront::GetNumberOfNodes{{{*/
395int Icefront::GetNumberOfNodes(void){
396
397 int type;
398 inputs->GetInputValue(&type,IceFrontTypeEnum);
399
400 switch(type){
401 case MacAyeal2dIceFrontEnum:
402 return NUMVERTICESSEG;
403#ifdef _HAVE_3D_
404 case MacAyeal3dIceFrontEnum:
405 return NUMVERTICESSEG;
406 case PattynIceFrontEnum:
407 return NUMVERTICESQUA;
408 case StokesIceFrontEnum:
409 return NUMVERTICESQUA;
410#endif
411 default:
412 _error_("Icefront type " << EnumToStringx(type) << " not supported yet");
413 }
414
415}
416/*}}}*/
417/*FUNCTION Icefront::IsPenalty{{{*/
418bool Icefront::IsPenalty(void){
419 return false;
420}
421/*}}}*/
422/*FUNCTION Icefront::PenaltyCreateKMatrix {{{*/
423void Icefront::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs, IssmDouble kmax){
424 /*do nothing: */
425 return;
426}
427/*}}}*/
428/*FUNCTION Icefront::PenaltyCreatePVector{{{*/
429void Icefront::PenaltyCreatePVector(Vector<IssmDouble>* pf,IssmDouble kmax){
430 /*do nothing: */
431 return;
432}
433/*}}}*/
434/*FUNCTION Icefront::PenaltyCreateJacobianMatrix{{{*/
435void Icefront::PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){
436 this->PenaltyCreateKMatrix(Jff,NULL,kmax);
437}
438/*}}}*/
439/*FUNCTION Icefront::SetwiseNodeConnectivity{{{*/
440void Icefront::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int set1_enum,int set2_enum){
441
442 /*Output */
443 int d_nz = 0;
444 int o_nz = 0;
445
446 /*Loop over all nodes*/
447 for(int i=0;i<this->GetNumberOfNodes();i++){
448
449 if(!flags[this->nodes[i]->Sid()]){
450
451 /*flag current node so that no other element processes it*/
452 flags[this->nodes[i]->Sid()]=true;
453
454 /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
455 switch(set2_enum){
456 case FsetEnum:
457 if(nodes[i]->indexing.fsize){
458 if(this->nodes[i]->IsClone())
459 o_nz += 1;
460 else
461 d_nz += 1;
462 }
463 break;
464 case GsetEnum:
465 if(nodes[i]->indexing.gsize){
466 if(this->nodes[i]->IsClone())
467 o_nz += 1;
468 else
469 d_nz += 1;
470 }
471 break;
472 case SsetEnum:
473 if(nodes[i]->indexing.ssize){
474 if(this->nodes[i]->IsClone())
475 o_nz += 1;
476 else
477 d_nz += 1;
478 }
479 break;
480 default: _error_("not supported");
481 }
482 }
483 }
484
485 /*Assign output pointers: */
486 *pd_nz=d_nz;
487 *po_nz=o_nz;
488}
489/*}}}*/
490/*FUNCTION Icefront::InAnalysis{{{*/
491bool Icefront::InAnalysis(int in_analysis_type){
492 if (in_analysis_type==this->analysis_type)return true;
493 else return false;
494}
495/*}}}*/
496
497/*Update virtual functions definitions:*/
498/*FUNCTION Icefront::InputUpdateFromVector(IssmDouble* vector, int name, int type) {{{*/
499void Icefront::InputUpdateFromVector(IssmDouble* vector, int name, int type){
500 /*Nothing updated yet*/
501}
502/*}}}*/
503/*FUNCTION Icefront::InputUpdateFromVector(int* vector, int name, int type) {{{*/
504void Icefront::InputUpdateFromVector(int* vector, int name, int type){
505 /*Nothing updated yet*/
506}
507/*}}}*/
508/*FUNCTION Icefront::InputUpdateFromVector(bool* vector, int name, int type) {{{*/
509void Icefront::InputUpdateFromVector(bool* vector, int name, int type){
510 /*Nothing updated yet*/
511}
512/*}}}*/
513/*FUNCTION Icefront::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type) {{{*/
514void Icefront::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){
515 /*Nothing updated yet*/
516}
517/*}}}*/
518/*FUNCTION Icefront::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type) {{{*/
519void Icefront::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){
520 /*Nothing updated yet*/
521}
522/*}}}*/
523/*FUNCTION Icefront::InputUpdateFromVectorDakota(int* vector, int name, int type) {{{*/
524void Icefront::InputUpdateFromVectorDakota(int* vector, int name, int type){
525 /*Nothing updated yet*/
526}
527/*}}}*/
528/*FUNCTION Icefront::InputUpdateFromVectorDakota(bool* vector, int name, int type) {{{*/
529void Icefront::InputUpdateFromVectorDakota(bool* vector, int name, int type){
530 /*Nothing updated yet*/
531}
532/*}}}*/
533/*FUNCTION Icefront::InputUpdateFromConstant(IssmDouble constant, int name) {{{*/
534void Icefront::InputUpdateFromConstant(IssmDouble constant, int name){
535 /*Nothing updated yet*/
536}
537/*}}}*/
538/*FUNCTION Icefront::InputUpdateFromConstant(int constant, int name) {{{*/
539void Icefront::InputUpdateFromConstant(int constant, int name){
540 /*Nothing updated yet*/
541}
542/*}}}*/
543/*FUNCTION Icefront::InputUpdateFromConstant(bool constant, int name) {{{*/
544void Icefront::InputUpdateFromConstant(bool constant, int name){
545 /*Nothing updated yet*/
546}
547/*}}}*/
548/*FUNCTION Icefront::InputUpdateFromSolution{{{*/
549void Icefront::InputUpdateFromSolution(IssmDouble* solution){
550 /*Nothing updated yet*/
551}
552/*}}}*/
553
554/*Icefront numerics: */
555#ifdef _HAVE_DIAGNOSTIC_
556/*FUNCTION Icefront::CreatePVectorDiagnosticHoriz {{{*/
557ElementVector* Icefront::CreatePVectorDiagnosticHoriz(void){
558
559 int type;
560 inputs->GetInputValue(&type,IceFrontTypeEnum);
561
562 switch(type){
563 case MacAyeal2dIceFrontEnum:
564 return CreatePVectorDiagnosticMacAyeal2d();
565 #ifdef _HAVE_3D_
566 case MacAyeal3dIceFrontEnum:
567 return CreatePVectorDiagnosticMacAyeal3d();
568 case PattynIceFrontEnum:
569 return CreatePVectorDiagnosticPattyn();
570 case StokesIceFrontEnum:
571 return CreatePVectorDiagnosticStokes();
572 #endif
573 default:
574 _error_("Icefront type " << EnumToStringx(type) << " not supported yet");
575 }
576}
577/*}}}*/
578/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal2d{{{*/
579ElementVector* Icefront::CreatePVectorDiagnosticMacAyeal2d(void){
580
581 /*Intermediary*/
582 int ig,index1,index2,fill;
583 IssmDouble Jdet;
584 IssmDouble thickness,bed,pressure,ice_pressure,rho_water,rho_ice,gravity;
585 IssmDouble water_pressure,air_pressure,surface_under_water,base_under_water;
586 IssmDouble xyz_list[NUMVERTICESSEG][3];
587 IssmDouble normal[2];
588 GaussTria *gauss;
589
590 /*return of element is on water*/
591 Tria* tria=((Tria*)element);
592 if(tria->IsOnWater()) return NULL;
593
594 /*Fetch number of nodes and dof for this finite element*/
595 //int numnodes = this->NumberofNodes();
596 int numnodes = 2;
597 int numdof = numnodes*NDOF2;
598
599 /*Initialize Element vector and vectors*/
600 ElementVector* pe=new ElementVector(nodes,NUMVERTICESSEG,this->parameters,MacAyealApproximationEnum);
601 IssmDouble* basis = xNew<IssmDouble>(numnodes);
602
603 /*Retrieve all inputs and parameters*/
604 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICESSEG);
605 Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
606 Input* bed_input =tria->inputs->GetInput(BedEnum); _assert_(bed_input);
607 inputs->GetInputValue(&fill,FillEnum);
608 rho_water=matpar->GetRhoWater();
609 rho_ice =matpar->GetRhoIce();
610 gravity =matpar->GetG();
611 GetSegmentNormal(&normal[0],xyz_list);
612
613 /*Start looping on Gaussian points*/
614 index1=tria->GetNodeIndex(nodes[0]);
615 index2=tria->GetNodeIndex(nodes[1]);
616 gauss=new GaussTria(index1,index2,3);
617
618 for(ig=gauss->begin();ig<gauss->end();ig++){
619
620 gauss->GaussPoint(ig);
621
622 thickness_input->GetInputValue(&thickness,gauss);
623 bed_input->GetInputValue(&bed,gauss);
624
625 switch(fill){
626 case WaterEnum:
627 surface_under_water=min(0.,thickness+bed); // 0 if the top of the glacier is above water level
628 base_under_water=min(0.,bed); // 0 if the bottom of the glacier is above water level
629 water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
630 break;
631 case AirEnum:
632 water_pressure=0;
633 break;
634 case IceEnum:
635 water_pressure=-1.0/2.0*gravity*rho_ice*pow(thickness,2); // we are facing a wall of ice. use water_pressure to cancel the lithostatic pressure.
636 break;
637 default:
638 _error_("fill type " << EnumToStringx(fill) << " not supported yet");
639 }
640 ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
641 air_pressure=0;
642 pressure = ice_pressure + water_pressure + air_pressure;
643
644 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
645 tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2);
646
647 for (int i=0;i<numnodes;i++){
648 pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
649 pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
650 }
651 }
652
653 /*Transform load vector*/
654 TransformLoadVectorCoord(pe,nodes,NUMVERTICESSEG,XYEnum);
655
656 /*Clean up and return*/
657 xDelete<IssmDouble>(basis);
658 delete gauss;
659 return pe;
660}
661/*}}}*/
662#endif
663
664#ifdef _HAVE_CONTROL_
665/*FUNCTION Icefront::CreatePVectorAdjointHoriz {{{*/
666ElementVector* Icefront::CreatePVectorAdjointHoriz(void){
667
668 /*No load vector applied to the adjoint*/
669 return NULL;
670}
671/*}}}*/
672#endif
673#ifdef _HAVE_3D_
674/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal3d{{{*/
675ElementVector* Icefront::CreatePVectorDiagnosticMacAyeal3d(void){
676
677 Icefront *icefront = NULL;
678 Penta *penta = NULL;
679 Tria *tria = NULL;
680
681 /*Cast element onto Penta*/
682 penta =(Penta*)this->element;
683
684 /*Return if not on bed*/
685 if(!penta->IsOnBed() || penta->IsOnWater()) return NULL;
686
687 /*Spawn Tria and call MacAyeal2d*/
688 tria =(Tria*)penta->SpawnTria(0,1,2);
689 icefront=(Icefront*)this->copy();
690 icefront->element=tria;
691 icefront->inputs->AddInput(new IntInput(IceFrontTypeEnum,MacAyeal2dIceFrontEnum));
692 ElementVector* pe=icefront->CreatePVectorDiagnosticMacAyeal2d();
693
694 /*clean-up and return*/
695 delete tria->material;
696 delete tria;
697 delete icefront;
698 return pe;
699}
700/*}}}*/
701/*FUNCTION Icefront::CreatePVectorDiagnosticPattyn{{{*/
702ElementVector* Icefront::CreatePVectorDiagnosticPattyn(void){
703
704 /*Constants*/
705 const int numdofs = NUMVERTICESQUA *NDOF2;
706
707 /*Intermediaries*/
708 int i,j,ig,index1,index2,index3,index4;
709 int fill;
710 IssmDouble surface,pressure,ice_pressure,rho_water,rho_ice,gravity;
711 IssmDouble water_pressure,air_pressure;
712 IssmDouble Jdet,z_g;
713 IssmDouble xyz_list[NUMVERTICESQUA][3];
714 IssmDouble normal[3];
715 IssmDouble l1l4[4];
716 GaussPenta *gauss = NULL;
717
718 Penta* penta=(Penta*)element;
719
720 /*Initialize Element vector and return if necessary*/
721 if(penta->IsOnWater()) return NULL;
722 ElementVector* pe=new ElementVector(nodes,NUMVERTICESQUA,this->parameters,PattynApproximationEnum);
723
724 /*Retrieve all inputs and parameters*/
725 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICESQUA);
726 Input* surface_input =penta->inputs->GetInput(SurfaceEnum); _assert_(surface_input);
727 inputs->GetInputValue(&fill,FillEnum);
728 rho_water=matpar->GetRhoWater();
729 rho_ice =matpar->GetRhoIce();
730 gravity =matpar->GetG();
731 GetQuadNormal(&normal[0],xyz_list);
732
733 /*Identify which nodes are in the quad: */
734 index1=element->GetNodeIndex(nodes[0]);
735 index2=element->GetNodeIndex(nodes[1]);
736 index3=element->GetNodeIndex(nodes[2]);
737 index4=element->GetNodeIndex(nodes[3]);
738
739 /* Start looping on the number of gaussian points: */
740 IssmDouble zmax=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
741 IssmDouble zmin=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
742 if(zmax>0 && zmin<0) gauss=new GaussPenta(index1,index2,index3,index4,3,10); //refined in vertical because of the sea level discontinuity
743 else gauss=new GaussPenta(index1,index2,index3,index4,3,3);
744 for(ig=gauss->begin();ig<gauss->end();ig++){
745
746 gauss->GaussPoint(ig);
747
748 penta->GetQuadNodalFunctions(l1l4,gauss,index1,index2,index3,index4);
749 penta->GetQuadJacobianDeterminant(&Jdet,xyz_list,gauss);
750 z_g=penta->GetZcoord(gauss);
751 surface_input->GetInputValue(&surface,gauss);
752
753 switch(fill){
754 case WaterEnum:
755 water_pressure=rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level
756 break;
757 case AirEnum:
758 water_pressure=0;
759 break;
760 default:
761 _error_("fill type " << EnumToStringx(fill) << " not supported yet");
762 }
763 ice_pressure=rho_ice*gravity*(surface-z_g);
764 air_pressure=0;
765 pressure = ice_pressure + water_pressure + air_pressure;
766
767 for(i=0;i<NUMVERTICESQUA;i++) for(j=0;j<NDOF2;j++) pe->values[i*NDOF2+j]+=Jdet*gauss->weight*pressure*l1l4[i]*normal[j];
768 }
769
770 /*Transform load vector*/
771 TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,XYEnum);
772
773 /*Clean up and return*/
774 delete gauss;
775 return pe;
776}
777/*}}}*/
778/*FUNCTION Icefront::CreatePVectorDiagnosticStokes{{{*/
779ElementVector* Icefront::CreatePVectorDiagnosticStokes(void){
780
781 /*Constants*/
782 const int numdofs = NUMVERTICESQUA *NDOF4;
783
784 /*Intermediaries*/
785 int i,j,ig,index1,index2,index3,index4;
786 int fill;
787 IssmDouble pressure,rho_water,gravity;
788 IssmDouble water_pressure,air_pressure;
789 IssmDouble Jdet,z_g;
790 IssmDouble xyz_list[NUMVERTICESQUA][3];
791 IssmDouble normal[3];
792 IssmDouble l1l4[4];
793 GaussPenta *gauss = NULL;
794
795 Penta* penta=(Penta*)element;
796
797 /*Initialize Element vector and return if necessary*/
798 if(penta->IsOnWater()) return NULL;
799 ElementVector* pe=new ElementVector(nodes,NUMVERTICESQUA,this->parameters,StokesApproximationEnum);
800
801 /*Retrieve all inputs and parameters*/
802 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICESQUA);
803 inputs->GetInputValue(&fill,FillEnum);
804 rho_water=matpar->GetRhoWater();
805 gravity =matpar->GetG();
806 GetQuadNormal(&normal[0],xyz_list);
807
808 /*Identify which nodes are in the quad: */
809 index1=element->GetNodeIndex(nodes[0]);
810 index2=element->GetNodeIndex(nodes[1]);
811 index3=element->GetNodeIndex(nodes[2]);
812 index4=element->GetNodeIndex(nodes[3]);
813
814 /* Start looping on the number of gaussian points: */
815 IssmDouble zmax=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
816 IssmDouble zmin=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
817 if(zmax>0 && zmin<0) gauss=new GaussPenta(index1,index2,index3,index4,3,30); //refined in vertical because of the sea level discontinuity
818 else gauss=new GaussPenta(index1,index2,index3,index4,3,3);
819 for(ig=gauss->begin();ig<gauss->end();ig++){
820
821 gauss->GaussPoint(ig);
822
823 penta->GetQuadNodalFunctions(l1l4,gauss,index1,index2,index3,index4);
824 penta->GetQuadJacobianDeterminant(&Jdet,xyz_list,gauss);
825 z_g=penta->GetZcoord(gauss);
826
827 switch(fill){
828 case WaterEnum:
829 water_pressure=rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level
830 break;
831 case AirEnum:
832 water_pressure=0;
833 break;
834 default:
835 _error_("fill type " << EnumToStringx(fill) << " not supported yet");
836 }
837 air_pressure=0;
838 pressure = water_pressure + air_pressure; //no ice pressure fore Stokes
839
840 for(i=0;i<NUMVERTICESQUA;i++){
841 for(j=0;j<NDOF4;j++){
842 if(j<3) pe->values[i*NDOF4+j]+=Jdet*gauss->weight*pressure*l1l4[i]*normal[j];
843 else pe->values[i*NDOF4+j]+=0; //pressure term
844 }
845 }
846 }
847
848 /*Transform load vector*/
849 TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,XYZPEnum);
850
851 /*Clean up and return*/
852 delete gauss;
853 return pe;
854}
855/*}}}*/
856#endif
857/*FUNCTION Icefront::GetDofList {{{*/
858void Icefront::GetDofList(int** pdoflist,int approximation_enum,int setenum){
859
860 int numberofdofs=0;
861 int count=0;
862 int type;
863 int numberofnodes=2;
864
865 /*output: */
866 int* doflist=NULL;
867
868 /*recover type: */
869 inputs->GetInputValue(&type,IceFrontTypeEnum);
870
871 /*Some checks for debugging*/
872 _assert_(nodes);
873
874 /*How many nodes? :*/
875 if(type==MacAyeal2dIceFrontEnum || type==MacAyeal3dIceFrontEnum)
876 numberofnodes=2;
877 else
878 numberofnodes=4;
879
880 /*Figure out size of doflist: */
881 for(int i=0;i<numberofnodes;i++){
882 numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
883 }
884
885 /*Allocate: */
886 doflist=xNew<int>(numberofdofs);
887
888 /*Populate: */
889 count=0;
890 for(int i=0;i<numberofnodes;i++){
891 nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
892 count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
893 }
894
895 /*Assign output pointers:*/
896 *pdoflist=doflist;
897}
898/*}}}*/
899/*FUNCTION Icefront::GetSegmentNormal {{{*/
900void Icefront:: GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){
901
902 /*Build unit outward pointing vector*/
903 const int numnodes=NUMVERTICESSEG;
904 IssmDouble vector[2];
905 IssmDouble norm;
906
907 vector[0]=xyz_list[1][0] - xyz_list[0][0];
908 vector[1]=xyz_list[1][1] - xyz_list[0][1];
909
910 norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0));
911
912 normal[0]= + vector[1]/norm;
913 normal[1]= - vector[0]/norm;
914}
915/*}}}*/
916/*FUNCTION Icefront::GetQuadNormal {{{*/
917void Icefront:: GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){
918
919 /*Build unit outward pointing vector*/
920 IssmDouble AB[3];
921 IssmDouble AC[3];
922 IssmDouble norm;
923
924 AB[0]=xyz_list[1][0] - xyz_list[0][0];
925 AB[1]=xyz_list[1][1] - xyz_list[0][1];
926 AB[2]=xyz_list[1][2] - xyz_list[0][2];
927 AC[0]=xyz_list[2][0] - xyz_list[0][0];
928 AC[1]=xyz_list[2][1] - xyz_list[0][1];
929 AC[2]=xyz_list[2][2] - xyz_list[0][2];
930
931 cross(normal,AB,AC);
932 norm=sqrt(pow(normal[0],2.0)+pow(normal[1],2.0)+pow(normal[2],2.0));
933
934 for(int i=0;i<3;i++) normal[i]=normal[i]/norm;
935}
936/*}}}*/
Note: See TracBrowser for help on using the repository browser.