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

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

Added DeepEcho function
File size: 16.6 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){
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
34 node_offset=UNDEF;
35 node=NULL;
36 matpar=NULL;
37 matpar_offset=UNDEF;
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);
50 printf(" mparid: %i\n",mparid);
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);
57 printf(" matpar_offset=%i\n",matpar_offset);
58
59 if(node)node->Echo();
60 if(matpar)matpar->Echo();
61 return;
62}
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}
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);
96 memcpy(marshalled_dataset,&mparid,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
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);
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);
105
106 *pmarshalled_dataset=marshalled_dataset;
107 return;
108}
109
110int Pengrid::MarshallSize(){
111
112 return sizeof(id)+
113 sizeof(mparid)+
114 sizeof(dof)+
115 sizeof(active)+
116 sizeof(penalty_offset)+
117 sizeof(thermal_steadystate)+
118 sizeof(node_id)+
119 sizeof(node_offset)+
120 sizeof(matpar)+
121 sizeof(matpar_offset)+
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);
141 memcpy(&mparid,marshalled_dataset,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
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);
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);
150
151
152 node=NULL;
153 matpar=NULL;
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}
170void Pengrid::DistributeNumDofs(int* numdofpernode,int analysis_type,int sub_analysis_type){return;}
171
172#undef __FUNCT__
173#define __FUNCT__ "Pengrid::Configure"
174
175void Pengrid::Configure(void* pelementsin,void* pnodesin,void* pmaterialsin){
176
177 DataSet* nodesin=NULL;
178 DataSet* materialsin=NULL;
179
180 /*Recover pointers :*/
181 nodesin=(DataSet*)pnodesin;
182 materialsin=(DataSet*)pmaterialsin;
183
184 /*Link this load with its nodes: */
185 ResolvePointers((Object**)&node,&node_id,&node_offset,1,nodesin);
186 ResolvePointers((Object**)&matpar,&mparid,&matpar_offset,1,materialsin);
187}
188
189
190#undef __FUNCT__
191#define __FUNCT__ "Pengrid::CreateKMatrix"
192
193void Pengrid::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
194
195 /*No loads applied, do nothing: */
196 return;
197
198}
199
200#undef __FUNCT__
201#define __FUNCT__ "Pengrid::CreatePVector"
202void Pengrid::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
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"
216void Pengrid::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
217
218 if ((analysis_type==DiagnosticAnalysisEnum()) && ((sub_analysis_type==StokesAnalysisEnum()))){
219
220 PenaltyCreateKMatrixDiagnosticStokes( Kgg,inputs,kmax,analysis_type,sub_analysis_type);
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);
230
231 }
232 else{
233 throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
234 }
235
236}
237
238#undef __FUNCT__
239#define __FUNCT__ "Pengrid::PenaltyCreateKMatrixDiagnosticStokes"
240void Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
241
242 const int numgrids=1;
243 const int NDOF4=4;
244 const int numdof=numgrids*NDOF4;
245 int doflist[numdof];
246 int numberofdofspernode;
247
248 int dofs1[1]={0};
249 int dofs2[1]={1};
250 double slope[2];
251 int found=0;
252 double Ke[4][4]={0.0};
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: */
263 found=inputs->Recover("bedslopex",&slope[0],1,dofs1,numgrids,(void**)&node);
264 if(!found)throw ErrorException(__FUNCT__," bedslopex needed in inputs!");
265 found=inputs->Recover("bedslopey",&slope[1],1,dofs2,numgrids,(void**)&node);
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)
269 Ke[2][0]=-slope[0]*kmax*pow(10.0,penalty_offset);
270 Ke[2][1]=-slope[1]*kmax*pow(10.0,penalty_offset);
271 Ke[2][2]=kmax*pow(10,penalty_offset);
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__
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
301 Ke[0][0]=kmax*pow(10,penalty_offset);
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
311 int found=0;
312 const int numgrids=1;
313 const int NDOF1=1;
314 const int numdof=numgrids*NDOF1;
315 double Ke[numdof][numdof]={0.0};
316 int dofs1[1]={0};
317 int doflist[numdof];
318 int numberofdofspernode;
319 double meltingpoint;
320
321 double pressure;
322 double temperature;
323 double beta,t_pmp;
324
325 ParameterInputs* inputs=NULL;
326
327 /*recover pointers: */
328 inputs=(ParameterInputs*)vinputs;
329
330 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
331 if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
332
333 found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
334 if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
335
336 /*Get dof list: */
337 GetDofList(&doflist[0],&numberofdofspernode);
338
339 //Compute pressure melting point
340 meltingpoint=matpar->GetMeltingPoint();
341 beta=matpar->GetBeta();
342 t_pmp=meltingpoint-beta*pressure;
343
344 //Add penalty load
345 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
346 Ke[0][0]=kmax*pow(10,penalty_offset);
347 }
348
349 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
350}
351
352#undef __FUNCT__
353#define __FUNCT__ "Pengrid::PenaltyCreatePVector"
354void Pengrid::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
355
356 if (analysis_type==ThermalAnalysisEnum()){
357
358 PenaltyCreatePVectorThermal( pg,inputs,kmax,analysis_type,sub_analysis_type);
359
360 }
361 else if (analysis_type==MeltingAnalysisEnum()){
362
363 PenaltyCreatePVectorMelting( pg,inputs,kmax,analysis_type,sub_analysis_type);
364
365 }
366 else if (analysis_type==DiagnosticAnalysisEnum()){
367
368 /*No loads applied, do nothing: */
369 return;
370
371 }
372 else{
373 throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
374 }
375
376}
377
378Object* Pengrid::copy() {
379 return new Pengrid(*this);
380}
381
382
383void Pengrid::GetDofList(int* doflist,int* pnumberofdofspernode){
384
385 int j;
386 int doflist_per_node[MAXDOFSPERNODE];
387 int numberofdofspernode;
388
389 node->GetDofList(&doflist_per_node[0],&numberofdofspernode);
390 for(j=0;j<numberofdofspernode;j++){
391 doflist[j]=doflist_per_node[j];
392 }
393
394 /*Assign output pointers:*/
395 *pnumberofdofspernode=numberofdofspernode;
396}
397
398void Pengrid::PenaltyCreatePVectorThermal(Vec pg, void* vinputs, double kmax,int analysis_type,int sub_analysis_type){
399
400 const int numgrids=1;
401 const int NDOF1=1;
402 const int numdof=numgrids*NDOF1;
403 int doflist[numdof];
404 double P_terms[numdof]={0.0};
405 int numberofdofspernode;
406 int found=0;
407 double pressure;
408 int dofs1[1]={0};
409 double meltingpoint;
410 double beta;
411 double t_pmp;
412
413 ParameterInputs* inputs=NULL;
414
415 /*recover pointers: */
416 inputs=(ParameterInputs*)vinputs;
417
418 if(!active)return;
419
420 /*Get dof list: */
421 GetDofList(&doflist[0],&numberofdofspernode);
422
423 //First recover pressure
424 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
425 if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
426
427 //Compute pressure melting point
428 meltingpoint=matpar->GetMeltingPoint();
429 beta=matpar->GetBeta();
430 t_pmp=meltingpoint-beta*pressure;
431
432 //Add penalty load
433 P_terms[0]=kmax*pow(10,penalty_offset)*t_pmp;
434
435 /*Add P_terms to global vector pg: */
436 VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
437}
438
439void Pengrid::PenaltyCreatePVectorMelting(Vec pg, void* vinputs, double kmax,int analysis_type,int sub_analysis_type){
440
441 const int numgrids=1;
442 const int NDOF1=1;
443 const int numdof=numgrids*NDOF1;
444 int doflist[numdof];
445 double P_terms[numdof]={0.0};
446 int numberofdofspernode;
447 int found=0;
448 int dofs1[1]={0};
449 double pressure;
450 double temperature;
451 double melting_offset;
452 double meltingpoint;
453 double beta, heatcapacity;
454 double latentheat;
455 double t_pmp;
456 double dt;
457
458 ParameterInputs* inputs=NULL;
459
460 /*recover pointers: */
461 inputs=(ParameterInputs*)vinputs;
462
463 /*Get dof list: */
464 GetDofList(&doflist[0],&numberofdofspernode);
465
466 //First recover pressure,melting offset and temperature vectors
467 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
468 if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
469
470 found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
471 if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
472
473 found=inputs->Recover("melting_offset",&melting_offset);
474 if(!found)throw ErrorException(__FUNCT__," could not find melting_offset in inputs!");
475
476 found=inputs->Recover("dt",&dt);
477 if((!found) && (sub_analysis_type==TransientAnalysisEnum()))throw ErrorException(__FUNCT__," could not find dt in inputs!");
478
479 meltingpoint=matpar->GetMeltingPoint();
480 beta=matpar->GetBeta();
481 heatcapacity=matpar->GetHeatCapacity();
482 latentheat=matpar->GetLatentHeat();
483
484 //Compute pressure melting point
485 t_pmp=meltingpoint-beta*pressure;
486
487 //Add penalty load
488 //This time, the penalty must have the same value as the one used for the thermal computation
489 //so that the corresponding melting can be computed correctly
490 //In the thermal computation, we used kmax=melting_offset, and the same penalty_offset
491 if (temperature<t_pmp){ //%no melting
492 P_terms[0]=0;
493 }
494 else{
495 if (sub_analysis_type==SteadyAnalysisEnum()){
496 P_terms[0]=melting_offset*pow(10,penalty_offset)*(temperature-t_pmp);
497 }
498 else{
499 P_terms[0]=melting_offset*pow(10,penalty_offset)*(temperature-t_pmp)/dt;
500 }
501 }
502
503 /*Add P_terms to global vector pg: */
504 VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
505}
506
507
508#undef __FUNCT__
509#define __FUNCT__ "Pengrid::PenaltyConstrain"
510void Pengrid::PenaltyConstrain(int* punstable,void* vinputs,int analysis_type,int sub_analysis_type){
511
512 if ((analysis_type==DiagnosticAnalysisEnum()) && ((sub_analysis_type==StokesAnalysisEnum()))){
513
514 /*No penalty to check*/
515 return;
516
517 }
518 else if (analysis_type==ThermalAnalysisEnum()){
519
520 PenaltyConstrainThermal(punstable,vinputs,analysis_type,sub_analysis_type);
521
522 }
523 else if (analysis_type==MeltingAnalysisEnum()){
524
525 /*No penalty to check*/
526 return;
527
528 }
529 else{
530 throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
531 }
532
533}
534
535#undef __FUNCT__
536#define __FUNCT__ "Pengrid::PenaltyConstrainThermal"
537void Pengrid::PenaltyConstrainThermal(int* punstable,void* vinputs,int analysis_type,int sub_analysis_type){
538
539 // The penalty is stable if it doesn't change during to successive iterations.
540
541 int found=0;
542 const int numgrids=1;
543
544
545 double pressure;
546 double temperature;
547 double beta,t_pmp;
548 double meltingpoint;
549 int new_active;
550 int dofs1[1]={0};
551 int unstable=0;
552
553 ParameterInputs* inputs=NULL;
554
555 /*recover pointers: */
556 inputs=(ParameterInputs*)vinputs;
557
558
559 //First recover beta, pressure and temperature vectors;
560 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
561 if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
562
563 found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
564 if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
565
566
567 //Compute pressure melting point
568 meltingpoint=matpar->GetMeltingPoint();
569 beta=matpar->GetBeta();
570
571 t_pmp=meltingpoint-beta*pressure;
572
573 //Figure out if temperature is over melting_point, in which case, this penalty needs to be activated.
574
575 if (temperature>t_pmp){
576 new_active=1;
577 }
578 else{
579 new_active=0;
580 }
581
582
583 //Figure out stability of this penalty
584 if (active==new_active){
585 unstable=0;
586 }
587 else{
588 unstable=1;
589 }
590
591 //Set penalty flag
592 active=new_active;
593
594 //*Assign output pointers:*/
595 *punstable=unstable;
596}
Note: See TracBrowser for help on using the repository browser.