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

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

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