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

Last change on this file since 20500 was 20500, checked in by Mathieu Morlighem, 9 years ago

merged trunk-jpl and trunk for revision 20497
File size: 21.4 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*/
22Riftfront::Riftfront(){/*{{{*/
23 this->parameters=NULL;
24 this->hnodes=NULL;
25 this->helements=NULL;
26 this->hmatpar=NULL;
27 this->nodes=NULL;
28 this->elements=NULL;
29 this->matpar=NULL;
30}
31/*}}}*/
32Riftfront::Riftfront(int riftfront_id,int i, IoModel* iomodel,int riftfront_analysis_type){/*{{{*/
33
34 /*data: */
35 const int RIFTINFOSIZE = 12;
36 int riftfront_node_ids[2];
37 int riftfront_elem_ids[2];
38 int riftfront_matpar_id;
39 IssmDouble riftfront_friction;
40 IssmDouble riftfront_fractionincrement;
41 int penalty_lock;
42
43 /*intermediary: */
44 int el1 ,el2;
45 int node1 ,node2;
46
47 /*Fetch parameters: */
48 iomodel->Constant(&penalty_lock,StressbalanceRiftPenaltyLockEnum);
49
50 /*Ok, retrieve all the data needed to add a penalty between the two nodes: */
51 el1=reCast<int,IssmDouble>(*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+2));
52 el2=reCast<int,IssmDouble>(*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+3)) ;
53
54 node1=reCast<int,IssmDouble>(*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+0));
55 node2=reCast<int,IssmDouble>(*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+1));
56
57 /*id: */
58 this->id=riftfront_id;
59 this->analysis_type=riftfront_analysis_type;
60
61 /*hooks: */
62 riftfront_node_ids[0]=iomodel->nodecounter+node1;
63 riftfront_node_ids[1]=iomodel->nodecounter+node2;
64 riftfront_elem_ids[0]=el1;
65 riftfront_elem_ids[1]=el2;
66 riftfront_matpar_id=iomodel->numberofelements+1; //matlab indexing
67
68 /*Hooks: */
69 this->hnodes=new Hook(riftfront_node_ids,2);
70 this->helements=new Hook(riftfront_elem_ids,2);
71 this->hmatpar=new Hook(&riftfront_matpar_id,1);
72
73 /*computational parameters: */
74 this->active=0;
75 this->frozen=0;
76 this->counter=0;
77 this->prestable=0;
78 this->penalty_lock=penalty_lock;
79 this->material_converged=0;
80 this->normal[0]=*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+4);
81 this->normal[1]=*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+5);
82 this->length=*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+6);
83 this->fraction=*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+9);
84 this->state=reCast<int,IssmDouble>(*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+11));
85
86 //intialize properties
87 this->type=SegmentRiftfrontEnum;
88 this->fill = reCast<int,IssmDouble>(*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+7));
89 this->friction=*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+8);
90 this->fractionincrement=*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+10);
91 this->shelf=reCast<bool,IssmDouble>(iomodel->Data(MaskGroundediceLevelsetEnum)[node1-1]<0.);
92
93 //parameters and hooked fields: we still can't point to them, they may not even exist. Configure will handle this.
94 this->parameters=NULL;
95 this->nodes= NULL;
96 this->elements= NULL;
97 this->matpar= NULL;
98
99}
100/*}}}*/
101Riftfront::~Riftfront(){/*{{{*/
102 this->parameters=NULL;
103 delete hnodes;
104 delete helements;
105 delete hmatpar;
106}
107/*}}}*/
108
109/*Object virtual functions definitions:*/
110Object* Riftfront::copy() {/*{{{*/
111
112 Riftfront* riftfront=NULL;
113
114 riftfront=new Riftfront();
115
116 /*copy fields: */
117 riftfront->id=this->id;
118 riftfront->analysis_type=this->analysis_type;
119 riftfront->type=this->type;
120 riftfront->fill=this->fill;
121 riftfront->friction=this->friction;
122 riftfront->fractionincrement=this->fractionincrement;
123 riftfront->shelf=this->shelf;
124
125 /*point parameters: */
126 riftfront->parameters=this->parameters;
127
128 /*now deal with hooks and objects: */
129 riftfront->hnodes=(Hook*)this->hnodes->copy();
130 riftfront->helements=(Hook*)this->helements->copy();
131 riftfront->hmatpar=(Hook*)this->hmatpar->copy();
132
133 /*corresponding fields*/
134 riftfront->nodes =(Node**)riftfront->hnodes->deliverp();
135 riftfront->elements=(Element**)riftfront->helements->deliverp();
136 riftfront->matpar =(Matpar*)riftfront->hmatpar->delivers();
137
138 /*internal data: */
139 riftfront->penalty_lock=this->penalty_lock;
140 riftfront->active=this->active;
141 riftfront->frozen=this->frozen;
142 riftfront->state=this->state;
143 riftfront->counter=this->counter;
144 riftfront->prestable=this->prestable;
145 riftfront->material_converged=this->material_converged;
146 riftfront->normal[0]=this->normal[0];
147 riftfront->normal[1]=this->normal[1];
148 riftfront->length=this->length;
149 riftfront->fraction=this->fraction;
150
151 return riftfront;
152
153}
154/*}}}*/
155void Riftfront::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
156
157 _assert_(this);
158
159 /*ok, marshall operations: */
160 MARSHALLING_ENUM(RiftfrontEnum);
161 MARSHALLING(id);
162 MARSHALLING(analysis_type);
163 MARSHALLING(type);
164 MARSHALLING(fill);
165 MARSHALLING(friction);
166 MARSHALLING(fractionincrement);
167 MARSHALLING(shelf);
168
169 if(marshall_direction==MARSHALLING_BACKWARD){
170 this->hnodes = new Hook();
171 this->hmatpar = new Hook();
172 this->helements = new Hook();
173 }
174
175 this->hnodes->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
176 this->hmatpar->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
177 this->helements->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
178
179 /*corresponding fields*/
180 nodes =(Node**)this->hnodes->deliverp();
181 matpar =(Matpar*)this->hmatpar->delivers();
182 elements =(Element**)this->helements->deliverp();
183
184 MARSHALLING(penalty_lock);
185 MARSHALLING(active);
186 MARSHALLING(frozen);
187 MARSHALLING(state);
188 MARSHALLING(counter);
189 MARSHALLING(prestable);
190 MARSHALLING(material_converged);
191 MARSHALLING(normal[0]);
192 MARSHALLING(normal[1]);
193 MARSHALLING(length);
194 MARSHALLING(fraction);
195
196}
197/*}}}*/
198void Riftfront::DeepEcho(void){/*{{{*/
199
200 _printf_("Riftfront:\n");
201 _printf_(" id: " << id << "\n");
202 _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n");
203 hnodes->DeepEcho();
204 helements->DeepEcho();
205 hmatpar->DeepEcho();
206 _printf_(" parameters\n");
207 if(parameters)parameters->DeepEcho();
208}
209/*}}}*/
210void Riftfront::Echo(void){/*{{{*/
211
212 _printf_("Riftfront:\n");
213 _printf_(" id: " << id << "\n");
214 _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n");
215 _printf_(" hnodes: " << hnodes << "\n");
216 _printf_(" helements: " << helements << "\n");
217 _printf_(" hmatpar: " << hmatpar << "\n");
218 _printf_(" parameters: " << parameters << "\n");
219 _printf_(" internal parameters: \n");
220 _printf_(" normal: " << normal[0] << "|" << normal[1] << "\n");
221 _printf_(" length: " << length << "\n");
222 _printf_(" penalty_lock: " << penalty_lock << "\n");
223 _printf_(" active: " <<(active ? "true":"false") << "\n");
224 _printf_(" counter: " << counter << "\n");
225 _printf_(" prestable: " << (prestable ? "true":"false") << "\n");
226 _printf_(" material_converged: " << (material_converged ? "true":"false") << "\n");
227 _printf_(" fill: " << fill << "\n");
228 _printf_(" friction: " << friction << "\n");
229 _printf_(" fraction: " << fraction << "\n");
230 _printf_(" fractionincrement: " << fractionincrement << "\n");
231 _printf_(" state: " << state << "\n");
232 _printf_(" frozen: " << (frozen ? "true":"false") << "\n");
233
234}
235/*}}}*/
236int Riftfront::Id(void){ return id; }/*{{{*/
237/*}}}*/
238int Riftfront::ObjectEnum(void){/*{{{*/
239
240 return RiftfrontEnum;
241
242}
243/*}}}*/
244
245/*Update virtual functions definitions:*/
246void Riftfront::InputUpdateFromConstant(bool constant,int name){/*{{{*/
247}
248/*}}}*/
249void Riftfront::InputUpdateFromConstant(IssmDouble constant,int name){/*{{{*/
250
251}
252/*}}}*/
253void Riftfront::InputUpdateFromVector(IssmDouble* vector, int name, int type){/*{{{*/
254
255 /*Nothing to update*/
256
257}
258/*}}}*/
259
260/*Load virtual functions definitions:*/
261void Riftfront::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
262
263 /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective
264 * datasets, using internal ids and offsets hidden in hooks: */
265 hnodes->configure(nodesin);
266 helements->configure(elementsin);
267 hmatpar->configure(materialsin);
268
269 /*Initialize hooked fields*/
270 this->nodes =(Node**)hnodes->deliverp();
271 this->elements=(Element**)helements->deliverp();
272 this->matpar =(Matpar*)hmatpar->delivers();
273
274 /*point parameters to real dataset: */
275 this->parameters=parametersin;
276
277}
278/*}}}*/
279void Riftfront::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){/*{{{*/
280 /*do nothing: */
281 return;
282}
283/*}}}*/
284void Riftfront::CreatePVector(Vector<IssmDouble>* pf){/*{{{*/
285 /*do nothing: */
286 return;
287}
288/*}}}*/
289void Riftfront::GetNodesLidList(int* lidlist){/*{{{*/
290
291 _assert_(lidlist);
292 _assert_(nodes);
293
294 for(int i=0;i<NUMVERTICES;i++) lidlist[i]=nodes[i]->Lid();
295}
296/*}}}*/
297void Riftfront::GetNodesSidList(int* sidlist){/*{{{*/
298
299 _assert_(sidlist);
300 _assert_(nodes);
301
302 for(int i=0;i<NUMVERTICES;i++) sidlist[i]=nodes[i]->Sid();
303}
304/*}}}*/
305int Riftfront::GetNumberOfNodes(void){/*{{{*/
306
307 return NUMVERTICES;
308}
309/*}}}*/
310bool Riftfront::InAnalysis(int in_analysis_type){/*{{{*/
311 if (in_analysis_type==this->analysis_type) return true;
312 else return false;
313}
314/*}}}*/
315bool Riftfront::IsPenalty(void){/*{{{*/
316 return true;
317}
318/*}}}*/
319void Riftfront::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,IssmDouble kmax){/*{{{*/
320
321 /*Retrieve parameters: */
322 ElementMatrix* Ke=NULL;
323 int analysis_type;
324 this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
325
326 switch(analysis_type){
327 case StressbalanceAnalysisEnum:
328 Ke=PenaltyCreateKMatrixStressbalanceHoriz(kmax);
329 break;
330 case AdjointHorizAnalysisEnum:
331 Ke=PenaltyCreateKMatrixStressbalanceHoriz(kmax);
332 break;
333 default:
334 _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
335 }
336
337 /*Add to global Vector*/
338 if(Ke){
339 Ke->AddToGlobal(Kff,Kfs);
340 delete Ke;
341 }
342}
343/*}}}*/
344void Riftfront::PenaltyCreatePVector(Vector<IssmDouble>* pf,IssmDouble kmax){/*{{{*/
345
346 /*Retrieve parameters: */
347 ElementVector* pe=NULL;
348 int analysis_type;
349 this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
350
351 switch(analysis_type){
352 case StressbalanceAnalysisEnum:
353 pe=PenaltyCreatePVectorStressbalanceHoriz(kmax);
354 break;
355 case AdjointHorizAnalysisEnum:
356 /*No penalty applied on load vector*/
357 break;
358 default:
359 _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
360 }
361
362 /*Add to global Vector*/
363 if(pe){
364 pe->AddToGlobal(pf);
365 delete pe;
366 }
367}
368/*}}}*/
369void Riftfront::ResetHooks(){/*{{{*/
370
371 this->nodes=NULL;
372 this->elements=NULL;
373 this->matpar=NULL;
374 this->parameters=NULL;
375
376 /*Get Element type*/
377 this->hnodes->reset();
378 this->helements->reset();
379 this->hmatpar->reset();
380
381}
382/*}}}*/
383void Riftfront::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
384
385}
386/*}}}*/
387void Riftfront::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/
388
389 /*Output */
390 int d_nz = 0;
391 int o_nz = 0;
392
393 /*Loop over all nodes*/
394 for(int i=0;i<NUMVERTICES;i++){
395
396 if(!flags[this->nodes[i]->Lid()]){
397
398 /*flag current node so that no other element processes it*/
399 flags[this->nodes[i]->Lid()]=true;
400
401 int counter=0;
402 while(flagsindices[counter]>=0) counter++;
403 flagsindices[counter]=this->nodes[i]->Lid();
404
405 /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
406 switch(set2_enum){
407 case FsetEnum:
408 if(nodes[i]->indexing.fsize){
409 if(this->nodes[i]->IsClone())
410 o_nz += 1;
411 else
412 d_nz += 1;
413 }
414 break;
415 case GsetEnum:
416 if(nodes[i]->indexing.gsize){
417 if(this->nodes[i]->IsClone())
418 o_nz += 1;
419 else
420 d_nz += 1;
421 }
422 break;
423 case SsetEnum:
424 if(nodes[i]->indexing.ssize){
425 if(this->nodes[i]->IsClone())
426 o_nz += 1;
427 else
428 d_nz += 1;
429 }
430 break;
431 default: _error_("not supported");
432 }
433 }
434 }
435
436 /*Assign output pointers: */
437 *pd_nz=d_nz;
438 *po_nz=o_nz;
439}
440/*}}}*/
441
442/*Riftfront numerics*/
443ElementMatrix* Riftfront::PenaltyCreateKMatrixStressbalanceHoriz(IssmDouble kmax){/*{{{*/
444
445 const int numdof = NDOF2*NUMVERTICES;
446 IssmDouble thickness;
447 IssmDouble h[2];
448 IssmDouble penalty_offset;
449
450 /*Objects: */
451 Tria *tria1 = NULL;
452 Tria *tria2 = NULL;
453
454 /*enum of element? */
455 if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
456 tria1=(Tria*)elements[0];
457 tria2=(Tria*)elements[1];
458
459 /*Initialize Element Matrix*/
460 if(!this->active) return NULL;
461 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters);
462
463 /*Get some parameters: */
464 this->parameters->FindParam(&penalty_offset,StressbalancePenaltyFactorEnum);
465 tria1->GetInputValue(&h[0],nodes[0],ThicknessEnum);
466 tria2->GetInputValue(&h[1],nodes[1],ThicknessEnum);
467 if (h[0]!=h[1])_error_("different thicknesses not supported for rift fronts");
468 thickness=h[0];
469
470 /*There is contact, we need to constrain the normal velocities (zero penetration), and the
471 *contact slip friction. */
472
473 /*From Peter Wriggers book (Computational Contact Mechanics, p191): */
474 Ke->values[0*numdof+0]+= +pow(normal[0],2)*kmax*pow(10,penalty_offset);
475 Ke->values[0*numdof+1]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset);
476 Ke->values[0*numdof+2]+= -pow(normal[0],2)*kmax*pow(10,penalty_offset);
477 Ke->values[0*numdof+3]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset);
478
479 Ke->values[1*numdof+0]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset);
480 Ke->values[1*numdof+1]+= +pow(normal[1],2)*kmax*pow(10,penalty_offset);
481 Ke->values[1*numdof+2]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset);
482 Ke->values[1*numdof+3]+= -pow(normal[1],2)*kmax*pow(10,penalty_offset);
483
484 Ke->values[2*numdof+0]+= -pow(normal[0],2)*kmax*pow(10,penalty_offset);
485 Ke->values[2*numdof+1]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset);
486 Ke->values[2*numdof+2]+= +pow(normal[0],2)*kmax*pow(10,penalty_offset);
487 Ke->values[2*numdof+3]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset);
488
489 Ke->values[3*numdof+0]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset);
490 Ke->values[3*numdof+1]+= -pow(normal[1],2)*kmax*pow(10,penalty_offset);
491 Ke->values[3*numdof+2]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset);
492 Ke->values[3*numdof+3]+= +pow(normal[1],2)*kmax*pow(10,penalty_offset);
493
494 /*Now take care of the friction: of type sigma=frictiontangent_velocity2-tangent_velocity1)*/
495
496 Ke->values[0*numdof+0]+= +pow(normal[1],2)*thickness*length*friction;
497 Ke->values[0*numdof+1]+= -normal[0]*normal[1]*thickness*length*friction;
498 Ke->values[0*numdof+2]+= -pow(normal[1],2)*thickness*length*friction;
499 Ke->values[0*numdof+3]+= +normal[0]*normal[1]*thickness*length*friction;
500
501 Ke->values[1*numdof+0]+= -normal[0]*normal[1]*thickness*length*friction;
502 Ke->values[1*numdof+1]+= +pow(normal[0],2)*thickness*length*friction;
503 Ke->values[1*numdof+2]+= +normal[0]*normal[1]*thickness*length*friction;
504 Ke->values[1*numdof+3]+= -pow(normal[0],2)*thickness*length*friction;
505
506 Ke->values[2*numdof+0]+= -pow(normal[1],2)*thickness*length*friction;
507 Ke->values[2*numdof+1]+= +normal[0]*normal[1]*thickness*length*friction;
508 Ke->values[2*numdof+2]+= +pow(normal[1],2)*thickness*length*friction;
509 Ke->values[2*numdof+3]+= -normal[0]*normal[1]*thickness*length*friction;
510
511 Ke->values[3*numdof+0]+= +normal[0]*normal[1]*thickness*length*friction;
512 Ke->values[3*numdof+1]+= -pow(normal[0],2)*thickness*length*friction;
513 Ke->values[3*numdof+2]+= -normal[0]*normal[1]*thickness*length*friction;
514 Ke->values[3*numdof+3]+= +pow(normal[0],2)*thickness*length*friction;
515
516 /*Clean up and return*/
517 return Ke;
518}
519/*}}}*/
520ElementVector* Riftfront::PenaltyCreatePVectorStressbalanceHoriz(IssmDouble kmax){/*{{{*/
521
522 int j;
523 IssmDouble rho_ice;
524 IssmDouble rho_water;
525 IssmDouble gravity;
526 IssmDouble thickness;
527 IssmDouble h[2];
528 IssmDouble bed;
529 IssmDouble b[2];
530 IssmDouble pressure;
531 IssmDouble pressure_litho;
532 IssmDouble pressure_air;
533 IssmDouble pressure_melange;
534 IssmDouble pressure_water;
535
536 /*Objects: */
537 Tria *tria1 = NULL;
538 Tria *tria2 = NULL;
539
540 /*enum of element? */
541 if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
542 tria1=(Tria*)elements[0];
543 tria2=(Tria*)elements[1];
544
545 /*Initialize Element Matrix*/
546 if(this->active) return NULL; /*The penalty is active. No loads implied here.*/
547 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
548
549 /*Get some inputs: */
550 rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
551 rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
552 gravity=matpar->GetMaterialParameter(ConstantsGEnum);
553 tria1->GetInputValue(&h[0],nodes[0],ThicknessEnum);
554 tria2->GetInputValue(&h[1],nodes[1],ThicknessEnum);
555 if (h[0]!=h[1])_error_("different thicknesses not supported for rift fronts");
556 thickness=h[0];
557 tria1->GetInputValue(&b[0],nodes[0],BaseEnum);
558 tria2->GetInputValue(&b[1],nodes[1],BaseEnum);
559 if (b[0]!=b[1])_error_("different beds not supported for rift fronts");
560 bed=b[0];
561
562 /*Ok, this rift is opening. We should put loads on both sides of the rift flanks. Because we are dealing with contact mechanics,
563 * and we want to avoid zigzagging of the loads, we want lump the loads onto nodes, not onto surfaces between nodes.:*/
564
565 /*Ok, to compute the pressure, we are going to need material properties, thickness and bed for the two nodes. We assume those properties to
566 * be the same across the rift.: */
567
568 /*Ok, now compute the pressure (in norm) that is being applied to the flanks, depending on the type of fill: */
569 if(fill==WaterEnum){
570 if(shelf){
571 /*We are on an ice shelf, hydrostatic equilibrium is used to determine the pressure for water fill: */
572 pressure=rho_ice*gravity*pow(thickness,2)/2.- rho_water*gravity*pow(bed,2)/2.;
573 }
574 else{
575 //We are on an icesheet, we assume the water column fills the entire front: */
576 pressure=rho_ice*gravity*pow(thickness,2)/2.- rho_water*gravity*pow(thickness,2)/2.;
577 }
578 }
579 else if(fill==AirEnum){
580 pressure=rho_ice*gravity*pow(thickness,2)/2.; //icefront on an ice sheet, pressure imbalance ice vs air.
581 }
582 else if(fill==IceEnum){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
583 pressure=0;
584 }
585 else if(fill==MelangeEnum){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
586
587 if(!shelf) _error_("fill type " << fill << " not supported on ice sheets yet.");
588
589 pressure_litho=rho_ice*gravity*pow(thickness,2)/2.;
590 pressure_air=0;
591 pressure_melange=rho_ice*gravity*pow(fraction*thickness,2)/2.;
592 pressure_water=1.0/2.0*rho_water*gravity* ( pow(bed,2.0)-pow(rho_ice/rho_water*fraction*thickness,2.0) );
593
594 pressure=pressure_litho-pressure_air-pressure_melange-pressure_water;
595 }
596 else{
597 _error_("fill type " << fill << " not supported yet.");
598 }
599
600 /*Ok, add contribution to first node, along the normal i==0: */
601 for (j=0;j<2;j++){
602 pe->values[j]+=pressure*normal[j]*length;
603 }
604
605 /*Add contribution to second node, along the opposite normal: i==1 */
606 for (j=0;j<2;j++){
607 pe->values[2+j]+= -pressure*normal[j]*length;
608 }
609
610 /*Clean up and return*/
611 return pe;
612}
613/*}}}*/
614#define _ZIGZAGCOUNTER_
615int Riftfront::Constrain(int* punstable){/*{{{*/
616
617 IssmDouble penetration;
618 bool activate;
619 int unstable;
620 IssmDouble vx1;
621 IssmDouble vy1;
622 IssmDouble vx2;
623 IssmDouble vy2;
624
625 /*Objects: */
626 Tria *tria1 = NULL;
627 Tria *tria2 = NULL;
628
629 /*enum of element? */
630 if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
631
632 /*recover elements on both side of rift: */
633 tria1=(Tria*)elements[0];
634 tria2=(Tria*)elements[1];
635
636 /*Is this constraint frozen? In which case we don't touch: */
637 if (this->frozen){
638 *punstable=0;
639 return 1;
640 }
641
642 /*Is this rift segment state specified by user input? :*/
643 if (this->state==OpenEnum || this->state==ClosedEnum){
644
645 if(this->state==OpenEnum)this->active=0;
646 if(this->state==ClosedEnum)this->active=1;
647
648 /*this segment is like frozen, no instability here: */
649 *punstable=0;
650 return 1;
651 }
652
653 /*First recover velocity: */
654 tria1->GetInputValue(&vx1,nodes[0],VxEnum);
655 tria2->GetInputValue(&vx2,nodes[1],VxEnum);
656 tria1->GetInputValue(&vy1,nodes[0],VyEnum);
657 tria2->GetInputValue(&vy2,nodes[1],VyEnum);
658
659 /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
660 penetration=(vx2-vx1)*normal[0]+(vy2-vy1)*normal[1];
661
662 /*activation: */
663 if(penetration<0)activate=true;
664 else activate=false;
665
666 /*Here, we try to avoid zigzaging. When a penalty activates and deactivates for more than penalty_lock times,
667 * we increase the fraction of melange:*/
668 if(this->counter>this->penalty_lock){
669 /*reset counter: */
670 this->counter=0;
671 /*increase melange fraction: */
672 this->fraction+=fractionincrement;
673 if (this->fraction>1)this->fraction=1.;
674 //_printf_("riftfront " << this->Id() << " fraction: " << this->fraction << "\n");
675 }
676
677 //Figure out stability of this penalty
678 if(this->active==activate){
679 unstable=0;
680 }
681 else{
682 unstable=1;
683 this->counter++;
684 }
685
686 //Set penalty flag
687 this->active=activate;
688
689 //if ((penetration>0) && (this->active==1))_printf_("Riftfront " << Id() << " wants to be released\n");
690
691 /*assign output pointer: */
692 *punstable=unstable;
693 return 1;
694}
695/*}}}*/
696void Riftfront::FreezeConstraints(void){/*{{{*/
697
698 /*Just set frozen flag to 1: */
699 this->frozen=1;
700
701}
702/*}}}*/
703bool Riftfront::IsFrozen(void){/*{{{*/
704
705 /*Just set frozen flag to 1: */
706 if(this->frozen)return 1;
707 else return 0;
708}
709/*}}}*/
Note: See TracBrowser for help on using the repository browser.