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

Last change on this file since 26744 was 26744, checked in by Mathieu Morlighem, 3 years ago

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