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

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

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