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

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

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