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

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

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