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

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

CHG: iomodel->numberofelements now always available
File size: 29.5 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==MacAyeal2dIceFrontEnum || in_icefront_type==MacAyeal3dIceFrontEnum){
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==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){
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(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=reCast<int,IssmDouble>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+segment_width-2)-1); //element is in the penultimate column (node1 node2 ... elem fill)
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==MacAyeal2dIceFrontEnum || in_icefront_type==MacAyeal3dIceFrontEnum){
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==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){
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==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){
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 MacAyeal2dIceFrontEnum:
378 case MacAyeal3dIceFrontEnum:
379 for(int i=0;i<NUMVERTICESSEG;i++) sidlist[i]=nodes[i]->Sid();
380 return;
381#ifdef _HAVE_3D_
382 case PattynIceFrontEnum:
383 case StokesIceFrontEnum:
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 MacAyeal2dIceFrontEnum:
400 return NUMVERTICESSEG;
401#ifdef _HAVE_3D_
402 case MacAyeal3dIceFrontEnum:
403 return NUMVERTICESSEG;
404 case PattynIceFrontEnum:
405 return NUMVERTICESQUA;
406 case StokesIceFrontEnum:
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 MacAyeal2dIceFrontEnum:
562 return CreatePVectorDiagnosticMacAyeal2d();
563 #ifdef _HAVE_3D_
564 case MacAyeal3dIceFrontEnum:
565 return CreatePVectorDiagnosticMacAyeal3d();
566 case PattynIceFrontEnum:
567 return CreatePVectorDiagnosticPattyn();
568 case StokesIceFrontEnum:
569 return CreatePVectorDiagnosticStokes();
570 #endif
571 default:
572 _error_("Icefront type " << EnumToStringx(type) << " not supported yet");
573 }
574}
575/*}}}*/
576/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal2d{{{*/
577ElementVector* Icefront::CreatePVectorDiagnosticMacAyeal2d(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 = this->NumberofNodes();
594 int numnodes = 2;
595 int numdof = numnodes*NDOF2;
596
597 /*Initialize Element vector and vectors*/
598 ElementVector* pe=new ElementVector(nodes,NUMVERTICESSEG,this->parameters,MacAyealApproximationEnum);
599 IssmDouble* basis = xNew<IssmDouble>(numnodes);
600
601 /*Retrieve all inputs and parameters*/
602 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICESSEG);
603 Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
604 Input* bed_input =tria->inputs->GetInput(BedEnum); _assert_(bed_input);
605 inputs->GetInputValue(&fill,FillEnum);
606 rho_water=matpar->GetRhoWater();
607 rho_ice =matpar->GetRhoIce();
608 gravity =matpar->GetG();
609 GetSegmentNormal(&normal[0],xyz_list);
610
611 /*Start looping on Gaussian points*/
612 index1=tria->GetNodeIndex(nodes[0]);
613 index2=tria->GetNodeIndex(nodes[1]);
614 gauss=new GaussTria(index1,index2,3);
615
616 for(ig=gauss->begin();ig<gauss->end();ig++){
617
618 gauss->GaussPoint(ig);
619
620 thickness_input->GetInputValue(&thickness,gauss);
621 bed_input->GetInputValue(&bed,gauss);
622
623 switch(fill){
624 case WaterEnum:
625 surface_under_water=min(0.,thickness+bed); // 0 if the top of the glacier is above water level
626 base_under_water=min(0.,bed); // 0 if the bottom of the glacier is above water level
627 water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
628 break;
629 case AirEnum:
630 water_pressure=0;
631 break;
632 case IceEnum:
633 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.
634 break;
635 default:
636 _error_("fill type " << EnumToStringx(fill) << " not supported yet");
637 }
638 ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
639 air_pressure=0;
640 pressure = ice_pressure + water_pressure + air_pressure;
641
642 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
643 tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2);
644
645 for (int i=0;i<numnodes;i++){
646 pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
647 pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
648 }
649 }
650
651 /*Transform load vector*/
652 TransformLoadVectorCoord(pe,nodes,NUMVERTICESSEG,XYEnum);
653
654 /*Clean up and return*/
655 xDelete<IssmDouble>(basis);
656 delete gauss;
657 return pe;
658}
659/*}}}*/
660#endif
661
662#ifdef _HAVE_CONTROL_
663/*FUNCTION Icefront::CreatePVectorAdjointHoriz {{{*/
664ElementVector* Icefront::CreatePVectorAdjointHoriz(void){
665
666 /*No load vector applied to the adjoint*/
667 return NULL;
668}
669/*}}}*/
670#endif
671#ifdef _HAVE_3D_
672/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal3d{{{*/
673ElementVector* Icefront::CreatePVectorDiagnosticMacAyeal3d(void){
674
675 Icefront *icefront = NULL;
676 Penta *penta = NULL;
677 Tria *tria = NULL;
678
679 /*Cast element onto Penta*/
680 penta =(Penta*)this->element;
681
682 /*Return if not on bed*/
683 if(!penta->IsOnBed() || penta->IsOnWater()) return NULL;
684
685 /*Spawn Tria and call MacAyeal2d*/
686 tria =(Tria*)penta->SpawnTria(0,1,2);
687 icefront=(Icefront*)this->copy();
688 icefront->element=tria;
689 icefront->inputs->AddInput(new IntInput(IceFrontTypeEnum,MacAyeal2dIceFrontEnum));
690 ElementVector* pe=icefront->CreatePVectorDiagnosticMacAyeal2d();
691
692 /*clean-up and return*/
693 delete tria->material;
694 delete tria;
695 delete icefront;
696 return pe;
697}
698/*}}}*/
699/*FUNCTION Icefront::CreatePVectorDiagnosticPattyn{{{*/
700ElementVector* Icefront::CreatePVectorDiagnosticPattyn(void){
701
702 /*Constants*/
703 const int numdofs = NUMVERTICESQUA *NDOF2;
704
705 /*Intermediaries*/
706 int i,j,ig,index1,index2,index3,index4;
707 int fill;
708 IssmDouble surface,pressure,ice_pressure,rho_water,rho_ice,gravity;
709 IssmDouble water_pressure,air_pressure;
710 IssmDouble Jdet,z_g;
711 IssmDouble xyz_list[NUMVERTICESQUA][3];
712 IssmDouble normal[3];
713 IssmDouble l1l4[4];
714 GaussPenta *gauss = NULL;
715
716 Penta* penta=(Penta*)element;
717
718 /*Initialize Element vector and return if necessary*/
719 if(penta->IsOnWater()) return NULL;
720 ElementVector* pe=new ElementVector(nodes,NUMVERTICESQUA,this->parameters,PattynApproximationEnum);
721
722 /*Retrieve all inputs and parameters*/
723 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICESQUA);
724 Input* surface_input =penta->inputs->GetInput(SurfaceEnum); _assert_(surface_input);
725 inputs->GetInputValue(&fill,FillEnum);
726 rho_water=matpar->GetRhoWater();
727 rho_ice =matpar->GetRhoIce();
728 gravity =matpar->GetG();
729 GetQuadNormal(&normal[0],xyz_list);
730
731 /*Identify which nodes are in the quad: */
732 index1=element->GetNodeIndex(nodes[0]);
733 index2=element->GetNodeIndex(nodes[1]);
734 index3=element->GetNodeIndex(nodes[2]);
735 index4=element->GetNodeIndex(nodes[3]);
736
737 /* Start looping on the number of gaussian points: */
738 IssmDouble zmax=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
739 IssmDouble zmin=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
740 if(zmax>0 && zmin<0) gauss=new GaussPenta(index1,index2,index3,index4,3,10); //refined in vertical because of the sea level discontinuity
741 else gauss=new GaussPenta(index1,index2,index3,index4,3,3);
742 for(ig=gauss->begin();ig<gauss->end();ig++){
743
744 gauss->GaussPoint(ig);
745
746 penta->GetQuadNodalFunctions(l1l4,gauss,index1,index2,index3,index4);
747 penta->GetQuadJacobianDeterminant(&Jdet,xyz_list,gauss);
748 z_g=penta->GetZcoord(gauss);
749 surface_input->GetInputValue(&surface,gauss);
750
751 switch(fill){
752 case WaterEnum:
753 water_pressure=rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level
754 break;
755 case AirEnum:
756 water_pressure=0;
757 break;
758 default:
759 _error_("fill type " << EnumToStringx(fill) << " not supported yet");
760 }
761 ice_pressure=rho_ice*gravity*(surface-z_g);
762 air_pressure=0;
763 pressure = ice_pressure + water_pressure + air_pressure;
764
765 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];
766 }
767
768 /*Transform load vector*/
769 TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,XYEnum);
770
771 /*Clean up and return*/
772 delete gauss;
773 return pe;
774}
775/*}}}*/
776/*FUNCTION Icefront::CreatePVectorDiagnosticStokes{{{*/
777ElementVector* Icefront::CreatePVectorDiagnosticStokes(void){
778
779 /*Constants*/
780 const int numdofs = NUMVERTICESQUA *NDOF4;
781
782 /*Intermediaries*/
783 int i,j,ig,index1,index2,index3,index4;
784 int fill;
785 IssmDouble pressure,rho_water,gravity;
786 IssmDouble water_pressure,air_pressure;
787 IssmDouble Jdet,z_g;
788 IssmDouble xyz_list[NUMVERTICESQUA][3];
789 IssmDouble normal[3];
790 IssmDouble l1l4[4];
791 GaussPenta *gauss = NULL;
792
793 Penta* penta=(Penta*)element;
794
795 /*Initialize Element vector and return if necessary*/
796 if(penta->IsOnWater()) return NULL;
797 ElementVector* pe=new ElementVector(nodes,NUMVERTICESQUA,this->parameters,StokesApproximationEnum);
798
799 /*Retrieve all inputs and parameters*/
800 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICESQUA);
801 inputs->GetInputValue(&fill,FillEnum);
802 rho_water=matpar->GetRhoWater();
803 gravity =matpar->GetG();
804 GetQuadNormal(&normal[0],xyz_list);
805
806 /*Identify which nodes are in the quad: */
807 index1=element->GetNodeIndex(nodes[0]);
808 index2=element->GetNodeIndex(nodes[1]);
809 index3=element->GetNodeIndex(nodes[2]);
810 index4=element->GetNodeIndex(nodes[3]);
811
812 /* Start looping on the number of gaussian points: */
813 IssmDouble zmax=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
814 IssmDouble zmin=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
815 if(zmax>0 && zmin<0) gauss=new GaussPenta(index1,index2,index3,index4,3,30); //refined in vertical because of the sea level discontinuity
816 else gauss=new GaussPenta(index1,index2,index3,index4,3,3);
817 for(ig=gauss->begin();ig<gauss->end();ig++){
818
819 gauss->GaussPoint(ig);
820
821 penta->GetQuadNodalFunctions(l1l4,gauss,index1,index2,index3,index4);
822 penta->GetQuadJacobianDeterminant(&Jdet,xyz_list,gauss);
823 z_g=penta->GetZcoord(gauss);
824
825 switch(fill){
826 case WaterEnum:
827 water_pressure=rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level
828 break;
829 case AirEnum:
830 water_pressure=0;
831 break;
832 default:
833 _error_("fill type " << EnumToStringx(fill) << " not supported yet");
834 }
835 air_pressure=0;
836 pressure = water_pressure + air_pressure; //no ice pressure fore Stokes
837
838 for(i=0;i<NUMVERTICESQUA;i++){
839 for(j=0;j<NDOF4;j++){
840 if(j<3) pe->values[i*NDOF4+j]+=Jdet*gauss->weight*pressure*l1l4[i]*normal[j];
841 else pe->values[i*NDOF4+j]+=0; //pressure term
842 }
843 }
844 }
845
846 /*Transform load vector*/
847 TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,XYZPEnum);
848
849 /*Clean up and return*/
850 delete gauss;
851 return pe;
852}
853/*}}}*/
854#endif
855/*FUNCTION Icefront::GetDofList {{{*/
856void Icefront::GetDofList(int** pdoflist,int approximation_enum,int setenum){
857
858 int numberofdofs=0;
859 int count=0;
860 int type;
861 int numberofnodes=2;
862
863 /*output: */
864 int* doflist=NULL;
865
866 /*recover type: */
867 inputs->GetInputValue(&type,IceFrontTypeEnum);
868
869 /*Some checks for debugging*/
870 _assert_(nodes);
871
872 /*How many nodes? :*/
873 if(type==MacAyeal2dIceFrontEnum || type==MacAyeal3dIceFrontEnum)
874 numberofnodes=2;
875 else
876 numberofnodes=4;
877
878 /*Figure out size of doflist: */
879 for(int i=0;i<numberofnodes;i++){
880 numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
881 }
882
883 /*Allocate: */
884 doflist=xNew<int>(numberofdofs);
885
886 /*Populate: */
887 count=0;
888 for(int i=0;i<numberofnodes;i++){
889 nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
890 count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
891 }
892
893 /*Assign output pointers:*/
894 *pdoflist=doflist;
895}
896/*}}}*/
897/*FUNCTION Icefront::GetSegmentNormal {{{*/
898void Icefront:: GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){
899
900 /*Build unit outward pointing vector*/
901 const int numnodes=NUMVERTICESSEG;
902 IssmDouble vector[2];
903 IssmDouble norm;
904
905 vector[0]=xyz_list[1][0] - xyz_list[0][0];
906 vector[1]=xyz_list[1][1] - xyz_list[0][1];
907
908 norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0));
909
910 normal[0]= + vector[1]/norm;
911 normal[1]= - vector[0]/norm;
912}
913/*}}}*/
914/*FUNCTION Icefront::GetQuadNormal {{{*/
915void Icefront:: GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){
916
917 /*Build unit outward pointing vector*/
918 IssmDouble AB[3];
919 IssmDouble AC[3];
920 IssmDouble norm;
921
922 AB[0]=xyz_list[1][0] - xyz_list[0][0];
923 AB[1]=xyz_list[1][1] - xyz_list[0][1];
924 AB[2]=xyz_list[1][2] - xyz_list[0][2];
925 AC[0]=xyz_list[2][0] - xyz_list[0][0];
926 AC[1]=xyz_list[2][1] - xyz_list[0][1];
927 AC[2]=xyz_list[2][2] - xyz_list[0][2];
928
929 cross(normal,AB,AC);
930 norm=sqrt(pow(normal[0],2.0)+pow(normal[1],2.0)+pow(normal[2],2.0));
931
932 for(int i=0;i<3;i++) normal[i]=normal[i]/norm;
933}
934/*}}}*/
Note: See TracBrowser for help on using the repository browser.