Changeset 5877

Timestamp:

09/17/10 17:18:54 (15 years ago)

Author:

Mathieu Morlighem

Message:

ElementMatrix in Vert

Location:

issm/trunk/src/c/objects/Elements

Files:

• Penta.cpp (modified) (5 diffs)
• Penta.h (modified) (1 diff)
• Tria.cpp (modified) (3 diffs)
• Tria.h (modified) (1 diff)

issm/trunk/src/c/objects/Elements/Penta.cpp

@@ -709 +709 @@
                         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 +2283 @@
 ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyeal3d(void){
 
-
         /*compute all stiffness matrices for this element*/
         ElementMatrix* Ke1=CreateKMatrixDiagnosticMacAyeal3dViscous();
@@ -2719 +2720 @@
 /*}}}*/
 /*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 +2747 @@
         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 +2760 @@
                 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; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
+                                        &Ke_gg[0][0],0);
+
+                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[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;
 }
 /*}}}*/

issm/trunk/src/c/objects/Elements/Penta.h

@@ -141 +141 @@
                 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);

issm/trunk/src/c/objects/Elements/Tria.cpp

@@ -3045 +3045 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreateKMatrixDiagnosticSurfaceVert {{{1*/
-void  Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg){
+/*FUNCTION Tria::CreateKMatrixDiagnosticVertSurface {{{1*/
+ElementMatrix* Tria::CreateKMatrixDiagnosticVertSurface(void){
 
         /*Constants*/
@@ -3059 +3059 @@
         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 +3085 @@
                                         &DL_scalar,1,1,0,
                                         L,1,3,0,
-                                        &Ke_gg_gaussian[0][0],0);
-
-                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
-        }
-
-        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
+                                        &Ke_g[0][0],0);
+
+                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
+        }
 
         /*Clean up and return*/
         delete gauss;
-        xfree((void**)&doflist);
+        return Ke;
 }
 /*}}}*/

issm/trunk/src/c/objects/Elements/Tria.h

@@ -132 +132 @@
                 ElementMatrix* CreateKMatrixDiagnosticPattynFriction(void);
                 ElementMatrix* CreateKMatrixDiagnosticHutter(void);
-                void      CreateKMatrixDiagnosticSurfaceVert(Mat Kgg);
+                ElementMatrix* CreateKMatrixDiagnosticVertSurface(void);
                 ElementMatrix* CreateKMatrixMelting(void);
                 ElementMatrix* CreateKMatrixPrognostic(void);