Ice Sheet System Model  4.18
Code documentation
Public Member Functions | Data Fields
Riftfront Class Reference

#include <Riftfront.h>

Inheritance diagram for Riftfront:
Load Object

Public Member Functions

 Riftfront ()
 
 Riftfront (int riftfront_id, int i, IoModel *iomodel)
 
 ~Riftfront ()
 
Object * copy ()
 
void DeepEcho ()
 
void Echo ()
 
int Id ()
 
void Marshall (char **pmarshalled_data, int *pmarshalled_data_size, int marshall_direction)
 
int ObjectEnum ()
 
void InputUpdateFromConstant (IssmDouble constant, int name)
 
void InputUpdateFromConstant (int constant, int name)
 
void InputUpdateFromConstant (bool constant, int name)
 
void InputUpdateFromIoModel (int index, IoModel *iomodel)
 
void InputUpdateFromMatrixDakota (IssmDouble *matrix, int nrows, int ncols, int name, int type)
 
void InputUpdateFromVector (IssmDouble *vector, int name, int type)
 
void InputUpdateFromVectorDakota (IssmDouble *vector, int name, int type)
 
void Configure (Elements *elements, Loads *loads, Nodes *nodes, Vertices *vertices, Materials *materials, Parameters *parameters)
 
void CreateJacobianMatrix (Matrix< IssmDouble > *Jff)
 
void CreateKMatrix (Matrix< IssmDouble > *Kff, Matrix< IssmDouble > *Kfs)
 
void CreatePVector (Vector< IssmDouble > *pf)
 
void GetNodesLidList (int *lidlist)
 
void GetNodesSidList (int *sidlist)
 
int GetNumberOfNodes (void)
 
bool IsPenalty (void)
 
void PenaltyCreateJacobianMatrix (Matrix< IssmDouble > *Jff, IssmDouble kmax)
 
void PenaltyCreateKMatrix (Matrix< IssmDouble > *Kff, Matrix< IssmDouble > *kfs, IssmDouble kmax)
 
void PenaltyCreatePVector (Vector< IssmDouble > *pf, IssmDouble kmax)
 
void ResetHooks ()
 
void SetCurrentConfiguration (Elements *elements, Loads *loads, Nodes *nodes, Vertices *vertices, Materials *materials, Parameters *parameters)
 
void SetwiseNodeConnectivity (int *d_nz, int *o_nz, Node *node, bool *flags, int *flagsindices, int set1_enum, int set2_enum)
 
int Constrain (int *punstable)
 
void FreezeConstraints (void)
 
bool IsFrozen (void)
 
ElementMatrix * PenaltyCreateKMatrixStressbalanceHoriz (IssmDouble kmax)
 
ElementVector * PenaltyCreatePVectorStressbalanceHoriz (IssmDouble kmax)
 
- Public Member Functions inherited from Load
virtual ~Load ()
 
- Public Member Functions inherited from Object
virtual ~Object ()
 

Data Fields

int id
 
int type
 
int fill
 
IssmDouble friction
 
IssmDouble fractionincrement
 
bool shelf
 
Hook * hnodes
 
Hook * hvertices
 
Hook * helements
 
Node ** nodes
 
Vertex ** vertices
 
Element ** elements
 
int penalty_lock
 
bool active
 
bool frozen
 
int counter
 
bool prestable
 
bool material_converged
 
IssmDouble normal [2]
 
IssmDouble length
 
IssmDouble fraction
 
int state
 
Parameters * parameters
 

Detailed Description

Definition at line 16 of file Riftfront.h.

Constructor & Destructor Documentation

◆ Riftfront() [1/2]

Riftfront::Riftfront ( )

Definition at line 22 of file Riftfront.cpp.

22  {/*{{{*/
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 }

◆ Riftfront() [2/2]

Riftfront::Riftfront ( int  riftfront_id,
int  i,
IoModel *  iomodel 
)

Definition at line 32 of file Riftfront.cpp.

32  {/*{{{*/
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 }

◆ ~Riftfront()

Riftfront::~Riftfront ( )

Definition at line 98 of file Riftfront.cpp.

98  {/*{{{*/
99  this->parameters=NULL;
100  delete hnodes;
101  delete hvertices;
102  delete helements;
103 }

Member Function Documentation

◆ copy()

Object * Riftfront::copy ( void  )
virtual

Implements Object.

Definition at line 107 of file Riftfront.cpp.

107  {/*{{{*/
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 }

◆ DeepEcho()

void Riftfront::DeepEcho ( void  )
virtual

Implements Object.

Definition at line 151 of file Riftfront.cpp.

151  {/*{{{*/
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 }

◆ Echo()

void Riftfront::Echo ( void  )
virtual

Implements Object.

