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

Last change on this file since 17806 was 17806, checked in by Mathieu Morlighem, 11 years ago

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