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

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

Model and FemModel are now classes in their own right.
This changes the cores quite a bit.

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