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

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

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