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

Last change on this file since 15099 was 15099, checked in by Eric.Larour, 12 years ago

CHG: greatly simplified the shared/io/Print routines. Replaced
_printf_ by _pprintString_ , then replaced all _printLine_ by _printString_
and _pprintLine_ by _pprintString_
We will then replace the _printString_ by _printf_ and _pprintString_ by _printf0_

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