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

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

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