Changeset 5863

Timestamp:

09/17/10 09:36:52 (15 years ago)

Author:

Mathieu Morlighem

Message:

Added Hutter 3d

Location:

issm/trunk/src/c/objects/Elements

Files:

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

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

@@ -702 +702 @@
                         break;
                 case DiagnosticHutterAnalysisEnum:
-                        CreateKMatrixDiagnosticHutter( Kgg);
+                        Ke=CreateKMatrixDiagnosticHutter();
+                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
+                        delete Ke;
                         break;
                 case DiagnosticVertAnalysisEnum:
@@ -2248 +2250 @@
 /*}}}*/
 /*FUNCTION Penta::CreateKMatrixDiagnosticHutter{{{1*/
-void  Penta::CreateKMatrixDiagnosticHutter(Mat Kgg){
-
-        /*Collapsed formulation: */
-        int       i;
-        const int numdofs=NDOF2*NUMVERTICES;
-        int*      doflist=NULL;
-        double    Ke_gg[numdofs][numdofs]={0.0};
+ElementMatrix* Penta::CreateKMatrixDiagnosticHutter(void){
+
+        /*Intermediaries*/
+        const int numdof=NDOF2*NUMVERTICES;
         int       connectivity[2];
         double    one0,one1;
-        int       i0,i1,j0,j1;
-        
-        /*If on water, skip: */
-        if(IsOnWater())return;
-
-        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+        int       i,i0,i1,j0,j1;
+
+        /*Initialize Element matrix and return if necessary*/
+        if(IsOnWater()) return NULL;
+        ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum);
 
         /*Spawn 3 beam elements: */
@@ -2281 +2279 @@
                 /*Create matrix for these two nodes*/
                 if (IsOnBed() && IsOnSurface()){
-                        Ke_gg[i0][i0]=one0;
-                        Ke_gg[i1][i1]=one0;
-                        Ke_gg[j0][i0]=-one1;
-                        Ke_gg[j0][j0]=one1;
-                        Ke_gg[j1][i1]=-one1;
-                        Ke_gg[j1][j1]=one1;
+                        Ke->values[i0*numdof+i0]=one0;
+                        Ke->values[i1*numdof+i1]=one0;
+                        Ke->values[j0*numdof+i0]=-one1;
+                        Ke->values[j0*numdof+j0]=one1;
+                        Ke->values[j1*numdof+i1]=-one1;
+                        Ke->values[j1*numdof+j1]=one1;
                 }
                 else if (IsOnBed()){
-                        Ke_gg[i0][i0]=one0;
-                        Ke_gg[i1][i1]=one0;
-                        Ke_gg[j0][i0]=-2*one1;
-                        Ke_gg[j0][j0]=2*one1;
-                        Ke_gg[j1][i1]=-2*one1;
-                        Ke_gg[j1][j1]=2*one1;
+                        Ke->values[i0*numdof+i0]=one0;
+                        Ke->values[i1*numdof+i1]=one0;
+                        Ke->values[j0*numdof+i0]=-2*one1;
+                        Ke->values[j0*numdof+j0]=2*one1;
+                        Ke->values[j1*numdof+i1]=-2*one1;
+                        Ke->values[j1*numdof+j1]=2*one1;
                 }
                 else if (IsOnSurface()){
-                        Ke_gg[j0][i0]=-one1;
-                        Ke_gg[j0][j0]=one1;
-                        Ke_gg[j1][i1]=-one1;
-                        Ke_gg[j1][j1]=one1;
+                        Ke->values[j0*numdof+i0]=-one1;
+                        Ke->values[j0*numdof+j0]=one1;
+                        Ke->values[j1*numdof+i1]=-one1;
+                        Ke->values[j1*numdof+j1]=one1;
                 }
                 else{ //node is on two horizontal layers and beams include the values only once, so the have to use half of the connectivity
-                        Ke_gg[j0][i0]=-2*one1;
-                        Ke_gg[j0][j0]=2*one1;
-                        Ke_gg[j1][i1]=-2*one1;
-                        Ke_gg[j1][j1]=2*one1;
-                }
-        }
-
-        /*Add Ke_gg to global matrix Kgg: */
-        MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
-
-        /*Free ressources:*/
-        xfree((void**)&doflist);
-
+                        Ke->values[j0*numdof+i0]=-2*one1;
+                        Ke->values[j0*numdof+j0]=2*one1;
+                        Ke->values[j1*numdof+i1]=-2*one1;
+                        Ke->values[j1*numdof+j1]=2*one1;
+                }
+        }
+
+        /*Clean up and return*/
+        return Ke;
 }
 /*FUNCTION Penta::CreateKMatrixDiagnosticMacAyeal2d{{{1*/

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

@@ -128 +128 @@
                 void      CreateKMatrixCouplingPattynStokes( Mat Kgg);
                 void      CreateKMatrixDiagnosticHoriz( Mat Kgg);
-                void      CreateKMatrixDiagnosticHutter( Mat Kgg);
+                ElementMatrix* CreateKMatrixDiagnosticHutter(void);
                 ElementMatrix* CreateKMatrixDiagnosticMacAyeal2d(void);
                 ElementMatrix* CreateKMatrixDiagnosticMacAyeal3d(void);