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

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

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