Definition at line 162 of file Riftfront.cpp.

162  {/*{{{*/
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 }

◆ Id()

int Riftfront::Id ( void  )
virtual

Implements Object.

Definition at line 187 of file Riftfront.cpp.

187 { return id; }/*{{{*/

◆ Marshall()

void Riftfront::Marshall ( char **  pmarshalled_data,
int *  pmarshalled_data_size,
int  marshall_direction 
)
virtual

Implements Object.

Definition at line 189 of file Riftfront.cpp.

189  { /*{{{*/
190 
191  _assert_(this);
192 
193  /*ok, marshall operations: */
194  MARSHALLING_ENUM(RiftfrontEnum);
195  MARSHALLING(id);
196  MARSHALLING(type);
197  MARSHALLING(fill);
198  MARSHALLING(friction);
199  MARSHALLING(fractionincrement);
200  MARSHALLING(shelf);
201 
202  if(marshall_direction==MARSHALLING_BACKWARD){
203  this->hnodes = new Hook();
204  this->hvertices = new Hook();
205  this->helements = new Hook();
206  }
207 
208  this->hnodes->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
209  this->hvertices->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
210  this->helements->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
211 
212  /*corresponding fields*/
213  nodes =(Node**)this->hnodes->deliverp();
214  vertices =(Vertex**)this->hvertices->deliverp();
215  elements =(Element**)this->helements->deliverp();
216 
217  MARSHALLING(penalty_lock);
218  MARSHALLING(active);
219  MARSHALLING(frozen);
220  MARSHALLING(state);
221  MARSHALLING(counter);
222  MARSHALLING(prestable);
223  MARSHALLING(material_converged);
224  MARSHALLING(normal[0]);
225  MARSHALLING(normal[1]);
226  MARSHALLING(length);
227  MARSHALLING(fraction);
228 
229 }

◆ ObjectEnum()

int Riftfront::ObjectEnum ( void  )
virtual

Implements Object.

Definition at line 231 of file Riftfront.cpp.

231  {/*{{{*/
232 
233  return RiftfrontEnum;
234 
235 }

◆ InputUpdateFromConstant() [1/3]

void Riftfront::InputUpdateFromConstant ( IssmDouble  constant,
int  name 
)

Definition at line 242 of file Riftfront.cpp.

242  {/*{{{*/
243 
244 }

◆ InputUpdateFromConstant() [2/3]

void Riftfront::InputUpdateFromConstant ( int  constant,
int  name 
)
inline

Definition at line 67 of file Riftfront.h.

67 {_error_("Not implemented yet!");}

◆ InputUpdateFromConstant() [3/3]

void Riftfront::InputUpdateFromConstant ( bool  constant,
int  name 
)

Definition at line 239 of file Riftfront.cpp.

239  {/*{{{*/
240 }

◆ InputUpdateFromIoModel()

void Riftfront::InputUpdateFromIoModel ( int  index,
IoModel *  iomodel 
)
inline

Definition at line 69 of file Riftfront.h.

69 {_error_("not implemented yet");};

◆ InputUpdateFromMatrixDakota()

void Riftfront::InputUpdateFromMatrixDakota ( IssmDouble *  matrix,
int  nrows,
int  ncols,
int  name,
int  type 
)
inline

Definition at line 70 of file Riftfront.h.

70 {_error_("Not implemented yet!");}

◆ InputUpdateFromVector()

void Riftfront::InputUpdateFromVector ( IssmDouble *  vector,
int  name,
int  type 
)

Definition at line 246 of file Riftfront.cpp.

246  {/*{{{*/
247 
248  /*Nothing to update*/
249 
250 }

◆ InputUpdateFromVectorDakota()

void Riftfront::InputUpdateFromVectorDakota ( IssmDouble *  vector,
int  name,
int  type 
)
inline

Definition at line 72 of file Riftfront.h.

72 {_error_("Not implemented yet!");}

◆ Configure()

void Riftfront::Configure ( Elements *  elements,
Loads *  loads,
Nodes *  nodes,
Vertices *  vertices,
Materials *  materials,
Parameters *  parameters 
)
virtual

Implements Load.

Definition at line 254 of file Riftfront.cpp.

254  {/*{{{*/
255 
256  /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective
257  * datasets, using internal ids and offsets hidden in hooks: */
258  hnodes->configure(nodesin);
259  hvertices->configure(verticesin);
260  helements->configure(elementsin);
261 
262  /*Initialize hooked fields*/
263  this->nodes =(Node**)hnodes->deliverp();
264  this->vertices=(Vertex**)hvertices->deliverp();
265  this->elements=(Element**)helements->deliverp();
266 
267  /*point parameters to real dataset: */
268  this->parameters=parametersin;
269 
270 }

