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

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

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