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

Last change on this file since 15401 was 15401, checked in by seroussi, 12 years ago

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