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

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

Added DofObject virtual abstract class. Redid all DistributeDofs,
so that it works for Node, Vertex, and anything that has dofs.
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::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 #ifdef _ISSM_DEBUG_
607 printf("Dealing with grid pair (%i,%i)\n",nodes[0]->GetId(),nodes[1]->GetId());
608 #endif
609
610 /*Recover input parameters: */
611 inputs->Recover("thickness",&h[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
612 if (h[0]!=h[1])ISSMERROR(" different thicknesses not supported for rift fronts");
613 thickness=h[0];
614
615 #ifdef _ISSM_DEBUG_
616 printf("Thickness at grid (%i,%i): %lg\n",nodes[0]->GetId(),nodes[1]->GetID(),thickness);
617 #endif
618
619 /*From Peter Wriggers book (Computational Contact Mechanics, p191): */
620 //First line:
621 Ke_gg[0][0]+=pow(normal[0],2)*kmax*pow(10,penalty_offset);
622 Ke_gg[0][1]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset);
623 Ke_gg[0][2]+=-pow(normal[0],2)*kmax*pow(10,penalty_offset);
624 Ke_gg[0][3]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset);
625 //Second line:
626 Ke_gg[1][0]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset);
627 Ke_gg[1][1]+=pow(normal[1],2)*kmax*pow(10,penalty_offset);
628 Ke_gg[1][2]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset);
629 Ke_gg[1][3]+=-pow(normal[1],2)*kmax*pow(10,penalty_offset);
630 //Third line:
631 Ke_gg[2][0]+=-pow(normal[0],2)*kmax*pow(10,penalty_offset);
632 Ke_gg[2][1]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset);
633 Ke_gg[2][2]+=pow(normal[0],2)*kmax*pow(10,penalty_offset);
634 Ke_gg[2][3]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset);
635 //Fourth line:
636 Ke_gg[3][0]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset);
637 Ke_gg[3][1]+=-pow(normal[1],2)*kmax*pow(10,penalty_offset);
638 Ke_gg[3][2]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset);
639 Ke_gg[3][3]+=pow(normal[1],2)*kmax*pow(10,penalty_offset);
640
641 /*Now take care of the friction: of type sigma=frictiontangent_velocity2-tangent_velocity1)*/
642
643 //First line:
644 Ke_gg[0][0]+=pow(normal[1],2)*thickness*length*friction;
645 Ke_gg[0][1]+=-normal[0]*normal[1]*thickness*length*friction;
646 Ke_gg[0][2]+=-pow(normal[1],2)*thickness*length*friction;
647 Ke_gg[0][3]+=normal[0]*normal[1]*thickness*length*friction;
648 //Second line:
649 Ke_gg[1][0]+=-normal[0]*normal[1]*thickness*length*friction;
650 Ke_gg[1][1]+=pow(normal[0],2)*thickness*length*friction;
651 Ke_gg[1][2]+=normal[0]*normal[1]*thickness*length*friction;
652 Ke_gg[1][3]+=-pow(normal[0],2)*thickness*length*friction;
653 //Third line:
654 Ke_gg[2][0]+=-pow(normal[1],2)*thickness*length*friction;
655 Ke_gg[2][1]+=normal[0]*normal[1]*thickness*length*friction;
656 Ke_gg[2][2]+=pow(normal[1],2)*thickness*length*friction;
657 Ke_gg[2][3]+=-normal[0]*normal[1]*thickness*length*friction;
658 //Fourth line:
659 Ke_gg[3][0]+=normal[0]*normal[1]*thickness*length*friction;
660 Ke_gg[3][1]+=-pow(normal[0],2)*thickness*length*friction;
661 Ke_gg[3][2]+=-normal[0]*normal[1]*thickness*length*friction;
662 Ke_gg[3][3]+=pow(normal[0],2)*thickness*length*friction;
663
664 /*Add Ke_gg to global matrix Kgg: */
665 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
666 }
667 else{
668 /*the grids on both sides of the rift do not penetrate. PenaltyCreatePVector will
669 *take care of adding point loads to simulate pressure on the rift flanks. But as far as stiffness,
670 there is none (0 stiffness implies decoupling of the flank rifts, which is exactly what we want): */
671 }
672
673}
674/*}}}1*/
675/*FUNCTION Riftfront::PenaltyCreatePVector {{{1*/
676void Riftfront::PenaltyCreatePVector(Vec pg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
677
678 int i,j;
679 const int numgrids=MAX_RIFTFRONT_GRIDS;
680 int dofs[1]={0};
681 double pe_g[4];
682 const int numdof=2*numgrids;
683 int doflist[numdof];
684 int numberofdofspernode;
685 ParameterInputs* inputs=NULL;
686 double rho_ice;
687 double rho_water;
688 double gravity;
689 double thickness;
690 double bed;
691 double pressure;
692 double pressure_litho;
693 double pressure_air;
694 double pressure_melange;
695 double pressure_water;
696
697 /*Some pointer intialization: */
698 inputs=(ParameterInputs*)vinputs;
699
700 /* Get node coordinates and dof list: */
701 GetDofList(&doflist[0],&numberofdofspernode);
702
703 /* Set pe_g to 0: */
704 for(i=0;i<numdof;i++) pe_g[i]=0;
705
706 if(!this->active){
707 /*Ok, this rift is opening. We should put loads on both sides of the rift flanks. Because we are dealing with contact mechanics,
708 * and we want to avoid zigzagging of the loads, we want lump the loads onto grids, not onto surfaces between grids.:*/
709
710 #ifdef _ISSM_DEBUG_
711 _printf_("Grids (%i,%i) are free of constraints\n",nodes[0]->GetId(),nodes[1]->GetID());
712 #endif
713
714 /*Ok, to compute the pressure, we are going to need material properties, thickness and bed for the two grids. We assume those properties to
715 * be the same across the rift.: */
716
717 rho_ice=matpar->GetRhoIce();
718 rho_water=matpar->GetRhoWater();
719 gravity=matpar->GetG();
720
721 /*get thickness: */
722 inputs->Recover("thickness",&h[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
723 if (h[0]!=h[1])ISSMERROR(" different thicknesses not supported for rift fronts");
724 thickness=h[0];
725
726 inputs->Recover("bed",&b[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
727 if (b[0]!=b[1])ISSMERROR(" different beds not supported for rift fronts");
728 bed=b[0];
729
730 /*Ok, now compute the pressure (in norm) that is being applied to the flanks, depending on the type of fill: */
731 if(fill==WaterEnum()){
732 if(shelf){
733 /*We are on an ice shelf, hydrostatic equilibrium is used to determine the pressure for water fill: */
734 pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2 - rho_water*gravity*pow(bed,(double)2)/(double)2;
735 }
736 else{
737 //We are on an icesheet, we assume the water column fills the entire front: */
738 pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2 - rho_water*gravity*pow(thickness,(double)2)/(double)2;
739 }
740 }
741 else if(fill==AirEnum()){
742 pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2; //icefront on an ice sheet, pressure imbalance ice vs air.
743 }
744 else if(fill==IceEnum()){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
745 pressure=0;
746 }
747 else if(fill==MelangeEnum()){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
748
749 if(!shelf) ISSMERROR(exprintf("%s%i%s","fill type ",fill," not supported on ice sheets yet."));
750
751 pressure_litho=rho_ice*gravity*pow(thickness,(double)2)/(double)2;
752 pressure_air=0;
753 pressure_melange=rho_ice*gravity*pow(fraction*thickness,(double)2)/(double)2;
754 pressure_water=1.0/2.0*rho_water*gravity* ( pow(bed,2.0)-pow(rho_ice/rho_water*fraction*thickness,2.0) );
755
756 pressure=pressure_litho-pressure_air-pressure_melange-pressure_water;
757 }
758 else{
759 ISSMERROR(exprintf("%s%i%s","fill type ",fill," not supported yet."));
760 }
761
762 /*Ok, add contribution to first grid, along the normal i==0: */
763 for (j=0;j<2;j++){
764 pe_g[j]+=pressure*normal[j]*length;
765 }
766
767 /*Add contribution to second grid, along the opposite normal: i==1 */
768 for (j=0;j<2;j++){
769 pe_g[2+j]+= -pressure*normal[j]*length;
770 }
771 /*Add pe_g to global vector pg; */
772 VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
773
774 }
775 else{
776 /*The penalty is active. No loads implied here.*/
777 }
778}
779/*}}}1*/
780/*FUNCTION Riftfront::Penetration {{{1*/
781int Riftfront::Penetration(double* ppenetration, void* vinputs, int analysis_type){
782
783 const int numgrids=2;
784 int dofs[2]={0,1};
785 double vxvy_list[2][2]; //velocities for all grids
786 double max_penetration;
787 double penetration;
788 int found;
789
790 ParameterInputs* inputs=NULL;
791
792 inputs=(ParameterInputs*)vinputs;
793
794
795 found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
796 if(!found)ISSMERROR(" could not find velocity in inputs!");
797
798 /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
799 penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1];
800
801 /*Now, we return penetration only if we are active!: */
802 if(this->active==0)penetration=0;
803
804 /*assign output pointer: */
805 *ppenetration=penetration;
806
807}
808/*}}}1*/
809/*FUNCTION Riftfront::PotentialUnstableConstraint {{{1*/
810int Riftfront::PotentialUnstableConstraint(int* punstable, void* vinputs, int analysis_type){
811
812
813 const int numgrids=2;
814 int dofs[2]={0,1};
815 double vxvy_list[2][2]; //velocities for all grids
816 double max_penetration;
817 double penetration;
818 int activate;
819 int unstable;
820 int found;
821
822 ParameterInputs* inputs=NULL;
823
824 inputs=(ParameterInputs*)vinputs;
825
826 found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
827 if(!found)ISSMERROR(" could not find velocity in inputs!");
828
829 /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
830 penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1];
831
832 /*Ok, we are looking for positive penetration in an active constraint: */
833 if(this->active){
834 if (penetration>=0){
835 unstable=1;
836 }
837 else{
838 unstable=0;
839 }
840 }
841 else{
842 unstable=0;
843 }
844
845 /*assign output pointer: */
846 *punstable=unstable;
847}
848/*}}}1*/
849/*FUNCTION Riftfront::PreConstrain {{{1*/
850int Riftfront::PreConstrain(int* punstable, void* vinputs, int analysis_type){
851
852 const int numgrids=2;
853 int dofs[2]={0,1};
854 double vxvy_list[2][2]; //velocities for all grids
855 double penetration;
856 int unstable;
857 ParameterInputs* inputs=NULL;
858 int found;
859
860 inputs=(ParameterInputs*)vinputs;
861
862 /*First recover velocity: */
863 found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
864 if(!found)ISSMERROR(" could not find velocity in inputs!");
865
866 /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
867 penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1];
868
869 /*Ok, we are preconstraining here. Ie, anything that penetrates is constrained until stability of the entire set
870 * of constraints is reached.: */
871 if(penetration<0){
872 if (!this->active){
873 /*This is the first time penetration happens: */
874 this->active=1;
875 unstable=1;
876 }
877 else{
878 /*This constraint was already active: */
879 this->active=1;
880 unstable=0;
881 }
882 }
883 else{
884 /*No penetration happening. : */
885 if (!this->active){
886 /*This penalty was not active, and no penetration happening. Do nonthing: */
887 this->active=0;
888 unstable=0;
889 }
890 else{
891 /*Ok, this penalty wants to get released. But not now, this is preconstraint, not constraint: */
892 this->active=1;
893 unstable=0;
894 }
895 }
896
897 /*assign output pointer: */
898 *punstable=unstable;
899}
900/*}}}1*/
901/*FUNCTION Riftfront::PreStable {{{1*/
902bool Riftfront::PreStable(){
903 return prestable;
904}
905/*}}}1*/
906/*FUNCTION Riftfront::SetPreStable {{{1*/
907void Riftfront::SetPreStable(){
908 prestable=1;
909}
910/*}}}1*/
911/*FUNCTION Riftfront::UpdateFromInputs {{{1*/
912void Riftfront::UpdateFromInputs(void* vinputs){
913
914 int dofs[1]={0};
915 ParameterInputs* inputs=NULL;
916
917 inputs=(ParameterInputs*)vinputs;
918
919 inputs->Recover("thickness",&h[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
920 inputs->Recover("bed",&b[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
921 inputs->Recover("surface",&s[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
922
923}
924/*}}}1*/
Note: See TracBrowser for help on using the repository browser.