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

Last change on this file since 539 was 539, checked in by Mathieu Morlighem, 17 years ago

get doflist missing
File size: 16.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_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 if (analysis_type==DiagnosticAnalysisEnum()){
351
352 /*No loads applied, do nothing: */
353 return;
354
355 }
356 else{
357 throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
358 }
359
360}
361
362Object* Pengrid::copy() {
363 return new Pengrid(*this);
364}
365
366
367void Pengrid::GetDofList(int* doflist,int* pnumberofdofspernode){
368
369 int j;
370 int doflist_per_node[MAXDOFSPERNODE];
371 int numberofdofspernode;
372
373 node->GetDofList(&doflist_per_node[0],&numberofdofspernode);
374 for(j=0;j<numberofdofspernode;j++){
375 doflist[j]=doflist_per_node[j];
376 }
377
378 /*Assign output pointers:*/
379 *pnumberofdofspernode=numberofdofspernode;
380}
381
382void Pengrid::PenaltyCreatePVectorThermal(Vec pg, void* vinputs, double kmax,int analysis_type,int sub_analysis_type){
383
384 const int numgrids=1;
385 const int NDOF1=1;
386 const int numdof=numgrids*NDOF1;
387 int doflist[numdof];
388 double P_terms[numdof]={0.0};
389 int numberofdofspernode;
390 int found=0;
391 double pressure;
392 int dofs1[1]={0};
393 double meltingpoint;
394 double beta;
395 double t_pmp;
396
397 ParameterInputs* inputs=NULL;
398
399 /*recover pointers: */
400 inputs=(ParameterInputs*)vinputs;
401
402 if(!active)return;
403
404 /*Get dof list: */
405 GetDofList(&doflist[0],&numberofdofspernode);
406
407 //First recover pressure
408 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
409 if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
410
411 //Compute pressure melting point
412 meltingpoint=matpar->GetMeltingPoint();
413 beta=matpar->GetBeta();
414 t_pmp=meltingpoint-beta*pressure;
415
416 //Add penalty load
417 P_terms[0]=kmax*pow(10,penalty_offset)*t_pmp;
418
419 /*Add P_terms to global vector pg: */
420 VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
421}
422
423void Pengrid::PenaltyCreatePVectorMelting(Vec pg, void* vinputs, double kmax,int analysis_type,int sub_analysis_type){
424
425 const int numgrids=1;
426 const int NDOF1=1;
427 const int numdof=numgrids*NDOF1;
428 int doflist[numdof];
429 double P_terms[numdof]={0.0};
430 int numberofdofspernode;
431 int found=0;
432 int dofs1[1]={0};
433 double pressure;
434 double temperature;
435 double melting_offset;
436 double meltingpoint;
437 double beta, heatcapacity;
438 double latentheat;
439 double t_pmp;
440 double dt;
441
442 ParameterInputs* inputs=NULL;
443
444 /*recover pointers: */
445 inputs=(ParameterInputs*)vinputs;
446
447 /*Get dof list: */
448 GetDofList(&doflist[0],&numberofdofspernode);
449
450 //First recover pressure,melting offset and temperature vectors
451 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
452 if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
453
454 found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
455 if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
456
457 found=inputs->Recover("melting_offset",&melting_offset);
458 if(!found)throw ErrorException(__FUNCT__," could not find melting_offset in inputs!");
459
460 found=inputs->Recover("dt",&dt);
461 if((!found) && (sub_analysis_type==TransientAnalysisEnum()))throw ErrorException(__FUNCT__," could not find dt in inputs!");
462
463 meltingpoint=matpar->GetMeltingPoint();
464 beta=matpar->GetBeta();
465 heatcapacity=matpar->GetHeatCapacity();
466 latentheat=matpar->GetLatentHeat();
467
468 //Compute pressure melting point
469 t_pmp=meltingpoint-beta*pressure;
470
471 //Add penalty load
472 //This time, the penalty must have the same value as the one used for the thermal computation
473 //so that the corresponding melting can be computed correctly
474 //In the thermal computation, we used kmax=melting_offset, and the same penalty_offset
475 if (temperature<t_pmp){ //%no melting
476 P_terms[0]=0;
477 }
478 else{
479 if (sub_analysis_type==SteadyAnalysisEnum()){
480 P_terms[0]=melting_offset*pow(10,penalty_offset)*(temperature-t_pmp);
481 }
482 else{
483 P_terms[0]=melting_offset*pow(10,penalty_offset)*(temperature-t_pmp)/dt;
484 }
485 }
486
487 /*Add P_terms to global vector pg: */
488 VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
489}
490
491
492#undef __FUNCT__
493#define __FUNCT__ "Pengrid::PenaltyConstrain"
494void Pengrid::PenaltyConstrain(int* punstable,void* vinputs,int analysis_type,int sub_analysis_type){
495
496 if ((analysis_type==DiagnosticAnalysisEnum()) && ((sub_analysis_type==StokesAnalysisEnum()))){
497
498 /*No penalty to check*/
499 return;
500
501 }
502 else if (analysis_type==ThermalAnalysisEnum()){
503
504 PenaltyConstrainThermal(punstable,vinputs,analysis_type,sub_analysis_type);
505
506 }
507 else if (analysis_type==MeltingAnalysisEnum()){
508
509 /*No penalty to check*/
510 return;
511
512 }
513 else{
514 throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
515 }
516
517}
518
519#undef __FUNCT__
520#define __FUNCT__ "Pengrid::PenaltyConstrainThermal"
521void Pengrid::PenaltyConstrainThermal(int* punstable,void* vinputs,int analysis_type,int sub_analysis_type){
522
523 // The penalty is stable if it doesn't change during to successive iterations.
524
525 int found=0;
526 const int numgrids=1;
527
528
529 double pressure;
530 double temperature;
531 double beta,t_pmp;
532 double meltingpoint;
533 int new_active;
534 int dofs1[1]={0};
535 int unstable=0;
536
537 ParameterInputs* inputs=NULL;
538
539 /*recover pointers: */
540 inputs=(ParameterInputs*)vinputs;
541
542
543 //First recover beta, pressure and temperature vectors;
544 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
545 if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
546
547 found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
548 if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
549
550
551 //Compute pressure melting point
552 meltingpoint=matpar->GetMeltingPoint();
553 beta=matpar->GetBeta();
554
555 t_pmp=meltingpoint-beta*pressure;
556
557 //Figure out if temperature is over melting_point, in which case, this penalty needs to be activated.
558
559 if (temperature>t_pmp){
560 new_active=1;
561 }
562 else{
563 new_active=0;
564 }
565
566
567 //Figure out stability of this penalty
568 if (active==new_active){
569 unstable=0;
570 }
571 else{
572 unstable=1;
573 }
574
575 //Set penalty flag
576 active=new_active;
577
578 //*Assign output pointers:*/
579 *punstable=unstable;
580}
Note: See TracBrowser for help on using the repository browser.