Changeset 5875

Timestamp:

09/17/10 17:00:10 (15 years ago)

Author:

seroussi

Message:

cleanup in tria

File:

• issm/trunk/src/c/objects/Elements/Tria.cpp (modified) (8 diffs)

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

@@ -2862 +2862 @@
 ElementMatrix* Tria::CreateKMatrixCouplingMacAyealPattynFriction(void){
 
-        const int numdof        = 2 *NUMVERTICES;
-        const int numdoftotal   = 4 *NUMVERTICES;
-        int       i,j;
-        int     ig;
-        int       analysis_type;
+        /*Constants*/
+        const int numdof        = NDOF2 *NUMVERTICES;
+        const int numdoftotal   = NDOF4 *NUMVERTICES;
+        
+        /*Intermediaries */
+        int       i,j,ig,analysis_type,drag_type;
+        double    Jdet,slope_magnitude,alpha2;
         double    xyz_list[NUMVERTICES][3];
-        int       numberofdofspernode=2;
+        double    slope[2]={0.0,0.0};
+        double    MAXSLOPE=.06; // 6 %
+        double    MOUNTAINKEXPONENT=10;
+        double    L[2][numdof];
+        double    DL[2][2]                  ={{ 0,0 },{0,0}}; //for basal drag
+        double    DL_scalar;
+        double    Ke_gg[numdof][numdof]     ={0.0};
+        double    Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
+        Friction  *friction = NULL;
         GaussTria *gauss=NULL;
-        double L[2][numdof];
-        double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
-        double DL_scalar;
-        double Ke_gg[numdof][numdof]={0.0};
-        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
-        double Jdet;
-        double  slope[2]={0.0,0.0};
-        double  slope_magnitude;
-        Friction *friction = NULL;
-        double    alpha2;
-        double MAXSLOPE=.06; // 6 %
-        double MOUNTAINKEXPONENT=10;
-        int  drag_type;
 
         /*Initialize Element matrix and return if necessary*/
@@ -2926 +2923 @@
 
                 /*Get L matrix: */
-                GetL(&L[0][0], &xyz_list[0][0], gauss,numberofdofspernode);
+                GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2);
 
                 
@@ -2955 +2952 @@
 ElementMatrix* Tria::CreateKMatrixDiagnosticPattynFriction(void){
 
+        /*Constants*/
         const int numdof   = NDOF2*NUMVERTICES;
-        int       i,j;
-        int     ig;
-        int       analysis_type;
+        
+        /*Intermediaries */
+        int       i,j,ig;
+        int       analysis_type,drag_type;
         double    xyz_list[NUMVERTICES][3];
-        int       numberofdofspernode=2;
+        double    slope_magnitude,alpha2,Jdet;
+        double    slope[2]={0.0,0.0};
+        double    MAXSLOPE=.06; // 6 %
+        double    MOUNTAINKEXPONENT=10;
+        double    L[2][numdof];
+        double    DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
+        double    DL_scalar;
+        double    Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
+        Friction  *friction = NULL;
         GaussTria *gauss=NULL;
-        double L[2][numdof];
-        double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
-        double DL_scalar;
-        double Ke_gg[numdof][numdof]={0.0};
-        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
-        double Jdet;
-        double  slope[2]={0.0,0.0};
-        double  slope_magnitude;
-        Friction *friction = NULL;
-        double    alpha2;
-        double MAXSLOPE=.06; // 6 %
-        double MOUNTAINKEXPONENT=10;
-        int  drag_type;
 
         /*Initialize Element matrix and return if necessary*/
@@ -2999 +2993 @@
                 gauss->GaussPoint(ig);
 
+                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+                GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2);
+
+                surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
                 friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
+                slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
 
                 // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case, 
                 //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
-                surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
-                slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
-
                 if (slope_magnitude>MAXSLOPE){
                         alpha2=pow((double)10,MOUNTAINKEXPONENT);
                 }
-
-                GetL(&L[0][0], &xyz_list[0][0], gauss,numberofdofspernode);
-                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
                 
                 DL_scalar=alpha2*gauss->weight*Jdet;
@@ -3055 +3048 @@
 void  Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg){
 
-        int i,j;
-
-        /* node data: */
-        const int    numdof=NDOF1*NUMVERTICES;
-        double       xyz_list[NUMVERTICES][3];
-        int*         doflist=NULL;
-
-        /* gaussian points: */
-        int     ig;
+        /*Constants*/
+        const int numdof=NDOF1*NUMVERTICES;
+
+        /*Intermediaries */
+        int       i,j,ig;
+        int*      doflist=NULL;
+        double    xyz_list[NUMVERTICES][3];
+        double    surface_normal[3];
+        double    Jdet,DL_scalar;
+        double    L[3];
         GaussTria *gauss=NULL;
-
-        /* surface normal: */
-        double x4,y4,z4;
-        double x5,y5,z5;
-        double x6,y6,z6;
-        double v46[3];
-        double v56[3];
-        double normal[3];
-        double norm_normal;
-        double nz;
-
-        /*Matrices: */
-        double DL_scalar;
-        double L[3];
-        double Jdet;
 
         /* local element matrices: */
@@ -3088 +3067 @@
         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-        /*Build normal vector to the surface:*/
-        x4=xyz_list[0][0];
-        y4=xyz_list[0][1];
-        z4=xyz_list[0][2];
-
-        x5=xyz_list[1][0];
-        y5=xyz_list[1][1];
-        z5=xyz_list[1][2];
-
-        x6=xyz_list[2][0];
-        y6=xyz_list[2][1];
-        z6=xyz_list[2][2];
-
-        v46[0]=x4-x6;
-        v46[1]=y4-y6;
-        v46[2]=z4-z6;
-
-        v56[0]=x5-x6;
-        v56[1]=y5-y6;
-        v56[2]=z5-z6;
-
-        normal[0]=(y4-y6)*(z5-z6)-(z4-z6)*(y5-y6);
-        normal[1]=(z4-z6)*(x5-x6)-(x4-x6)*(z5-z6);
-        normal[2]=(x4-x6)*(y5-y6)-(y4-y6)*(x5-x6);
-
-        norm_normal=sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2));
-        nz=1.0/norm_normal*normal[2];
+        SurfaceNormal(&surface_normal[0],xyz_list);
 
         /* Start  looping on the number of gaussian points: */
@@ -3123 +3075 @@
                 gauss->GaussPoint(ig);
 
-                /* Get Jacobian determinant: */
                 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss);
-
-                //Get L matrix if viscous basal drag present:
                 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
 
                 /**********************Do not forget the sign**********************************/
-                DL_scalar=- gauss->weight*Jdet*nz; 
+                DL_scalar=- gauss->weight*Jdet*surface_normal[2]; 
                 /******************************************************************************/
 
-                /*  Do the triple producte tL*D*L: */
                 TripleMultiply( L,1,3,1,
                                         &DL_scalar,1,1,0,
@@ -3139 +3087 @@
                                         &Ke_gg_gaussian[0][0],0);
 
-                /* Add the 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];
-
-
-        }
-
-        /*Add Ke_gg to global matrix Kgg: */
+        }
+
         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);