Changeset 5850

Timestamp:

09/16/10 14:46:34 (15 years ago)

Author:

Mathieu Morlighem

Message:

MacAyeal3d is now using ElemMatrix

Location:

issm/trunk/src/c/objects/Elements

Files:

• Penta.cpp (modified) (6 diffs)
• Penta.h (modified) (1 diff)

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

@@ -2212 +2212 @@
                         break;
                 case MacAyealPattynApproximationEnum:
-                        CreateKMatrixDiagnosticMacAyeal3d( Kgg);
+                        Ke=CreateKMatrixDiagnosticMacAyeal3d();
+                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
+                        delete Ke;
                         Ke=CreateKMatrixDiagnosticPattyn();
                         if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
@@ -2324 +2326 @@
 /*}}}*/
 /*FUNCTION Penta::CreateKMatrixDiagnosticMacAyeal3d{{{1*/
-void Penta::CreateKMatrixDiagnosticMacAyeal3d( Mat Kgg){
+ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyeal3d(void){
+
+
+        /*compute all stiffness matrices for this element*/
+        ElementMatrix* Ke1=CreateKMatrixDiagnosticMacAyeal3dViscous();
+        ElementMatrix* Ke2=CreateKMatrixDiagnosticMacAyeal3dFriction();
+        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
+
+        /*clean-up and return*/
+        delete Ke1;
+        delete Ke2;
+        return Ke;
+}
+/*}}}*/
+/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyeal3dViscous{{{1*/
+ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyeal3dViscous(void){
 
         const int    numdof2d=2*NUMVERTICES2D;
         int             i,j,ig;
         double       xyz_list[NUMVERTICES][3];
-        int*         doflist=NULL;
         GaussPenta *gauss=NULL;
         GaussTria  *gauss_tria=NULL;
@@ -2342 +2358 @@
         double DL[2][2]={0.0}; //for basal drag
         double DL_scalar;
-        double Ke_gg[numdof2d][numdof2d]={0.0}; //local element stiffness matrix 
         double Ke_gg_gaussian[numdof2d][numdof2d]; //stiffness matrix evaluated at the gaussian point.
         double Jdet;
@@ -2354 +2369 @@
         Penta* pentabase=NULL;
 
-        /*If on water, skip stiffness: */
-        if(IsOnWater())return;
-
         /*Find penta on bed as this is a macayeal elements: */
         pentabase=GetBasalElement();
         tria=pentabase->SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 
-        /* Get dof list: */
-        tria->GetDofList(&doflist,MacAyealApproximationEnum,GsetEnum);  //Pattyn dof list
+        /*Initialize Element matrix and return if necessary*/
+        if(IsOnWater()) return NULL;
+        ElementMatrix* Ke=tria->NewElementMatrix(MacAyealApproximationEnum);
 
         /*Retrieve all inputs and parameters*/
@@ -2398 +2411 @@
                                         &Ke_gg_gaussian[0][0],0);
 
-                for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
-        }
-
-        MatSetValues(Kgg,numdof2d,doflist,numdof2d,doflist,(const double*)Ke_gg,ADD_VALUES);
-
-        //Deal with 2d friction at the bedrock interface
-        if(IsOnBed() && !IsOnShelf()){
-                /*Build a tria element using the 3 grids of the base of the penta. Then use 
-                 * the tria functionality to build a friction stiffness matrix on these 3
-                 * grids: */
-
-                tria->CreateKMatrixDiagnosticMacAyealFriction(Kgg,NULL,NULL);
+                for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof2d+j]+=Ke_gg_gaussian[i][j];
         }
 
@@ -2417 +2419 @@
         delete gauss_tria;
         delete gauss;
-        xfree((void**)&doflist);
+        return Ke;
+}
+/*}}}*/
+/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyeal3dFriction{{{1*/
+ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyeal3dFriction(void){
+
+        /*Initialize Element matrix and return if necessary*/
+        if(IsOnWater() || IsOnShelf() || !IsOnBed()) return NULL;
+
+        /*Build a tria element using the 3 grids of the base of the penta. Then use 
+         * the tria functionality to build a friction stiffness matrix on these 3
+         * grids: */
+        Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+        ElementMatrix* Ke=tria->CreateKMatrixDiagnosticMacAyealFriction();
+        delete tria->matice; delete tria;
+
+        /*clean-up and return*/
+        return Ke;
 }
 /*}}}*/

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

@@ -128 +128 @@
                 void      CreateKMatrixDiagnosticHutter( Mat Kgg);
                 ElementMatrix* CreateKMatrixDiagnosticMacAyeal2d(void);
-                void      CreateKMatrixDiagnosticMacAyeal3d(Mat Kgg);
+                ElementMatrix* CreateKMatrixDiagnosticMacAyeal3d(void);
+                ElementMatrix* CreateKMatrixDiagnosticMacAyeal3dViscous(void);
+                ElementMatrix* CreateKMatrixDiagnosticMacAyeal3dFriction(void);
                 ElementMatrix* CreateKMatrixDiagnosticPattyn(void);
                 ElementMatrix* CreateKMatrixDiagnosticPattynViscous(void);