◆ CreateJacobianMatrix()

void Riftfront::CreateJacobianMatrix ( Matrix< IssmDouble > *  Jff)
inlinevirtual

Implements Load.

Definition at line 76 of file Riftfront.h.

76 {_error_("Not implemented yet");};

◆ CreateKMatrix()

void Riftfront::CreateKMatrix ( Matrix< IssmDouble > *  Kff,
Matrix< IssmDouble > *  Kfs 
)
virtual

Implements Load.

Definition at line 272 of file Riftfront.cpp.

272  {/*{{{*/
273  /*do nothing: */
274  return;
275 }

◆ CreatePVector()

void Riftfront::CreatePVector ( Vector< IssmDouble > *  pf)
virtual

Implements Load.

Definition at line 277 of file Riftfront.cpp.

277  {/*{{{*/
278  /*do nothing: */
279  return;
280 }

◆ GetNodesLidList()

void Riftfront::GetNodesLidList ( int *  lidlist)
virtual

Implements Load.

Definition at line 282 of file Riftfront.cpp.

282  {/*{{{*/
283 
284  _assert_(lidlist);
285  _assert_(nodes);
286 
287  for(int i=0;i<NUMVERTICES;i++) lidlist[i]=nodes[i]->Lid();
288 }

◆ GetNodesSidList()

void Riftfront::GetNodesSidList ( int *  sidlist)
virtual

Implements Load.

Definition at line 290 of file Riftfront.cpp.

290  {/*{{{*/
291 
292  _assert_(sidlist);
293  _assert_(nodes);
294 
295  for(int i=0;i<NUMVERTICES;i++) sidlist[i]=nodes[i]->Sid();
296 }

◆ GetNumberOfNodes()

int Riftfront::GetNumberOfNodes ( void  )
virtual

Implements Load.

Definition at line 298 of file Riftfront.cpp.

298  {/*{{{*/
299 
300  return NUMVERTICES;
301 }

◆ IsPenalty()

bool Riftfront::IsPenalty ( void  )
virtual

Implements Load.

Definition at line 303 of file Riftfront.cpp.

303  {/*{{{*/
304  return true;
305 }

◆ PenaltyCreateJacobianMatrix()

void Riftfront::PenaltyCreateJacobianMatrix ( Matrix< IssmDouble > *  Jff,
IssmDouble  kmax 
)
inlinevirtual

Implements Load.

Definition at line 83 of file Riftfront.h.

83 {_error_("Not implemented yet");};

◆ PenaltyCreateKMatrix()

void Riftfront::PenaltyCreateKMatrix ( Matrix< IssmDouble > *  Kff,
Matrix< IssmDouble > *  kfs,
IssmDouble  kmax 
)
virtual

Implements Load.

Definition at line 307 of file Riftfront.cpp.

307  {/*{{{*/
308 
309  /*Retrieve parameters: */
310  ElementMatrix* Ke=NULL;
311  int analysis_type;
312  this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
313 
314  switch(analysis_type){
315  case StressbalanceAnalysisEnum:
316  Ke=PenaltyCreateKMatrixStressbalanceHoriz(kmax);
317  break;
318  case AdjointHorizAnalysisEnum:
319  Ke=PenaltyCreateKMatrixStressbalanceHoriz(kmax);
320  break;
321  default:
322  _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
323  }
324 
325  /*Add to global Vector*/
326  if(Ke){
327  Ke->AddToGlobal(Kff,Kfs);
328  delete Ke;
329  }
330 }

◆ PenaltyCreatePVector()

void Riftfront::PenaltyCreatePVector ( Vector< IssmDouble > *  pf,
IssmDouble  kmax 
)
virtual

Implements Load.

Definition at line 332 of file Riftfront.cpp.

332  {/*{{{*/
333 
334  /*Retrieve parameters: */
335  ElementVector* pe=NULL;
336  int analysis_type;
337  this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
338 
339  switch(analysis_type){
340  case StressbalanceAnalysisEnum:
341  pe=PenaltyCreatePVectorStressbalanceHoriz(kmax);
342  break;
343  case AdjointHorizAnalysisEnum:
344  /*No penalty applied on load vector*/
345  break;
346  default:
347  _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
348  }
349 
350  /*Add to global Vector*/
351  if(pe){
352  pe->AddToGlobal(pf);
353  delete pe;
354  }
355 }

◆ ResetHooks()

void Riftfront::ResetHooks ( )
virtual

Implements Load.

Definition at line 357 of file Riftfront.cpp.

357  {/*{{{*/
358 
359  this->nodes=NULL;
360  this->vertices=NULL;
361  this->elements=NULL;
362  this->parameters=NULL;
363 
364  /*Get Element type*/
365  this->hnodes->reset();
366  this->hvertices->reset();
367  this->helements->reset();
368 
369 }

