Changeset 5934

Timestamp:

09/21/10 16:58:02 (15 years ago)

Author:

seroussi

Message:

some cleaning in KMAtrix penta

Location:

issm/trunk/src/c/objects/Elements

Files:

• Penta.cpp (modified) (8 diffs)
• Tria.cpp (modified) (1 diff)

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

@@ -2371 +2371 @@
 ElementMatrix* Penta::CreateKMatrixDiagnosticPattynViscous(void){
 
-        /* local declarations */
-        const int    numdof=2*NUMVERTICES;
-        int             i,j,ig;
-        double       xyz_list[NUMVERTICES][3];
+        /*Constants*/
+        const int    numdof=NDOF2*NUMVERTICES;
+
+        /*Intermediaries */
+        int        i,j,ig;
+        double     xyz_list[NUMVERTICES][3];
+        double     Jdet;
+        double     viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity
+        double     epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
+        double     D_scalar;
+        double     D[5][5]={0.0};            // material matrix, simple scalar matrix.
+        double     B[5][numdof];
+        double     Bprime[5][numdof];
+        double     Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
+        Tria*      tria=NULL;
         GaussPenta *gauss=NULL;
-        double viscosity; //viscosity
-        double oldviscosity; //viscosity
-        double newviscosity; //viscosity
-        double viscosity_overshoot;
-        double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
-        double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
-        double B[5][numdof];
-        double Bprime[5][numdof];
-        double L[2][numdof];
-        double D[5][5]={0.0};            // material matrix, simple scalar matrix.
-        double D_scalar;
-        double DL[2][2]={0.0}; //for basal drag
-        double DL_scalar;
-        double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix 
-        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
-        double Jdet;
-        double  slope[2]={0.0};
-        double  slope_magnitude;
-        double  alpha2_list[3];
-        double  alpha2;
-        double MAXSLOPE=.06; // 6 %
-        double MOUNTAINKEXPONENT=10;
-        Tria*  tria=NULL;
 
         /*Initialize Element matrix and return if necessary*/
@@ -2418 +2406 @@
                 gauss->GaussPoint(ig);
 
+                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+                GetBPattyn(&B[0][0], &xyz_list[0][0], gauss);
+                GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss);
+
                 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
                 this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
                 matice->GetViscosity3d(&viscosity, &epsilon[0]);
                 matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
-
-                GetBPattyn(&B[0][0], &xyz_list[0][0], gauss);
-                GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss);
-
-                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 
                 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
@@ -2495 +2482 @@
 ElementMatrix* Penta::CreateKMatrixDiagnosticStokesViscous(void){
 
-        const int numdof=NUMVERTICES*NDOF4;
-        int i,j;
-        int     ig;
+        /*Intermediaries */
+        int        i,j,ig,approximation;
+        double     Jdet,viscosity,stokesreconditioning;
         double     xyz_list[NUMVERTICES][3];
-        double    xyz_list_tria[NUMVERTICES2D][3];
-        double    bed_normal[3];
-        double     Ke_temp[27][27]={0.0}; //for the six nodes and the bubble 
-        double     Ke_reduced[numdof][numdof]; //for the six nodes only
-        double     Ke_gaussian[27][27];
+        double     epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
         double     B[8][27];
         double     B_prime[8][27];
-        double     Jdet;
+        double     D_scalar;
         double     D[8][8]={0.0};
-        double     D_scalar;
-        double     DLStokes[14][14]={0.0};
+        double     Ke_temp[27][27]={0.0}; //for the six nodes and the bubble 
+        double     Ke_gaussian[27][27];
         GaussPenta *gauss=NULL;
-        double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
-        double  viscosity;
-        double  alpha2_gauss;
-        Friction* friction=NULL;
-        double stokesreconditioning;
-        int analysis_type;
-        int approximation;
 
         /*If on water or not Stokes, skip stiffness: */
@@ -2526 +2502 @@
         /*Retrieve all inputs and parameters*/
         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
         parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
         Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
@@ -2538 +2513 @@
                 gauss->GaussPoint(ig);
 
+                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+                GetBStokes(&B[0][0],&xyz_list[0][0],gauss); 
+                GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0],gauss); 
+
                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
                 matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
 
-                GetBStokes(&B[0][0],&xyz_list[0][0],gauss); 
-                GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0],gauss); 
-
-                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-
-                /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant 
-                 * onto this scalar matrix, so that we win some computational time: */
                 D_scalar=gauss->weight*Jdet;
                 for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity;
@@ -2571 +2543 @@
 ElementMatrix* Penta::CreateKMatrixDiagnosticStokesFriction(void){
 
+        /*Constants*/
         const int numdof=NUMVERTICES*NDOF4;
         const int numdof2d=NUMVERTICES2D*NDOF4;
-        int i,j;
-        int     ig;
+
+        /*Intermediaries */
+        int        i,j,ig;
+        int        analysis_type,approximation;
+        double     stokesreconditioning;
+        double     viscosity,alpha2_gauss,Jdet2d;
+        double    bed_normal[3];
+        double     epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
         double     xyz_list[NUMVERTICES][3];
         double    xyz_list_tria[NUMVERTICES2D][3];
-        double    bed_normal[3];
-        double     Ke_temp[27][27]={0.0}; //for the six nodes and the bubble 
-        double     Ke_reduced[numdof][numdof]; //for the six nodes only
-        double     Ke_gaussian[27][27];
-        double     Ke_drag_gaussian[numdof2d][numdof2d];
-        double     B[8][27];
-        double     B_prime[8][27];
         double     LStokes[14][numdof2d];
         double     LprimeStokes[14][numdof2d];
-        double     Jdet;
-        double     Jdet2d;
-        double     D[8][8]={0.0};
-        double     D_scalar;
         double     DLStokes[14][14]={0.0};
+        double     Ke_temp[27][27]={0.0}; //for the six nodes and the bubble 
+        double     Ke_drag_gaussian[numdof2d][numdof2d];
+        Friction*  friction=NULL;
         GaussPenta *gauss=NULL;
-        double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
-        double  viscosity;
-        double  alpha2_gauss;
-        Friction* friction=NULL;
-        double stokesreconditioning;
-        int analysis_type;
-        int approximation;
 
         /*If on water or not Stokes, skip stiffness: */
@@ -2623 +2587 @@
                 gauss->GaussPoint(ig);
 
+                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
                 GetLStokes(&LStokes[0][0], gauss);
                 GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss);
@@ -2630 +2595 @@
 
                 BedNormal(&bed_normal[0],xyz_list_tria);
-                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
                 friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
 

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

@@ -2401 +2401 @@
 ElementMatrix* Tria::CreateKMatrixBalancedthickness_CG(void){
 
-
         /*Constants*/
         const int    numdof=NDOF1*NUMVERTICES;