source: issm/trunk/src/c/objects/Riftfront.cpp@ 3595

Last change on this file since 3595 was 3595, checked in by Mathieu Morlighem, 15 years ago

ISSMASSERT is now compiled only when ISSM_DEBUG is defined. Removed old ISSM_DEBUG outputs
File size: 30.7 KB
Line 
1/*!\file Riftfront.cpp
2 * \brief: implementation of the Riftfront object
3 */
4
5
6#ifdef HAVE_CONFIG_H
7 #include "config.h"
8#else
9#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
10#endif
11
12#include "stdio.h"
13#include <string.h>
14#include "../EnumDefinitions/EnumDefinitions.h"
15#include "../shared/shared.h"
16#include "../include/typedefs.h"
17#include "../include/macros.h"
18#include "./objects.h"
19
20
21/*Object constructors and destructor*/
22/*FUNCTION Riftfront::default constructor {{{1*/
23Riftfront::Riftfront(){
24 /*in case :*/
25 material_converged=0;
26 return;
27}
28/*}}}1*/
29/*FUNCTION Riftfront::constructor {{{1*/
30Riftfront::Riftfront(char riftfront_type[RIFTFRONTSTRING],int riftfront_id, int riftfront_node_ids[MAX_RIFTFRONT_GRIDS], int riftfront_mparid, double riftfront_h[MAX_RIFTFRONT_GRIDS],double riftfront_b[MAX_RIFTFRONT_GRIDS],double riftfront_s[MAX_RIFTFRONT_GRIDS],double riftfront_normal[2],double riftfront_length,int riftfront_fill,double riftfront_friction, double riftfront_fraction,double riftfront_fractionincrement, double riftfront_penalty_offset, int riftfront_penalty_lock, bool riftfront_active,bool riftfront_frozen, int riftfront_counter,bool riftfront_prestable,bool riftfront_shelf){
31
32 this->Init(riftfront_type,riftfront_id, riftfront_node_ids, riftfront_mparid, riftfront_h,riftfront_b,riftfront_s,riftfront_normal,riftfront_length,riftfront_fill,riftfront_friction, riftfront_fraction,riftfront_fractionincrement, riftfront_penalty_offset, riftfront_penalty_lock, riftfront_active,riftfront_frozen, riftfront_counter,riftfront_prestable,riftfront_shelf);
33
34}
35/*}}}1*/
36/*FUNCTION Riftfront::constructor from iomodel{{{1*/
37Riftfront::Riftfront(int i, IoModel* iomodel){
38
39 /*rifts: */
40 char riftfront_type[RIFTFRONTSTRING];
41 int riftfront_id;
42 int riftfront_node_ids[2];
43 int riftfront_mparid;
44 double riftfront_h[2];
45 double riftfront_b[2];
46 double riftfront_s[2];
47 double riftfront_normal[2];
48 double riftfront_length;
49 int riftfront_fill;
50 double riftfront_friction;
51 double riftfront_fraction;
52 double riftfront_fractionincrement;
53 bool riftfront_shelf;
54 double riftfront_penalty_offset;
55 int riftfront_penalty_lock;
56 bool riftfront_active;
57 bool riftfront_frozen;
58 int riftfront_counter;
59 bool riftfront_prestable;
60 int el1,el2;
61 int grid1,grid2;
62 double normal[2];
63 double length;
64 int fill;
65 double friction;
66 double fraction;
67 double fractionincrement;
68
69 /*Ok, retrieve all the data needed to add a penalty between the two grids: */
70 el1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2);
71 el2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+3);
72
73 grid1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+0);
74 grid2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+1);
75
76 normal[0]=*(iomodel->riftinfo+RIFTINFOSIZE*i+4);
77 normal[1]=*(iomodel->riftinfo+RIFTINFOSIZE*i+5);
78 length=*(iomodel->riftinfo+RIFTINFOSIZE*i+6);
79
80 fill = (int)*(iomodel->riftinfo+RIFTINFOSIZE*i+7);
81 friction=*(iomodel->riftinfo+RIFTINFOSIZE*i+8);
82 fraction=*(iomodel->riftinfo+RIFTINFOSIZE*i+9);
83 fractionincrement=*(iomodel->riftinfo+RIFTINFOSIZE*i+10);
84
85 strcpy(riftfront_type,"2d");
86 riftfront_id=i+1; //matlab indexing
87 riftfront_node_ids[0]=grid1;
88 riftfront_node_ids[1]=grid2;
89 riftfront_mparid=iomodel->numberofelements+1; //matlab indexing
90
91 riftfront_h[0]=iomodel->thickness[grid1-1];
92 riftfront_h[1]=iomodel->thickness[grid2-1];
93
94 riftfront_b[0]=iomodel->bed[grid1-1];
95 riftfront_b[1]=iomodel->bed[grid2-1];
96
97 riftfront_s[0]=iomodel->surface[grid1-1];
98 riftfront_s[1]=iomodel->surface[grid2-1];
99
100 riftfront_normal[0]=normal[0];
101 riftfront_normal[1]=normal[1];
102 riftfront_length=length;
103
104 riftfront_fill=fill;
105 riftfront_friction=friction;
106 riftfront_fraction=fraction;
107 riftfront_fractionincrement=fractionincrement;
108 riftfront_shelf=(bool)iomodel->gridoniceshelf[grid1-1];
109
110 riftfront_penalty_offset=iomodel->penalty_offset;
111 riftfront_penalty_lock=iomodel->penalty_lock;
112
113 riftfront_active=0;
114 riftfront_frozen=0;
115 riftfront_counter=0;
116 riftfront_prestable=0;
117
118 this->Init(riftfront_type,riftfront_id, riftfront_node_ids, riftfront_mparid, riftfront_h,riftfront_b,riftfront_s,riftfront_normal,riftfront_length,riftfront_fill,riftfront_friction, riftfront_fraction,riftfront_fractionincrement, riftfront_penalty_offset, riftfront_penalty_lock, riftfront_active,riftfront_frozen, riftfront_counter,riftfront_prestable,riftfront_shelf);
119
120}
121/*}}}1*/
122/*FUNCTION Riftfront::Init, used by constructor {{{1*/
123void Riftfront::Init(char riftfront_type[RIFTFRONTSTRING],int riftfront_id, int riftfront_node_ids[MAX_RIFTFRONT_GRIDS], int riftfront_mparid, double riftfront_h[MAX_RIFTFRONT_GRIDS],double riftfront_b[MAX_RIFTFRONT_GRIDS],double riftfront_s[MAX_RIFTFRONT_GRIDS],double riftfront_normal[2],double riftfront_length,int riftfront_fill,double riftfront_friction, double riftfront_fraction,double riftfront_fractionincrement, double riftfront_penalty_offset, int riftfront_penalty_lock, bool riftfront_active,bool riftfront_frozen, int riftfront_counter,bool riftfront_prestable,bool riftfront_shelf){
124
125 int i;
126
127 strcpy(type,riftfront_type);
128 id=riftfront_id;
129
130 for(i=0;i<MAX_RIFTFRONT_GRIDS;i++){
131 node_ids[i]=riftfront_node_ids[i];
132 node_offsets[i]=UNDEF;
133 nodes[i]=NULL;
134 }
135
136 mparid=riftfront_mparid;
137 matpar=NULL;
138 matpar_offset=UNDEF;
139
140 for(i=0;i<MAX_RIFTFRONT_GRIDS;i++){
141 h[i]=riftfront_h[i];
142 b[i]=riftfront_b[i];
143 s[i]=riftfront_s[i];
144 }
145
146 normal[0]=riftfront_normal[0];
147 normal[1]=riftfront_normal[1];
148 length=riftfront_length;
149 fill=riftfront_fill;
150 friction=riftfront_friction;
151 fraction=riftfront_fraction;
152 fractionincrement=riftfront_fractionincrement;
153 penalty_offset=riftfront_penalty_offset;
154 penalty_lock=riftfront_penalty_lock;
155 active=riftfront_active;
156 frozen=riftfront_frozen;
157 counter=riftfront_counter;
158 prestable=riftfront_prestable;
159 shelf=riftfront_shelf;
160
161 /*not in the constructor, but still needed: */
162 material_converged=0;
163
164 return;
165}
166/*}}}1*/
167/*FUNCTION Riftfront::destructor {{{1*/
168Riftfront::~Riftfront(){
169 return;
170}
171/*}}}1*/
172
173/*Object marshall*/
174/*FUNCTION Riftfront::Marshall {{{1*/
175void Riftfront::Marshall(char** pmarshalled_dataset){
176
177 char* marshalled_dataset=NULL;
178 int enum_type=0;
179
180 /*recover marshalled_dataset: */
181 marshalled_dataset=*pmarshalled_dataset;
182
183 /*get enum type of Riftfront: */
184 enum_type=RiftfrontEnum;
185
186 /*marshall enum: */
187 memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
188
189 /*marshall Riftfront data: */
190 memcpy(marshalled_dataset,&type,sizeof(type));marshalled_dataset+=sizeof(type);
191 memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
192 memcpy(marshalled_dataset,&node_ids,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
193 memcpy(marshalled_dataset,&node_offsets,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets);
194 memcpy(marshalled_dataset,&mparid,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
195 memcpy(marshalled_dataset,&matpar_offset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
196
197 memcpy(marshalled_dataset,&h,sizeof(h));marshalled_dataset+=sizeof(h);
198 memcpy(marshalled_dataset,&b,sizeof(b));marshalled_dataset+=sizeof(b);
199 memcpy(marshalled_dataset,&s,sizeof(s));marshalled_dataset+=sizeof(s);
200
201 memcpy(marshalled_dataset,&normal,sizeof(normal));marshalled_dataset+=sizeof(normal);
202 memcpy(marshalled_dataset,&length,sizeof(length));marshalled_dataset+=sizeof(length);
203 memcpy(marshalled_dataset,&fill,sizeof(fill));marshalled_dataset+=sizeof(fill);
204 memcpy(marshalled_dataset,&friction,sizeof(friction));marshalled_dataset+=sizeof(friction);
205 memcpy(marshalled_dataset,&fraction,sizeof(fraction));marshalled_dataset+=sizeof(fraction);
206 memcpy(marshalled_dataset,&fractionincrement,sizeof(fractionincrement));marshalled_dataset+=sizeof(fractionincrement);
207 memcpy(marshalled_dataset,&penalty_offset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
208 memcpy(marshalled_dataset,&penalty_lock,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock);
209 memcpy(marshalled_dataset,&active,sizeof(active));marshalled_dataset+=sizeof(active);
210 memcpy(marshalled_dataset,&frozen,sizeof(frozen));marshalled_dataset+=sizeof(frozen);
211 memcpy(marshalled_dataset,&counter,sizeof(counter));marshalled_dataset+=sizeof(counter);
212 memcpy(marshalled_dataset,&prestable,sizeof(prestable));marshalled_dataset+=sizeof(prestable);
213 memcpy(marshalled_dataset,&material_converged,sizeof(material_converged));marshalled_dataset+=sizeof(material_converged);
214 memcpy(marshalled_dataset,&shelf,sizeof(shelf));marshalled_dataset+=sizeof(shelf);
215
216 *pmarshalled_dataset=marshalled_dataset;
217 return;
218}
219/*}}}1*/
220/*FUNCTION Riftfront::MarshallSize {{{1*/
221int Riftfront::MarshallSize(){
222
223 return sizeof(type)+
224 sizeof(id)+
225 sizeof(node_ids)+
226 sizeof(node_offsets)+
227 sizeof(mparid)+
228 sizeof(matpar_offset)+
229 sizeof(h)+
230 sizeof(b)+
231 sizeof(s)+
232 sizeof(normal)+
233 sizeof(length)+
234 sizeof(fill)+
235 sizeof(friction)+
236 sizeof(fraction)+
237 sizeof(fractionincrement)+
238 sizeof(penalty_offset)+
239 sizeof(penalty_lock)+
240 sizeof(active)+
241 sizeof(frozen)+
242 sizeof(counter)+
243 sizeof(prestable)+
244 sizeof(material_converged)+
245 sizeof(shelf)+
246 sizeof(int); //sizeof(int) for enum type
247}
248/*}}}1*/
249/*FUNCTION Riftfront::Demarshall {{{1*/
250void Riftfront::Demarshall(char** pmarshalled_dataset){
251
252 int i;
253 char* marshalled_dataset=NULL;
254
255 /*recover marshalled_dataset: */
256 marshalled_dataset=*pmarshalled_dataset;
257
258 /*this time, no need to get enum type, the pointer directly points to the beginning of the
259 *object data (thanks to DataSet::Demarshall):*/
260
261 memcpy(&type,marshalled_dataset,sizeof(type));marshalled_dataset+=sizeof(type);
262 memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
263 memcpy(&node_ids,marshalled_dataset,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
264 memcpy(&node_offsets,marshalled_dataset,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets);
265 memcpy(&mparid,marshalled_dataset,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
266 memcpy(&matpar_offset,marshalled_dataset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
267
268 memcpy(&h,marshalled_dataset,sizeof(h));marshalled_dataset+=sizeof(h);
269 memcpy(&b,marshalled_dataset,sizeof(b));marshalled_dataset+=sizeof(b);
270 memcpy(&s,marshalled_dataset,sizeof(s));marshalled_dataset+=sizeof(s);
271 memcpy(&normal,marshalled_dataset,sizeof(normal));marshalled_dataset+=sizeof(normal);
272 memcpy(&length,marshalled_dataset,sizeof(length));marshalled_dataset+=sizeof(length);
273 memcpy(&fill,marshalled_dataset,sizeof(fill));marshalled_dataset+=sizeof(fill);
274 memcpy(&friction,marshalled_dataset,sizeof(friction));marshalled_dataset+=sizeof(friction);
275 memcpy(&fraction,marshalled_dataset,sizeof(fraction));marshalled_dataset+=sizeof(fraction);
276 memcpy(&fractionincrement,marshalled_dataset,sizeof(fractionincrement));marshalled_dataset+=sizeof(fractionincrement);
277 memcpy(&penalty_offset,marshalled_dataset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
278 memcpy(&penalty_lock,marshalled_dataset,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock);
279 memcpy(&active,marshalled_dataset,sizeof(active));marshalled_dataset+=sizeof(active);
280 memcpy(&frozen,marshalled_dataset,sizeof(frozen));marshalled_dataset+=sizeof(frozen);
281 memcpy(&counter,marshalled_dataset,sizeof(counter));marshalled_dataset+=sizeof(counter);
282 memcpy(&prestable,marshalled_dataset,sizeof(prestable));marshalled_dataset+=sizeof(prestable);
283 memcpy(&material_converged,marshalled_dataset,sizeof(material_converged));marshalled_dataset+=sizeof(material_converged);
284 memcpy(&shelf,marshalled_dataset,sizeof(shelf));marshalled_dataset+=sizeof(shelf);
285
286 for(i=0;i<MAX_RIFTFRONT_GRIDS;i++)nodes[i]=NULL;
287 matpar=NULL;
288
289 /*return: */
290 *pmarshalled_dataset=marshalled_dataset;
291 return;
292}
293/*}}}1*/
294
295/*Object functions*/
296/*FUNCTION Riftfront::Configure {{{1*/
297void Riftfront::Configure(void* pelementsin,void* pnodesin,void* pmaterialsin){
298
299 DataSet* nodesin=NULL;
300 DataSet* materialsin=NULL;
301
302 /*Recover pointers :*/
303 nodesin=(DataSet*)pnodesin;
304 materialsin=(DataSet*)pmaterialsin;
305
306 /*Link this load with its nodes: */
307 ResolvePointers((Object**)nodes,node_ids,node_offsets,MAX_RIFTFRONT_GRIDS,nodesin);
308
309 /*Same for materials: */
310 ResolvePointers((Object**)&matpar,&mparid,&matpar_offset,1,materialsin);
311
312}
313/*}}}1*/
314/*FUNCTION Riftfront::Constrain {{{1*/
315#define _ZIGZAGCOUNTER_
316
317int Riftfront::Constrain(int* punstable, void* vinputs, int analysis_type){
318
319 const int numgrids=2;
320 int dofs[2]={0,1};
321 double vxvy_list[2][2]; //velocities for all grids
322 double max_penetration;
323 double penetration;
324 int activate;
325 int found;
326 int unstable;
327
328 ParameterInputs* inputs=NULL;
329
330 inputs=(ParameterInputs*)vinputs;
331
332 /*Is this constraint frozen? In which case we don't touch: */
333 if (this->frozen){
334 *punstable=0;
335 return 1;
336 }
337
338 /*First recover parameter inputs: */
339 found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
340 if(!found)ISSMERROR(" could not find velocity in inputs!");
341
342
343 /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
344 penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1];
345
346 /*activation: */
347 if(penetration<0)activate=1;
348 else activate=0;
349
350 /*Here, we try to avoid zigzaging. When a penalty activates and deactivates for more than penalty_lock times,
351 * we increase the fraction of melange:*/
352 if(this->counter>this->penalty_lock){
353 /*reset counter: */
354 this->counter=0;
355 /*increase melange fraction: */
356 this->fraction+=this->fractionincrement;
357 if (this->fraction>1)this->fraction=(double)1.0;
358 //printf("riftfront %i fraction: %g\n",this->GetId(),this->fraction);
359 }
360
361 //Figure out stability of this penalty
362 if(this->active==activate){
363 unstable=0;
364 }
365 else{
366 unstable=1;
367 this->counter++;
368 }
369
370 //Set penalty flag
371 this->active=activate;
372
373 //if ((penetration>0) & (this->active==1))printf("Riftfront %i wants to be released\n",GetId());
374
375 /*assign output pointer: */
376 *punstable=unstable;
377}
378/*}}}1*/
379/*FUNCTION Riftfront::copy {{{1*/
380Object* Riftfront::copy() {
381 return new Riftfront(*this);
382}
383/*}}}1*/
384/*FUNCTION Riftfront::CreateKMatrix {{{1*/
385void Riftfront::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
386 /*do nothing: */
387}
388/*}}}1*/
389/*FUNCTION Riftfront::CreatePVector {{{1*/
390void Riftfront::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
391 /*do nothing: */
392}
393/*}}}1*/
394/*FUNCTION Riftfront::DeepEcho {{{1*/
395void Riftfront::DeepEcho(void){
396
397 int i;
398
399 printf("Riftfront:\n");
400 printf(" type: %s\n",type);
401 printf(" id: %i\n",id);
402
403 printf(" node_ids: ");
404 for(i=0;i<MAX_RIFTFRONT_GRIDS;i++)printf("%i ",node_ids[i]);
405 for(i=0;i<MAX_RIFTFRONT_GRIDS;i++){
406 if (nodes[i])nodes[i]->Echo();
407 }
408 printf("\n");
409
410 printf(" mparid: %i\n",mparid);
411 if(matpar)matpar->Echo();
412
413 printf("normal [%g,%g], length %g\n",normal[0],normal[1],normal[2]);
414 printf("fill: %i friction %g fraction %g fractionincrement %g \n",fill,friction,fraction,fractionincrement);
415 printf("penalty_offset %g\n",penalty_offset);
416 printf("penalty_lock %i\n",penalty_lock);
417 printf("active %i\n",active);
418 printf("frozen %i\n",frozen);
419 printf("counter %i\n",counter);
420 printf("prestable %i\n",prestable);
421 printf("material_converged %i\n",material_converged);
422 printf("shelf %i\n",shelf);
423
424}
425/*}}}1*/
426/*FUNCTION Riftfront::Echo {{{1*/
427void Riftfront::Echo(void){
428
429 int i;
430
431 printf("Riftfront:\n");
432 printf(" type: %s\n",type);
433 printf(" id: %i\n",id);
434 printf(" mparid: %i\n",mparid);
435
436 printf(" node_ids: ");
437 for(i=0;i<MAX_RIFTFRONT_GRIDS;i++)printf("%i ",node_ids[i]);
438 printf("\n");
439
440 printf("normal [%g,%g], length %g\n",normal[0],normal[1],normal[2]);
441 printf("fill: %i friction %g fraction %g fractionincrement %g \n",fill,friction,fraction,fractionincrement);
442 printf("penalty_offset %g\n",penalty_offset);
443 printf("penalty_lock %i\n",penalty_lock);
444 printf("active %i\n",active);
445 printf("frozen %i\n",frozen);
446 printf("counter %i\n",counter);
447 printf("prestable %i\n",prestable);
448 printf("material_converged %i\n",material_converged);
449 printf("shelf %i\n",shelf);
450}
451/*}}}1*/
452/*FUNCTION Riftfront::Enum {{{1*/
453int Riftfront::Enum(void){
454
455 return RiftfrontEnum;
456
457}
458/*}}}1*/
459/*FUNCTION Riftfront::FreezeConstraints{{{1*/
460void Riftfront::FreezeConstraints(void* vinputs, int analysis_type){
461
462 /*Just set frozen flag to 1: */
463 this->frozen=1;
464
465}
466/*}}}1*/
467/*FUNCTION Riftfront::GetDofList {{{1*/
468
469void Riftfront::GetDofList(int* doflist,int* pnumberofdofspernode){
470
471 int i,j;
472 int doflist_per_node[MAXDOFSPERNODE];
473 int numberofdofspernode;
474
475 for(i=0;i<MAX_RIFTFRONT_GRIDS;i++){
476 nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
477 for(j=0;j<numberofdofspernode;j++){
478 doflist[i*numberofdofspernode+j]=doflist_per_node[j];
479 }
480 }
481
482 /*Assign output pointers:*/
483 *pnumberofdofspernode=numberofdofspernode;
484}
485/*}}}1*/
486/*FUNCTION Riftfront::GetId {{{1*/
487int Riftfront::GetId(void){ return id; }
488/*}}}1*/
489/*FUNCTION Riftfront::GetName {{{1*/
490char* Riftfront::GetName(void){
491 return "riftfront";
492}
493/*}}}1*/
494/*FUNCTION Riftfront::IsFrozen{{{1*/
495bool Riftfront::IsFrozen(void){
496
497 /*Just set frozen flag to 1: */
498 if(this->frozen)return 1;
499 else return 0;
500}
501/*}}}1*/
502/*FUNCTION Riftfront::IsMaterialStable {{{1*/
503int Riftfront::IsMaterialStable(void* vinputs, int analysis_type){
504
505 int found=0;
506 ParameterInputs* inputs=NULL;
507 double converged=0;
508
509 inputs=(ParameterInputs*)vinputs;
510
511 found=inputs->Recover("converged",&converged);
512 if(!found)ISSMERROR(" could not find converged flag in inputs!");
513
514 if(converged){
515 /*ok, material non-linearity has converged. If that was already the case, we keep
516 * constraining the rift front. If it was not, and this is the first time the material
517 * has converged, we start constraining now!: */
518 this->material_converged=1;
519 }
520
521 return this->material_converged;
522}
523/*}}}1*/
524/*FUNCTION Riftfront::MaxPenetration {{{1*/
525int Riftfront::MaxPenetration(double* ppenetration, void* vinputs, int analysis_type){
526
527 const int numgrids=2;
528 int dofs[2]={0,1};
529 double vxvy_list[2][2]; //velocities for all grids
530 double max_penetration;
531 double penetration=0;
532 int found;
533
534 ParameterInputs* inputs=NULL;
535
536 inputs=(ParameterInputs*)vinputs;
537
538 //initialize:
539 penetration=-1;
540
541 found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
542 if(!found)ISSMERROR(" could not find velocity in inputs!");
543
544 /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
545 penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1];
546
547 /*Now, we return penetration only if we are active!: */
548 if(this->active==0)penetration=-1;
549
550 /*If we are zigzag locked, same thing: */
551 if(this->counter>this->penalty_lock)penetration=-1;
552
553 /*assign output pointer: */
554 *ppenetration=penetration;
555
556}
557/*}}}1*/
558/*FUNCTION Riftfront::MyRank {{{1*/
559int Riftfront::MyRank(void){
560 extern int my_rank;
561 return my_rank;
562}
563/*}}}1*/
564/*FUNCTION Riftfront::OutputProperties {{{1*/
565void Riftfront::OutputProperties(Vec riftproperties){
566
567 int row_id=0;
568 double value;
569
570 /*recover id of penalty: */
571 row_id=this->GetId()-1; //c indexing, ids were matlab indexed
572 value=(double)this->fraction;
573
574 /*Plug id and fraction into riftproperties matrix: */
575 VecSetValues(riftproperties,1,&row_id,&value,INSERT_VALUES);
576}
577/*}}}1*/
578/*FUNCTION Riftfront::PenaltyCreateKMatrix {{{1*/
579void Riftfront::PenaltyCreateKMatrix(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
580
581 int i,j;
582 const int numgrids=MAX_RIFTFRONT_GRIDS;
583 int dofs[1]={0};
584 double Ke_gg[4][4];
585 const int numdof=2*numgrids;
586 int doflist[numdof];
587 int numberofdofspernode;
588 double thickness;
589 ParameterInputs* inputs=NULL;
590
591 /*Some pointer intialization: */
592 inputs=(ParameterInputs*)vinputs;
593
594 /* Get node coordinates and dof list: */
595 GetDofList(&doflist[0],&numberofdofspernode);
596
597 /* Set Ke_gg to 0: */
598 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
599
600
601 if(this->active){
602
603 /*There is contact, we need to constrain the normal velocities (zero penetration), and the
604 *contact slip friction. */
605
606 /*Recover input parameters: */
607 inputs->Recover("thickness",&h[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
608 if (h[0]!=h[1])ISSMERROR(" different thicknesses not supported for rift fronts");
609 thickness=h[0];
610
611 /*From Peter Wriggers book (Computational Contact Mechanics, p191): */
612 //First line:
613 Ke_gg[0][0]+=pow(normal[0],2)*kmax*pow(10,penalty_offset);
614 Ke_gg[0][1]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset);
615 Ke_gg[0][2]+=-pow(normal[0],2)*kmax*pow(10,penalty_offset);
616 Ke_gg[0][3]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset);
617 //Second line:
618 Ke_gg[1][0]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset);
619 Ke_gg[1][1]+=pow(normal[1],2)*kmax*pow(10,penalty_offset);
620 Ke_gg[1][2]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset);
621 Ke_gg[1][3]+=-pow(normal[1],2)*kmax*pow(10,penalty_offset);
622 //Third line:
623 Ke_gg[2][0]+=-pow(normal[0],2)*kmax*pow(10,penalty_offset);
624 Ke_gg[2][1]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset);
625 Ke_gg[2][2]+=pow(normal[0],2)*kmax*pow(10,penalty_offset);
626 Ke_gg[2][3]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset);
627 //Fourth line:
628 Ke_gg[3][0]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset);
629 Ke_gg[3][1]+=-pow(normal[1],2)*kmax*pow(10,penalty_offset);
630 Ke_gg[3][2]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset);
631 Ke_gg[3][3]+=pow(normal[1],2)*kmax*pow(10,penalty_offset);
632
633 /*Now take care of the friction: of type sigma=frictiontangent_velocity2-tangent_velocity1)*/
634
635 //First line:
636 Ke_gg[0][0]+=pow(normal[1],2)*thickness*length*friction;
637 Ke_gg[0][1]+=-normal[0]*normal[1]*thickness*length*friction;
638 Ke_gg[0][2]+=-pow(normal[1],2)*thickness*length*friction;
639 Ke_gg[0][3]+=normal[0]*normal[1]*thickness*length*friction;
640 //Second line:
641 Ke_gg[1][0]+=-normal[0]*normal[1]*thickness*length*friction;
642 Ke_gg[1][1]+=pow(normal[0],2)*thickness*length*friction;
643 Ke_gg[1][2]+=normal[0]*normal[1]*thickness*length*friction;
644 Ke_gg[1][3]+=-pow(normal[0],2)*thickness*length*friction;
645 //Third line:
646 Ke_gg[2][0]+=-pow(normal[1],2)*thickness*length*friction;
647 Ke_gg[2][1]+=normal[0]*normal[1]*thickness*length*friction;
648 Ke_gg[2][2]+=pow(normal[1],2)*thickness*length*friction;
649 Ke_gg[2][3]+=-normal[0]*normal[1]*thickness*length*friction;
650 //Fourth line:
651 Ke_gg[3][0]+=normal[0]*normal[1]*thickness*length*friction;
652 Ke_gg[3][1]+=-pow(normal[0],2)*thickness*length*friction;
653 Ke_gg[3][2]+=-normal[0]*normal[1]*thickness*length*friction;
654 Ke_gg[3][3]+=pow(normal[0],2)*thickness*length*friction;
655
656 /*Add Ke_gg to global matrix Kgg: */
657 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
658 }
659 else{
660 /*the grids on both sides of the rift do not penetrate. PenaltyCreatePVector will
661 *take care of adding point loads to simulate pressure on the rift flanks. But as far as stiffness,
662 there is none (0 stiffness implies decoupling of the flank rifts, which is exactly what we want): */
663 }
664
665}
666/*}}}1*/
667/*FUNCTION Riftfront::PenaltyCreatePVector {{{1*/
668void Riftfront::PenaltyCreatePVector(Vec pg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
669
670 int i,j;
671 const int numgrids=MAX_RIFTFRONT_GRIDS;
672 int dofs[1]={0};
673 double pe_g[4];
674 const int numdof=2*numgrids;
675 int doflist[numdof];
676 int numberofdofspernode;
677 ParameterInputs* inputs=NULL;
678 double rho_ice;
679 double rho_water;
680 double gravity;
681 double thickness;
682 double bed;
683 double pressure;
684 double pressure_litho;
685 double pressure_air;
686 double pressure_melange;
687 double pressure_water;
688
689 /*Some pointer intialization: */
690 inputs=(ParameterInputs*)vinputs;
691
692 /* Get node coordinates and dof list: */
693 GetDofList(&doflist[0],&numberofdofspernode);
694
695 /* Set pe_g to 0: */
696 for(i=0;i<numdof;i++) pe_g[i]=0;
697
698 if(!this->active){
699 /*Ok, this rift is opening. We should put loads on both sides of the rift flanks. Because we are dealing with contact mechanics,
700 * and we want to avoid zigzagging of the loads, we want lump the loads onto grids, not onto surfaces between grids.:*/
701
702 /*Ok, to compute the pressure, we are going to need material properties, thickness and bed for the two grids. We assume those properties to
703 * be the same across the rift.: */
704
705 rho_ice=matpar->GetRhoIce();
706 rho_water=matpar->GetRhoWater();
707 gravity=matpar->GetG();
708
709 /*get thickness: */
710 inputs->Recover("thickness",&h[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
711 if (h[0]!=h[1])ISSMERROR(" different thicknesses not supported for rift fronts");
712 thickness=h[0];
713
714 inputs->Recover("bed",&b[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
715 if (b[0]!=b[1])ISSMERROR(" different beds not supported for rift fronts");
716 bed=b[0];
717
718 /*Ok, now compute the pressure (in norm) that is being applied to the flanks, depending on the type of fill: */
719 if(fill==WaterEnum){
720 if(shelf){
721 /*We are on an ice shelf, hydrostatic equilibrium is used to determine the pressure for water fill: */
722 pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2 - rho_water*gravity*pow(bed,(double)2)/(double)2;
723 }
724 else{
725 //We are on an icesheet, we assume the water column fills the entire front: */
726 pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2 - rho_water*gravity*pow(thickness,(double)2)/(double)2;
727 }
728 }
729 else if(fill==AirEnum){
730 pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2; //icefront on an ice sheet, pressure imbalance ice vs air.
731 }
732 else if(fill==IceEnum){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
733 pressure=0;
734 }
735 else if(fill==MelangeEnum){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
736
737 if(!shelf) ISSMERROR("%s%i%s","fill type ",fill," not supported on ice sheets yet.");
738
739 pressure_litho=rho_ice*gravity*pow(thickness,(double)2)/(double)2;
740 pressure_air=0;
741 pressure_melange=rho_ice*gravity*pow(fraction*thickness,(double)2)/(double)2;
742 pressure_water=1.0/2.0*rho_water*gravity* ( pow(bed,2.0)-pow(rho_ice/rho_water*fraction*thickness,2.0) );
743
744 pressure=pressure_litho-pressure_air-pressure_melange-pressure_water;
745 }
746 else{
747 ISSMERROR("%s%i%s","fill type ",fill," not supported yet.");
748 }
749
750 /*Ok, add contribution to first grid, along the normal i==0: */
751 for (j=0;j<2;j++){
752 pe_g[j]+=pressure*normal[j]*length;
753 }
754
755 /*Add contribution to second grid, along the opposite normal: i==1 */
756 for (j=0;j<2;j++){
757 pe_g[2+j]+= -pressure*normal[j]*length;
758 }
759 /*Add pe_g to global vector pg; */
760 VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
761
762 }
763 else{
764 /*The penalty is active. No loads implied here.*/
765 }
766}
767/*}}}1*/
768/*FUNCTION Riftfront::Penetration {{{1*/
769int Riftfront::Penetration(double* ppenetration, void* vinputs, int analysis_type){
770
771 const int numgrids=2;
772 int dofs[2]={0,1};
773 double vxvy_list[2][2]; //velocities for all grids
774 double max_penetration;
775 double penetration;
776 int found;
777
778 ParameterInputs* inputs=NULL;
779
780 inputs=(ParameterInputs*)vinputs;
781
782
783 found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
784 if(!found)ISSMERROR(" could not find velocity in inputs!");
785
786 /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
787 penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1];
788
789 /*Now, we return penetration only if we are active!: */
790 if(this->active==0)penetration=0;
791
792 /*assign output pointer: */
793 *ppenetration=penetration;
794
795}
796/*}}}1*/
797/*FUNCTION Riftfront::PotentialUnstableConstraint {{{1*/
798int Riftfront::PotentialUnstableConstraint(int* punstable, void* vinputs, int analysis_type){
799
800
801 const int numgrids=2;
802 int dofs[2]={0,1};
803 double vxvy_list[2][2]; //velocities for all grids
804 double max_penetration;
805 double penetration;
806 int activate;
807 int unstable;
808 int found;
809
810 ParameterInputs* inputs=NULL;
811
812 inputs=(ParameterInputs*)vinputs;
813
814 found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
815 if(!found)ISSMERROR(" could not find velocity in inputs!");
816
817 /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
818 penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1];
819
820 /*Ok, we are looking for positive penetration in an active constraint: */
821 if(this->active){
822 if (penetration>=0){
823 unstable=1;
824 }
825 else{
826 unstable=0;
827 }
828 }
829 else{
830 unstable=0;
831 }
832
833 /*assign output pointer: */
834 *punstable=unstable;
835}
836/*}}}1*/
837/*FUNCTION Riftfront::PreConstrain {{{1*/
838int Riftfront::PreConstrain(int* punstable, void* vinputs, int analysis_type){
839
840 const int numgrids=2;
841 int dofs[2]={0,1};
842 double vxvy_list[2][2]; //velocities for all grids
843 double penetration;
844 int unstable;
845 ParameterInputs* inputs=NULL;
846 int found;
847
848 inputs=(ParameterInputs*)vinputs;
849
850 /*First recover velocity: */
851 found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
852 if(!found)ISSMERROR(" could not find velocity in inputs!");
853
854 /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
855 penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1];
856
857 /*Ok, we are preconstraining here. Ie, anything that penetrates is constrained until stability of the entire set
858 * of constraints is reached.: */
859 if(penetration<0){
860 if (!this->active){
861 /*This is the first time penetration happens: */
862 this->active=1;
863 unstable=1;
864 }
865 else{
866 /*This constraint was already active: */
867 this->active=1;
868 unstable=0;
869 }
870 }
871 else{
872 /*No penetration happening. : */
873 if (!this->active){
874 /*This penalty was not active, and no penetration happening. Do nonthing: */
875 this->active=0;
876 unstable=0;
877 }
878 else{
879 /*Ok, this penalty wants to get released. But not now, this is preconstraint, not constraint: */
880 this->active=1;
881 unstable=0;
882 }
883 }
884
885 /*assign output pointer: */
886 *punstable=unstable;
887}
888/*}}}1*/
889/*FUNCTION Riftfront::PreStable {{{1*/
890bool Riftfront::PreStable(){
891 return prestable;
892}
893/*}}}1*/
894/*FUNCTION Riftfront::SetPreStable {{{1*/
895void Riftfront::SetPreStable(){
896 prestable=1;
897}
898/*}}}1*/
899/*FUNCTION Riftfront::UpdateFromInputs {{{1*/
900void Riftfront::UpdateFromInputs(void* vinputs){
901
902 int dofs[1]={0};
903 ParameterInputs* inputs=NULL;
904
905 inputs=(ParameterInputs*)vinputs;
906
907 inputs->Recover("thickness",&h[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
908 inputs->Recover("bed",&b[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
909 inputs->Recover("surface",&s[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
910
911}
912/*}}}1*/
Note: See TracBrowser for help on using the repository browser.