◆ SetCurrentConfiguration()

void Riftfront::SetCurrentConfiguration ( Elements *  elements,
Loads *  loads,
Nodes *  nodes,
Vertices *  vertices,
Materials *  materials,
Parameters *  parameters 
)
virtual

Implements Load.

Definition at line 371 of file Riftfront.cpp.

371  {/*{{{*/
372 
373 }

◆ SetwiseNodeConnectivity()

void Riftfront::SetwiseNodeConnectivity ( int *  d_nz,
int *  o_nz,
Node *  node,
bool *  flags,
int *  flagsindices,
int  set1_enum,
int  set2_enum 
)
virtual

Implements Load.

Definition at line 375 of file Riftfront.cpp.

375  {/*{{{*/
376 
377  /*Output */
378  int d_nz = 0;
379  int o_nz = 0;
380 
381  /*Loop over all nodes*/
382  for(int i=0;i<NUMVERTICES;i++){
383 
384  if(!flags[this->nodes[i]->Lid()]){
385 
386  /*flag current node so that no other element processes it*/
387  flags[this->nodes[i]->Lid()]=true;
388 
389  int counter=0;
390  while(flagsindices[counter]>=0) counter++;
391  flagsindices[counter]=this->nodes[i]->Lid();
392 
393  /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
394  switch(set2_enum){
395  case FsetEnum:
396  if(nodes[i]->fsize){
397  if(this->nodes[i]->IsClone())
398  o_nz += 1;
399  else
400  d_nz += 1;
401  }
402  break;
403  case GsetEnum:
404  if(nodes[i]->gsize){
405  if(this->nodes[i]->IsClone())
406  o_nz += 1;
407  else
408  d_nz += 1;
409  }
410  break;
411  case SsetEnum:
412  if(nodes[i]->ssize){
413  if(this->nodes[i]->IsClone())
414  o_nz += 1;
415  else
416  d_nz += 1;
417  }
418  break;
419  default: _error_("not supported");
420  }
421  }
422  }
423 
424  /*Assign output pointers: */
425  *pd_nz=d_nz;
426  *po_nz=o_nz;
427 }

◆ Constrain()

int Riftfront::Constrain ( int *  punstable)

Definition at line 594 of file Riftfront.cpp.

594  {/*{{{*/
595 
596  IssmDouble penetration;
597  bool activate;
598  int unstable;
599  IssmDouble vx1;
600  IssmDouble vy1;
601  IssmDouble vx2;
602  IssmDouble vy2;
603 
604  /*Objects: */
605  Tria *tria1 = NULL;
606  Tria *tria2 = NULL;
607 
608  /*enum of element? */
609  if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
610 
611  /*recover elements on both side of rift: */
612  tria1=(Tria*)elements[0];
613  tria2=(Tria*)elements[1];
614 
615  /*Is this constraint frozen? In which case we don't touch: */
616  if (this->frozen){
617  *punstable=0;
618  return 1;
619  }
620 
621  /*Is this rift segment state specified by user input? :*/
622  if (this->state==OpenEnum || this->state==ClosedEnum){
623 
624  if(this->state==OpenEnum)this->active=0;
625  if(this->state==ClosedEnum)this->active=1;
626 
627  /*this segment is like frozen, no instability here: */
628  *punstable=0;
629  return 1;
630  }
631 
632  /*First recover velocity: */
633  tria1->GetInputValue(&vx1,vertices[0],VxEnum);
634  tria2->GetInputValue(&vx2,vertices[1],VxEnum);
635  tria1->GetInputValue(&vy1,vertices[0],VyEnum);
636  tria2->GetInputValue(&vy2,vertices[1],VyEnum);
637 
638  /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
639  penetration=(vx2-vx1)*normal[0]+(vy2-vy1)*normal[1];
640 
641  /*activation: */
642  if(penetration<0)activate=true;
643  else activate=false;
644 
645  /*Here, we try to avoid zigzaging. When a penalty activates and deactivates for more than penalty_lock times,
646  * we increase the fraction of melange:*/
647  if(this->counter>this->penalty_lock){
648  /*reset counter: */
649  this->counter=0;
650  /*increase melange fraction: */
651  this->fraction+=fractionincrement;
652  if (this->fraction>1)this->fraction=1.;
653  //_printf_("riftfront " << this->Id() << " fraction: " << this->fraction << "\n");
654  }
655 
656  //Figure out stability of this penalty
657  if(this->active==activate){
658  unstable=0;
659  }
660  else{
661  unstable=1;
662  this->counter++;
663  }
664 
665  //Set penalty flag
666  this->active=activate;
667 
668  //if ((penetration>0) && (this->active==1))_printf_("Riftfront " << Id() << " wants to be released\n");
669 
670  /*assign output pointer: */
671  *punstable=unstable;
672  return 1;
673 }

