Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15694) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15695) @@ -391,68 +391,12 @@ } /*}}}*/ -/*FUNCTION Penta::CreateKMatrix {{{*/ +/*FUNCTION Penta::CreateKMatrix(Matrix* Kff, Matrix* Kfs) {{{*/ void Penta::CreateKMatrix(Matrix* Kff, Matrix* Kfs){ - - /*retrieve parameters: */ - ElementMatrix* Ke=NULL; - int analysis_type; - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*Checks in debugging {{{*/ - _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs); - /*}}}*/ /*Skip if water element*/ if(IsOnWater()) return; - /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - switch(analysis_type){ - #ifdef _HAVE_DIAGNOSTIC_ - case DiagnosticHorizAnalysisEnum: - Ke=CreateKMatrixDiagnosticHoriz(); - break; - case AdjointHorizAnalysisEnum: - Ke=CreateKMatrixAdjointHoriz(); - break; - case DiagnosticSIAAnalysisEnum: - Ke=CreateKMatrixDiagnosticSIA(); - break; - case DiagnosticVertAnalysisEnum: - Ke=CreateKMatrixDiagnosticVert(); - break; - #endif - case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: - Ke=CreateBasalMassMatrix(); - break; - case PrognosticAnalysisEnum: - Ke=CreateKMatrixPrognostic(); - break; - #ifdef _HAVE_BALANCED_ - case BalancethicknessAnalysisEnum: - Ke=CreateKMatrixBalancethickness(); - break; - #endif - #ifdef _HAVE_THERMAL_ - case ThermalAnalysisEnum: - Ke=CreateKMatrixThermal(); - break; - case EnthalpyAnalysisEnum: - Ke=CreateKMatrixEnthalpy(); - break; - case MeltingAnalysisEnum: - Ke=CreateKMatrixMelting(); - break; - #endif - #ifdef _HAVE_HYDROLOGY_ - case HydrologyDCInefficientAnalysisEnum: - Ke=CreateKMatrixHydrologyDCInefficient(); - break; - case HydrologyDCEfficientAnalysisEnum: - Ke=CreateKMatrixHydrologyDCEfficient(); - break; - #endif - default: - _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); - } + /*Create element stiffness matrix*/ + ElementMatrix* Ke=CreateKMatrix(); if(Ke){ @@ -462,4 +406,13 @@ Ke->StaticCondensation(3,&indices[0]); } + else if(this->element_type==P1bubblecondensedEnum){ + int size = nodes[6]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum); + int offset = 0; + for(int i=0;i<6;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum); + int* indices=xNew(size); + for(int i=0;iStaticCondensation(size,indices); + xDelete(indices); + } /*Add to global matrix*/ @@ -467,4 +420,73 @@ delete Ke; } +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrix(void){{{*/ +ElementMatrix* Penta::CreateKMatrix(void){ + + /*retrieve parameters: */ + ElementMatrix* Ke=NULL; + int analysis_type; + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*Checks in debugging {{{*/ + _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs); + /*}}}*/ + + /*Skip if water element*/ + if(IsOnWater()) return NULL; + + /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ + switch(analysis_type){ + #ifdef _HAVE_DIAGNOSTIC_ + case DiagnosticHorizAnalysisEnum: + return CreateKMatrixDiagnosticHoriz(); + break; + case AdjointHorizAnalysisEnum: + return CreateKMatrixAdjointHoriz(); + break; + case DiagnosticSIAAnalysisEnum: + return CreateKMatrixDiagnosticSIA(); + break; + case DiagnosticVertAnalysisEnum: + return CreateKMatrixDiagnosticVert(); + break; + #endif + case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: + return CreateBasalMassMatrix(); + break; + case PrognosticAnalysisEnum: + return CreateKMatrixPrognostic(); + break; + #ifdef _HAVE_BALANCED_ + case BalancethicknessAnalysisEnum: + return CreateKMatrixBalancethickness(); + break; + #endif + #ifdef _HAVE_THERMAL_ + case ThermalAnalysisEnum: + return CreateKMatrixThermal(); + break; + case EnthalpyAnalysisEnum: + return CreateKMatrixEnthalpy(); + break; + case MeltingAnalysisEnum: + return CreateKMatrixMelting(); + break; + #endif + #ifdef _HAVE_HYDROLOGY_ + case HydrologyDCInefficientAnalysisEnum: + return CreateKMatrixHydrologyDCInefficient(); + break; + case HydrologyDCEfficientAnalysisEnum: + return CreateKMatrixHydrologyDCEfficient(); + break; + #endif + default: + _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); + } + + /*Make compiler happy*/ + return NULL; } /*}}}*/ @@ -531,9 +553,50 @@ } /*}}}*/ -/*FUNCTION Penta::CreatePVector {{{*/ +/*FUNCTION Penta::CreatePVector(Vector* pf) {{{*/ void Penta::CreatePVector(Vector* pf){ - /*retrive parameters: */ - ElementVector* pe=NULL; + /*Skip if water element*/ + if(IsOnWater()) return; + + /*Create element load vector*/ + ElementVector* pe = CreatePVector(); + + if(pe){ + /*StaticCondensation if requested*/ + if(this->element_type==MINIcondensedEnum){ + int indices[3]={18,19,20}; + + this->element_type=MINIEnum; + ElementMatrix* Ke = CreateKMatrixDiagnosticFS(); + this->element_type=MINIcondensedEnum; + + pe->StaticCondensation(Ke,3,&indices[0]); + delete Ke; + } + else if(this->element_type==P1bubblecondensedEnum){ + int size = nodes[6]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum); + int offset = 0; + for(int i=0;i<6;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum); + int* indices=xNew(size); + for(int i=0;ielement_type=P1bubbleEnum; + ElementMatrix* Ke = CreateKMatrix(); + this->element_type=P1bubblecondensedEnum; + pe->StaticCondensation(Ke,size,indices); + xDelete(indices); + delete Ke; + } + + /*Add to global Vector*/ + pe->AddToGlobal(pf); + delete pe; + } +} +/*}}}*/ +/*FUNCTION Penta::CreatePVector(void) {{{*/ +ElementVector* Penta::CreatePVector(void){ + + /*retrieve parameters: */ int analysis_type; parameters->FindParam(&analysis_type,AnalysisTypeEnum); @@ -544,5 +607,5 @@ /*Skip if water element*/ - if(IsOnWater()) return; + if(IsOnWater()) return NULL; /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ @@ -550,46 +613,46 @@ #ifdef _HAVE_DIAGNOSTIC_ case DiagnosticHorizAnalysisEnum: - pe=CreatePVectorDiagnosticHoriz(); + return CreatePVectorDiagnosticHoriz(); break; case DiagnosticSIAAnalysisEnum: - pe=CreatePVectorDiagnosticSIA(); + return CreatePVectorDiagnosticSIA(); break; case DiagnosticVertAnalysisEnum: - pe=CreatePVectorDiagnosticVert(); + return CreatePVectorDiagnosticVert(); break; #endif #ifdef _HAVE_CONTROL_ case AdjointHorizAnalysisEnum: - pe=CreatePVectorAdjointHoriz(); + return CreatePVectorAdjointHoriz(); break; #endif #ifdef _HAVE_THERMAL_ case ThermalAnalysisEnum: - pe=CreatePVectorThermal(); + return CreatePVectorThermal(); break; case EnthalpyAnalysisEnum: - pe=CreatePVectorEnthalpy(); + return CreatePVectorEnthalpy(); break; case MeltingAnalysisEnum: - pe=CreatePVectorMelting(); + return CreatePVectorMelting(); break; #endif case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: - pe=CreatePVectorSlope(); + return CreatePVectorSlope(); break; case PrognosticAnalysisEnum: - pe=CreatePVectorPrognostic(); + return CreatePVectorPrognostic(); break; #ifdef _HAVE_BALANCED_ case BalancethicknessAnalysisEnum: - pe=CreatePVectorBalancethickness(); + return CreatePVectorBalancethickness(); break; #endif #ifdef _HAVE_HYDROLOGY_ case HydrologyDCInefficientAnalysisEnum: - pe=CreatePVectorHydrologyDCInefficient(); + return CreatePVectorHydrologyDCInefficient(); break; case HydrologyDCEfficientAnalysisEnum: - pe=CreatePVectorHydrologyDCEfficient(); + return CreatePVectorHydrologyDCEfficient(); break; #endif @@ -598,22 +661,6 @@ } - if(pe){ - - /*StaticCondensation if requested*/ - if(this->element_type==MINIcondensedEnum){ - int indices[3]={18,19,20}; - - this->element_type=MINIEnum; - ElementMatrix* Ke = CreateKMatrixDiagnosticFS(); - this->element_type=MINIcondensedEnum; - - pe->StaticCondensation(Ke,3,&indices[0]); - delete Ke; - } - - /*Add to global Vector*/ - pe->AddToGlobal(pf); - delete pe; - } + /*make compiler happy*/ + return NULL; } /*}}}*/ @@ -8328,7 +8375,4 @@ ElementVector* Penta::CreatePVectorDiagnosticHODrivingStress(void){ - /*Constants*/ - const int numdof=NDOF2*NUMVERTICES; - /*Intermediaries*/ int i,j; @@ -8337,9 +8381,13 @@ IssmDouble driving_stress_baseline,thickness; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble basis[6]; GaussPenta *gauss=NULL; - /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,HOApproximationEnum); + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF2; + + /*Initialize vector*/ + ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters,HOApproximationEnum); + IssmDouble* basis = xNewZeroInit(numnodes); /*Retrieve all inputs and parameters*/ @@ -8355,5 +8403,5 @@ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsP1(basis, gauss); + GetNodalFunctions(basis,gauss); thickness_input->GetInputValue(&thickness, gauss); @@ -8362,12 +8410,16 @@ driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG(); - for(i=0;ivalues[i*NDOF2+j]+= -driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i]; + for(i=0;ivalues[i*NDOF2+0]+= -driving_stress_baseline*slope[0]*Jdet*gauss->weight*basis[i]; + pe->values[i*NDOF2+1]+= -driving_stress_baseline*slope[1]*Jdet*gauss->weight*basis[i]; + } } /*Transform coordinate system*/ - TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYEnum); + TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum); /*Clean up and return*/ delete gauss; + xDelete(basis); return pe; } Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15694) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15695) @@ -179,5 +179,7 @@ void BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]); ElementMatrix* CreateBasalMassMatrix(void); + ElementMatrix* CreateKMatrix(void); ElementMatrix* CreateKMatrixPrognostic(void); + ElementVector* CreatePVector(void); ElementVector* CreatePVectorPrognostic(void); ElementVector* CreatePVectorSlope(void); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15694) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15695) @@ -172,62 +172,12 @@ } /*}}}*/ -/*FUNCTION Tria::CreateKMatrix {{{*/ +/*FUNCTION Tria::CreateKMatrix(Matrix* Kff, Matrix* Kfs) {{{*/ void Tria::CreateKMatrix(Matrix* Kff, Matrix* Kfs){ - - /*retreive parameters: */ - ElementMatrix* Ke=NULL; - int analysis_type; - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*Checks in debugging mode{{{*/ - _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs); - /*}}}*/ /*Skip if water element*/ if(IsOnWater()) return; - /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - switch(analysis_type){ - #ifdef _HAVE_DIAGNOSTIC_ - case DiagnosticHorizAnalysisEnum: - Ke=CreateKMatrixDiagnosticSSA(); - break; - case DiagnosticSIAAnalysisEnum: - Ke=CreateKMatrixDiagnosticSIA(); - break; - #endif - case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: - Ke=CreateMassMatrix(); - break; - case PrognosticAnalysisEnum: - Ke=CreateKMatrixPrognostic(); - break; - #ifdef _HAVE_HYDROLOGY_ - case HydrologyShreveAnalysisEnum: - Ke=CreateKMatrixHydrologyShreve(); - break; - case HydrologyDCInefficientAnalysisEnum: - Ke=CreateKMatrixHydrologyDCInefficient(); - break; - case HydrologyDCEfficientAnalysisEnum: - Ke=CreateKMatrixHydrologyDCEfficient(); - break; - #endif - #ifdef _HAVE_BALANCED_ - case BalancethicknessAnalysisEnum: - Ke=CreateKMatrixBalancethickness(); - break; - #endif - #ifdef _HAVE_CONTROL_ - case AdjointBalancethicknessAnalysisEnum: - Ke=CreateKMatrixAdjointBalancethickness(); - break; - case AdjointHorizAnalysisEnum: - Ke=CreateKMatrixAdjointSSA(); - break; - #endif - default: - _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); - } + /*Create element stiffness matrix*/ + ElementMatrix* Ke=CreateKMatrix(); if(Ke){ @@ -239,5 +189,4 @@ int* indices=xNew(size); for(int i=0;iStaticCondensation(size,indices); xDelete(indices); @@ -248,4 +197,66 @@ delete Ke; } +} +/*}}}*/ +/*FUNCTION Tria::CreateKMatrix(void) {{{*/ +ElementMatrix* Tria::CreateKMatrix(void){ + + /*retreive parameters: */ + int analysis_type; + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*Checks in debugging mode{{{*/ + _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs); + /*}}}*/ + + /*Skip if water element*/ + if(IsOnWater()) return NULL; + + /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ + switch(analysis_type){ + #ifdef _HAVE_DIAGNOSTIC_ + case DiagnosticHorizAnalysisEnum: + return CreateKMatrixDiagnosticSSA(); + break; + case DiagnosticSIAAnalysisEnum: + return CreateKMatrixDiagnosticSIA(); + break; + #endif + case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: + return CreateMassMatrix(); + break; + case PrognosticAnalysisEnum: + return CreateKMatrixPrognostic(); + break; + #ifdef _HAVE_HYDROLOGY_ + case HydrologyShreveAnalysisEnum: + return CreateKMatrixHydrologyShreve(); + break; + case HydrologyDCInefficientAnalysisEnum: + return CreateKMatrixHydrologyDCInefficient(); + break; + case HydrologyDCEfficientAnalysisEnum: + return CreateKMatrixHydrologyDCEfficient(); + break; + #endif + #ifdef _HAVE_BALANCED_ + case BalancethicknessAnalysisEnum: + return CreateKMatrixBalancethickness(); + break; + #endif + #ifdef _HAVE_CONTROL_ + case AdjointBalancethicknessAnalysisEnum: + return CreateKMatrixAdjointBalancethickness(); + break; + case AdjointHorizAnalysisEnum: + return CreateKMatrixAdjointSSA(); + break; + #endif + default: + _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); + } + + /*Make compiler happy*/ + return NULL; } /*}}}*/ @@ -296,63 +307,12 @@ } /*}}}*/ -/*FUNCTION Tria::CreatePVector {{{*/ +/*FUNCTION Tria::CreatePVector(Vector* pf) {{{*/ void Tria::CreatePVector(Vector* pf){ - - /*retrive parameters: */ - ElementVector* pe=NULL; - int analysis_type; - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*asserts: {{{*/ - /*if debugging mode, check that all pointers exist*/ - _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs); - /*}}}*/ /*Skip if water element*/ if(IsOnWater()) return; - /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ - switch(analysis_type){ - #ifdef _HAVE_DIAGNOSTIC_ - case DiagnosticHorizAnalysisEnum: - pe=CreatePVectorDiagnosticSSA(); - break; - case DiagnosticSIAAnalysisEnum: - pe=CreatePVectorDiagnosticSIA(); - break; - #endif - case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: - pe=CreatePVectorSlope(); - break; - case PrognosticAnalysisEnum: - pe=CreatePVectorPrognostic(); - break; - #ifdef _HAVE_HYDROLOGY_ - case HydrologyShreveAnalysisEnum: - pe=CreatePVectorHydrologyShreve(); - break; - case HydrologyDCInefficientAnalysisEnum: - pe=CreatePVectorHydrologyDCInefficient(); - break; - case HydrologyDCEfficientAnalysisEnum: - pe=CreatePVectorHydrologyDCEfficient(); - break; - #endif - #ifdef _HAVE_BALANCED_ - case BalancethicknessAnalysisEnum: - pe=CreatePVectorBalancethickness(); - break; - #endif - #ifdef _HAVE_CONTROL_ - case AdjointBalancethicknessAnalysisEnum: - pe=CreatePVectorAdjointBalancethickness(); - break; - case AdjointHorizAnalysisEnum: - pe=CreatePVectorAdjointHoriz(); - break; - #endif - default: - _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); - } + /*Create element load vector*/ + ElementVector* pe = CreatePVector(); /*Add to global Vector*/ @@ -367,5 +327,5 @@ this->element_type=P1bubbleEnum; - ElementMatrix* Ke = CreateKMatrixDiagnosticSSA(); + ElementMatrix* Ke = CreateKMatrix(); this->element_type=P1bubblecondensedEnum; @@ -378,4 +338,68 @@ delete pe; } +} +/*}}}*/ +/*FUNCTION Tria::CreatePVector(void){{{*/ +ElementVector* Tria::CreatePVector(void){ + + /*retrive parameters: */ + ElementVector* pe=NULL; + int analysis_type; + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*asserts: {{{*/ + /*if debugging mode, check that all pointers exist*/ + _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs); + /*}}}*/ + + /*Skip if water element*/ + if(IsOnWater()) return NULL; + + /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ + switch(analysis_type){ +#ifdef _HAVE_DIAGNOSTIC_ + case DiagnosticHorizAnalysisEnum: + return CreatePVectorDiagnosticSSA(); + break; + case DiagnosticSIAAnalysisEnum: + return CreatePVectorDiagnosticSIA(); + break; +#endif + case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: + return CreatePVectorSlope(); + break; + case PrognosticAnalysisEnum: + return CreatePVectorPrognostic(); + break; +#ifdef _HAVE_HYDROLOGY_ + case HydrologyShreveAnalysisEnum: + return CreatePVectorHydrologyShreve(); + break; + case HydrologyDCInefficientAnalysisEnum: + return CreatePVectorHydrologyDCInefficient(); + break; + case HydrologyDCEfficientAnalysisEnum: + return CreatePVectorHydrologyDCEfficient(); + break; +#endif +#ifdef _HAVE_BALANCED_ + case BalancethicknessAnalysisEnum: + return CreatePVectorBalancethickness(); + break; +#endif +#ifdef _HAVE_CONTROL_ + case AdjointBalancethicknessAnalysisEnum: + return CreatePVectorAdjointBalancethickness(); + break; + case AdjointHorizAnalysisEnum: + return CreatePVectorAdjointHoriz(); + break; +#endif + default: + _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); + } + + /*make compiler happy*/ + return NULL; } /*}}}*/ @@ -3153,8 +3177,7 @@ /*Build load vector: */ - for (i=0;ivalues[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i]; - } + for(i=0;ivalues[i*NDOF2+0]+=-driving_stress_baseline*slope[0]*Jdet*gauss->weight*basis[i]; + pe->values[i*NDOF2+1]+=-driving_stress_baseline*slope[1]*Jdet*gauss->weight*basis[i]; } } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15694) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15695) @@ -179,4 +179,5 @@ /*}}}*/ /*Tria specific routines:{{{*/ + ElementMatrix* CreateKMatrix(void); ElementMatrix* CreateKMatrixBalancethickness(void); ElementMatrix* CreateKMatrixBalancethickness_DG(void); @@ -187,4 +188,5 @@ ElementMatrix* CreateKMatrixPrognostic_DG(void); ElementMatrix* CreateMassMatrix(void); + ElementVector* CreatePVector(void); ElementVector* CreatePVectorBalancethickness(void); ElementVector* CreatePVectorBalancethickness_DG(void);