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

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

Almost there
File size: 22.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 "../include/macros.h"
17#include "../shared/shared.h"
18#include "../include/typedefs.h"
19#include "../include/macros.h"
20#include "../DataSet/DataSet.h"
21#include "../DataSet/Inputs.h"
22
23/*Object constructors and destructor*/
24/*FUNCTION Pengrid::constructor {{{1*/
25Pengrid::Pengrid(){
26 this->inputs=NULL;
27 this->parameters=NULL;
28
29 /*not active, not zigzagging: */
30 active=0;
31 zigzag_counter=0;
32
33}
34/*}}}1*/
35/*FUNCTION Pengrid::Pengrid(int id, int node_ids int matpar_id){{{1*/
36Pengrid::Pengrid(int pengrid_id,int pengrid_node_id, int pengrid_matpar_id):
37 hnode(pengrid_node_ids,1),
38 hmatice(&pengrid_matice_id,1),
39 hmatpar(&pengrid_matpar_id,1)
40{
41
42 /*all the initialization has been done by the initializer, just fill in the id: */
43 this->id=pengrid_id;
44 this->parameters=NULL;
45 this->inputs=new Inputs();
46
47 /*not active, not zigzagging: */
48 active=0;
49 zigzag_counter=0;
50
51}
52/*}}}*/
53/*FUNCTION Pengrid::Pengrid(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, DataSet* parameters, Inputs* pengrid_inputs) {{{1*/
54Pengrid::Pengrid(int pengrid_id,Hook* pengrid_hnode, Hook* pengrid_hmatpar, Parameters* pengrid_parameters, Inputs* pengrid_inputs):
55 hnode(pengrid_hnode),
56 hmatpar(pengrid_hmatpar)
57{
58
59 /*all the initialization has been done by the initializer, just fill in the id: */
60 this->id=pengrid_id;
61 if(pengrid_inputs){
62 this->inputs=(Inputs*)pengrid_inputs->Copy();
63 }
64 else{
65 this->inputs=new Inputs();
66 }
67 /*point parameters: */
68 this->parameters=pengrid_parameters;
69
70 /*not active, not zigzagging: */
71 active=0;
72 zigzag_counter=0;
73
74}
75/*}}}*/
76Pengrid::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,int pengrid_stabilize_constraints){
77
78 id=pengrid_id;
79 node_id=pengrid_node_id;
80 mparid=pengrid_mparid;
81 dof=pengrid_dof;
82 active=pengrid_active;
83 penalty_offset =pengrid_penalty_offset;
84 thermal_steadystate=pengrid_thermal_steadystate;
85 stabilize_constraints=pengrid_stabilize_constraints;
86
87 node_offset=UNDEF;
88 node=NULL;
89 matpar=NULL;
90 matpar_offset=UNDEF;
91
92
93 return;
94}
95/*}}}1*/
96/*FUNCTION Pengrid::Pengrid(int i, IoModel* iomodel){{{1*/
97Pengrid::Pengrid(int index, IoModel* iomodel){ //i is the element index
98
99 int i,j;
100 int tria_node_ids[3];
101 int tria_matice_id;
102 int tria_matpar_id;
103 double nodeinputs[3];
104
105 /*id: */
106 this->id=index+1;
107
108 /*hooks: */
109 //go recover node ids, needed to initialize the node hook.
110 if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==Balancedthickness2AnalysisEnum){
111 /*Discontinuous Galerkin*/
112 tria_node_ids[0]=3*index+1;
113 tria_node_ids[1]=3*index+2;
114 tria_node_ids[2]=3*index+3;
115 }
116 else{
117 /*Continuous Galerkin*/
118 for(i=0;i<3;i++){
119 tria_node_ids[i]=(int)*(iomodel->elements+3*index+i); //ids for vertices are in the elements array from Matlab
120 }
121 }
122 tria_matice_id=index+1; //refers to the corresponding ice material object
123 tria_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object
124
125 this->hnodes.Init(tria_node_ids,3);
126 this->hmatice.Init(&tria_matice_id,1);
127 this->hmatpar.Init(&tria_matpar_id,1);
128
129 //intialize inputs, and add as many inputs per element as requested:
130 this->inputs=new Inputs();
131
132 if (iomodel->thickness) {
133 for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness[tria_node_ids[i]-1];
134 this->inputs->AddInput(new PengridVertexInput(ThicknessEnum,nodeinputs));
135 }
136 if (iomodel->surface) {
137 for(i=0;i<3;i++)nodeinputs[i]=iomodel->surface[tria_node_ids[i]-1];
138 this->inputs->AddInput(new PengridVertexInput(SurfaceEnum,nodeinputs));
139 }
140 if (iomodel->bed) {
141 for(i=0;i<3;i++)nodeinputs[i]=iomodel->bed[tria_node_ids[i]-1];
142 this->inputs->AddInput(new PengridVertexInput(BedEnum,nodeinputs));
143 }
144 if (iomodel->drag_coefficient) {
145 for(i=0;i<3;i++)nodeinputs[i]=iomodel->drag_coefficient[tria_node_ids[i]-1];
146 this->inputs->AddInput(new PengridVertexInput(DragCoefficientEnum,nodeinputs));
147
148 if (iomodel->drag_p) this->inputs->AddInput(new DoubleInput(DragPEnum,iomodel->drag_p[index]));
149 if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index]));
150 this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type));
151
152 }
153 if (iomodel->melting_rate) {
154 for(i=0;i<3;i++)nodeinputs[i]=iomodel->melting_rate[tria_node_ids[i]-1];
155 this->inputs->AddInput(new PengridVertexInput(MeltingRateEnum,nodeinputs));
156 }
157 if (iomodel->accumulation_rate) {
158 for(i=0;i<3;i++)nodeinputs[i]=iomodel->accumulation_rate[tria_node_ids[i]-1];
159 this->inputs->AddInput(new PengridVertexInput(AccumulationRateEnum,nodeinputs));
160 }
161 if (iomodel->geothermalflux) {
162 for(i=0;i<3;i++)nodeinputs[i]=iomodel->geothermalflux[tria_node_ids[i]-1];
163 this->inputs->AddInput(new PengridVertexInput(GeothermalFluxEnum,nodeinputs));
164 }
165
166 if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index]));
167 if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index]));
168 if (iomodel->elementonwater) this->inputs->AddInput(new BoolInput(ElementOnWaterEnum,(IssmBool)iomodel->elementonwater[index]));
169 if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index]));
170
171 //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
172 this->parameters=NULL;
173
174
175}
176/*}}}*/
177/*FUNCTION Pengrid::destructor {{{1*/
178Pengrid::~Pengrid(){
179 return;
180}
181/*}}}1*/
182
183/*Object marshall*/
184/*FUNCTION Pengrid::Marshall {{{1*/
185void Pengrid::Marshall(char** pmarshalled_dataset){
186
187 char* marshalled_dataset=NULL;
188 int enum_type=0;
189
190 /*recover marshalled_dataset: */
191 marshalled_dataset=*pmarshalled_dataset;
192
193 /*get enum type of Pengrid: */
194 enum_type=PengridEnum;
195
196 /*marshall enum: */
197 memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
198
199 /*marshall Pengrid data: */
200 memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
201 memcpy(marshalled_dataset,&mparid,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
202 memcpy(marshalled_dataset,&dof,sizeof(dof));marshalled_dataset+=sizeof(dof);
203 memcpy(marshalled_dataset,&active,sizeof(active));marshalled_dataset+=sizeof(active);
204 memcpy(marshalled_dataset,&penalty_offset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
205 memcpy(marshalled_dataset,&thermal_steadystate,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
206 memcpy(marshalled_dataset,&node_id,sizeof(node_id));marshalled_dataset+=sizeof(node_id);
207 memcpy(marshalled_dataset,&node_offset,sizeof(node_offset));marshalled_dataset+=sizeof(node_offset);
208 memcpy(marshalled_dataset,&matpar,sizeof(matpar));marshalled_dataset+=sizeof(matpar);
209 memcpy(marshalled_dataset,&matpar_offset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
210 memcpy(marshalled_dataset,&stabilize_constraints,sizeof(stabilize_constraints));marshalled_dataset+=sizeof(stabilize_constraints);
211 memcpy(marshalled_dataset,&zigzag_counter,sizeof(zigzag_counter));marshalled_dataset+=sizeof(zigzag_counter);
212
213 *pmarshalled_dataset=marshalled_dataset;
214 return;
215}
216/*}}}1*/
217/*FUNCTION Pengrid::MarshallSize {{{1*/
218int Pengrid::MarshallSize(){
219
220 return sizeof(id)+
221 sizeof(mparid)+
222 sizeof(dof)+
223 sizeof(active)+
224 sizeof(penalty_offset)+
225 sizeof(thermal_steadystate)+
226 sizeof(node_id)+
227 sizeof(node_offset)+
228 sizeof(matpar)+
229 sizeof(matpar_offset)+
230 sizeof(stabilize_constraints)+
231 sizeof(zigzag_counter)+
232 sizeof(int); //sizeof(int) for enum type
233}
234/*}}}1*/
235/*FUNCTION Pengrid::Demarshall {{{1*/
236void Pengrid::Demarshall(char** pmarshalled_dataset){
237
238 char* marshalled_dataset=NULL;
239
240 /*recover marshalled_dataset: */
241 marshalled_dataset=*pmarshalled_dataset;
242
243 /*this time, no need to get enum type, the pointer directly points to the beginning of the
244 *object data (thanks to DataSet::Demarshall):*/
245
246 memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
247 memcpy(&mparid,marshalled_dataset,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
248 memcpy(&dof,marshalled_dataset,sizeof(dof));marshalled_dataset+=sizeof(dof);
249 memcpy(&active,marshalled_dataset,sizeof(active));marshalled_dataset+=sizeof(active);
250 memcpy(&penalty_offset,marshalled_dataset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
251 memcpy(&thermal_steadystate,marshalled_dataset,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
252 memcpy(&node_id,marshalled_dataset,sizeof(node_id));marshalled_dataset+=sizeof(node_id);
253 memcpy(&node_offset,marshalled_dataset,sizeof(node_offset));marshalled_dataset+=sizeof(node_offset);
254 memcpy(&matpar,marshalled_dataset,sizeof(matpar));marshalled_dataset+=sizeof(matpar);
255 memcpy(&matpar_offset,marshalled_dataset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
256 memcpy(&stabilize_constraints,marshalled_dataset,sizeof(stabilize_constraints));marshalled_dataset+=sizeof(stabilize_constraints);
257 memcpy(&zigzag_counter,marshalled_dataset,sizeof(zigzag_counter));marshalled_dataset+=sizeof(zigzag_counter);
258
259 node=NULL;
260 matpar=NULL;
261
262 /*return: */
263 *pmarshalled_dataset=marshalled_dataset;
264 return;
265}
266/*}}}1*/
267
268/*Object functions*/
269/*FUNCTION Pengrid::copy {{{1*/
270Object* Pengrid::copy() {
271 return new Pengrid(*this);
272}
273/*}}}1*/
274/*FUNCTION Pengrid::Configure {{{1*/
275
276void Pengrid::Configure(void* pelementsin,void* pnodesin,void* pmaterialsin){
277
278 DataSet* nodesin=NULL;
279 DataSet* materialsin=NULL;
280
281 /*Recover pointers :*/
282 nodesin=(DataSet*)pnodesin;
283 materialsin=(DataSet*)pmaterialsin;
284
285 /*Link this load with its nodes: */
286 ResolvePointers((Object**)&node,&node_id,&node_offset,1,nodesin);
287 ResolvePointers((Object**)&matpar,&mparid,&matpar_offset,1,materialsin);
288}
289/*}}}1*/
290/*FUNCTION Pengrid::CreateKMatrix {{{1*/
291
292void Pengrid::CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type){
293
294 /*No loads applied, do nothing: */
295 return;
296
297}
298/*}}}1*/
299/*FUNCTION Pengrid::CreatePVector {{{1*/
300void Pengrid::CreatePVector(Vec pg, int analysis_type,int sub_analysis_type){
301
302 /*No loads applied, do nothing: */
303 return;
304
305}
306/*}}}1*/
307/*FUNCTION Pengrid::DeepEcho {{{1*/
308void Pengrid::DeepEcho(void){
309
310 printf("Pengrid:\n");
311 printf(" id: %i\n",id);
312 printf(" mparid: %i\n",mparid);
313 printf(" dof: %i\n",dof);
314 printf(" active: %i\n",active);
315 printf(" penalty_offset: %g\n",penalty_offset);
316 printf(" thermal_steadystate: %i\n",thermal_steadystate);
317 printf(" node_id: [%i]\n",node_id);
318 printf(" node_offset: [%i]\n",node_offset);
319 printf(" matpar_offset=%i\n",matpar_offset);
320
321 if(node)node->Echo();
322 if(matpar)matpar->Echo();
323 return;
324}
325/*}}}1*/
326/*FUNCTION Pengrid::DistributenumDofs {{{1*/
327void Pengrid::DistributeNumDofs(int* numdofpernode,int analysis_type,int sub_analysis_type){return;}
328/*}}}1*/
329/*FUNCTION Pengrid::Echo {{{1*/
330void Pengrid::Echo(void){
331
332 printf("Pengrid:\n");
333 printf(" id: %i\n",id);
334 printf(" mparid: %i\n",mparid);
335 printf(" dof: %i\n",dof);
336 printf(" active: %i\n",active);
337 printf(" penalty_offset: %g\n",penalty_offset);
338 printf(" thermal_steadystate: %i\n",thermal_steadystate);
339 printf(" node_id: [%i]\n",node_id);
340 printf(" node_offset: [%i]\n",node_offset);
341 printf(" matpar_offset=%i\n",matpar_offset);
342
343 return;
344}
345/*}}}1*/
346/*FUNCTION Pengrid::Enum {{{1*/
347int Pengrid::Enum(void){
348
349 return PengridEnum;
350}
351/*}}}1*/
352/*FUNCTION Pengrid::GetDofList {{{1*/
353void Pengrid::GetDofList(int* doflist,int* pnumberofdofspernode){
354
355 int j;
356 int doflist_per_node[MAXDOFSPERNODE];
357 int numberofdofspernode;
358
359 node->GetDofList(&doflist_per_node[0],&numberofdofspernode);
360 for(j=0;j<numberofdofspernode;j++){
361 doflist[j]=doflist_per_node[j];
362 }
363
364 /*Assign output pointers:*/
365 *pnumberofdofspernode=numberofdofspernode;
366}
367/*}}}1*/
368/*FUNCTION Pengrid::GetId {{{1*/
369int Pengrid::GetId(void){ return id; }
370/*}}}1*/
371/*FUNCTION Pengrid::GetName {{{1*/
372char* Pengrid::GetName(void){
373 return "pengrid";
374}
375/*}}}1*/
376/*FUNCTION Pengrid::MyRank {{{1*/
377int Pengrid::MyRank(void){
378 extern int my_rank;
379 return my_rank;
380}
381/*}}}1*/
382/*FUNCTION Pengrid::PenaltyConstrain {{{1*/
383void Pengrid::PenaltyConstrain(int* punstable,int analysis_type,int sub_analysis_type){
384
385 if ((analysis_type==DiagnosticAnalysisEnum) && ((sub_analysis_type==StokesAnalysisEnum))){
386
387 /*No penalty to check*/
388 return;
389
390 }
391 else if (analysis_type==ThermalAnalysisEnum){
392
393 PenaltyConstrainThermal(punstable,analysis_type,sub_analysis_type);
394
395 }
396 else if (analysis_type==MeltingAnalysisEnum){
397
398 /*No penalty to check*/
399 return;
400
401 }
402 else{
403 ISSMERROR("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet");
404 }
405
406}
407/*}}}1*/
408/*FUNCTION Pengrid::PenaltyConstrainThermal {{{1*/
409void Pengrid::PenaltyConstrainThermal(int* punstable,int analysis_type,int sub_analysis_type){
410
411 // The penalty is stable if it doesn't change during to successive iterations.
412
413 int found=0;
414 const int numgrids=1;
415
416
417 double pressure;
418 double temperature;
419 double beta,t_pmp;
420 double meltingpoint;
421 int new_active;
422 int dofs1[1]={0};
423 int unstable=0;
424 int reset_penalties=0;
425
426 ParameterInputs* inputs=NULL;
427
428 /*check that pengrid is not a clone (penalty to be added only once)*/
429 if (node->IsClone()){
430 unstable=0;
431 *punstable=unstable;
432 return;
433 }
434
435 /*recover pointers: */
436 inputs=(ParameterInputs*)vinputs;
437
438 //First recover beta, pressure and temperature vectors;
439 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
440 if(!found)ISSMERROR(" could not find pressure in inputs!");
441
442 found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
443 if(!found)ISSMERROR(" could not find temperature in inputs!");
444
445 found=inputs->Recover("reset_penalties",&reset_penalties);
446
447 if(reset_penalties)zigzag_counter=0;
448
449 //Compute pressure melting point
450 meltingpoint=matpar->GetMeltingPoint();
451 beta=matpar->GetBeta();
452
453 t_pmp=meltingpoint-beta*pressure;
454
455 //Figure out if temperature is over melting_point, in which case, this penalty needs to be activated.
456
457 if (temperature>t_pmp){
458 new_active=1;
459 }
460 else{
461 new_active=0;
462 }
463
464
465 //Figure out stability of this penalty
466 if (active==new_active){
467 unstable=0;
468 }
469 else{
470 unstable=1;
471 if(stabilize_constraints)zigzag_counter++;
472 }
473
474 /*If penalty keeps zigzagging more than 5 times: */
475 if(stabilize_constraints){
476 if(zigzag_counter>stabilize_constraints){
477 unstable=0;
478 active=1;
479 }
480 }
481
482 //Set penalty flag
483 active=new_active;
484
485 //*Assign output pointers:*/
486 *punstable=unstable;
487}
488/*}}}1*/
489/*FUNCTION Pengrid::PenaltyCreateMatrix {{{1*/
490void Pengrid::PenaltyCreateKMatrix(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type){
491
492 if ((analysis_type==DiagnosticAnalysisEnum) && ((sub_analysis_type==StokesAnalysisEnum))){
493
494 PenaltyCreateKMatrixDiagnosticStokes( Kgg,kmax,analysis_type,sub_analysis_type);
495 }
496 else if (analysis_type==ThermalAnalysisEnum){
497
498 PenaltyCreateKMatrixThermal( Kgg,kmax,analysis_type,sub_analysis_type);
499
500 }
501 else if (analysis_type==MeltingAnalysisEnum){
502
503 PenaltyCreateKMatrixMelting( Kgg,kmax,analysis_type,sub_analysis_type);
504
505 }
506 else{
507 ISSMERROR("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet");
508 }
509
510}
511/*}}}1*/
512/*FUNCTION Pengrid::PenaltyCreateKMatrixDiagnosticStokes {{{1*/
513void Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type){
514
515 const int numgrids=1;
516 const int NDOF4=4;
517 const int numdof=numgrids*NDOF4;
518 int doflist[numdof];
519 int numberofdofspernode;
520
521 int dofs1[1]={0};
522 int dofs2[1]={1};
523 double slope[2];
524 int found=0;
525 double Ke[4][4]={0.0};
526
527 ParameterInputs* inputs=NULL;
528
529 /*recover pointers: */
530 inputs=(ParameterInputs*)vinputs;
531
532 /*Get dof list: */
533 GetDofList(&doflist[0],&numberofdofspernode);
534
535 /*recover slope: */
536 found=inputs->Recover("bedslopex",&slope[0],1,dofs1,numgrids,(void**)&node);
537 if(!found)ISSMERROR(" bedslopex needed in inputs!");
538 found=inputs->Recover("bedslopey",&slope[1],1,dofs2,numgrids,(void**)&node);
539 if(!found)ISSMERROR(" bedslopey needed in inputs!");
540
541 //Create elementary matrix: add penalty to contrain wb (wb=ub*db/dx+vb*db/dy)
542 Ke[2][0]=-slope[0]*kmax*pow((double)10.0,penalty_offset);
543 Ke[2][1]=-slope[1]*kmax*pow((double)10.0,penalty_offset);
544 Ke[2][2]=kmax*pow((double)10,penalty_offset);
545
546 /*Add Ke to global matrix Kgg: */
547 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
548}
549/*}}}1*/
550/*FUNCTION Pengrid::PenaltyCreateKMatrixMelting {{{1*/
551void Pengrid::PenaltyCreateKMatrixMelting(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type){
552
553
554 int found=0;
555 const int numgrids=1;
556 const int NDOF1=1;
557 const int numdof=numgrids*NDOF1;
558 double Ke[numdof][numdof]={0.0};
559 int dofs1[1]={0};
560 int doflist[numdof];
561 int numberofdofspernode;
562 double meltingpoint;
563
564 double pressure;
565 double temperature;
566 double beta,t_pmp;
567
568 ParameterInputs* inputs=NULL;
569
570 /*check that pengrid is not a clone (penalty to be added only once)*/
571 if (node->IsClone()) return;
572
573 /*recover pointers: */
574 inputs=(ParameterInputs*)vinputs;
575
576 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
577 if(!found)ISSMERROR(" could not find pressure in inputs!");
578
579 found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
580 if(!found)ISSMERROR(" could not find temperature in inputs!");
581
582 /*Get dof list: */
583 GetDofList(&doflist[0],&numberofdofspernode);
584
585 //Compute pressure melting point
586 meltingpoint=matpar->GetMeltingPoint();
587 beta=matpar->GetBeta();
588 t_pmp=meltingpoint-beta*pressure;
589
590 //Add penalty load
591 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
592 Ke[0][0]=kmax*pow((double)10,penalty_offset);
593 }
594
595 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
596}
597/*}}}1*/
598/*FUNCTION Pengrid::PenaltyCreateKMatrixThermal {{{1*/
599void Pengrid::PenaltyCreateKMatrixThermal(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type){
600
601 int found=0;
602
603 const int numgrids=1;
604 const int NDOF1=1;
605 const int numdof=numgrids*NDOF1;
606 double Ke[numdof][numdof];
607 int doflist[numdof];
608 int numberofdofspernode;
609
610 ParameterInputs* inputs=NULL;
611
612 /*recover pointers: */
613 inputs=(ParameterInputs*)vinputs;
614
615
616 if(!active)return;
617
618 /*Get dof list: */
619 GetDofList(&doflist[0],&numberofdofspernode);
620
621 Ke[0][0]=kmax*pow((double)10,penalty_offset);
622
623 /*Add Ke to global matrix Kgg: */
624 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
625}
626/*}}}1*/
627/*FUNCTION Pengrid::PenaltyCreatePVector {{{1*/
628void Pengrid::PenaltyCreatePVector(Vec pg,double kmax,int analysis_type,int sub_analysis_type){
629
630 if (analysis_type==ThermalAnalysisEnum){
631
632 PenaltyCreatePVectorThermal( pg,kmax,analysis_type,sub_analysis_type);
633
634 }
635 else if (analysis_type==MeltingAnalysisEnum){
636
637 PenaltyCreatePVectorMelting( pg,kmax,analysis_type,sub_analysis_type);
638
639 }
640 else if (analysis_type==DiagnosticAnalysisEnum){
641
642 /*No loads applied, do nothing: */
643 return;
644
645 }
646 else{
647 ISSMERROR("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet");
648 }
649
650}
651/*}}}1*/
652/*FUNCTION Pengrid::PenaltyCreatePVectorMelting {{{1*/
653void Pengrid::PenaltyCreatePVectorMelting(Vec pg, double kmax,int analysis_type,int sub_analysis_type){
654
655 const int numgrids=1;
656 const int NDOF1=1;
657 const int numdof=numgrids*NDOF1;
658 int doflist[numdof];
659 double P_terms[numdof]={0.0};
660 int numberofdofspernode;
661 int found=0;
662 int dofs1[1]={0};
663 double pressure;
664 double temperature;
665 double melting_offset;
666 double meltingpoint;
667 double beta, heatcapacity;
668 double latentheat;
669 double t_pmp;
670 double dt;
671
672 ParameterInputs* inputs=NULL;
673
674 /*check that pengrid is not a clone (penalty to be added only once)*/
675 if (node->IsClone()) return;
676
677 /*recover pointers: */
678 inputs=(ParameterInputs*)vinputs;
679
680 /*Get dof list: */
681 GetDofList(&doflist[0],&numberofdofspernode);
682
683 //First recover pressure,melting offset and temperature vectors
684 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
685 if(!found)ISSMERROR(" could not find pressure in inputs!");
686
687 found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
688 if(!found)ISSMERROR(" could not find temperature in inputs!");
689
690 found=inputs->Recover("melting_offset",&melting_offset);
691 if(!found)ISSMERROR(" could not find melting_offset in inputs!");
692
693 found=inputs->Recover("dt",&dt);
694 if(!found)ISSMERROR(" could not find dt in inputs!");
695
696 meltingpoint=matpar->GetMeltingPoint();
697 beta=matpar->GetBeta();
698 heatcapacity=matpar->GetHeatCapacity();
699 latentheat=matpar->GetLatentHeat();
700
701 //Compute pressure melting point
702 t_pmp=meltingpoint-beta*pressure;
703
704 //Add penalty load
705 //This time, the penalty must have the same value as the one used for the thermal computation
706 //so that the corresponding melting can be computed correctly
707 //In the thermal computation, we used kmax=melting_offset, and the same penalty_offset
708 if (temperature<t_pmp){ //%no melting
709 P_terms[0]=0;
710 }
711 else{
712 if (dt){
713 P_terms[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp)/dt;
714 }
715 else{
716 P_terms[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp);
717 }
718 }
719
720 /*Add P_terms to global vector pg: */
721 VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
722}
723/*}}}1*/
724/*FUNCTION Pengrid::PenaltyCreatePVectorThermal {{{1*/
725void Pengrid::PenaltyCreatePVectorThermal(Vec pg, double kmax,int analysis_type,int sub_analysis_type){
726
727 const int numgrids=1;
728 const int NDOF1=1;
729 const int numdof=numgrids*NDOF1;
730 int doflist[numdof];
731 double P_terms[numdof]={0.0};
732 int numberofdofspernode;
733 int found=0;
734 double pressure;
735 int dofs1[1]={0};
736 double meltingpoint;
737 double beta;
738 double t_pmp;
739
740 ParameterInputs* inputs=NULL;
741
742 /*recover pointers: */
743 inputs=(ParameterInputs*)vinputs;
744
745 if(!active)return;
746
747 /*Get dof list: */
748 GetDofList(&doflist[0],&numberofdofspernode);
749
750 //First recover pressure
751 found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
752 if(!found)ISSMERROR(" could not find pressure in inputs!");
753
754 //Compute pressure melting point
755 meltingpoint=matpar->GetMeltingPoint();
756 beta=matpar->GetBeta();
757 t_pmp=meltingpoint-beta*pressure;
758
759 //Add penalty load
760 P_terms[0]=kmax*pow((double)10,penalty_offset)*t_pmp;
761
762 /*Add P_terms to global vector pg: */
763 VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
764}
765/*}}}1*/
766/*FUNCTION Pengrid::UpdateFromInputs {{{1*/
767void Pengrid::UpdateFromInputs(void* inputs){
768
769}
770/*}}}1*/
Note: See TracBrowser for help on using the repository browser.