◆ FreezeConstraints()

void Riftfront::FreezeConstraints ( void  )

Definition at line 675 of file Riftfront.cpp.

675  {/*{{{*/
676 
677  /*Just set frozen flag to 1: */
678  this->frozen=1;
679 
680 }

◆ IsFrozen()

bool Riftfront::IsFrozen ( void  )

Definition at line 682 of file Riftfront.cpp.

682  {/*{{{*/
683 
684  /*Just set frozen flag to 1: */
685  if(this->frozen)return 1;
686  else return 0;
687 }

◆ PenaltyCreateKMatrixStressbalanceHoriz()

ElementMatrix * Riftfront::PenaltyCreateKMatrixStressbalanceHoriz ( IssmDouble  kmax)

Definition at line 431 of file Riftfront.cpp.

431  {/*{{{*/
432 
433  const int numdof = 2*NUMVERTICES;
434  IssmDouble thickness;
435  IssmDouble h[2];
436  IssmDouble penalty_offset;
437 
438  /*enum of element? */
439  if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
440  Tria* tria1=(Tria*)elements[0];
441  Tria* tria2=(Tria*)elements[1];
442 
443  /*Initialize Element Matrix*/
444  if(!this->active) return NULL;
445  ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters);
446 
447  /*Get some parameters: */
448  this->parameters->FindParam(&penalty_offset,StressbalancePenaltyFactorEnum);
449  tria1->GetInputValue(&h[0],vertices[0],ThicknessEnum);
450  tria2->GetInputValue(&h[1],vertices[1],ThicknessEnum);
451  if (h[0]!=h[1])_error_("different thicknesses not supported for rift fronts");
452  thickness=h[0];
453 
454  /*There is contact, we need to constrain the normal velocities (zero penetration), and the
455  *contact slip friction. */
456 
457  /*From Peter Wriggers book (Computational Contact Mechanics, p191): */
458  Ke->values[0*numdof+0]+= +pow(normal[0],2)*kmax*pow(10,penalty_offset);
459  Ke->values[0*numdof+1]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset);
460  Ke->values[0*numdof+2]+= -pow(normal[0],2)*kmax*pow(10,penalty_offset);
461  Ke->values[0*numdof+3]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset);
462 
463  Ke->values[1*numdof+0]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset);
464  Ke->values[1*numdof+1]+= +pow(normal[1],2)*kmax*pow(10,penalty_offset);
465  Ke->values[1*numdof+2]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset);
466  Ke->values[1*numdof+3]+= -pow(normal[1],2)*kmax*pow(10,penalty_offset);
467 
468  Ke->values[2*numdof+0]+= -pow(normal[0],2)*kmax*pow(10,penalty_offset);
469  Ke->values[2*numdof+1]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset);
470  Ke->values[2*numdof+2]+= +pow(normal[0],2)*kmax*pow(10,penalty_offset);
471  Ke->values[2*numdof+3]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset);
472 
473  Ke->values[3*numdof+0]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset);
474  Ke->values[3*numdof+1]+= -pow(normal[1],2)*kmax*pow(10,penalty_offset);
475  Ke->values[3*numdof+2]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset);
476  Ke->values[3*numdof+3]+= +pow(normal[1],2)*kmax*pow(10,penalty_offset);
477 
478  /*Now take care of the friction: of type sigma=frictiontangent_velocity2-tangent_velocity1)*/
479 
480  Ke->values[0*numdof+0]+= +pow(normal[1],2)*thickness*length*friction;
481  Ke->values[0*numdof+1]+= -normal[0]*normal[1]*thickness*length*friction;
482  Ke->values[0*numdof+2]+= -pow(normal[1],2)*thickness*length*friction;
483  Ke->values[0*numdof+3]+= +normal[0]*normal[1]*thickness*length*friction;
484 
485  Ke->values[1*numdof+0]+= -normal[0]*normal[1]*thickness*length*friction;
486  Ke->values[1*numdof+1]+= +pow(normal[0],2)*thickness*length*friction;
487  Ke->values[1*numdof+2]+= +normal[0]*normal[1]*thickness*length*friction;
488  Ke->values[1*numdof+3]+= -pow(normal[0],2)*thickness*length*friction;
489 
490  Ke->values[2*numdof+0]+= -pow(normal[1],2)*thickness*length*friction;
491  Ke->values[2*numdof+1]+= +normal[0]*normal[1]*thickness*length*friction;
492  Ke->values[2*numdof+2]+= +pow(normal[1],2)*thickness*length*friction;
493  Ke->values[2*numdof+3]+= -normal[0]*normal[1]*thickness*length*friction;
494 
495  Ke->values[3*numdof+0]+= +normal[0]*normal[1]*thickness*length*friction;
496  Ke->values[3*numdof+1]+= -pow(normal[0],2)*thickness*length*friction;
497  Ke->values[3*numdof+2]+= -normal[0]*normal[1]*thickness*length*friction;
498  Ke->values[3*numdof+3]+= +pow(normal[0],2)*thickness*length*friction;
499 
500  /*Clean up and return*/
501  return Ke;
502 }

