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

Last change on this file since 3433 was 3433, checked in by Eric.Larour, 15 years ago

Debug
File size: 31.1 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::GetDofList {{{1*/
460
461void Riftfront::GetDofList(int* doflist,int* pnumberofdofspernode){
462
463 int i,j;
464 int doflist_per_node[MAXDOFSPERNODE];
465 int numberofdofspernode;
466
467 for(i=0;i<MAX_RIFTFRONT_GRIDS;i++){
468 nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
469 for(j=0;j<numberofdofspernode;j++){
470 doflist[i*numberofdofspernode+j]=doflist_per_node[j];
471 }
472 }
473
474 /*Assign output pointers:*/
475 *pnumberofdofspernode=numberofdofspernode;
476}
477/*}}}1*/
478/*FUNCTION Riftfront::GetId {{{1*/
479int Riftfront::GetId(void){ return id; }
480/*}}}1*/
481/*FUNCTION Riftfront::GetName {{{1*/
482char* Riftfront::GetName(void){
483 return "riftfront";
484}
485/*}}}1*/
486/*FUNCTION Riftfront::IsMaterialStable {{{1*/
487int Riftfront::IsMaterialStable(void* vinputs, int analysis_type){
488
489 int found=0;
490 ParameterInputs* inputs=NULL;
491 double converged=0;
492
493 inputs=(ParameterInputs*)vinputs;
494
495 found=inputs->Recover("converged",&converged);
496 if(!found)ISSMERROR(" could not find converged flag in inputs!");
497
498 if(converged){
499 /*ok, material non-linearity has converged. If that was already the case, we keep
500 * constraining the rift front. If it was not, and this is the first time the material
501 * has converged, we start constraining now!: */
502 this->material_converged=1;
503 }
504
505 return this->material_converged;
506}
507/*}}}1*/
508/*FUNCTION Riftfront::MaxPenetration {{{1*/
509int Riftfront::MaxPenetration(double* ppenetration, void* vinputs, int analysis_type){
510
511 const int numgrids=2;
512 int dofs[2]={0,1};
513 double vxvy_list[2][2]; //velocities for all grids
514 double max_penetration;
515 double penetration=0;
516 int found;
517
518 ParameterInputs* inputs=NULL;
519
520 inputs=(ParameterInputs*)vinputs;
521
522 //initialize:
523 penetration=-1;
524
525 found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
526 if(!found)ISSMERROR(" could not find velocity in inputs!");
527
528 /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
529 penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1];
530
531 /*Now, we return penetration only if we are active!: */
532 if(this->active==0)penetration=-1;
533
534 /*If we are zigzag locked, same thing: */
535 if(this->counter>this->penalty_lock)penetration=-1;
536
537 /*assign output pointer: */
538 *ppenetration=penetration;
539
540}
541/*}}}1*/
542/*FUNCTION Riftfront::MyRank {{{1*/
543int Riftfront::MyRank(void){
544 extern int my_rank;
545 return my_rank;
546}
547/*}}}1*/
548/*FUNCTION Riftfront::OutputProperties {{{1*/
549void Riftfront::OutputProperties(Vec riftproperties){
550
551 int row_id=0;
552 double value;
553
554 /*recover id of penalty: */
555 row_id=this->GetId()-1; //c indexing, ids were matlab indexed
556 value=(double)this->fraction;
557
558 /*Plug id and fraction into riftproperties matrix: */
559 VecSetValues(riftproperties,1,&row_id,&value,INSERT_VALUES);
560}
561/*}}}1*/
562/*FUNCTION Riftfront::PenaltyCreateKMatrix {{{1*/
563void Riftfront::PenaltyCreateKMatrix(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
564
565 int i,j;
566 const int numgrids=MAX_RIFTFRONT_GRIDS;
567 int dofs[1]={0};
568 double Ke_gg[4][4];
569 const int numdof=2*numgrids;
570 int doflist[numdof];
571 int numberofdofspernode;
572 double thickness;
573 ParameterInputs* inputs=NULL;
574
575 /*Some pointer intialization: */
576 inputs=(ParameterInputs*)vinputs;
577
578 /* Get node coordinates and dof list: */
579 GetDofList(&doflist[0],&numberofdofspernode);
580
581 /* Set Ke_gg to 0: */
582 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
583
584
585 if(this->active){
586
587 /*There is contact, we need to constrain the normal velocities (zero penetration), and the
588 *contact slip friction. */
589
590 #ifdef _ISSM_DEBUG_
591 printf("Dealing with grid pair (%i,%i)\n",nodes[0]->GetId(),nodes[1]->GetId());
592 #endif
593
594 /*Recover input parameters: */
595 inputs->Recover("thickness",&h[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
596 if (h[0]!=h[1])ISSMERROR(" different thicknesses not supported for rift fronts");
597 thickness=h[0];
598
599 #ifdef _ISSM_DEBUG_
600 printf("Thickness at grid (%i,%i): %lg\n",nodes[0]->GetId(),nodes[1]->GetID(),thickness);
601 #endif
602
603 /*From Peter Wriggers book (Computational Contact Mechanics, p191): */
604 //First line:
605 Ke_gg[0][0]+=pow(normal[0],2)*kmax*pow(10,penalty_offset);
606 Ke_gg[0][1]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset);
607 Ke_gg[0][2]+=-pow(normal[0],2)*kmax*pow(10,penalty_offset);
608 Ke_gg[0][3]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset);
609 //Second line:
610 Ke_gg[1][0]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset);
611 Ke_gg[1][1]+=pow(normal[1],2)*kmax*pow(10,penalty_offset);
612 Ke_gg[1][2]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset);
613 Ke_gg[1][3]+=-pow(normal[1],2)*kmax*pow(10,penalty_offset);
614 //Third line:
615 Ke_gg[2][0]+=-pow(normal[0],2)*kmax*pow(10,penalty_offset);
616 Ke_gg[2][1]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset);
617 Ke_gg[2][2]+=pow(normal[0],2)*kmax*pow(10,penalty_offset);
618 Ke_gg[2][3]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset);
619 //Fourth line:
620 Ke_gg[3][0]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset);
621 Ke_gg[3][1]+=-pow(normal[1],2)*kmax*pow(10,penalty_offset);
622 Ke_gg[3][2]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset);
623 Ke_gg[3][3]+=pow(normal[1],2)*kmax*pow(10,penalty_offset);
624
625 /*Now take care of the friction: of type sigma=frictiontangent_velocity2-tangent_velocity1)*/
626
627 //First line:
628 Ke_gg[0][0]+=pow(normal[1],2)*thickness*length*friction;
629 Ke_gg[0][1]+=-normal[0]*normal[1]*thickness*length*friction;
630 Ke_gg[0][2]+=-pow(normal[1],2)*thickness*length*friction;
631 Ke_gg[0][3]+=normal[0]*normal[1]*thickness*length*friction;
632 //Second line:
633 Ke_gg[1][0]+=-normal[0]*normal[1]*thickness*length*friction;
634 Ke_gg[1][1]+=pow(normal[0],2)*thickness*length*friction;
635 Ke_gg[1][2]+=normal[0]*normal[1]*thickness*length*friction;
636 Ke_gg[1][3]+=-pow(normal[0],2)*thickness*length*friction;
637 //Third line:
638 Ke_gg[2][0]+=-pow(normal[1],2)*thickness*length*friction;
639 Ke_gg[2][1]+=normal[0]*normal[1]*thickness*length*friction;
640 Ke_gg[2][2]+=pow(normal[1],2)*thickness*length*friction;
641 Ke_gg[2][3]+=-normal[0]*normal[1]*thickness*length*friction;
642 //Fourth line:
643 Ke_gg[3][0]+=normal[0]*normal[1]*thickness*length*friction;
644 Ke_gg[3][1]+=-pow(normal[0],2)*thickness*length*friction;
645 Ke_gg[3][2]+=-normal[0]*normal[1]*thickness*length*friction;
646 Ke_gg[3][3]+=pow(normal[0],2)*thickness*length*friction;
647
648 /*Add Ke_gg to global matrix Kgg: */
649 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
650 }
651 else{
652 /*the grids on both sides of the rift do not penetrate. PenaltyCreatePVector will
653 *take care of adding point loads to simulate pressure on the rift flanks. But as far as stiffness,
654 there is none (0 stiffness implies decoupling of the flank rifts, which is exactly what we want): */
655 }
656
657}
658/*}}}1*/
659/*FUNCTION Riftfront::PenaltyCreatePVector {{{1*/
660void Riftfront::PenaltyCreatePVector(Vec pg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
661
662 int i,j;
663 const int numgrids=MAX_RIFTFRONT_GRIDS;
664 int dofs[1]={0};
665 double pe_g[4];
666 const int numdof=2*numgrids;
667 int doflist[numdof];
668 int numberofdofspernode;
669 ParameterInputs* inputs=NULL;
670 double rho_ice;
671 double rho_water;
672 double gravity;
673 double thickness;
674 double bed;
675 double pressure;
676 double pressure_litho;
677 double pressure_air;
678 double pressure_melange;
679 double pressure_water;
680
681 /*Some pointer intialization: */
682 inputs=(ParameterInputs*)vinputs;
683
684 /* Get node coordinates and dof list: */
685 GetDofList(&doflist[0],&numberofdofspernode);
686
687 /* Set pe_g to 0: */
688 for(i=0;i<numdof;i++) pe_g[i]=0;
689
690 if(!this->active){
691 /*Ok, this rift is opening. We should put loads on both sides of the rift flanks. Because we are dealing with contact mechanics,
692 * and we want to avoid zigzagging of the loads, we want lump the loads onto grids, not onto surfaces between grids.:*/
693
694 #ifdef _ISSM_DEBUG_
695 _printf_("Grids (%i,%i) are free of constraints\n",nodes[0]->GetId(),nodes[1]->GetID());
696 #endif
697
698 /*Ok, to compute the pressure, we are going to need material properties, thickness and bed for the two grids. We assume those properties to
699 * be the same across the rift.: */
700
701 rho_ice=matpar->GetRhoIce();
702 rho_water=matpar->GetRhoWater();
703 gravity=matpar->GetG();
704
705 /*get thickness: */
706 inputs->Recover("thickness",&h[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
707 if (h[0]!=h[1])ISSMERROR(" different thicknesses not supported for rift fronts");
708 thickness=h[0];
709
710 inputs->Recover("bed",&b[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
711 if (b[0]!=b[1])ISSMERROR(" different beds not supported for rift fronts");
712 bed=b[0];
713
714 /*Ok, now compute the pressure (in norm) that is being applied to the flanks, depending on the type of fill: */
715 if(fill==WaterEnum()){
716 if(shelf){
717 /*We are on an ice shelf, hydrostatic equilibrium is used to determine the pressure for water fill: */
718 pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2 - rho_water*gravity*pow(bed,(double)2)/(double)2;
719 }
720 else{
721 //We are on an icesheet, we assume the water column fills the entire front: */
722 pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2 - rho_water*gravity*pow(thickness,(double)2)/(double)2;
723 }
724 }
725 else if(fill==AirEnum()){
726 pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2; //icefront on an ice sheet, pressure imbalance ice vs air.
727 }
728 else if(fill==IceEnum()){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
729 pressure=0;
730 }
731 else if(fill==MelangeEnum()){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
732
733 if(!shelf) ISSMERROR(exprintf("%s%i%s","fill type ",fill," not supported on ice sheets yet."));
734
735 pressure_litho=rho_ice*gravity*pow(thickness,(double)2)/(double)2;
736 pressure_air=0;
737 pressure_melange=rho_ice*gravity*pow(fraction*thickness,(double)2)/(double)2;
738 pressure_water=1.0/2.0*rho_water*gravity* ( pow(bed,2.0)-pow(rho_ice/rho_water*fraction*thickness,2.0) );
739
740 pressure=pressure_litho-pressure_air-pressure_melange-pressure_water;
741 }
742 else{
743 ISSMERROR(exprintf("%s%i%s","fill type ",fill," not supported yet."));
744 }
745
746 /*Ok, add contribution to first grid, along the normal i==0: */
747 for (j=0;j<2;j++){
748 pe_g[j]+=pressure*normal[j]*length;
749 }
750
751 /*Add contribution to second grid, along the opposite normal: i==1 */
752 for (j=0;j<2;j++){
753 pe_g[2+j]+= -pressure*normal[j]*length;
754 }
755 /*Add pe_g to global vector pg; */
756 VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
757
758 }
759 else{
760 /*The penalty is active. No loads implied here.*/
761 }
762}
763/*}}}1*/
764/*FUNCTION Riftfront::Penetration {{{1*/
765int Riftfront::Penetration(double* ppenetration, void* vinputs, int analysis_type){
766
767 const int numgrids=2;
768 int dofs[2]={0,1};
769 double vxvy_list[2][2]; //velocities for all grids
770 double max_penetration;
771 double penetration;
772 int found;
773
774 ParameterInputs* inputs=NULL;
775
776 inputs=(ParameterInputs*)vinputs;
777
778
779 found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
780 if(!found)ISSMERROR(" could not find velocity in inputs!");
781
782 /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
783 penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1];
784
785 /*Now, we return penetration only if we are active!: */
786 if(this->active==0)penetration=0;
787
788 /*assign output pointer: */
789 *ppenetration=penetration;
790
791}
792/*}}}1*/
793/*FUNCTION Riftfront::PotentialUnstableC
794 * nstraint {{{1*/
795int Riftfront::PotentialUnstableConstraint(int* punstable, void* vinputs, int analysis_type){
796
797
798 const int numgrids=2;
799 int dofs[2]={0,1};
800 double vxvy_list[2][2]; //velocities for all grids
801 double max_penetration;
802 double penetration;
803 int activate;
804 int unstable;
805 int found;
806
807 ParameterInputs* inputs=NULL;
808
809 inputs=(ParameterInputs*)vinputs;
810
811 found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
812 if(!found)ISSMERROR(" could not find velocity in inputs!");
813
814 /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
815 penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1];
816
817 /*Ok, we are looking for positive penetration in an active constraint: */
818 if(this->active){
819 if (penetration>=0){
820 unstable=1;
821 }
822 else{
823 unstable=0;
824 }
825 }
826 else{
827 unstable=0;
828 }
829
830 /*assign output pointer: */
831 *punstable=unstable;
832}
833/*}}}1*/
834/*FUNCTION Riftfront::PreConstrain {{{1*/
835int Riftfront::PreConstrain(int* punstable, void* vinputs, int analysis_type){
836
837 const int numgrids=2;
838 int dofs[2]={0,1};
839 double vxvy_list[2][2]; //velocities for all grids
840 double penetration;
841 int unstable;
842 ParameterInputs* inputs=NULL;
843 int found;
844
845 inputs=(ParameterInputs*)vinputs;
846
847 /*First recover velocity: */
848 found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
849 if(!found)ISSMERROR(" could not find velocity in inputs!");
850
851 /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
852 penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1];
853
854 /*Ok, we are preconstraining here. Ie, anything that penetrates is constrained until stability of the entire set
855 * of constraints is reached.: */
856 if(penetration<0){
857 if (!this->active){
858 /*This is the first time penetration happens: */
859 this->active=1;
860 unstable=1;
861 }
862 else{
863 /*This constraint was already active: */
864 this->active=1;
865 unstable=0;
866 }
867 }
868 else{
869 /*No penetration happening. : */
870 if (!this->active){
871 /*This penalty was not active, and no penetration happening. Do nonthing: */
872 this->active=0;
873 unstable=0;
874 }
875 else{
876 /*Ok, this penalty wants to get released. But not now, this is preconstraint, not constraint: */
877 this->active=1;
878 unstable=0;
879 }
880 }
881
882 /*assign output pointer: */
883 *punstable=unstable;
884}
885/*}}}1*/
886/*FUNCTION Riftfront::PreStable {{{1*/
887bool Riftfront::PreStable(){
888 return prestable;
889}
890/*}}}1*/
891/*FUNCTION Riftfront::SetPreStable {{{1*/
892void Riftfront::SetPreStable(){
893 prestable=1;
894}
895/*}}}1*/
896/*FUNCTION Riftfront::UpdateFromInputs {{{1*/
897void Riftfront::UpdateFromInputs(void* vinputs){
898
899 int dofs[1]={0};
900 ParameterInputs* inputs=NULL;
901
902 inputs=(ParameterInputs*)vinputs;
903
904 inputs->Recover("thickness",&h[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
905 inputs->Recover("bed",&b[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
906 inputs->Recover("surface",&s[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
907
908}
909/*}}}1*/
910/*FUNCTION Riftfront::FreezeConstraints{{{1*/
911void Riftfront::FreezeConstraints(void* vinputs, int analysis_type){
912
913 /*Just set frozen flag to 1: */
914 this->frozen=1;
915
916}
917/*}}}1*/
918/*FUNCTION Riftfront::IsFrozen{{{1*/
919bool Riftfront::IsFrozen(void){
920
921 /*Just set frozen flag to 1: */
922 if(this->frozen)return 1;
923 else return 0;
924}
925/*}}}1*/
Note: See TracBrowser for help on using the repository browser.