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

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

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