◆ PenaltyCreatePVectorStressbalanceHoriz()

ElementVector * Riftfront::PenaltyCreatePVectorStressbalanceHoriz ( IssmDouble  kmax)

Definition at line 504 of file Riftfront.cpp.

504  {/*{{{*/
505 
506  int j;
507  IssmDouble rho_ice;
508  IssmDouble rho_water;
509  IssmDouble gravity;
510  IssmDouble thickness;
511  IssmDouble h[2];
512  IssmDouble bed;
513  IssmDouble b[2];
514  IssmDouble pressure;
515  IssmDouble pressure_litho;
516  IssmDouble pressure_air;
517  IssmDouble pressure_melange;
518  IssmDouble pressure_water;
519 
520  /*enum of element? */
521  if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
522  Tria* tria1=(Tria*)elements[0];
523  Tria* tria2=(Tria*)elements[1];
524 
525  /*Initialize Element Matrix*/
526  if(this->active) return NULL; /*The penalty is active. No loads implied here.*/
527  ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
528 
529  /*Get some inputs: */
530  rho_ice=tria1->FindParam(MaterialsRhoIceEnum);
531  rho_water=tria1->FindParam(MaterialsRhoSeawaterEnum);
532  gravity=tria1->FindParam(ConstantsGEnum);
533  tria1->GetInputValue(&h[0],vertices[0],ThicknessEnum);
534  tria2->GetInputValue(&h[1],vertices[1],ThicknessEnum);
535  if (h[0]!=h[1])_error_("different thicknesses not supported for rift fronts");
536  thickness=h[0];
537  tria1->GetInputValue(&b[0],vertices[0],BaseEnum);
538  tria2->GetInputValue(&b[1],vertices[1],BaseEnum);
539  if (b[0]!=b[1])_error_("different beds not supported for rift fronts");
540  bed=b[0];
541 
542  /*Ok, this rift is opening. We should put loads on both sides of the rift flanks. Because we are dealing with contact mechanics,
543  * and we want to avoid zigzagging of the loads, we want lump the loads onto nodes, not onto surfaces between nodes.:*/
544 
545  /*Ok, to compute the pressure, we are going to need material properties, thickness and bed for the two nodes. We assume those properties to
546  * be the same across the rift.: */
547 
548  /*Ok, now compute the pressure (in norm) that is being applied to the flanks, depending on the type of fill: */
549  if(fill==WaterEnum){
550  if(shelf){
551  /*We are on an ice shelf, hydrostatic equilibrium is used to determine the pressure for water fill: */
552  pressure=rho_ice*gravity*pow(thickness,2)/2.- rho_water*gravity*pow(bed,2)/2.;
553  }
554  else{
555  //We are on an icesheet, we assume the water column fills the entire front: */
556  pressure=rho_ice*gravity*pow(thickness,2)/2.- rho_water*gravity*pow(thickness,2)/2.;
557  }
558  }
559  else if(fill==AirEnum){
560  pressure=rho_ice*gravity*pow(thickness,2)/2.; //icefront on an ice sheet, pressure imbalance ice vs air.
561  }
562  else if(fill==IceEnum){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
563  pressure=0;
564  }
565  else if(fill==MelangeEnum){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
566 
567  if(!shelf) _error_("fill type " << fill << " not supported on ice sheets yet.");
568 
569  pressure_litho=rho_ice*gravity*pow(thickness,2)/2.;
570  pressure_air=0;
571  pressure_melange=rho_ice*gravity*pow(fraction*thickness,2)/2.;
572  pressure_water=1.0/2.0*rho_water*gravity* ( pow(bed,2.0)-pow(rho_ice/rho_water*fraction*thickness,2.0) );
573 
574  pressure=pressure_litho-pressure_air-pressure_melange-pressure_water;
575  }
576  else{
577  _error_("fill type " << fill << " not supported yet.");
578  }
579 
580  /*Ok, add contribution to first node, along the normal i==0: */
581  for(int j=0;j<2;j++){
582  pe->values[j]+=pressure*normal[j]*length;
583  }
584 
585  /*Add contribution to second node, along the opposite normal: i==1 */
586  for(int j=0;j<2;j++){
587  pe->values[2+j]+= -pressure*normal[j]*length;
588  }
589 
590  /*Clean up and return*/
591  return pe;
592 }

