Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5935) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5936) @@ -688,6 +688,6 @@ /*retrieve parameters: */ + ElementMatrix* Ke=NULL; int analysis_type; - ElementMatrix* Ke=NULL; parameters->FindParam(&analysis_type,AnalysisTypeEnum); Index: /issm/trunk/src/c/objects/Loads/Pengrid.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Pengrid.cpp (revision 5935) +++ /issm/trunk/src/c/objects/Loads/Pengrid.cpp (revision 5936) @@ -3,5 +3,6 @@ */ - +/*Headers*/ +/*{{{1*/ #ifdef HAVE_CONFIG_H #include "config.h" @@ -17,5 +18,13 @@ #include "../../shared/shared.h" #include "../../Container/Container.h" - +/*}}}*/ + +/*Element macros*/ +#define NUMVERTICES 1 +#define NDOF1 1 +#define NDOF2 2 +#define NDOF3 3 +#define NDOF4 4 + /*Pengrid constructors and destructor*/ /*FUNCTION Pengrid::Pengrid(){{{1*/ @@ -299,22 +308,28 @@ void Pengrid::PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs,double kmax){ + /*Retrieve parameters: */ + ElementMatrix* Ke=NULL; int analysis_type; - - /*Retrieve parameters: */ this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); - if (analysis_type==DiagnosticHorizAnalysisEnum || analysis_type==AdjointHorizAnalysisEnum){ - PenaltyCreateKMatrixDiagnosticStokes( Kgg,kmax); - } - else if (analysis_type==ThermalAnalysisEnum){ - PenaltyCreateKMatrixThermal( Kgg,kmax); - } - else if (analysis_type==MeltingAnalysisEnum){ - PenaltyCreateKMatrixMelting( Kgg,kmax); - } - else{ - ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); - } - + switch(analysis_type){ + case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum: + Ke=PenaltyCreateKMatrixDiagnosticStokes(kmax); + break; + case ThermalAnalysisEnum: + Ke=PenaltyCreateKMatrixThermal(kmax); + break; + case MeltingAnalysisEnum: + Ke=PenaltyCreateKMatrixMelting(kmax); + break; + default: + ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); + } + + /*Add to global matrix*/ + if(Ke){ + Ke->AddToGlobal(Kgg,Kff,Kfs); + delete Ke; + } } /*}}}1*/ @@ -322,25 +337,27 @@ void Pengrid::PenaltyCreatePVector(Vec pg,Vec pf,double kmax){ + /*Retrieve parameters: */ + ElementVector* pe=NULL; int analysis_type; - - /*Retrieve parameters: */ this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); - if (analysis_type==ThermalAnalysisEnum){ - PenaltyCreatePVectorThermal( pg,kmax); - } - else if (analysis_type==MeltingAnalysisEnum){ - PenaltyCreatePVectorMelting( pg,kmax); - } - else if (analysis_type==DiagnosticHorizAnalysisEnum || analysis_type==AdjointHorizAnalysisEnum){ - - /*No loads applied, do nothing: */ - return; - - } - else{ - ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); - } - + switch(analysis_type){ + case ThermalAnalysisEnum: + pe=PenaltyCreatePVectorThermal(kmax); + break; + case MeltingAnalysisEnum: + pe=PenaltyCreatePVectorMelting(kmax); + break; + case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum: + break; + default: + ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); + } + + /*Add to global Vector*/ + if(pe){ + pe->AddToGlobal(pg,pf); + delete pe; + } } /*}}}1*/ @@ -541,146 +558,88 @@ /*}}}1*/ /*FUNCTION Pengrid::PenaltyCreateKMatrixDiagnosticStokes {{{1*/ -void Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax){ - - const int numgrids=1; - const int NDOF4=4; - const int numdof=numgrids*NDOF4; - int* doflist=NULL; - - int dofs1[1]={0}; - int dofs2[1]={1}; - double slope[2]; - int found=0; - double Ke[4][4]={0.0}; - double penalty_offset; - int approximation; - - /*pointers: */ - Penta* penta=NULL; - - /*recover pointers: */ - penta=(Penta*)element; - - /*If not Stokes, return*/ +ElementMatrix* Pengrid::PenaltyCreateKMatrixDiagnosticStokes(double kmax){ + + const int numdof = NUMVERTICES *NDOF4; + double slope[2]; + double penalty_offset; + int approximation; + + Penta* penta=(Penta*)element; + + /*Initialize Element vector and return if necessary*/ penta->inputs->GetParameterValue(&approximation,ApproximationEnum); - if(approximation!=StokesApproximationEnum) return; - - /*Get dof list: */ - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - //recover slope: */ + if(approximation!=StokesApproximationEnum) return NULL; + ElementMatrix* Ke=NewElementMatrix(&node,1,this->parameters,StokesApproximationEnum); + + /*Retrieve all inputs and parameters*/ + parameters->FindParam(&penalty_offset,PenaltyOffsetEnum); penta->GetParameterValue(&slope[0],node,BedSlopeXEnum); penta->GetParameterValue(&slope[1],node,BedSlopeYEnum); - + + /*Create elementary matrix: add penalty to constrain wb (wb=ub*db/dx+vb*db/dy)*/ + Ke->values[2*NDOF4+0]=-slope[0]*kmax*pow((double)10.0,penalty_offset); + Ke->values[2*NDOF4+1]=-slope[1]*kmax*pow((double)10.0,penalty_offset); + Ke->values[2*NDOF4+2]= kmax*pow((double)10,penalty_offset); + + /*Clean up and return*/ + return Ke; +} +/*}}}1*/ +/*FUNCTION Pengrid::PenaltyCreateKMatrixMelting {{{1*/ +ElementMatrix* Pengrid::PenaltyCreateKMatrixMelting(double kmax){ + + const int numdof=NUMVERTICES*NDOF1; + double pressure,temperature,t_pmp; + double penalty_offset; + + Penta* penta=(Penta*)element; + + /*check that pengrid is not a clone (penalty to be added only once)*/ + if (node->IsClone()) return NULL; + ElementMatrix* Ke=NewElementMatrix(&node,1,this->parameters); + + /*Retrieve all inputs and parameters*/ + penta->GetParameterValue(&pressure,node,PressureEnum); + penta->GetParameterValue(&temperature,node,TemperatureEnum); + parameters->FindParam(&penalty_offset,PenaltyOffsetEnum); + + /*Compute pressure melting point*/ + t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure; + + /*Add penalty load*/ + if (temperaturevalues[0]=kmax*pow((double)10,penalty_offset); + } + + /*Clean up and return*/ + return Ke; +} +/*}}}1*/ +/*FUNCTION Pengrid::PenaltyCreateKMatrixThermal {{{1*/ +ElementMatrix* Pengrid::PenaltyCreateKMatrixThermal(double kmax){ + + const int numdof=NUMVERTICES*NDOF1; + double penalty_offset; + + /*Initialize Element matrix and return if necessary*/ + if(!this->active) return NULL; + ElementMatrix* Ke=NewElementMatrix(&node,NUMVERTICES,this->parameters); + /*recover parameters: */ parameters->FindParam(&penalty_offset,PenaltyOffsetEnum); - //Create elementary matrix: add penalty to constrain wb (wb=ub*db/dx+vb*db/dy) - Ke[2][0]=-slope[0]*kmax*pow((double)10.0,penalty_offset); - Ke[2][1]=-slope[1]*kmax*pow((double)10.0,penalty_offset); - Ke[2][2]=kmax*pow((double)10,penalty_offset); - - /*Add Ke to global matrix Kgg: */ - MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES); - - /*Free ressources:*/ - xfree((void**)&doflist); - -} -/*}}}1*/ -/*FUNCTION Pengrid::PenaltyCreateKMatrixMelting {{{1*/ -void Pengrid::PenaltyCreateKMatrixMelting(Mat Kgg,double kmax){ - - int found=0; - const int numgrids=1; - const int NDOF1=1; - const int numdof=numgrids*NDOF1; - double Ke[numdof][numdof]={0.0}; - int dofs1[1]={0}; - int* doflist=NULL; - double meltingpoint; - - double pressure; - double temperature; - double beta,t_pmp; - double penalty_offset; - - /*recover pointers: */ - Penta* penta=(Penta*)element; - - /*check that pengrid is not a clone (penalty to be added only once)*/ - if (node->IsClone()) return; - - //First recover pressure and temperature values, using the element: */ - penta->GetParameterValue(&pressure,node,PressureEnum); - penta->GetParameterValue(&temperature,node,TemperatureEnum); - - /*recover parameters: */ - parameters->FindParam(&penalty_offset,PenaltyOffsetEnum); - - /*Get dof list: */ - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - //Compute pressure melting point - meltingpoint=matpar->GetMeltingPoint(); - beta=matpar->GetBeta(); - t_pmp=meltingpoint-beta*pressure; - - //Add penalty load - if (temperatureactive)return; - - /*recover parameters: */ - parameters->FindParam(&penalty_offset,PenaltyOffsetEnum); - - /*Get dof list: */ - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - Ke[0][0]=kmax*pow((double)10,penalty_offset); - - /*Add Ke to global matrix Kgg: */ - MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES); - - /*Free ressources:*/ - xfree((void**)&doflist); + Ke->values[0]=kmax*pow((double)10,penalty_offset); + + /*Clean up and return*/ + return Ke; } /*}}}1*/ /*FUNCTION Pengrid::PenaltyCreatePVectorMelting {{{1*/ -void Pengrid::PenaltyCreatePVectorMelting(Vec pg, double kmax){ - - const int numgrids=1; - const int NDOF1=1; - const int numdof=numgrids*NDOF1; - int* doflist=NULL; - double P_terms[numdof]={0.0}; - int found=0; - int dofs1[1]={0}; +ElementVector* Pengrid::PenaltyCreatePVectorMelting(double kmax){ + + const int numdof=NUMVERTICES*NDOF1; double pressure; double temperature; double melting_offset; - double meltingpoint; - double beta, heatcapacity; - double latentheat; double t_pmp; double dt,penalty_offset; @@ -690,10 +649,8 @@ /*check that pengrid is not a clone (penalty to be added only once)*/ - if (node->IsClone()) return; - - /*Get dof list: */ - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - //First recover pressure and temperature values, using the element: */ + if (node->IsClone()) return NULL; + ElementVector* pe=NewElementVector(&node,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ penta->GetParameterValue(&pressure,node,PressureEnum); penta->GetParameterValue(&temperature,node,TemperatureEnum); @@ -702,78 +659,48 @@ parameters->FindParam(&penalty_offset,PenaltyOffsetEnum); - meltingpoint=matpar->GetMeltingPoint(); - beta=matpar->GetBeta(); - heatcapacity=matpar->GetHeatCapacity(); - latentheat=matpar->GetLatentHeat(); - - //Compute pressure melting point - t_pmp=meltingpoint-beta*pressure; - - //Add penalty load - //This time, the penalty must have the same value as the one used for the thermal computation - //so that the corresponding melting can be computed correctly - //In the thermal computation, we used kmax=melting_offset, and the same penalty_offset + /*Compute pressure melting point*/ + t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure; + + /*Add penalty load + This time, the penalty must have the same value as the one used for the thermal computation + so that the corresponding melting can be computed correctly + In the thermal computation, we used kmax=melting_offset, and the same penalty_offset*/ if (temperaturevalues[0]=0; } else{ - if (dt){ - P_terms[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp)/dt; - } - else{ - P_terms[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp); - } - } - - /*Add P_terms to global vector pg: */ - VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES); - - /*Clean up*/ - xfree((void**)&doflist); - + if (dt) pe->values[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp)/dt; + else pe->values[0]=melting_offset*pow((double)10,penalty_offset)*(temperature-t_pmp); + } + + /*Clean up and return*/ + return pe; } /*}}}1*/ /*FUNCTION Pengrid::PenaltyCreatePVectorThermal {{{1*/ -void Pengrid::PenaltyCreatePVectorThermal(Vec pg, double kmax){ - - const int numgrids=1; - const int NDOF1=1; - const int numdof=numgrids*NDOF1; - int* doflist=NULL; - double P_terms[numdof]={0.0}; - int found=0; +ElementVector* Pengrid::PenaltyCreatePVectorThermal(double kmax){ + + const int numdof=NUMVERTICES*NDOF1; double pressure; - int dofs1[1]={0}; - double meltingpoint; - double beta; double t_pmp; double penalty_offset; - /*recover pointers: */ Penta* penta=(Penta*)element; - if(!this->active)return; - - /*Get dof list: */ - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - //First recover pressure and penalty_offset + /*Initialize Element matrix and return if necessary*/ + if(!this->active) return NULL; + ElementVector* pe=NewElementVector(&node,1,this->parameters); + + /*Retrieve all inputs and parameters*/ penta->GetParameterValue(&pressure,node,PressureEnum); parameters->FindParam(&penalty_offset,PenaltyOffsetEnum); - //Compute pressure melting point - meltingpoint=matpar->GetMeltingPoint(); - beta=matpar->GetBeta(); - t_pmp=meltingpoint-beta*pressure; - - //Add penalty load - P_terms[0]=kmax*pow((double)10,penalty_offset)*t_pmp; - - /*Add P_terms to global vector pg: */ - VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES); - - /*Clean up*/ - xfree((void**)&doflist); - + /*Compute pressure melting point*/ + t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure; + + pe->values[0]=kmax*pow((double)10,penalty_offset)*t_pmp; + + /*Clean up and return*/ + return pe; } /*}}}1*/ Index: /issm/trunk/src/c/objects/Loads/Pengrid.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Pengrid.h (revision 5935) +++ /issm/trunk/src/c/objects/Loads/Pengrid.h (revision 5936) @@ -78,10 +78,10 @@ /*}}}*/ /*Pengrid management {{{1*/ - void PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax); void GetDofList(int** pdoflist,int approximation_enum,int setenum); - void PenaltyCreateKMatrixThermal(Mat Kgg,double kmax); - void PenaltyCreateKMatrixMelting(Mat Kgg,double kmax); - void PenaltyCreatePVectorThermal(Vec pg, double kmax); - void PenaltyCreatePVectorMelting(Vec pg, double kmax); + ElementMatrix* PenaltyCreateKMatrixDiagnosticStokes(double kmax); + ElementMatrix* PenaltyCreateKMatrixThermal(double kmax); + ElementMatrix* PenaltyCreateKMatrixMelting(double kmax); + ElementVector* PenaltyCreatePVectorThermal(double kmax); + ElementVector* PenaltyCreatePVectorMelting(double kmax); void PenaltyConstrain(int* punstable); void PenaltyConstrainThermal(int* punstable);