Index: /issm/trunk/src/c/objects/Loads/Pengrid.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Pengrid.cpp (revision 5936) +++ /issm/trunk/src/c/objects/Loads/Pengrid.cpp (revision 5937) @@ -435,29 +435,4 @@ /*Pengrid management:*/ -/*FUNCTION Pengrid::GetDofList {{{1*/ -void Pengrid::GetDofList(int** pdoflist,int approximation,int setenum){ - - int i,j; - int numberofdofs=0; - - /*output: */ - int* doflist=NULL; - - /*Some checks for debugging*/ - ISSMASSERT(node); - - /*First, figure out size of doflist: */ - numberofdofs=node->GetNumberOfDofs(approximation,setenum); - - /*Allocate: */ - doflist=(int*)xmalloc(numberofdofs*sizeof(int)); - - /*Populate: */ - node->GetDofList(doflist,approximation,setenum); - - /*Assign output pointers:*/ - *pdoflist=doflist; -} -/*}}}*/ /*FUNCTION Pengrid::PenaltyConstrain {{{1*/ void Pengrid::PenaltyConstrain(int* punstable){ Index: /issm/trunk/src/c/objects/Loads/Pengrid.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Pengrid.h (revision 5936) +++ /issm/trunk/src/c/objects/Loads/Pengrid.h (revision 5937) @@ -78,5 +78,4 @@ /*}}}*/ /*Pengrid management {{{1*/ - void GetDofList(int** pdoflist,int approximation_enum,int setenum); ElementMatrix* PenaltyCreateKMatrixDiagnosticStokes(double kmax); ElementMatrix* PenaltyCreateKMatrixThermal(double kmax); Index: /issm/trunk/src/c/objects/Loads/Penpair.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Penpair.cpp (revision 5936) +++ /issm/trunk/src/c/objects/Loads/Penpair.cpp (revision 5937) @@ -3,7 +3,8 @@ */ - +/*Headers*/ +/*{{{1*/ #ifdef HAVE_CONFIG_H - #include "config.h" +#include "config.h" #else #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" @@ -16,5 +17,13 @@ #include "../../include/include.h" #include "../../shared/shared.h" - +/*}}}*/ + +/*Element macros*/ +#define NUMVERTICES 2 +#define NDOF1 1 +#define NDOF2 2 +#define NDOF3 3 +#define NDOF4 4 + /*Penpair constructors and destructor*/ /*FUNCTION Penpair::constructor {{{1*/ @@ -201,50 +210,22 @@ /*FUNCTION Penpair::PenaltyCreateKMatrix {{{1*/ void Penpair::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){ - Node** nodes=NULL; - int approximation0,approximation1; - nodes=(Node**)hnodes->deliverp(); - approximation0=nodes[0]->GetApproximation(); - approximation1=nodes[1]->GetApproximation(); - - switch(approximation0){ - case MacAyealApproximationEnum: - switch(approximation1){ - case MacAyealApproximationEnum: PenaltyCreateKMatrixDiagnosticHoriz(Kgg,kmax); break; - case PattynApproximationEnum: PenaltyCreateKMatrixDiagnosticHoriz(Kgg,kmax); break; - default: ISSMERROR("not supported yet"); - } - break; - case PattynApproximationEnum: - switch(approximation1){ - case MacAyealApproximationEnum: PenaltyCreateKMatrixDiagnosticHoriz(Kgg,kmax); break; - case PattynApproximationEnum: PenaltyCreateKMatrixDiagnosticHoriz(Kgg,kmax); break; - default: ISSMERROR("not supported yet"); - } - break; - case StokesApproximationEnum: - switch(approximation1){ - case StokesApproximationEnum: PenaltyCreateKMatrixDiagnosticStokes(Kgg,kmax); break; - case NoneApproximationEnum: PenaltyCreateKMatrixDiagnosticStokes(Kgg,kmax); break; - default: ISSMERROR("not supported yet"); - } - break; - case NoneApproximationEnum: - switch(approximation1){ - case StokesApproximationEnum: PenaltyCreateKMatrixDiagnosticStokes(Kgg,kmax); break; - case NoneApproximationEnum: PenaltyCreateKMatrixDiagnosticStokes(Kgg,kmax); break; - default: ISSMERROR("not supported yet"); - } - break; - default: ISSMERROR("not supported yet"); - } + switch(analysis_type){ + case DiagnosticHorizAnalysisEnum: + Ke=PenaltyCreateKMatrixDiagnosticHoriz(kmax); + break; + default: + ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); } - else{ - ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); + + /*Add to global Vector*/ + if(Ke){ + Ke->AddToGlobal(Kgg,Kff,Kfs); + delete Ke; } } @@ -281,24 +262,51 @@ /*Penpair management:*/ -/*FUNCTION Penpair::PenaltyCreateKMatrixDiagnosticHoriz {{{1*/ -void Penpair::PenaltyCreateKMatrixDiagnosticHoriz(Mat Kgg,double kmax){ - - const int numgrids=2; - const int NDOF2=2; - const int numdof=numgrids*NDOF2; - int* doflist=NULL; - int numberofdofspernode; - - double Ke[4][4]={0.0}; +/*FUNCTION Penpair::PenaltyCreateKMatrixDiagnosticHoriz{{{1*/ +ElementMatrix* Penpair::PenaltyCreateKMatrixDiagnosticHoriz(double kmax){ + + Node** nodes=(Node**)hnodes->deliverp(); + int approximation0=nodes[0]->GetApproximation(); + int approximation1=nodes[1]->GetApproximation(); + + switch(approximation0){ + case MacAyealApproximationEnum: + switch(approximation1){ + case MacAyealApproximationEnum: return PenaltyCreateKMatrixDiagnosticMacAyealPattyn(kmax); + case PattynApproximationEnum: return PenaltyCreateKMatrixDiagnosticMacAyealPattyn(kmax); + default: ISSMERROR("not supported yet"); + } + case PattynApproximationEnum: + switch(approximation1){ + case MacAyealApproximationEnum: return PenaltyCreateKMatrixDiagnosticMacAyealPattyn(kmax); + case PattynApproximationEnum: return PenaltyCreateKMatrixDiagnosticMacAyealPattyn(kmax); + default: ISSMERROR("not supported yet"); + } + case StokesApproximationEnum: + switch(approximation1){ + case StokesApproximationEnum: return PenaltyCreateKMatrixDiagnosticStokes(kmax); + case NoneApproximationEnum: return PenaltyCreateKMatrixDiagnosticStokes(kmax); + default: ISSMERROR("not supported yet"); + } + case NoneApproximationEnum: + switch(approximation1){ + case StokesApproximationEnum: return PenaltyCreateKMatrixDiagnosticStokes(kmax); + case NoneApproximationEnum: return PenaltyCreateKMatrixDiagnosticStokes(kmax); + default: ISSMERROR("not supported yet"); + } + default: ISSMERROR("not supported yet"); + } +} +/*}}}1*/ +/*FUNCTION Penpair::PenaltyCreateKMatrixDiagnosticMacAyealPattyn {{{1*/ +ElementMatrix* Penpair::PenaltyCreateKMatrixDiagnosticMacAyealPattyn(double kmax){ + + const int numdof=NUMVERTICES*NDOF2; double penalty_offset; - /*pointers: */ - Node** nodes=NULL; - - /*Get dof list: */ - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - /*recover pointers: */ - nodes=(Node**)hnodes->deliverp(); + Node** nodes=(Node**)hnodes->deliverp(); + + /*Initialize Element vector and return if necessary*/ + ElementMatrix* Ke=NewElementMatrix(nodes,NUMVERTICES,this->parameters); /*recover parameters: */ @@ -306,42 +314,29 @@ //Create elementary matrix: add penalty to - Ke[0][0]=kmax*pow((double)10.0,penalty_offset); - Ke[0][2]=-kmax*pow((double)10.0,penalty_offset); - Ke[2][0]=-kmax*pow((double)10.0,penalty_offset); - Ke[2][2]=kmax*pow((double)10.0,penalty_offset); - - Ke[1][1]=kmax*pow((double)10.0,penalty_offset); - Ke[1][3]=-kmax*pow((double)10.0,penalty_offset); - Ke[3][1]=-kmax*pow((double)10.0,penalty_offset); - Ke[3][3]=kmax*pow((double)10.0,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*numdof+0]=+kmax*pow((double)10.0,penalty_offset); + Ke->values[0*numdof+2]=-kmax*pow((double)10.0,penalty_offset); + Ke->values[2*numdof+0]=-kmax*pow((double)10.0,penalty_offset); + Ke->values[2*numdof+2]=+kmax*pow((double)10.0,penalty_offset); + + Ke->values[1*numdof+1]=+kmax*pow((double)10.0,penalty_offset); + Ke->values[1*numdof+3]=-kmax*pow((double)10.0,penalty_offset); + Ke->values[3*numdof+1]=-kmax*pow((double)10.0,penalty_offset); + Ke->values[3*numdof+3]=+kmax*pow((double)10.0,penalty_offset); + + /*Clean up and return*/ + return Ke; } /*}}}1*/ /*FUNCTION Penpair::PenaltyCreateKMatrixDiagnosticStokes {{{1*/ -void Penpair::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax){ - - const int numgrids=2; - const int NDOF3=4; - const int numdof=numgrids*NDOF3; - int* doflist=NULL; - int numberofdofspernode; - - double Ke[8][8]={0.0}; +ElementMatrix* Penpair::PenaltyCreateKMatrixDiagnosticStokes(double kmax){ + + const int numdof=NUMVERTICES*NDOF4; double penalty_offset; - /*pointers: */ - Node** nodes=NULL; - - /*Get dof list: */ - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - /*recover pointers: */ - nodes=(Node**)hnodes->deliverp(); + Node** nodes=(Node**)hnodes->deliverp(); + + /*Initialize Element vector and return if necessary*/ + ElementMatrix* Ke=NewElementMatrix(nodes,NUMVERTICES,this->parameters); /*recover parameters: */ @@ -349,68 +344,26 @@ //Create elementary matrix: add penalty to - Ke[0][0]=kmax*pow((double)10.0,penalty_offset); - Ke[0][4]=-kmax*pow((double)10.0,penalty_offset); - Ke[4][0]=-kmax*pow((double)10.0,penalty_offset); - Ke[4][4]=kmax*pow((double)10.0,penalty_offset); - - Ke[1][1]=kmax*pow((double)10.0,penalty_offset); - Ke[1][5]=-kmax*pow((double)10.0,penalty_offset); - Ke[5][1]=-kmax*pow((double)10.0,penalty_offset); - Ke[5][5]=kmax*pow((double)10.0,penalty_offset); - - Ke[2][2]=kmax*pow((double)10.0,penalty_offset); - Ke[2][6]=-kmax*pow((double)10.0,penalty_offset); - Ke[6][2]=-kmax*pow((double)10.0,penalty_offset); - Ke[6][6]=kmax*pow((double)10.0,penalty_offset); - - Ke[3][3]=kmax*pow((double)10.0,penalty_offset); - Ke[3][7]=-kmax*pow((double)10.0,penalty_offset); - Ke[7][3]=-kmax*pow((double)10.0,penalty_offset); - Ke[7][7]=kmax*pow((double)10.0,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 Penpair::GetDofList {{{1*/ -void Penpair::GetDofList(int** pdoflist,int approximation,int setenum){ - - int i,j; - int numberofdofs=0; - int count=0; - - /*output: */ - int* doflist=NULL; - - /*pointers: */ - Node** nodes=NULL; - - /*recover pointers: */ - nodes=(Node**)hnodes->deliverp(); - - /*Some checks for debugging*/ - ISSMASSERT(nodes); - - /*First, figure out size of doflist: */ - for(i=0;i<2;i++){ - numberofdofs+=nodes[i]->GetNumberOfDofs(approximation,setenum); - } - - /*Allocate: */ - doflist=(int*)xmalloc(numberofdofs*sizeof(int)); - - /*Populate: */ - count=0; - for(i=0;i<2;i++){ - nodes[i]->GetDofList(doflist+count,approximation,setenum); - count+=nodes[i]->GetNumberOfDofs(approximation,setenum); - } - - /*Assign output pointers:*/ - *pdoflist=doflist; -} -/*}}}*/ + Ke->values[0*numdof+0]=+kmax*pow((double)10.0,penalty_offset); + Ke->values[0*numdof+4]=-kmax*pow((double)10.0,penalty_offset); + Ke->values[4*numdof+0]=-kmax*pow((double)10.0,penalty_offset); + Ke->values[4*numdof+4]=+kmax*pow((double)10.0,penalty_offset); + + Ke->values[1*numdof+1]=+kmax*pow((double)10.0,penalty_offset); + Ke->values[1*numdof+5]=-kmax*pow((double)10.0,penalty_offset); + Ke->values[5*numdof+1]=-kmax*pow((double)10.0,penalty_offset); + Ke->values[5*numdof+5]=+kmax*pow((double)10.0,penalty_offset); + + Ke->values[2*numdof+2]=+kmax*pow((double)10.0,penalty_offset); + Ke->values[2*numdof+6]=-kmax*pow((double)10.0,penalty_offset); + Ke->values[6*numdof+2]=-kmax*pow((double)10.0,penalty_offset); + Ke->values[6*numdof+6]=+kmax*pow((double)10.0,penalty_offset); + + Ke->values[3*numdof+3]=+kmax*pow((double)10.0,penalty_offset); + Ke->values[3*numdof+7]=-kmax*pow((double)10.0,penalty_offset); + Ke->values[7*numdof+3]=-kmax*pow((double)10.0,penalty_offset); + Ke->values[7*numdof+7]=+kmax*pow((double)10.0,penalty_offset); + + /*Clean up and return*/ + return Ke; +} +/*}}}1*/ Index: /issm/trunk/src/c/objects/Loads/Penpair.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Penpair.h (revision 5936) +++ /issm/trunk/src/c/objects/Loads/Penpair.h (revision 5937) @@ -64,7 +64,7 @@ /*}}}*/ /*Penpair management: {{{1*/ - void PenaltyCreateKMatrixDiagnosticHoriz(Mat Kgg,double kmax); - void PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax); - void GetDofList(int** pdoflist,int approximation_enum,int setenum); + ElementMatrix* PenaltyCreateKMatrixDiagnosticHoriz(double kmax); + ElementMatrix* PenaltyCreateKMatrixDiagnosticMacAyealPattyn(double kmax); + ElementMatrix* PenaltyCreateKMatrixDiagnosticStokes(double kmax); /*}}}*/ };