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

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

Temporary
File size: 19.1 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_element_id,int pengrid_matpar_id):
37 hnode(&pengrid_node_id,1),
38 helement(&pengrid_element_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_helement,Hook* pengrid_hmatpar, Parameters* pengrid_parameters, Inputs* pengrid_inputs):
55 hnode(pengrid_hnode),
56 helement(pengrid_helement),
57 hmatpar(pengrid_hmatpar)
58{
59
60 /*all the initialization has been done by the initializer, just fill in the id: */
61 this->id=pengrid_id;
62 if(pengrid_inputs){
63 this->inputs=(Inputs*)pengrid_inputs->Copy();
64 }
65 else{
66 this->inputs=new Inputs();
67 }
68 /*point parameters: */
69 this->parameters=pengrid_parameters;
70
71 /*not active, not zigzagging: */
72 active=0;
73 zigzag_counter=0;
74
75}
76/*}}}*/
77/*FUNCTION Pengrid::Pengrid(int index, int id, IoModel* iomodel){{{1*/
78Pengrid::Pengrid(int id, int index, IoModel* iomodel){ //i is the element index
79
80 int i,j;
81 int pengrid_node_id;
82 int pengrid_matpar_id;
83 int pengrid_element_id;
84
85 /*id: */
86 this->id=id;
87
88 /*hooks: */
89 pengrid_node_id=index+1;
90 pengrid_element_id=iomodel->singlenodetoelementconnectivity[index];
91 pengrid_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object
92
93 this->hnode.Init(&pengrid_node_id,1);
94 this->helement.Init(&pengrid_element_id,1);
95 this->hmatpar.Init(&pengrid_matpar_id,1);
96
97 //initialize inputs: none needed
98 this->inputs=new Inputs();
99
100 //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
101 this->parameters=NULL;
102
103 //let's not forget internals
104 this->active=0;
105 this->zigzag_counter=0;
106
107}
108/*}}}*/
109/*FUNCTION Pengrid::destructor {{{1*/
110Pengrid::~Pengrid(){
111 return;
112}
113/*}}}1*/
114
115/*Object management*/
116/*FUNCTION Pengrid::Configure {{{1*/
117void Pengrid::Configure(DataSet* elementsin,DataSet* loadsin,DataSet* nodesin,DataSet* verticesin,DataSet* materialsin,Parameters* parametersin){
118
119 /*Take care of hooking up all objects for this load, ie links the objects in the hooks to their respective
120 * datasets, using internal ids and offsets hidden in hooks: */
121 hnode.configure(nodesin);
122 helement.configure(elementsin);
123 hmatpar.configure(materialsin);
124
125 /*point parameters to real dataset: */
126 this->parameters=parametersin;
127}
128/*}}}1*/
129/*FUNCTION Pengrid::copy {{{1*/
130Object* Pengrid::copy() {
131 return new Pengrid(this->id,&this->hnode,&this->helement,&this->hmatpar,this->parameters,this->inputs);
132}
133/*}}}1*/
134/*FUNCTION Pengrid::DeepEcho{{{1*/
135void Pengrid::DeepEcho(void){
136
137 printf("Pengrid:\n");
138 printf(" id: %i\n",id);
139 hnode.DeepEcho();
140 helement.DeepEcho();
141 hmatpar.DeepEcho();
142 printf(" active %i\n",this->active);
143 printf(" zigzag_counter %i\n",this->zigzag_counter);
144 printf(" parameters\n");
145 parameters->DeepEcho();
146 printf(" inputs\n");
147 inputs->DeepEcho();
148}
149/*}}}*/
150/*FUNCTION Pengrid::Demarshall {{{1*/
151void Pengrid::Demarshall(char** pmarshalled_dataset){
152
153 char* marshalled_dataset=NULL;
154 int i;
155
156 /*recover marshalled_dataset: */
157 marshalled_dataset=*pmarshalled_dataset;
158
159 /*this time, no need to get enum type, the pointer directly points to the beginning of the
160 *object data (thanks to DataSet::Demarshall):*/
161
162 memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
163 memcpy(&active,marshalled_dataset,sizeof(active));marshalled_dataset+=sizeof(active);
164 memcpy(&zigzag_counter,marshalled_dataset,sizeof(zigzag_counter));marshalled_dataset+=sizeof(zigzag_counter);
165
166 /*demarshall hooks: */
167 hnode.Demarshall(&marshalled_dataset);
168 helement.Demarshall(&marshalled_dataset);
169 hmatpar.Demarshall(&marshalled_dataset);
170
171 /*demarshall inputs: */
172 inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset);
173
174 /*parameters: may not exist even yet, so let Configure handle it: */
175 this->parameters=NULL;
176
177 /*return: */
178 *pmarshalled_dataset=marshalled_dataset;
179 return;
180}
181/*}}}*/
182/*FUNCTION Pengrid::Echo {{{1*/
183void Pengrid::Echo(void){
184 this->DeepEcho();
185}
186/*}}}1*/
187/*FUNCTION Pengrid::Enum {{{1*/
188int Pengrid::Enum(void){
189
190 return PengridEnum;
191}
192/*}}}1*/
193/*FUNCTION Pengrid::GetId {{{1*/
194int Pengrid::GetId(void){ return id; }
195/*}}}1*/
196/*FUNCTION Pengrid::GetName {{{1*/
197char* Pengrid::GetName(void){
198 return "pengrid";
199}
200/*}}}1*/
201/*FUNCTION Pengrid::Marshall {{{1*/
202void Pengrid::Marshall(char** pmarshalled_dataset){
203
204 char* marshalled_dataset=NULL;
205 int enum_type=0;
206 char* marshalled_inputs=NULL;
207 int marshalled_inputs_size;
208
209 /*recover marshalled_dataset: */
210 marshalled_dataset=*pmarshalled_dataset;
211
212 /*get enum type of Tria: */
213 enum_type=PengridEnum;
214
215 /*marshall enum: */
216 memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
217
218 /*marshall Tria data: */
219 memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
220 memcpy(marshalled_dataset,&active,sizeof(active));marshalled_dataset+=sizeof(active);
221 memcpy(marshalled_dataset,&zigzag_counter,sizeof(zigzag_counter));marshalled_dataset+=sizeof(zigzag_counter);
222
223 /*Marshall hooks: */
224 hnode.Marshall(&marshalled_dataset);
225 helement.Marshall(&marshalled_dataset);
226 hmatpar.Marshall(&marshalled_dataset);
227
228 /*Marshall inputs: */
229 marshalled_inputs_size=inputs->MarshallSize();
230 marshalled_inputs=inputs->Marshall();
231 memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char));
232 marshalled_dataset+=marshalled_inputs_size;
233
234 /*parameters: don't do anything about it. parameters are marshalled somewhere else!*/
235
236 xfree((void**)&marshalled_inputs);
237
238 *pmarshalled_dataset=marshalled_dataset;
239 return;
240}
241/*}}}*/
242/*FUNCTION Pengrid::MarshallSize {{{1*/
243int Pengrid::MarshallSize(){
244
245 return sizeof(id)
246 +sizeof(active)
247 +sizeof(zigzag_counter)
248 +hnode.MarshallSize()
249 +helement.MarshallSize()
250 +hmatpar.MarshallSize()
251 +inputs->MarshallSize()
252 +sizeof(int); //sizeof(int) for enum type
253}
254/*}}}*/
255/*FUNCTION Pengrid::MyRank {{{1*/
256int Pengrid::MyRank(void){
257 extern int my_rank;
258 return my_rank;
259}
260/*}}}1*/
261
262/*Object functions*/
263/*FUNCTION Pengrid::CreateKMatrix {{{1*/
264
265void Pengrid::CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type){
266
267 /*No loads applied, do nothing: */
268 return;
269
270}
271/*}}}1*/
272/*FUNCTION Pengrid::CreatePVector {{{1*/
273void Pengrid::CreatePVector(Vec pg, int analysis_type,int sub_analysis_type){
274
275 /*No loads applied, do nothing: */
276 return;
277
278}
279/*}}}1*/
280/*FUNCTION Pengrid::DistributenumDofs {{{1*/
281void Pengrid::DistributeNumDofs(int* numdofpernode,int analysis_type,int sub_analysis_type){return;}
282/*}}}1*/
283/*FUNCTION Pengrid::GetDofList {{{1*/
284void Pengrid::GetDofList(int* doflist,int* pnumberofdofspernode){
285
286 int i,j;
287 int doflist_per_node[MAXDOFSPERNODE];
288 int numberofdofspernode;
289
290 /*dynamic objects pointed to by hooks: */
291 Node* node=NULL;
292
293 /*recover objects from hooks: */
294 node=(Node*)hnode.delivers();
295
296 node->GetDofList(&doflist_per_node[0],&numberofdofspernode);
297 for(j=0;j<numberofdofspernode;j++){
298 doflist[j]=doflist_per_node[j];
299 }
300
301 /*Assign output pointers:*/
302 *pnumberofdofspernode=numberofdofspernode;
303
304}
305/*}}}*/
306/*FUNCTION Pengrid::PenaltyConstrain {{{1*/
307void Pengrid::PenaltyConstrain(int* punstable,int analysis_type,int sub_analysis_type){
308
309 if ((analysis_type==DiagnosticAnalysisEnum) && ((sub_analysis_type==StokesAnalysisEnum))){
310
311 /*No penalty to check*/
312 return;
313
314 }
315 else if (analysis_type==ThermalAnalysisEnum){
316
317 PenaltyConstrainThermal(punstable,analysis_type,sub_analysis_type);
318
319 }
320 else if (analysis_type==MeltingAnalysisEnum){
321
322 /*No penalty to check*/
323 return;
324
325 }
326 else{
327 ISSMERROR("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet");
328 }
329
330}
331/*}}}1*/
332/*FUNCTION Pengrid::PenaltyConstrainThermal {{{1*/
333void Pengrid::PenaltyConstrainThermal(int* punstable,int analysis_type,int sub_analysis_type){
334
335 // The penalty is stable if it doesn't change during to successive iterations.
336
337 int found=0;
338 const int numgrids=1;
339 double pressure;
340 double temperature;
341 double beta,t_pmp;
342 double meltingpoint;
343 int new_active;
344 int unstable=0;
345 int reset_penalties=0;
346 int stabilize_constraints;
347
348 /*pointers: */
349 Node* node=NULL;
350 Penta* penta=NULL;
351 Matpar* matpar=NULL;
352
353 /*check that pengrid is not a clone (penalty to be added only once)*/
354 if (node->IsClone()){
355 unstable=0;
356 *punstable=unstable;
357 return;
358 }
359
360 /*recover pointers: */
361 node=(Node*)hnode.delivers();
362 penta=(Penta*)helement.delivers();
363 matpar=(Matpar*)hmatpar.delivers();
364
365 //First recover pressure and temperature values, using the element: */
366 penta->inputs->GetParameterValueAtNode(&pressure,node,PressureEnum);
367 penta->inputs->GetParameterValueAtNode(&temperature,node,TemperatureEnum);
368
369 //Recover our data:
370 inputs->GetParameterValue(&reset_penalties,ResetPenaltiesEnum);
371 meltingpoint=matpar->GetMeltingPoint();
372 beta=matpar->GetBeta();
373 parameters->FindParam(&stabilize_constraints,"stabilize_constraints");
374
375 if(reset_penalties)zigzag_counter=0;
376
377 //Compute pressure melting point
378 t_pmp=meltingpoint-beta*pressure;
379
380 //Figure out if temperature is over melting_point, in which case, this penalty needs to be activated.
381
382 if (temperature>t_pmp){
383 new_active=1;
384 }
385 else{
386 new_active=0;
387 }
388
389
390 //Figure out stability of this penalty
391 if (active==new_active){
392 unstable=0;
393 }
394 else{
395 unstable=1;
396 if(stabilize_constraints)zigzag_counter++;
397 }
398
399 /*If penalty keeps zigzagging more than 5 times: */
400 if(stabilize_constraints){
401 if(zigzag_counter>stabilize_constraints){
402 unstable=0;
403 active=1;
404 }
405 }
406
407 //Set penalty flag
408 active=new_active;
409
410 //*Assign output pointers:*/
411 *punstable=unstable;
412}
413/*}}}1*/
414/*FUNCTION Pengrid::PenaltyCreateMatrix {{{1*/
415void Pengrid::PenaltyCreateKMatrix(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type){
416
417 if ((analysis_type==DiagnosticAnalysisEnum) && ((sub_analysis_type==StokesAnalysisEnum))){
418
419 PenaltyCreateKMatrixDiagnosticStokes( Kgg,kmax,analysis_type,sub_analysis_type);
420 }
421 else if (analysis_type==ThermalAnalysisEnum){
422
423 PenaltyCreateKMatrixThermal( Kgg,kmax,analysis_type,sub_analysis_type);
424
425 }
426 else if (analysis_type==MeltingAnalysisEnum){
427
428 PenaltyCreateKMatrixMelting( Kgg,kmax,analysis_type,sub_analysis_type);
429
430 }
431 else{
432 ISSMERROR("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet");
433 }
434
435}
436/*}}}1*/
437/*FUNCTION Pengrid::PenaltyCreateKMatrixDiagnosticStokes {{{1*/
438void Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type){
439
440 const int numgrids=1;
441 const int NDOF4=4;
442 const int numdof=numgrids*NDOF4;
443 int doflist[numdof];
444 int numberofdofspernode;
445
446 int dofs1[1]={0};
447 int dofs2[1]={1};
448 double slope[2];
449 int found=0;
450 double Ke[4][4]={0.0};
451 double penalty_offset;
452
453 /*pointers: */
454 Node* node=NULL;
455 Penta* penta=NULL;
456
457 /*Get dof list: */
458 GetDofList(&doflist[0],&numberofdofspernode);
459
460 /*recover pointers: */
461 node=(Node*)hnode.delivers();
462 penta=(Penta*)helement.delivers();
463
464 //recover slope: */
465 penta->inputs->GetParameterValueAtNode(&slope[0],node,BedSlopexEnum);
466 penta->inputs->GetParameterValueAtNode(&slope[1],node,BedSlopeyEnum);
467
468 /*recover parameters: */
469 parameters->FindParam(&penalty_offset,"penalty_offset");
470
471 //Create elementary matrix: add penalty to contrain wb (wb=ub*db/dx+vb*db/dy)
472 Ke[2][0]=-slope[0]*kmax*pow((double)10.0,penalty_offset);
473 Ke[2][1]=-slope[1]*kmax*pow((double)10.0,penalty_offset);
474 Ke[2][2]=kmax*pow((double)10,penalty_offset);
475
476 /*Add Ke to global matrix Kgg: */
477 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
478}
479/*}}}1*/
480/*FUNCTION Pengrid::PenaltyCreateKMatrixMelting {{{1*/
481void Pengrid::PenaltyCreateKMatrixMelting(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type){
482
483
484 int found=0;
485 const int numgrids=1;
486 const int NDOF1=1;
487 const int numdof=numgrids*NDOF1;
488 double Ke[numdof][numdof]={0.0};
489 int dofs1[1]={0};
490 int doflist[numdof];
491 int numberofdofspernode;
492 double meltingpoint;
493
494 double pressure;
495 double temperature;
496 double beta,t_pmp;
497 double penalty_offset;
498
499 /*pointers: */
500 Node* node=NULL;
501 Penta* penta=NULL;
502 Matpar* matpar=NULL;
503
504 /*recover pointers: */
505 node=(Node*)hnode.delivers();
506 penta=(Penta*)helement.delivers();
507 matpar=(Matpar*)hmatpar.delivers();
508
509 /*check that pengrid is not a clone (penalty to be added only once)*/
510 if (node->IsClone()) return;
511
512 //First recover pressure and temperature values, using the element: */
513 penta->inputs->GetParameterValueAtNode(&pressure,node,PressureEnum);
514 penta->inputs->GetParameterValueAtNode(&temperature,node,TemperatureEnum);
515
516 /*recover parameters: */
517 parameters->FindParam(&penalty_offset,"penalty_offset");
518
519 /*Get dof list: */
520 GetDofList(&doflist[0],&numberofdofspernode);
521
522 //Compute pressure melting point
523 meltingpoint=matpar->GetMeltingPoint();
524 beta=matpar->GetBeta();
525 t_pmp=meltingpoint-beta*pressure;
526
527 //Add penalty load
528 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
529 Ke[0][0]=kmax*pow((double)10,penalty_offset);
530 }
531
532 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
533}
534/*}}}1*/
535/*FUNCTION Pengrid::PenaltyCreateKMatrixThermal {{{1*/
536void Pengrid::PenaltyCreateKMatrixThermal(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type){
537
538 int found=0;
539
540 const int numgrids=1;
541 const int NDOF1=1;
542 const int numdof=numgrids*NDOF1;
543 double Ke[numdof][numdof];
544 int doflist[numdof];
545 int numberofdofspernode;
546 double penalty_offset;
547
548 if(!this->active)return;
549
550 /*recover parameters: */
551 parameters->FindParam(&penalty_offset,"penalty_offset");
552
553 /*Get dof list: */
554 GetDofList(&doflist[0],&numberofdofspernode);
555
556 Ke[0][0]=kmax*pow((double)10,penalty_offset);
557
558 /*Add Ke to global matrix Kgg: */
559 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
560}
561/*}}}1*/
562/*FUNCTION Pengrid::PenaltyCreatePVector {{{1*/
563void Pengrid::PenaltyCreatePVector(Vec pg,double kmax,int analysis_type,int sub_analysis_type){
564
565 if (analysis_type==ThermalAnalysisEnum){
566
567 PenaltyCreatePVectorThermal( pg,kmax,analysis_type,sub_analysis_type);
568
569 }
570 else if (analysis_type==MeltingAnalysisEnum){
571
572 PenaltyCreatePVectorMelting( pg,kmax,analysis_type,sub_analysis_type);
573
574 }
575 else if (analysis_type==DiagnosticAnalysisEnum){
576
577 /*No loads applied, do nothing: */
578 return;
579
580 }
581 else{
582 ISSMERROR("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet");
583 }
584
585}
586/*}}}1*/
587/*FUNCTION Pengrid::PenaltyCreatePVectorMelting {{{1*/
588void Pengrid::PenaltyCreatePVectorMelting(Vec pg, double kmax,int analysis_type,int sub_analysis_type){
589
590 const int numgrids=1;
591 const int NDOF1=1;
592 const int numdof=numgrids*NDOF1;
593 int doflist[numdof];
594 double P_terms[numdof]={0.0};
595 int numberofdofspernode;
596 int found=0;
597 int dofs1[1]={0};
598 double pressure;
599 double temperature;
600 double melting_offset;
601 double meltingpoint;
602 double beta, heatcapacity;
603 double latentheat;
604 double t_pmp;
605 double dt,penalty_offset;
606
607 /*pointers: */
608 Node* node=NULL;
609 Penta* penta=NULL;
610 Matpar* matpar=NULL;
611
612 /*recover pointers: */
613 node=(Node*)hnode.delivers();
614 penta=(Penta*)helement.delivers();
615 matpar=(Matpar*)hmatpar.delivers();
616
617 /*check that pengrid is not a clone (penalty to be added only once)*/
618 if (node->IsClone()) return;
619
620 /*Get dof list: */
621 GetDofList(&doflist[0],&numberofdofspernode);
622
623 //First recover pressure and temperature values, using the element: */
624 penta->inputs->GetParameterValueAtNode(&pressure,node,PressureEnum);
625 penta->inputs->GetParameterValueAtNode(&temperature,node,TemperatureEnum);
626 inputs->GetParameterValue(&melting_offset,MeltingOffsetEnum);
627 parameters->FindParam(&dt,"dt");
628 parameters->FindParam(&penalty_offset,"penalty_offset");
629
630 meltingpoint=matpar->GetMeltingPoint();
631 beta=matpar->GetBeta();
632 heatcapacity=matpar->GetHeatCapacity();
633 latentheat=matpar->GetLatentHeat();
634
635 //Compute pressure melting point
636 t_pmp=meltingpoint-beta*pressure;
637
638 //Add penalty load
639 //This time, the penalty must have the same value as the one used for the thermal computation
640 //so that the corresponding melting can be computed correctly
641 //In the thermal computation, we used kmax=melting_offset, and the same penalty_offset
642 if (temperature<t_pmp){ //%no melting
643 P_terms[0]=0;
644 }
645 else{
646 if (dt){
647 P_terms[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp)/dt;
648 }
649 else{
650 P_terms[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp);
651 }
652 }
653
654 /*Add P_terms to global vector pg: */
655 VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
656}
657/*}}}1*/
658/*FUNCTION Pengrid::PenaltyCreatePVectorThermal {{{1*/
659void Pengrid::PenaltyCreatePVectorThermal(Vec pg, double kmax,int analysis_type,int sub_analysis_type){
660
661 const int numgrids=1;
662 const int NDOF1=1;
663 const int numdof=numgrids*NDOF1;
664 int doflist[numdof];
665 double P_terms[numdof]={0.0};
666 int numberofdofspernode;
667 int found=0;
668 double pressure;
669 int dofs1[1]={0};
670 double meltingpoint;
671 double beta;
672 double t_pmp;
673 double penalty_offset;
674
675 /*pointers: */
676 Node* node=NULL;
677 Penta* penta=NULL;
678 Matpar* matpar=NULL;
679
680 /*recover pointers: */
681 node=(Node*)hnode.delivers();
682 penta=(Penta*)helement.delivers();
683 matpar=(Matpar*)hmatpar.delivers();
684
685 if(!this->active)return;
686
687 /*Get dof list: */
688 GetDofList(&doflist[0],&numberofdofspernode);
689
690 //First recover pressure and penalty_offset
691 penta->inputs->GetParameterValueAtNode(&pressure,node,PressureEnum);
692 parameters->FindParam(&penalty_offset,"penalty_offset");
693
694 //Compute pressure melting point
695 meltingpoint=matpar->GetMeltingPoint();
696 beta=matpar->GetBeta();
697 t_pmp=meltingpoint-beta*pressure;
698
699 //Add penalty load
700 P_terms[0]=kmax*pow((double)10,penalty_offset)*t_pmp;
701
702 /*Add P_terms to global vector pg: */
703 VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
704}
705/*}}}1*/
706
707/*Updates: */
708/*FUNCTION Pengrid::UpdateFromDakota {{{1*/
709void Pengrid::UpdateFromDakota(void* inputs){
710 ISSMERROR("not supported yet!");
711}
712/*}}}1*/
713/*FUNCTION Pengrid::UpdateInputs {{{1*/
714void Pengrid::UpdateInputs(double* solution, int analysis_type, int sub_analysis_type){
715 ISSMERROR("not supported yet!");
716}
717/*}}}1*/
Note: See TracBrowser for help on using the repository browser.