Changeset 8416

Timestamp:

05/24/11 13:41:18 (14 years ago)

Author:

seroussi

Message:

moved KMatrix DiagnosticVertSurface from Tria to Penta

Location:

issm/trunk/src/c/objects

Files:

• Elements/Penta.cpp (modified) (1 diff)
• Elements/Tria.cpp (modified) (1 diff)
• Elements/Tria.h (modified) (1 diff)
• Gauss/GaussPenta.cpp (modified) (1 diff)

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

@@ -1628 +1628 @@
         if (!IsOnSurface()) 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*/
+        /*Constants*/
+        const int numdof=NDOF1*NUMVERTICES;
+
+        /*Intermediaries */
+        int       i,j,ig;
+        double    xyz_list[NUMVERTICES][3];
+        double    xyz_list_tria[NUMVERTICES2D][3];
+        double    surface_normal[3];
+        double    Jdet2d,DL_scalar;
+        double    l1l6[NUMVERTICES];
+        GaussPenta *gauss=NULL;
+        double Ke_g[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
+
+        /*Initialize Element matrix*/
+        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
+
+        /*Retrieve all inputs and parameters*/
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i+3][j];
+        SurfaceNormal(&surface_normal[0],xyz_list_tria);
+
+        /* Start  looping on the number of gaussian points: */
+        gauss=new GaussPenta(3,4,5,2);
+        for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
+                GetNodalFunctionsP1(&l1l6[0], gauss);
+
+                /**********************Do not forget the sign**********************************/
+                DL_scalar=- gauss->weight*Jdet2d*surface_normal[2]; 
+                /******************************************************************************/
+
+                TripleMultiply( l1l6,1,numdof,1,
+                                        &DL_scalar,1,1,0,
+                                        l1l6,1,numdof,0,
+                                        &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;
         return Ke;
 }

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

@@ -1087 +1087 @@
 
         /*Clean up and return*/
-        return Ke;
-}
-/*}}}*/
-/*FUNCTION Tria::CreateKMatrixDiagnosticVertSurface {{{1*/
-ElementMatrix* Tria::CreateKMatrixDiagnosticVertSurface(void){
-
-        /*Constants*/
-        const int numdof=NDOF1*NUMVERTICES;
-
-        /*Intermediaries */
-        int       i,j,ig;
-        double    xyz_list[NUMVERTICES][3];
-        double    surface_normal[3];
-        double    Jdet,DL_scalar;
-        double    L[3];
-        GaussTria *gauss=NULL;
-        double Ke_g[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
-
-        /*Initialize Element matrix*/
-        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
-
-        /*Retrieve all inputs and parameters*/
-        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-        SurfaceNormal(&surface_normal[0],xyz_list);
-
-        /* Start  looping on the number of gaussian points: */
-        gauss=new GaussTria(2);
-        for (ig=gauss->begin();ig<gauss->end();ig++){
-
-                gauss->GaussPoint(ig);
-
-                GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss);
-                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
-
-                /**********************Do not forget the sign**********************************/
-                DL_scalar=- gauss->weight*Jdet*surface_normal[2]; 
-                /******************************************************************************/
-
-                TripleMultiply( L,1,3,1,
-                                        &DL_scalar,1,1,0,
-                                        L,1,3,0,
-                                        &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;
         return Ke;
 }

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

@@ -153 +153 @@
                 ElementMatrix* CreateKMatrixDiagnosticPattynFriction(void);
                 ElementMatrix* CreateKMatrixDiagnosticHutter(void);
-                ElementMatrix* CreateKMatrixDiagnosticVertSurface(void);
                 ElementMatrix* CreateKMatrixMelting(void);
                 ElementMatrix* CreateKMatrixHydrology(void);

issm/trunk/src/c/objects/Gauss/GaussPenta.cpp ¶

@@ -149 +149 @@
                 coords4=(double*)xmalloc(numgauss*sizeof(double));
                 for(int i=0;i<numgauss;i++) coords4[i]=-1.0;
+        }
+        /*Upper surface Tria*/
+        else if(index1==3 && index2==4 && index3==5){
+
+                /*Get GaussTria*/
+                GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
+
+                /*compute z coordinate*/
+                coords4=(double*)xmalloc(numgauss*sizeof(double));
+                for(int i=0;i<numgauss;i++) coords4[i]=1.0;
         }
         else{