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

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

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