Field Documentation

◆ id

int Riftfront::id

Definition at line 19 of file Riftfront.h.

◆ type

int Riftfront::type

Definition at line 22 of file Riftfront.h.

◆ fill

int Riftfront::fill

Definition at line 23 of file Riftfront.h.

◆ friction

IssmDouble Riftfront::friction

Definition at line 24 of file Riftfront.h.

◆ fractionincrement

IssmDouble Riftfront::fractionincrement

Definition at line 25 of file Riftfront.h.

◆ shelf

bool Riftfront::shelf

Definition at line 26 of file Riftfront.h.

◆ hnodes

Hook* Riftfront::hnodes

Definition at line 29 of file Riftfront.h.

◆ hvertices

Hook* Riftfront::hvertices

Definition at line 30 of file Riftfront.h.

◆ helements

Hook* Riftfront::helements

Definition at line 31 of file Riftfront.h.

◆ nodes

Node** Riftfront::nodes

Definition at line 34 of file Riftfront.h.

◆ vertices

Vertex** Riftfront::vertices

Definition at line 35 of file Riftfront.h.

◆ elements

Element** Riftfront::elements

Definition at line 36 of file Riftfront.h.

◆ penalty_lock

int Riftfront::penalty_lock

Definition at line 39 of file Riftfront.h.

◆ active

bool Riftfront::active

Definition at line 40 of file Riftfront.h.

◆ frozen

bool Riftfront::frozen

Definition at line 41 of file Riftfront.h.

◆ counter

int Riftfront::counter

Definition at line 42 of file Riftfront.h.

◆ prestable

bool Riftfront::prestable

Definition at line 43 of file Riftfront.h.

◆ material_converged

bool Riftfront::material_converged

Definition at line 44 of file Riftfront.h.

◆ normal

IssmDouble Riftfront::normal[2]

Definition at line 45 of file Riftfront.h.

◆ length

IssmDouble Riftfront::length

Definition at line 46 of file Riftfront.h.

◆ fraction

IssmDouble Riftfront::fraction

Definition at line 47 of file Riftfront.h.

◆ state

int Riftfront::state

Definition at line 48 of file Riftfront.h.

◆ parameters

Parameters* Riftfront::parameters

Definition at line 50 of file Riftfront.h.

