Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5876) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5877) @@ -709,5 +709,7 @@ break; case DiagnosticVertAnalysisEnum: - CreateKMatrixDiagnosticVert( Kgg); + Ke=CreateKMatrixDiagnosticVert(); + if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); + delete Ke; break; case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: @@ -2281,5 +2283,4 @@ ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyeal3d(void){ - /*compute all stiffness matrices for this element*/ ElementMatrix* Ke1=CreateKMatrixDiagnosticMacAyeal3dViscous(); @@ -2719,21 +2720,26 @@ /*}}}*/ /*FUNCTION Penta::CreateKMatrixDiagnosticVert {{{1*/ -void Penta::CreateKMatrixDiagnosticVert( Mat Kgg){ +ElementMatrix* Penta::CreateKMatrixDiagnosticVert(void){ + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixDiagnosticVertVolume(); + ElementMatrix* Ke2=CreateKMatrixDiagnosticVertVolume(); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticVertVolume {{{1*/ +ElementMatrix* Penta::CreateKMatrixDiagnosticVertVolume(void){ /* local declarations */ - int i,j; - - /* node data: */ const int numdof=NDOF1*NUMVERTICES; + int i,j,ig; double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; - - /* 3d gaussian points: */ - int ig; GaussPenta *gauss=NULL; - - /* matrices: */ double Ke_gg[numdof][numdof]={0.0}; - double Ke_gg_gaussian[numdof][numdof]; double Jdet; double B[NDOF1][NUMVERTICES]; @@ -2741,23 +2747,10 @@ double DL_scalar; - /*Collapsed formulation: */ - Tria* tria=NULL; - - /*If on water, skip stiffness: */ - if(IsOnWater())return; - - /*If this element is on the surface, we have a dynamic boundary condition that applies, as a stiffness - * matrix: */ - if(IsOnSurface()){ - tria=(Tria*)SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface - tria->CreateKMatrixDiagnosticSurfaceVert(Kgg); - delete tria->matice; delete tria; - } - - /*Now, onto the formulation for the vertical velocity: */ - - /* Get node coordinates and dof list: */ + /*Initialize Element matrix and return if necessary*/ + if(IsOnWater()) return NULL; + ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum); + + /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /* Start looping on the number of gaussian points: */ @@ -2767,27 +2760,35 @@ gauss->GaussPoint(ig); - /*Get B and Bprime matrices: */ GetBVert(&B[0][0], &xyz_list[0][0], gauss); GetBprimeVert(&Bprime[0][0], &xyz_list[0][0], gauss); - /* Get Jacobian determinant: */ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); DL_scalar=gauss->weight*Jdet; - /* Do the triple product tB*D*Bprime: */ TripleMultiply( &B[0][0],1,NUMVERTICES,1, &DL_scalar,1,1,0, &Bprime[0][0],1,NUMVERTICES,0, - &Ke_gg_gaussian[0][0],0); - - /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ - for( i=0; ivalues[i*numdof+j]+=Ke_gg[i][j]; } - /*Add Ke_gg to global matrix Kgg: */ - MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); - - xfree((void**)&doflist); + /*Clean up and return*/ delete gauss; + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticVertSurface {{{1*/ +ElementMatrix* Penta::CreateKMatrixDiagnosticVertSurface(void){ + + if (!IsOnSurface() || IsOnWater()) return NULL; + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface + ElementMatrix* Ke=tria->CreateKMatrixDiagnosticVertSurface(); + delete tria->matice; delete tria; + + /*clean up and return*/ + return Ke; } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 5876) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 5877) @@ -141,5 +141,7 @@ ElementMatrix* CreateKMatrixDiagnosticStokesViscous(void); ElementMatrix* CreateKMatrixDiagnosticStokesFriction(void); - void CreateKMatrixDiagnosticVert( Mat Kgg); + ElementMatrix* CreateKMatrixDiagnosticVert(void); + ElementMatrix* CreateKMatrixDiagnosticVertVolume(void); + ElementMatrix* CreateKMatrixDiagnosticVertSurface(void); ElementMatrix* CreateKMatrixMelting(void); ElementMatrix* CreateKMatrixPrognostic(void); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5876) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5877) @@ -3045,6 +3045,6 @@ } /*}}}*/ -/*FUNCTION Tria::CreateKMatrixDiagnosticSurfaceVert {{{1*/ -void Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg){ +/*FUNCTION Tria::CreateKMatrixDiagnosticVertSurface {{{1*/ +ElementMatrix* Tria::CreateKMatrixDiagnosticVertSurface(void){ /*Constants*/ @@ -3059,12 +3059,12 @@ double L[3]; GaussTria *gauss=NULL; - - /* local element matrices: */ - double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix - double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. - - /* Get node coordinates and dof list: */ + double Ke_g[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. + + /*Initialize Element matrix and return if necessary*/ + if(IsOnWater()) return NULL; + ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum); + + /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); SurfaceNormal(&surface_normal[0],xyz_list); @@ -3085,14 +3085,12 @@ &DL_scalar,1,1,0, L,1,3,0, - &Ke_gg_gaussian[0][0],0); - - for( i=0; ivalues[i*numdof+j]+=Ke_g[i][j]; + } /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return Ke; } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 5876) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 5877) @@ -132,5 +132,5 @@ ElementMatrix* CreateKMatrixDiagnosticPattynFriction(void); ElementMatrix* CreateKMatrixDiagnosticHutter(void); - void CreateKMatrixDiagnosticSurfaceVert(Mat Kgg); + ElementMatrix* CreateKMatrixDiagnosticVertSurface(void); ElementMatrix* CreateKMatrixMelting(void); ElementMatrix* CreateKMatrixPrognostic(void);