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

Last change on this file since 488 was 488, checked in by seroussi, 16 years ago

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