source: issm/trunk/src/c/objects/Pengrid.cpp@ 2342

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

Forgot to marshall and demarshall zigzag_counter and stabilize_constraints -> holy schmuck!
File size: 17.9 KB
Line 
1/*!\file Pengrid.c
2 * \brief: implementation of the Pengrid 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 "./Pengrid.h"
14#include <string.h>
15#include "../EnumDefinitions/EnumDefinitions.h"
16#include "../shared/shared.h"
17#include "../include/typedefs.h"
18
19
20Pengrid::Pengrid(){
21 return;
22}
23
24Pengrid::Pengrid(int pengrid_id, int pengrid_node_id,int pengrid_mparid, int pengrid_dof, int pengrid_active, double pengrid_penalty_offset,int pengrid_thermal_steadystate,int pengrid_stabilize_constraints){
25
26 id=pengrid_id;
27 node_id=pengrid_node_id;
28 mparid=pengrid_mparid;
29 dof=pengrid_dof;
30 active=pengrid_active;
31 penalty_offset =pengrid_penalty_offset;
32 thermal_steadystate=pengrid_thermal_steadystate;
33 stabilize_constraints=pengrid_stabilize_constraints;
34
35 node_offset=UNDEF;
36 node=NULL;
37 matpar=NULL;
38 matpar_offset=UNDEF;
39
40 zigzag_counter=0;
41
42 return;
43}
44
45Pengrid::~Pengrid(){
46 return;
47}
48
49void Pengrid::Echo(void){
50
51 printf("Pengrid:\n");
52 printf(" id: %i\n",id);
53 printf(" mparid: %i\n",mparid);
54 printf(" dof: %i\n",dof);
55 printf(" active: %i\n",active);
56 printf(" penalty_offset: %g\n",penalty_offset);
57 printf(" thermal_steadystate: %i\n",thermal_steadystate);
58 printf(" node_id: [%i]\n",node_id);
59 printf(" node_offset: [%i]\n",node_offset);
60 printf(" matpar_offset=%i\n",matpar_offset);
61
62 return;
63}
64void Pengrid::DeepEcho(void){
65
66 printf("Pengrid:\n");
67 printf(" id: %i\n",id);
68 printf(" mparid: %i\n",mparid);
69 printf(" dof: %i\n",dof);
70 printf(" active: %i\n",active);
71 printf(" penalty_offset: %g\n",penalty_offset);
72 printf(" thermal_steadystate: %i\n",thermal_steadystate);
73 printf(" node_id: [%i]\n",node_id);
74 printf(" node_offset: [%i]\n",node_offset);
75 printf(" matpar_offset=%i\n",matpar_offset);
76
77 if(node)node->Echo();
78 if(matpar)matpar->Echo();
79 return;
80}
81void Pengrid::Marshall(char** pmarshalled_dataset){
82
83 char* marshalled_dataset=NULL;
84 int enum_type=0;
85
86 /*recover marshalled_dataset: */
87 marshalled_dataset=*pmarshalled_dataset;
88
89 /*get enum type of Pengrid: */
90 enum_type=PengridEnum();
91
92 /*marshall enum: */
93 memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
94
95 /*marshall Pengrid data: */
96 memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
97 memcpy(marshalled_dataset,&mparid,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
98 memcpy(marshalled_dataset,&dof,sizeof(dof));marshalled_dataset+=sizeof(dof);
99 memcpy(marshalled_dataset,&active,sizeof(active));marshalled_dataset+=sizeof(active);
100 memcpy(marshalled_dataset,&penalty_offset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
101 memcpy(marshalled_dataset,&thermal_steadystate,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
102 memcpy(marshalled_dataset,&node_id,sizeof(node_id));marshalled_dataset+=sizeof(node_id);
103 memcpy(marshalled_dataset,&node_offset,sizeof(node_offset));marshalled_dataset+=sizeof(node_offset);
104 memcpy(marshalled_dataset,&matpar,sizeof(matpar));marshalled_dataset+=sizeof(matpar);
105 memcpy(marshalled_dataset,&matpar_offset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
106 memcpy(marshalled_dataset,&stabilize_constraints,sizeof(stabilize_constraints));marshalled_dataset+=sizeof(stabilize_constraints);
107 memcpy(marshalled_dataset,&zigzag_counter,sizeof(zigzag_counter));marshalled_dataset+=sizeof(zigzag_counter);
108
109 *pmarshalled_dataset=marshalled_dataset;
110 return;
111}
112
113int Pengrid::MarshallSize(){
114
115 return sizeof(id)+
116 sizeof(mparid)+
117 sizeof(dof)+
118 sizeof(active)+
119 sizeof(penalty_offset)+
120 sizeof(thermal_steadystate)+
121 sizeof(node_id)+
122 sizeof(node_offset)+
123 sizeof(matpar)+
124 sizeof(matpar_offset)+
125 sizeof(stabilize_constraints)+
126 sizeof(zigzag_counter)+
127 sizeof(int); //sizeof(int) for enum type
128}
129
130char* Pengrid::GetName(void){
131 return "pengrid";
132}
133
134
135void Pengrid::Demarshall(char** pmarshalled_dataset){
136
137 char* marshalled_dataset=NULL;
138
139 /*recover marshalled_dataset: */
140 marshalled_dataset=*pmarshalled_dataset;
141
142 /*this time, no need to get enum type, the pointer directly points to the beginning of the
143 *object data (thanks to DataSet::Demarshall):*/
144
145 memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
146 memcpy(&mparid,marshalled_dataset,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
147 memcpy(&dof,marshalled_dataset,sizeof(dof));marshalled_dataset+=sizeof(dof);
148 memcpy(&active,marshalled_dataset,sizeof(active));marshalled_dataset+=sizeof(active);
149 memcpy(&penalty_offset,marshalled_dataset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
150 memcpy(&thermal_steadystate,marshalled_dataset,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
151 memcpy(&node_id,marshalled_dataset,sizeof(node_id));marshalled_dataset+=sizeof(node_id);
152 memcpy(&node_offset,marshalled_dataset,sizeof(node_offset));marshalled_dataset+=sizeof(node_offset);
153 memcpy(&matpar,marshalled_dataset,sizeof(matpar));marshalled_dataset+=sizeof(matpar);
154 memcpy(&matpar_offset,marshalled_dataset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
155 memcpy(&stabilize_constraints,marshalled_dataset,sizeof(stabilize_constraints));marshalled_dataset+=sizeof(stabilize_constraints);
156 memcpy(&zigzag_counter,marshalled_dataset,sizeof(zigzag_counter));marshalled_dataset+=sizeof(zigzag_counter);
157
158 node=NULL;
159 matpar=NULL;
160
161 /*return: */
162 *pmarshalled_dataset=marshalled_dataset;
163 return;
164}
165int Pengrid::Enum(void){
166
167 return PengridEnum();
168}
169
170int Pengrid::GetId(void){ return id; }
171
172int Pengrid::MyRank(void){
173 extern int my_rank;
174 return my_rank;
175}
176void Pengrid::DistributeNumDofs(int* numdofpernode,int analysis_type,int sub_analysis_type){return;}
177
178#undef __FUNCT__
179#define __FUNCT__ "Pengrid::Configure"
180
181void Pengrid::Configure(void* pelementsin,void* pnodesin,void* pmaterialsin){
182
183 DataSet* nodesin=NULL;
184 DataSet* materialsin=NULL;
185
186 /*Recover pointers :*/
187 nodesin=(DataSet*)pnodesin;
188 materialsin=(DataSet*)pmaterialsin;
189
190 /*Link this load with its nodes: */
191 ResolvePointers((Object**)&node,&node_id,&node_offset,1,nodesin);
192 ResolvePointers((Object**)&matpar,&mparid,&matpar_offset,1,materialsin);
193}
194
195
196#undef __FUNCT__
197#define __FUNCT__ "Pengrid::CreateKMatrix"
198
199void Pengrid::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
200
201 /*No loads applied, do nothing: */
202 return;
203
204}
205
206#undef __FUNCT__
207#define __FUNCT__ "Pengrid::CreatePVector"
208void Pengrid::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
209
210 /*No loads applied, do nothing: */
211 return;
212
213}
214#undef __FUNCT__
215#define __FUNCT__ "Pengrid::UpdateFromInputs"
216void Pengrid::UpdateFromInputs(void* inputs){
217
218}
219
220#undef __FUNCT__
221#define __FUNCT__ "Pengrid::PenaltyCreateKMatrix"
222void Pengrid::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
223
224 if ((analysis_type==DiagnosticAnalysisEnum()) && ((sub_analysis_type==StokesAnalysisEnum()))){
225
226 PenaltyCreateKMatrixDiagnosticStokes( Kgg,inputs,kmax,analysis_type,sub_analysis_type);
227 }
228 else if (analysis_type==ThermalAnalysisEnum()){
229
230 PenaltyCreateKMatrixThermal( Kgg,inputs,kmax,analysis_type,sub_analysis_type);
231
232 }
233 else if (analysis_type==MeltingAnalysisEnum()){
234
235 PenaltyCreateKMatrixMelting( Kgg,inputs,kmax,analysis_type,sub_analysis_type);
236
237 }
238 else{
239 throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
240 }
241
242}
243
244#undef __FUNCT__
245#define __FUNCT__ "Pengrid::PenaltyCreateKMatrixDiagnosticStokes"
246void Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
247
248 const int numgrids=1;
249 const int NDOF4=4;
250 const int numdof=numgrids*NDOF4;
251 int doflist[numdof];
252 int numberofdofspernode;
253
254 int dofs1[1]={0};
255 int dofs2[1]={1};
256 double slope[2];
257 int found=0;
258 double Ke[4][4]={0.0};
259
260 ParameterInputs* inputs=NULL;
261
262 /*recover pointers: */
263 inputs=(ParameterInputs*)vinputs;
264
265 /*Get dof list: */
266 GetDofList(&doflist[0],&numberofdofspernode);
267
268 /*recover slope: */
269 found=inputs->Recover("bedslopex",&slope[0],1,dofs1,numgrids,(void**)&node);
270 if(!found)throw ErrorException(__FUNCT__," bedslopex needed in inputs!");
271 found=inputs->Recover("bedslopey",&slope[1],1,dofs2,numgrids,(void**)&node);
272 if(!found)throw ErrorException(__FUNCT__," bedslopey needed in inputs!");
273
274 //Create elementary matrix: add penalty to contrain wb (wb=ub*db/dx+vb*db/dy)
275 Ke[2][0]=-slope[0]*kmax*pow((double)10.0,penalty_offset);
276 Ke[2][1]=-slope[1]*kmax*pow((double)10.0,penalty_offset);
277 Ke[2][2]=kmax*pow((double)10,penalty_offset);
278
279 /*Add Ke to global matrix Kgg: */
280 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
281}
282
283#undef __FUNCT__
284#define __FUNCT__ "Pengrid::PenaltyCreateKMatrixThermal"
285void Pengrid::PenaltyCreateKMatrixThermal(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
286
287 int found=0;
288
289 const int numgrids=1;
290 const int NDOF1=1;
291 const int numdof=numgrids*NDOF1;
292 double Ke[numdof][numdof];
293 int doflist[numdof];
294 int numberofdofspernode;
295
296 ParameterInputs* inputs=NULL;
297
298 /*recover pointers: */
299 inputs=(ParameterInputs*)vinputs;
300
301
302 if(!active)return;
303
304 /*Get dof list: */
305 GetDofList(&doflist[0],&numberofdofspernode);
306
307 Ke[0][0]=kmax*pow((double)10,penalty_offset);
308
309 /*Add Ke to global matrix Kgg: */
310 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
311}
312
313#undef __FUNCT__
314#define __FUNCT__ "Pengrid::PenaltyCreateKMatrixMelting"
315void Pengrid::PenaltyCreateKMatrixMelting(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
316
317
318 int found=0;
319 const int numgrids=1;
320 const int NDOF1=1;
321 const int numdof=numgrids*NDOF1;
322 double Ke[numdof][numdof]={0.0};
323 int dofs1[1]={0};
324 int doflist[numdof];
325 int numberofdofspernode;
326 double meltingpoint;
327
328 double pressure;
329 double temperature;
330 double beta,t_pmp;
331
332 ParameterInputs* inputs=NULL;
333
334 /*check that pengrid is not a clone (penalty to be added only once)*/
335 if (node->IsClone()) return;
336
337 /*recover pointers: */
338 inputs=(ParameterInputs*)vinputs;
339
340 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
341 if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
342
343 found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
344 if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
345
346 /*Get dof list: */
347 GetDofList(&doflist[0],&numberofdofspernode);
348
349 //Compute pressure melting point
350 meltingpoint=matpar->GetMeltingPoint();
351 beta=matpar->GetBeta();
352 t_pmp=meltingpoint-beta*pressure;
353
354 //Add penalty load
355 if (temperature<t_pmp){ //If T<Tpmp, there must be no melting. Therefore, melting should be constrained to 0 when T<Tpmp, instead of using spcs, use penalties
356 Ke[0][0]=kmax*pow((double)10,penalty_offset);
357 }
358
359 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
360}
361
362#undef __FUNCT__
363#define __FUNCT__ "Pengrid::PenaltyCreatePVector"
364void Pengrid::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
365
366 if (analysis_type==ThermalAnalysisEnum()){
367
368 PenaltyCreatePVectorThermal( pg,inputs,kmax,analysis_type,sub_analysis_type);
369
370 }
371 else if (analysis_type==MeltingAnalysisEnum()){
372
373 PenaltyCreatePVectorMelting( pg,inputs,kmax,analysis_type,sub_analysis_type);
374
375 }
376 else if (analysis_type==DiagnosticAnalysisEnum()){
377
378 /*No loads applied, do nothing: */
379 return;
380
381 }
382 else{
383 throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
384 }
385
386}
387
388Object* Pengrid::copy() {
389 return new Pengrid(*this);
390}
391
392
393void Pengrid::GetDofList(int* doflist,int* pnumberofdofspernode){
394
395 int j;
396 int doflist_per_node[MAXDOFSPERNODE];
397 int numberofdofspernode;
398
399 node->GetDofList(&doflist_per_node[0],&numberofdofspernode);
400 for(j=0;j<numberofdofspernode;j++){
401 doflist[j]=doflist_per_node[j];
402 }
403
404 /*Assign output pointers:*/
405 *pnumberofdofspernode=numberofdofspernode;
406}
407
408void Pengrid::PenaltyCreatePVectorThermal(Vec pg, void* vinputs, double kmax,int analysis_type,int sub_analysis_type){
409
410 const int numgrids=1;
411 const int NDOF1=1;
412 const int numdof=numgrids*NDOF1;
413 int doflist[numdof];
414 double P_terms[numdof]={0.0};
415 int numberofdofspernode;
416 int found=0;
417 double pressure;
418 int dofs1[1]={0};
419 double meltingpoint;
420 double beta;
421 double t_pmp;
422
423 ParameterInputs* inputs=NULL;
424
425 /*recover pointers: */
426 inputs=(ParameterInputs*)vinputs;
427
428 if(!active)return;
429
430 /*Get dof list: */
431 GetDofList(&doflist[0],&numberofdofspernode);
432
433 //First recover pressure
434 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
435 if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
436
437 //Compute pressure melting point
438 meltingpoint=matpar->GetMeltingPoint();
439 beta=matpar->GetBeta();
440 t_pmp=meltingpoint-beta*pressure;
441
442 //Add penalty load
443 P_terms[0]=kmax*pow((double)10,penalty_offset)*t_pmp;
444
445 /*Add P_terms to global vector pg: */
446 VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
447}
448
449void Pengrid::PenaltyCreatePVectorMelting(Vec pg, void* vinputs, double kmax,int analysis_type,int sub_analysis_type){
450
451 const int numgrids=1;
452 const int NDOF1=1;
453 const int numdof=numgrids*NDOF1;
454 int doflist[numdof];
455 double P_terms[numdof]={0.0};
456 int numberofdofspernode;
457 int found=0;
458 int dofs1[1]={0};
459 double pressure;
460 double temperature;
461 double melting_offset;
462 double meltingpoint;
463 double beta, heatcapacity;
464 double latentheat;
465 double t_pmp;
466 double dt;
467
468 ParameterInputs* inputs=NULL;
469
470 /*check that pengrid is not a clone (penalty to be added only once)*/
471 if (node->IsClone()) return;
472
473 /*recover pointers: */
474 inputs=(ParameterInputs*)vinputs;
475
476 /*Get dof list: */
477 GetDofList(&doflist[0],&numberofdofspernode);
478
479 //First recover pressure,melting offset and temperature vectors
480 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
481 if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
482
483 found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
484 if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
485
486 found=inputs->Recover("melting_offset",&melting_offset);
487 if(!found)throw ErrorException(__FUNCT__," could not find melting_offset in inputs!");
488
489 found=inputs->Recover("dt",&dt);
490 if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!");
491
492 meltingpoint=matpar->GetMeltingPoint();
493 beta=matpar->GetBeta();
494 heatcapacity=matpar->GetHeatCapacity();
495 latentheat=matpar->GetLatentHeat();
496
497 //Compute pressure melting point
498 t_pmp=meltingpoint-beta*pressure;
499
500 //Add penalty load
501 //This time, the penalty must have the same value as the one used for the thermal computation
502 //so that the corresponding melting can be computed correctly
503 //In the thermal computation, we used kmax=melting_offset, and the same penalty_offset
504 if (temperature<t_pmp){ //%no melting
505 P_terms[0]=0;
506 }
507 else{
508 if (dt){
509 P_terms[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp)/dt;
510 }
511 else{
512 P_terms[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp);
513 }
514 }
515
516 /*Add P_terms to global vector pg: */
517 VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
518}
519
520
521#undef __FUNCT__
522#define __FUNCT__ "Pengrid::PenaltyConstrain"
523void Pengrid::PenaltyConstrain(int* punstable,void* vinputs,int analysis_type,int sub_analysis_type){
524
525 if ((analysis_type==DiagnosticAnalysisEnum()) && ((sub_analysis_type==StokesAnalysisEnum()))){
526
527 /*No penalty to check*/
528 return;
529
530 }
531 else if (analysis_type==ThermalAnalysisEnum()){
532
533 PenaltyConstrainThermal(punstable,vinputs,analysis_type,sub_analysis_type);
534
535 }
536 else if (analysis_type==MeltingAnalysisEnum()){
537
538 /*No penalty to check*/
539 return;
540
541 }
542 else{
543 throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
544 }
545
546}
547
548#undef __FUNCT__
549#define __FUNCT__ "Pengrid::PenaltyConstrainThermal"
550void Pengrid::PenaltyConstrainThermal(int* punstable,void* vinputs,int analysis_type,int sub_analysis_type){
551
552 // The penalty is stable if it doesn't change during to successive iterations.
553
554 int found=0;
555 const int numgrids=1;
556
557
558 double pressure;
559 double temperature;
560 double beta,t_pmp;
561 double meltingpoint;
562 int new_active;
563 int dofs1[1]={0};
564 int unstable=0;
565 int reset_penalties=0;
566
567 ParameterInputs* inputs=NULL;
568
569 /*check that pengrid is not a clone (penalty to be added only once)*/
570 if (node->IsClone()){
571 unstable=0;
572 *punstable=unstable;
573 return;
574 }
575
576 /*recover pointers: */
577 inputs=(ParameterInputs*)vinputs;
578
579 //First recover beta, pressure and temperature vectors;
580 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
581 if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
582
583 found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
584 if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
585
586 found=inputs->Recover("reset_penalties",&reset_penalties);
587
588 if(reset_penalties)zigzag_counter=0;
589
590 //Compute pressure melting point
591 meltingpoint=matpar->GetMeltingPoint();
592 beta=matpar->GetBeta();
593
594 t_pmp=meltingpoint-beta*pressure;
595
596 //Figure out if temperature is over melting_point, in which case, this penalty needs to be activated.
597
598 if (temperature>t_pmp){
599 new_active=1;
600 }
601 else{
602 new_active=0;
603 }
604
605
606 //Figure out stability of this penalty
607 if (active==new_active){
608 unstable=0;
609 }
610 else{
611 unstable=1;
612 if(stabilize_constraints)zigzag_counter++;
613 }
614
615 /*If penalty keeps zigzagging more than 5 times: */
616 if(stabilize_constraints){
617 if(zigzag_counter>stabilize_constraints){
618 unstable=0;
619 active=1;
620 }
621 }
622
623 //Set penalty flag
624 active=new_active;
625
626 //*Assign output pointers:*/
627 *punstable=unstable;
628}
Note: See TracBrowser for help on using the repository browser.