source: issm/trunk/src/c/classes/Loads/Riftfront.cpp@ 16137

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

merged trunk-jpl and trunk for revision 16135
File size: 28.9 KB
Line 
1/*!\file Riftfront.cpp
2 * \brief: implementation of the Riftfront 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 "shared/shared.h"
14#include "modules/ModelProcessorx/ModelProcessorx.h"
15#include "../classes.h"
16/*}}}*/
17
18/*Element macros*/
19#define NUMVERTICES 2
20
21/*Riftfront constructors and destructor*/
22/*FUNCTION Riftfront::Riftfront(){{{*/
23Riftfront::Riftfront(){
24 this->inputs=NULL;
25 this->parameters=NULL;
26 this->hnodes=NULL;
27 this->helements=NULL;
28 this->hmatpar=NULL;
29 this->nodes=NULL;
30 this->elements=NULL;
31 this->matpar=NULL;
32}
33/*}}}*/
34/*FUNCTION Riftfront::Riftfront(int id, int i, IoModel* iomodel,int analysis_type){{{*/
35Riftfront::Riftfront(int riftfront_id,int i, IoModel* iomodel,int riftfront_analysis_type){
36
37 /*data: */
38 int riftfront_node_ids[2];
39 int riftfront_elem_ids[2];
40 int riftfront_matpar_id;
41 int riftfront_type;
42 int riftfront_fill;
43 IssmDouble riftfront_friction;
44 IssmDouble riftfront_fractionincrement;
45 bool riftfront_shelf;
46 int penalty_lock;
47
48 /*intermediary: */
49 int el1 ,el2;
50 int node1 ,node2;
51
52 /*Fetch parameters: */
53 iomodel->Constant(&penalty_lock,StressbalanceRiftPenaltyLockEnum);
54
55 /*Ok, retrieve all the data needed to add a penalty between the two nodes: */
56 el1=reCast<int,IssmDouble>(*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+2));
57 el2=reCast<int,IssmDouble>(*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+3)) ;
58
59 node1=reCast<int,IssmDouble>(*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+0));
60 node2=reCast<int,IssmDouble>(*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+1));
61
62 /*id: */
63 this->id=riftfront_id;
64 this->analysis_type=riftfront_analysis_type;
65
66 /*hooks: */
67 riftfront_node_ids[0]=iomodel->nodecounter+node1;
68 riftfront_node_ids[1]=iomodel->nodecounter+node2;
69 riftfront_elem_ids[0]=el1;
70 riftfront_elem_ids[1]=el2;
71 riftfront_matpar_id=iomodel->numberofelements+1; //matlab indexing
72
73 /*Hooks: */
74 this->hnodes=new Hook(riftfront_node_ids,2);
75 this->helements=new Hook(riftfront_elem_ids,2);
76 this->hmatpar=new Hook(&riftfront_matpar_id,1);
77
78 /*computational parameters: */
79 this->active=0;
80 this->frozen=0;
81 this->counter=0;
82 this->prestable=0;
83 this->penalty_lock=penalty_lock;
84 this->material_converged=0;
85 this->normal[0]=*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+4);
86 this->normal[1]=*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+5);
87 this->length=*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+6);
88 this->fraction=*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+9);
89 this->state=reCast<int,IssmDouble>(*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+11));
90
91 //intialize inputs, and add as many inputs per element as requested:
92 this->inputs=new Inputs();
93
94 riftfront_type=SegmentRiftfrontEnum;
95 riftfront_fill = reCast<int,IssmDouble>(*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+7));
96 riftfront_friction=*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+8);
97 riftfront_fractionincrement=*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+10);
98 riftfront_shelf=reCast<bool,IssmDouble>(iomodel->Data(MaskGroundediceLevelsetEnum)[node1-1]<0.);
99
100 this->inputs->AddInput(new IntInput(RiftfrontTypeEnum,riftfront_type));
101 this->inputs->AddInput(new IntInput(FillEnum,riftfront_fill));
102 this->inputs->AddInput(new DoubleInput(FrictionEnum,riftfront_friction));
103 this->inputs->AddInput(new DoubleInput(FractionIncrementEnum,riftfront_fractionincrement));
104 this->inputs->AddInput(new BoolInput(SegmentOnIceShelfEnum,riftfront_shelf));
105
106 //parameters and hooked fields: we still can't point to them, they may not even exist. Configure will handle this.
107 this->parameters=NULL;
108 this->nodes= NULL;
109 this->elements= NULL;
110 this->matpar= NULL;
111
112}
113/*}}}*/
114/*FUNCTION Riftfront::~Riftfront(){{{*/
115Riftfront::~Riftfront(){
116 delete inputs;
117 this->parameters=NULL;
118
119 delete hnodes;
120 delete helements;
121 delete hmatpar;
122}
123/*}}}*/
124
125/*Object virtual functions definitions:*/
126/*FUNCTION Riftfront::Echo {{{*/
127void Riftfront::Echo(void){
128
129 Input* input=NULL;
130 int fill;
131 IssmDouble friction,fractionincrement;
132
133 /*recover some inputs first: */
134 input=(Input*)this->inputs->GetInput(FillEnum); input->GetInputValue(&fill);
135 input=(Input*)this->inputs->GetInput(FrictionEnum); input->GetInputValue(&friction);
136 input=(Input*)this->inputs->GetInput(FractionIncrementEnum); input->GetInputValue(&fractionincrement);
137
138 _printf_("Riftfront:\n");
139 _printf_(" id: " << id << "\n");
140 _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n");
141 _printf_(" hnodes: " << hnodes << "\n");
142 _printf_(" helements: " << helements << "\n");
143 _printf_(" hmatpar: " << hmatpar << "\n");
144 _printf_(" parameters: " << parameters << "\n");
145 _printf_(" inputs: " << inputs << "\n");
146 _printf_(" internal parameters: \n");
147 _printf_(" normal: " << normal[0] << "|" << normal[1] << "\n");
148 _printf_(" length: " << length << "\n");
149 _printf_(" penalty_lock: " << penalty_lock << "\n");
150 _printf_(" active: " <<(active ? "true":"false") << "\n");
151 _printf_(" counter: " << counter << "\n");
152 _printf_(" prestable: " << (prestable ? "true":"false") << "\n");
153 _printf_(" material_converged: " << (material_converged ? "true":"false") << "\n");
154 _printf_(" fill: " << fill << "\n");
155 _printf_(" friction: " << friction << "\n");
156 _printf_(" fraction: " << fraction << "\n");
157 _printf_(" fractionincrement: " << fractionincrement << "\n");
158 _printf_(" state: " << state << "\n");
159 _printf_(" frozen: " << (frozen ? "true":"false") << "\n");
160
161}
162/*}}}*/
163/*FUNCTION Riftfront::DeepEcho{{{*/
164void Riftfront::DeepEcho(void){
165
166 _printf_("Riftfront:\n");
167 _printf_(" id: " << id << "\n");
168 _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n");
169 hnodes->DeepEcho();
170 helements->DeepEcho();
171 hmatpar->DeepEcho();
172 _printf_(" parameters\n");
173 if(parameters)parameters->DeepEcho();
174 _printf_(" inputs\n");
175 if(inputs)inputs->DeepEcho();
176}
177/*}}}*/
178/*FUNCTION Riftfront::Id {{{*/
179int Riftfront::Id(void){ return id; }
180/*}}}*/
181/*FUNCTION Riftfront::ObjectEnum{{{*/
182int Riftfront::ObjectEnum(void){
183
184 return RiftfrontEnum;
185
186}
187/*}}}*/
188/*FUNCTION Riftfront::copy {{{*/
189Object* Riftfront::copy() {
190
191 Riftfront* riftfront=NULL;
192
193 riftfront=new Riftfront();
194
195 /*copy fields: */
196 riftfront->id=this->id;
197 riftfront->analysis_type=this->analysis_type;
198 if(this->inputs){
199 riftfront->inputs=(Inputs*)this->inputs->Copy();
200 }
201 else{
202 riftfront->inputs=new Inputs();
203 }
204 /*point parameters: */
205 riftfront->parameters=this->parameters;
206
207 /*now deal with hooks and objects: */
208 riftfront->hnodes=(Hook*)this->hnodes->copy();
209 riftfront->helements=(Hook*)this->helements->copy();
210 riftfront->hmatpar=(Hook*)this->hmatpar->copy();
211
212 /*corresponding fields*/
213 riftfront->nodes =(Node**)riftfront->hnodes->deliverp();
214 riftfront->elements=(Element**)riftfront->helements->deliverp();
215 riftfront->matpar =(Matpar*)riftfront->hmatpar->delivers();
216
217 /*internal data: */
218 riftfront->penalty_lock=this->penalty_lock;
219 riftfront->active=this->active;
220 riftfront->frozen=this->frozen;
221 riftfront->state=this->state;
222 riftfront->counter=this->counter;
223 riftfront->prestable=this->prestable;
224 riftfront->material_converged=this->material_converged;
225 riftfront->normal[0]=this->normal[0];
226 riftfront->normal[1]=this->normal[1];
227 riftfront->length=this->length;
228 riftfront->fraction=this->fraction;
229
230 return riftfront;
231
232}
233/*}}}*/
234
235/*Update virtual functions definitions:*/
236/*FUNCTION Riftfront::InputUpdateFromConstant(bool constant,int name) {{{*/
237void Riftfront::InputUpdateFromConstant(bool constant,int name){
238
239 /*Check that name is a Riftfront input*/
240 if (!IsInput(name)) return;
241
242 /*update input*/
243 this->inputs->AddInput(new BoolInput(name,constant));
244
245}
246/*}}}*/
247/*FUNCTION Riftfront::InputUpdateFromConstant(IssmDouble constant,int name) {{{*/
248void Riftfront::InputUpdateFromConstant(IssmDouble constant,int name){
249
250 /*Check that name is a Riftfront input*/
251 if (!IsInput(name)) return;
252
253 /*update input*/
254 this->inputs->AddInput(new DoubleInput(name,constant));
255
256}
257/*}}}*/
258/*FUNCTION Riftfront::InputUpdateFromVector(IssmDouble* constant,int name) {{{*/
259void Riftfront::InputUpdateFromVector(IssmDouble* vector, int name, int type){
260
261 /*Check that name is a Riftfront input*/
262 if (!IsInput(name)) return;
263
264 /*update input*/
265 _error_("not implemented yet");
266 //this->inputs->AddInput(new DoubleInput(name,constant));
267
268}
269/*}}}*/
270
271/*Load virtual functions definitions:*/
272/*FUNCTION Riftfront::Configure {{{*/
273void Riftfront::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){
274
275 /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective
276 * datasets, using internal ids and offsets hidden in hooks: */
277 hnodes->configure(nodesin);
278 helements->configure(elementsin);
279 hmatpar->configure(materialsin);
280
281 /*Initialize hooked fields*/
282 this->nodes =(Node**)hnodes->deliverp();
283 this->elements=(Element**)helements->deliverp();
284 this->matpar =(Matpar*)hmatpar->delivers();
285
286 /*point parameters to real dataset: */
287 this->parameters=parametersin;
288
289}
290/*}}}*/
291/*FUNCTION Riftfront::IsPenalty{{{*/
292bool Riftfront::IsPenalty(void){
293 return true;
294}
295/*}}}*/
296/*FUNCTION Riftfront::SetCurrentConfiguration {{{*/
297void Riftfront::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){
298
299}
300/*}}}*/
301/*FUNCTION Riftfront::PenaltyCreateKMatrix {{{*/
302void Riftfront::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,IssmDouble kmax){
303
304 /*Retrieve parameters: */
305 ElementMatrix* Ke=NULL;
306 int analysis_type;
307 this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
308
309 switch(analysis_type){
310 case StressbalanceAnalysisEnum:
311 Ke=PenaltyCreateKMatrixStressbalanceHoriz(kmax);
312 break;
313 case AdjointHorizAnalysisEnum:
314 Ke=PenaltyCreateKMatrixStressbalanceHoriz(kmax);
315 break;
316 default:
317 _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
318 }
319
320 /*Add to global Vector*/
321 if(Ke){
322 Ke->AddToGlobal(Kff,Kfs);
323 delete Ke;
324 }
325}
326/*}}}*/
327/*FUNCTION Riftfront::PenaltyCreatePVector {{{*/
328void Riftfront::PenaltyCreatePVector(Vector<IssmDouble>* pf,IssmDouble kmax){
329
330 /*Retrieve parameters: */
331 ElementVector* pe=NULL;
332 int analysis_type;
333 this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
334
335 switch(analysis_type){
336 case StressbalanceAnalysisEnum:
337 pe=PenaltyCreatePVectorStressbalanceHoriz(kmax);
338 break;
339 case AdjointHorizAnalysisEnum:
340 /*No penalty applied on load vector*/
341 break;
342 default:
343 _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
344 }
345
346 /*Add to global Vector*/
347 if(pe){
348 pe->AddToGlobal(pf);
349 delete pe;
350 }
351}
352/*}}}*/
353/*FUNCTION Riftfront::CreateKMatrix {{{*/
354void Riftfront::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){
355 /*do nothing: */
356 return;
357}
358/*}}}*/
359/*FUNCTION Riftfront::CreatePVector {{{*/
360void Riftfront::CreatePVector(Vector<IssmDouble>* pf){
361 /*do nothing: */
362 return;
363}
364/*}}}*/
365/*FUNCTION Riftfront::GetNodesSidList{{{*/
366void Riftfront::GetNodesSidList(int* sidlist){
367
368 _assert_(sidlist);
369 _assert_(nodes);
370
371 for(int i=0;i<NUMVERTICES;i++) sidlist[i]=nodes[i]->Sid();
372}
373/*}}}*/
374/*FUNCTION Riftfront::GetNodesLidList{{{*/
375void Riftfront::GetNodesLidList(int* lidlist){
376
377 _assert_(lidlist);
378 _assert_(nodes);
379
380 for(int i=0;i<NUMVERTICES;i++) lidlist[i]=nodes[i]->Lid();
381}
382/*}}}*/
383/*FUNCTION Riftfront::GetNumberOfNodes{{{*/
384int Riftfront::GetNumberOfNodes(void){
385
386 return NUMVERTICES;
387}
388/*}}}*/
389/*FUNCTION Riftfront::InAnalysis{{{*/
390bool Riftfront::InAnalysis(int in_analysis_type){
391 if (in_analysis_type==this->analysis_type) return true;
392 else return false;
393}
394/*}}}*/
395/*FUNCTION Riftfront::SetwiseNodeConnectivity{{{*/
396void Riftfront::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){
397
398 /*Output */
399 int d_nz = 0;
400 int o_nz = 0;
401
402 /*Loop over all nodes*/
403 for(int i=0;i<NUMVERTICES;i++){
404
405 if(!flags[this->nodes[i]->Lid()]){
406
407 /*flag current node so that no other element processes it*/
408 flags[this->nodes[i]->Lid()]=true;
409
410 int counter=0;
411 while(flagsindices[counter]>=0) counter++;
412 flagsindices[counter]=this->nodes[i]->Lid();
413
414 /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
415 switch(set2_enum){
416 case FsetEnum:
417 if(nodes[i]->indexing.fsize){
418 if(this->nodes[i]->IsClone())
419 o_nz += 1;
420 else
421 d_nz += 1;
422 }
423 break;
424 case GsetEnum:
425 if(nodes[i]->indexing.gsize){
426 if(this->nodes[i]->IsClone())
427 o_nz += 1;
428 else
429 d_nz += 1;
430 }
431 break;
432 case SsetEnum:
433 if(nodes[i]->indexing.ssize){
434 if(this->nodes[i]->IsClone())
435 o_nz += 1;
436 else
437 d_nz += 1;
438 }
439 break;
440 default: _error_("not supported");
441 }
442 }
443 }
444
445 /*Assign output pointers: */
446 *pd_nz=d_nz;
447 *po_nz=o_nz;
448}
449/*}}}*/
450
451/*Riftfront numerics*/
452/*FUNCTION Riftfront::PenaltyCreateKMatrixStressbalanceHoriz {{{*/
453ElementMatrix* Riftfront::PenaltyCreateKMatrixStressbalanceHoriz(IssmDouble kmax){
454
455 const int numdof = NDOF2*NUMVERTICES;
456 int dofs[1] = {0};
457 IssmDouble thickness;
458 IssmDouble h[2];
459 IssmDouble penalty_offset;
460 IssmDouble friction;
461
462 /*Objects: */
463 Tria *tria1 = NULL;
464 Tria *tria2 = NULL;
465
466 /*enum of element? */
467 if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
468 tria1=(Tria*)elements[0];
469 tria2=(Tria*)elements[1];
470
471 /*Initialize Element Matrix*/
472 if(!this->active) return NULL;
473 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters);
474
475 /*Get some parameters: */
476 this->parameters->FindParam(&penalty_offset,StressbalancePenaltyFactorEnum);
477 this->inputs->GetInputValue(&friction,FrictionEnum);
478 tria1->GetInputValue(&h[0],nodes[0],ThicknessEnum);
479 tria2->GetInputValue(&h[1],nodes[1],ThicknessEnum);
480 if (h[0]!=h[1])_error_("different thicknesses not supported for rift fronts");
481 thickness=h[0];
482
483 /*There is contact, we need to constrain the normal velocities (zero penetration), and the
484 *contact slip friction. */
485
486 /*From Peter Wriggers book (Computational Contact Mechanics, p191): */
487 Ke->values[0*numdof+0]+= +pow(normal[0],2)*kmax*pow(10,penalty_offset);
488 Ke->values[0*numdof+1]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset);
489 Ke->values[0*numdof+2]+= -pow(normal[0],2)*kmax*pow(10,penalty_offset);
490 Ke->values[0*numdof+3]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset);
491
492 Ke->values[1*numdof+0]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset);
493 Ke->values[1*numdof+1]+= +pow(normal[1],2)*kmax*pow(10,penalty_offset);
494 Ke->values[1*numdof+2]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset);
495 Ke->values[1*numdof+3]+= -pow(normal[1],2)*kmax*pow(10,penalty_offset);
496
497 Ke->values[2*numdof+0]+= -pow(normal[0],2)*kmax*pow(10,penalty_offset);
498 Ke->values[2*numdof+1]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset);
499 Ke->values[2*numdof+2]+= +pow(normal[0],2)*kmax*pow(10,penalty_offset);
500 Ke->values[2*numdof+3]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset);
501
502 Ke->values[3*numdof+0]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset);
503 Ke->values[3*numdof+1]+= -pow(normal[1],2)*kmax*pow(10,penalty_offset);
504 Ke->values[3*numdof+2]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset);
505 Ke->values[3*numdof+3]+= +pow(normal[1],2)*kmax*pow(10,penalty_offset);
506
507 /*Now take care of the friction: of type sigma=frictiontangent_velocity2-tangent_velocity1)*/
508
509 Ke->values[0*numdof+0]+= +pow(normal[1],2)*thickness*length*friction;
510 Ke->values[0*numdof+1]+= -normal[0]*normal[1]*thickness*length*friction;
511 Ke->values[0*numdof+2]+= -pow(normal[1],2)*thickness*length*friction;
512 Ke->values[0*numdof+3]+= +normal[0]*normal[1]*thickness*length*friction;
513
514 Ke->values[1*numdof+0]+= -normal[0]*normal[1]*thickness*length*friction;
515 Ke->values[1*numdof+1]+= +pow(normal[0],2)*thickness*length*friction;
516 Ke->values[1*numdof+2]+= +normal[0]*normal[1]*thickness*length*friction;
517 Ke->values[1*numdof+3]+= -pow(normal[0],2)*thickness*length*friction;
518
519 Ke->values[2*numdof+0]+= -pow(normal[1],2)*thickness*length*friction;
520 Ke->values[2*numdof+1]+= +normal[0]*normal[1]*thickness*length*friction;
521 Ke->values[2*numdof+2]+= +pow(normal[1],2)*thickness*length*friction;
522 Ke->values[2*numdof+3]+= -normal[0]*normal[1]*thickness*length*friction;
523
524 Ke->values[3*numdof+0]+= +normal[0]*normal[1]*thickness*length*friction;
525 Ke->values[3*numdof+1]+= -pow(normal[0],2)*thickness*length*friction;
526 Ke->values[3*numdof+2]+= -normal[0]*normal[1]*thickness*length*friction;
527 Ke->values[3*numdof+3]+= +pow(normal[0],2)*thickness*length*friction;
528
529 /*Clean up and return*/
530 return Ke;
531}
532/*}}}*/
533/*FUNCTION Riftfront::PenaltyCreatePVectorStressbalanceHoriz {{{*/
534ElementVector* Riftfront::PenaltyCreatePVectorStressbalanceHoriz(IssmDouble kmax){
535
536 const int numdof = NDOF2*NUMVERTICES;
537 int j;
538 IssmDouble rho_ice;
539 IssmDouble rho_water;
540 IssmDouble gravity;
541 IssmDouble thickness;
542 IssmDouble h[2];
543 IssmDouble bed;
544 IssmDouble b[2];
545 IssmDouble pressure;
546 IssmDouble pressure_litho;
547 IssmDouble pressure_air;
548 IssmDouble pressure_melange;
549 IssmDouble pressure_water;
550 int fill;
551 bool shelf;
552
553 /*Objects: */
554 Tria *tria1 = NULL;
555 Tria *tria2 = NULL;
556
557 /*enum of element? */
558 if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
559 tria1=(Tria*)elements[0];
560 tria2=(Tria*)elements[1];
561
562 /*Initialize Element Matrix*/
563 if(this->active) return NULL; /*The penalty is active. No loads implied here.*/
564 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
565
566 /*Get some inputs: */
567 this->inputs->GetInputValue(&fill,FillEnum);
568 this->inputs->GetInputValue(&shelf,SegmentOnIceShelfEnum);
569 rho_ice=matpar->GetRhoIce();
570 rho_water=matpar->GetRhoWater();
571 gravity=matpar->GetG();
572 tria1->GetInputValue(&h[0],nodes[0],ThicknessEnum);
573 tria2->GetInputValue(&h[1],nodes[1],ThicknessEnum);
574 if (h[0]!=h[1])_error_("different thicknesses not supported for rift fronts");
575 thickness=h[0];
576 tria1->GetInputValue(&b[0],nodes[0],BedEnum);
577 tria2->GetInputValue(&b[1],nodes[1],BedEnum);
578 if (b[0]!=b[1])_error_("different beds not supported for rift fronts");
579 bed=b[0];
580
581 /*Ok, this rift is opening. We should put loads on both sides of the rift flanks. Because we are dealing with contact mechanics,
582 * and we want to avoid zigzagging of the loads, we want lump the loads onto nodes, not onto surfaces between nodes.:*/
583
584 /*Ok, to compute the pressure, we are going to need material properties, thickness and bed for the two nodes. We assume those properties to
585 * be the same across the rift.: */
586
587 /*Ok, now compute the pressure (in norm) that is being applied to the flanks, depending on the type of fill: */
588 if(fill==WaterEnum){
589 if(shelf){
590 /*We are on an ice shelf, hydrostatic equilibrium is used to determine the pressure for water fill: */
591 pressure=rho_ice*gravity*pow(thickness,2)/2.- rho_water*gravity*pow(bed,2)/2.;
592 }
593 else{
594 //We are on an icesheet, we assume the water column fills the entire front: */
595 pressure=rho_ice*gravity*pow(thickness,2)/2.- rho_water*gravity*pow(thickness,2)/2.;
596 }
597 }
598 else if(fill==AirEnum){
599 pressure=rho_ice*gravity*pow(thickness,2)/2.; //icefront on an ice sheet, pressure imbalance ice vs air.
600 }
601 else if(fill==IceEnum){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
602 pressure=0;
603 }
604 else if(fill==MelangeEnum){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
605
606 if(!shelf) _error_("fill type " << fill << " not supported on ice sheets yet.");
607
608 pressure_litho=rho_ice*gravity*pow(thickness,2)/2.;
609 pressure_air=0;
610 pressure_melange=rho_ice*gravity*pow(fraction*thickness,2)/2.;
611 pressure_water=1.0/2.0*rho_water*gravity* ( pow(bed,2.0)-pow(rho_ice/rho_water*fraction*thickness,2.0) );
612
613 pressure=pressure_litho-pressure_air-pressure_melange-pressure_water;
614 }
615 else{
616 _error_("fill type " << fill << " not supported yet.");
617 }
618
619 /*Ok, add contribution to first node, along the normal i==0: */
620 for (j=0;j<2;j++){
621 pe->values[j]+=pressure*normal[j]*length;
622 }
623
624 /*Add contribution to second node, along the opposite normal: i==1 */
625 for (j=0;j<2;j++){
626 pe->values[2+j]+= -pressure*normal[j]*length;
627 }
628
629 /*Clean up and return*/
630 return pe;
631}
632/*}}}*/
633/*FUNCTION Riftfront::Constrain {{{*/
634#define _ZIGZAGCOUNTER_
635
636int Riftfront::Constrain(int* punstable){
637
638 const int numnodes = 2;
639 IssmDouble penetration;
640 int activate;
641 int unstable;
642 IssmDouble vx1;
643 IssmDouble vy1;
644 IssmDouble vx2;
645 IssmDouble vy2;
646 IssmDouble fractionincrement;
647
648 /*Objects: */
649 Tria *tria1 = NULL;
650 Tria *tria2 = NULL;
651
652 /*enum of element? */
653 if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
654
655 /*recover elements on both side of rift: */
656 tria1=(Tria*)elements[0];
657 tria2=(Tria*)elements[1];
658
659 /*Is this constraint frozen? In which case we don't touch: */
660 if (this->frozen){
661 *punstable=0;
662 return 1;
663 }
664
665 /*Is this rift segment state specified by user input? :*/
666 if (this->state==OpenEnum || this->state==ClosedEnum){
667
668 if(this->state==OpenEnum)this->active=0;
669 if(this->state==ClosedEnum)this->active=1;
670
671 /*this segment is like frozen, no instability here: */
672 *punstable=0;
673 return 1;
674 }
675
676 /*recover parameters: */
677 this->inputs->GetInputValue(&fractionincrement,FractionIncrementEnum);
678
679 /*First recover velocity: */
680 tria1->GetInputValue(&vx1,nodes[0],VxEnum);
681 tria2->GetInputValue(&vx2,nodes[1],VxEnum);
682 tria1->GetInputValue(&vy1,nodes[0],VyEnum);
683 tria2->GetInputValue(&vy2,nodes[1],VyEnum);
684
685 /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
686 penetration=(vx2-vx1)*normal[0]+(vy2-vy1)*normal[1];
687
688 /*activation: */
689 if(penetration<0)activate=1;
690 else activate=0;
691
692 /*Here, we try to avoid zigzaging. When a penalty activates and deactivates for more than penalty_lock times,
693 * we increase the fraction of melange:*/
694 if(this->counter>this->penalty_lock){
695 /*reset counter: */
696 this->counter=0;
697 /*increase melange fraction: */
698 this->fraction+=fractionincrement;
699 if (this->fraction>1)this->fraction=1.;
700 //_printf_("riftfront " << this->Id() << " fraction: " << this->fraction << "\n");
701 }
702
703 //Figure out stability of this penalty
704 if(this->active==activate){
705 unstable=0;
706 }
707 else{
708 unstable=1;
709 this->counter++;
710 }
711
712 //Set penalty flag
713 this->active=activate;
714
715 //if ((penetration>0) && (this->active==1))_printf_("Riftfront " << Id() << " wants to be released\n");
716
717 /*assign output pointer: */
718 *punstable=unstable;
719 return 1;
720}
721/*}}}*/
722/*FUNCTION Riftfront::FreezeConstraints{{{*/
723void Riftfront::FreezeConstraints(void){
724
725 /*Just set frozen flag to 1: */
726 this->frozen=1;
727
728}
729/*}}}*/
730/*FUNCTION Riftfront::IsFrozen{{{*/
731bool Riftfront::IsFrozen(void){
732
733 /*Just set frozen flag to 1: */
734 if(this->frozen)return 1;
735 else return 0;
736}
737/*}}}*/
738/*FUNCTION Riftfront::IsMaterialStable {{{*/
739int Riftfront::IsMaterialStable(void){
740
741 IssmDouble converged=0;
742
743 this->inputs->GetInputValue(&converged,ConvergedEnum);
744
745 if(reCast<int,IssmDouble>(converged)){
746 /*ok, material non-linearity has converged. If that was already the case, we keep
747 * constraining the rift front. If it was not, and this is the first time the material
748 * has converged, we start constraining now!: */
749 this->material_converged=1;
750 }
751
752 return this->material_converged;
753}
754/*}}}*/
755/*FUNCTION Riftfront::MaxPenetration {{{*/
756int Riftfront::MaxPenetration(IssmDouble* ppenetration){
757
758 const int numnodes=2;
759 IssmDouble penetration=0;
760 IssmDouble vx1;
761 IssmDouble vy1;
762 IssmDouble vx2;
763 IssmDouble vy2;
764
765 /*Objects: */
766 Tria *tria1 = NULL;
767 Tria *tria2 = NULL;
768
769 /*enum of element? */
770 if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
771
772 /*recover elements on both side of rift: */
773 tria1=(Tria*)elements[0];
774 tria2=(Tria*)elements[1];
775
776 //initialize:
777 penetration=-1;
778
779 /*recover velocity: */
780 tria1->GetInputValue(&vx1,nodes[0],VxEnum);
781 tria2->GetInputValue(&vx2,nodes[1],VxEnum);
782 tria1->GetInputValue(&vy1,nodes[0],VyEnum);
783 tria2->GetInputValue(&vy2,nodes[1],VyEnum);
784
785 /*Node1 faces node2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
786 penetration=(vx2-vx1)*normal[0]+(vy2-vy1)*normal[1];
787
788 /*Now, we return penetration only if we are active!: */
789 if(this->active==0)penetration=-1;
790
791 /*If we are zigzag locked, same thing: */
792 if(this->counter>this->penalty_lock)penetration=-1;
793
794 /*assign output pointer: */
795 *ppenetration=penetration;
796 return 1;
797}
798/*}}}*/
799/*FUNCTION Riftfront::Penetration {{{*/
800int Riftfront::Penetration(IssmDouble* ppenetration){
801
802 IssmDouble vx1;
803 IssmDouble vy1;
804 IssmDouble vx2;
805 IssmDouble vy2;
806
807 IssmDouble penetration;
808
809 /*Objects: */
810 Tria *tria1 = NULL;
811 Tria *tria2 = NULL;
812
813 /*enum of element? */
814 if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
815
816 /*recover elements on both side of rift: */
817 tria1=(Tria*)elements[0];
818 tria2=(Tria*)elements[1];
819
820 /*First recover velocity: */
821 tria1->GetInputValue(&vx1,nodes[0],VxEnum);
822 tria2->GetInputValue(&vx2,nodes[1],VxEnum);
823 tria1->GetInputValue(&vy1,nodes[0],VyEnum);
824 tria2->GetInputValue(&vy2,nodes[1],VyEnum);
825
826 /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
827 penetration=(vx2-vx1)*normal[0]+(vy2-vy1)*normal[1];
828
829 /*Now, we return penetration only if we are active!: */
830 if(this->active==0)penetration=0;
831
832 /*assign output pointer: */
833 *ppenetration=penetration;
834 return 1;
835}
836/*}}}*/
837/*FUNCTION Riftfront::PotentialUnstableConstraint {{{*/
838int Riftfront::PotentialUnstableConstraint(int* punstable){
839
840 const int numnodes = 2;
841 IssmDouble penetration;
842 int unstable;
843 IssmDouble vx1;
844 IssmDouble vy1;
845 IssmDouble vx2;
846 IssmDouble vy2;
847
848 /*Objects: */
849 Tria *tria1 = NULL;
850 Tria *tria2 = NULL;
851
852 /*enum of element? */
853 if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
854
855 /*recover elements on both side of rift: */
856 tria1=(Tria*)elements[0];
857 tria2=(Tria*)elements[1];
858
859 /*First recover velocity: */
860 tria1->GetInputValue(&vx1,nodes[0],VxEnum);
861 tria2->GetInputValue(&vx2,nodes[1],VxEnum);
862 tria1->GetInputValue(&vy1,nodes[0],VyEnum);
863 tria2->GetInputValue(&vy2,nodes[1],VyEnum);
864
865 /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
866 penetration=(vx2-vx1)*normal[0]+(vy2-vy1)*normal[1];
867
868 /*Ok, we are looking for positive penetration in an active constraint: */
869 if(this->active){
870 if (penetration>=0){
871 unstable=1;
872 }
873 else{
874 unstable=0;
875 }
876 }
877 else{
878 unstable=0;
879 }
880
881 /*assign output pointer: */
882 *punstable=unstable;
883 return 1;
884}
885/*}}}*/
886/*FUNCTION Riftfront::PreConstrain {{{*/
887int Riftfront::PreConstrain(int* punstable){
888
889 const int numnodes = 2;
890 IssmDouble penetration;
891 int unstable;
892 IssmDouble vx1;
893 IssmDouble vy1;
894 IssmDouble vx2;
895 IssmDouble vy2;
896
897 /*Objects: */
898 Tria *tria1 = NULL;
899 Tria *tria2 = NULL;
900
901 /*enum of element? */
902 if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
903
904 /*recover elements on both side of rift: */
905 tria1=(Tria*)elements[0];
906 tria2=(Tria*)elements[1];
907
908 /*First recover velocity: */
909 tria1->GetInputValue(&vx1,nodes[0],VxEnum);
910 tria2->GetInputValue(&vx2,nodes[1],VxEnum);
911 tria1->GetInputValue(&vy1,nodes[0],VyEnum);
912 tria2->GetInputValue(&vy2,nodes[1],VyEnum);
913
914 /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
915 penetration=(vx2-vx1)*normal[0]+(vy2-vy1)*normal[1];
916
917 /*Ok, we are preconstraining here. Ie, anything that penetrates is constrained until stability of the entire set
918 * of constraints is reached.: */
919 if(penetration<0){
920 if (!this->active){
921 /*This is the first time penetration happens: */
922 this->active=1;
923 unstable=1;
924 }
925 else{
926 /*This constraint was already active: */
927 this->active=1;
928 unstable=0;
929 }
930 }
931 else{
932 /*No penetration happening. : */
933 if (!this->active){
934 /*This penalty was not active, and no penetration happening. Do nonthing: */
935 this->active=0;
936 unstable=0;
937 }
938 else{
939 /*Ok, this penalty wants to get released. But not now, this is preconstraint, not constraint: */
940 this->active=1;
941 unstable=0;
942 }
943 }
944
945 /*assign output pointer: */
946 *punstable=unstable;
947 return 1;
948}
949/*}}}*/
950/*FUNCTION Riftfront::PreStable {{{*/
951bool Riftfront::PreStable(){
952 return prestable;
953}
954/*}}}*/
955/*FUNCTION Riftfront::SetPreStable {{{*/
956void Riftfront::SetPreStable(){
957 prestable=1;
958}
959/*}}}*/
960/*FUNCTION Riftfront::IsInput{{{*/
961bool Riftfront::IsInput(int name){
962 if (
963 name==ConvergedEnum ||
964 name==ThicknessEnum ||
965 name==SurfaceEnum ||
966 name==BedEnum
967 ){
968 return true;
969 }
970 else return false;
971}
972/*}}}*/
Note: See TracBrowser for help on using the repository browser.