The documentation for this class was generated from the following files:
BaseEnum
@ BaseEnum
Definition: EnumDefinitions.h:495
Riftfront::normal
IssmDouble normal[2]
Definition: Riftfront.h:45
StressbalancePenaltyFactorEnum
@ StressbalancePenaltyFactorEnum
Definition: EnumDefinitions.h:409
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
Riftfront::material_converged
bool material_converged
Definition: Riftfront.h:44
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
StressbalanceAnalysisEnum
@ StressbalanceAnalysisEnum
Definition: EnumDefinitions.h:1285
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
Hook::deliverp
Object ** deliverp(void)
Definition: Hook.cpp:187
MARSHALLING_ENUM
#define MARSHALLING_ENUM(EN)
Definition: Marshalling.h:14
Hook::DeepEcho
void DeepEcho(void)
Definition: Hook.cpp:77
Riftfront::shelf
bool shelf
Definition: Riftfront.h:26
Riftfront::penalty_lock
int penalty_lock
Definition: Riftfront.h:39
SsetEnum
@ SsetEnum
Definition: EnumDefinitions.h:1282
MaterialsRhoIceEnum
@ MaterialsRhoIceEnum
Definition: EnumDefinitions.h:264
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
Riftfront::hvertices
Hook * hvertices
Definition: Riftfront.h:30
Riftfront::Riftfront
Riftfront()
Definition: Riftfront.cpp:22
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
Riftfront::length
IssmDouble length
Definition: Riftfront.h:46
Hook::reset
void reset(void)
Definition: Hook.cpp:211
Riftfront::elements
Element ** elements
Definition: Riftfront.h:36
NUMVERTICES
#define NUMVERTICES
Definition: Riftfront.cpp:19
Element
Definition: Element.h:41
Riftfront::parameters
Parameters * parameters
Definition: Riftfront.h:50
Riftfront::ObjectEnum
int ObjectEnum()
Definition: Riftfront.cpp:231
Riftfront::hnodes
Hook * hnodes
Definition: Riftfront.h:29
Parameters::DeepEcho
void DeepEcho()
Definition: Parameters.cpp:99
ConstantsGEnum
@ ConstantsGEnum
Definition: EnumDefinitions.h:102
AdjointHorizAnalysisEnum
@ AdjointHorizAnalysisEnum
Definition: EnumDefinitions.h:972
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
Hook
Definition: Hook.h:16
Hook::configure
void configure(DataSet *dataset)
Definition: Hook.cpp:145
GsetEnum
@ GsetEnum
Definition: EnumDefinitions.h:1093
IoModel::FindConstant
void FindConstant(bool *pvalue, const char *constant_name)
Definition: IoModel.cpp:2362
AirEnum
@ AirEnum
Definition: EnumDefinitions.h:469
Riftfront::state
int state
Definition: Riftfront.h:48
MARSHALLING
#define MARSHALLING(FIELD)
Definition: Marshalling.h:29
MelangeEnum
@ MelangeEnum
Definition: EnumDefinitions.h:1181
Riftfront::active
bool active
Definition: Riftfront.h:40
Riftfront::friction
IssmDouble friction
Definition: Riftfront.h:24
Riftfront::helements
Hook * helements
Definition: Riftfront.h:31
MaterialsRhoSeawaterEnum
@ MaterialsRhoSeawaterEnum
Definition: EnumDefinitions.h:265
RiftfrontEnum
@ RiftfrontEnum
Definition: EnumDefinitions.h:1240
Hook::copy
Object * copy(void)
Definition: Hook.cpp:61
Riftfront::fraction
IssmDouble fraction
Definition: Riftfront.h:47
Riftfront::fractionincrement
IssmDouble fractionincrement
Definition: Riftfront.h:25
MARSHALLING_BACKWARD
@ MARSHALLING_BACKWARD
Definition: Marshalling.h:10
SegmentRiftfrontEnum
@ SegmentRiftfrontEnum
Definition: EnumDefinitions.h:1271
Riftfront::PenaltyCreateKMatrixStressbalanceHoriz
ElementMatrix * PenaltyCreateKMatrixStressbalanceHoriz(IssmDouble kmax)
Definition: Riftfront.cpp:431
Tria::GetInputValue
void GetInputValue(IssmDouble *pvalue, Node *node, int enumtype)
Definition: Tria.cpp:2172
IoModel::Data
IssmDouble * Data(const char *data_name)
Definition: IoModel.cpp:437
Node
Definition: Node.h:23
Node::Lid
int Lid(void)
Definition: Node.cpp:618
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Tria
Definition: Tria.h:29
Riftfront::vertices
Vertex ** vertices
Definition: Riftfront.h:35
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
Riftfront::nodes
Node ** nodes
Definition: Riftfront.h:34
IceEnum
@ IceEnum
Definition: EnumDefinitions.h:626
Riftfront::PenaltyCreatePVectorStressbalanceHoriz
ElementVector * PenaltyCreatePVectorStressbalanceHoriz(IssmDouble kmax)
Definition: Riftfront.cpp:504
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
ThicknessEnum
@ ThicknessEnum
Definition: EnumDefinitions.h:840
ElementVector::AddToGlobal
void AddToGlobal(Vector< IssmDouble > *pf)
Definition: ElementVector.cpp:155
Riftfront::counter
int counter
Definition: Riftfront.h:42
Riftfront::frozen
bool frozen
Definition: Riftfront.h:41
AnalysisTypeEnum
@ AnalysisTypeEnum
Definition: EnumDefinitions.h:36
ElementVector
Definition: ElementVector.h:20
ElementMatrix::AddToGlobal
void AddToGlobal(Matrix< IssmDouble > *Kff, Matrix< IssmDouble > *Kfs)
Definition: ElementMatrix.cpp:271
Riftfront::fill
int fill
Definition: Riftfront.h:23
Riftfront::type
int type
Definition: Riftfront.h:22
Vertex
Definition: Vertex.h:19
IoRiftfillToEnum
int IoRiftfillToEnum(int enum_in)
Definition: IoCodeConversions.cpp:327
FsetEnum
@ FsetEnum
Definition: EnumDefinitions.h:1075
ElementMatrix
Definition: ElementMatrix.h:19
ClosedEnum
@ ClosedEnum
Definition: EnumDefinitions.h:1011
TriaEnum
@ TriaEnum
Definition: EnumDefinitions.h:1318
OpenEnum
@ OpenEnum
Definition: EnumDefinitions.h:1209
Riftfront::prestable
bool prestable
Definition: Riftfront.h:43
Riftfront
Definition: Riftfront.h:16
Hook::Marshall
void Marshall(char **pmarshalled_data, int *pmarshalled_data_size, int marshall_direction)
Definition: Hook.cpp:122
WaterEnum
@ WaterEnum
Definition: EnumDefinitions.h:1328
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
Riftfront::id
int id
Definition: Riftfront.h:19
Generated on Thu Jul 2 2020 08:09:30 for Ice Sheet System Model by   doxygen 1.8.19