Changeset 5678

Timestamp:

09/03/10 14:10:50 (15 years ago)

Author:

seroussi

Message:

remove ugly gauss in PVector Stokes

File:

• issm/trunk/src/c/objects/Elements/Penta.cpp (modified) (7 diffs)

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

@@ -2807 +2807 @@
 
         /* Start  looping on the number of gaussian points: */
-        /* Start  looping on the number of gaussian points: */
         gauss=new GaussPenta(5,5);
         for (ig=gauss->begin();ig<gauss->end();ig++){
@@ -3728 +3727 @@
 
         /* gaussian points: */
-        int     num_area_gauss;
-        int     igarea,igvert;
-        double* first_gauss_area_coord  =  NULL;
-        double* second_gauss_area_coord =  NULL;
-        double* third_gauss_area_coord  =  NULL;
-        double* vert_gauss_coord = NULL;
-        double* area_gauss_weights  =  NULL;
-        double* vert_gauss_weights  =  NULL;
-
-        /* specific gaussian point: */
-        double  gauss_weight,area_gauss_weight,vert_gauss_weight;
-        double  gauss_coord[4];
-        double  gauss_coord_tria[3];
-
-        int     area_order=5;
-        int       num_vert_gauss=5;
+        int     ig;
+        GaussPenta *gauss=NULL;
+        GaussTria  *gauss_tria=NULL;
 
         double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
@@ -3805 +3791 @@
         GetDofList(&doflist,StokesApproximationEnum);
 
-        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
-                get tria gaussian points as well as segment gaussian points. For tria gaussian 
-                points, order of integration is 2, because we need to integrate the product tB*D*B' 
-                which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian 
-                points, same deal, which yields 3 gaussian points.*/
-
-        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
-        gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
-
         /*Retrieve all inputs we will be needing: */
         vx_input=inputs->GetInput(VxEnum);
@@ -3821 +3798 @@
 
         /* Start  looping on the number of gaussian points: */
-        for (igarea=0; igarea<num_area_gauss; igarea++){
-                for (igvert=0; igvert<num_vert_gauss; igvert++){
-                        /*Pick up the gaussian point: */
-                        area_gauss_weight=*(area_gauss_weights+igarea);
-                        vert_gauss_weight=*(vert_gauss_weights+igvert);
-                        gauss_weight=area_gauss_weight*vert_gauss_weight;
-                        gauss_coord[0]=*(first_gauss_area_coord+igarea); 
-                        gauss_coord[1]=*(second_gauss_area_coord+igarea);
-                        gauss_coord[2]=*(third_gauss_area_coord+igarea);
-                        gauss_coord[3]=*(vert_gauss_coord+igvert);
-
-                        /*Compute strain rate and viscosity: */
-                        this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input);
-                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
-
-                        /* Get Jacobian determinant: */
-                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
-
-                        /* Get nodal functions */
-                        GetNodalFunctionsMINI(&l1l7[0], gauss_coord);
-
-                        /* Build gaussian vector */
-                        for(i=0;i<NUMVERTICES+1;i++){
-                                Pe_gaussian[i*NDOF4+2]=-rho_ice*gravity*Jdet*gauss_weight*l1l7[i];
-                        }
-
-                        /*Add Pe_gaussian to terms in Pe_temp. Watch out for column orientation from matlab: */
-                        for(i=0;i<27;i++){
-                                Pe_temp[i]+=Pe_gaussian[i];
-                        }
-
-                        /*Get B and Bprime matrices: */
-                        GetBStokes(&B[0][0],&xyz_list[0][0],gauss_coord); 
-                        GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss_coord); 
-
-                        /*Get bubble part of Bprime */
-                        for(i=0;i<8;i++){
-                                for(j=0;j<3;j++){
-                                        B_prime_bubble[i][j]=B_prime[i][j+24];
-                                }
-                        }
-
-                        /* 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;
-                        }
-                        for (i=6;i<8;i++){
-                                D[i][i]=-D_scalar*stokesreconditioning;
-                        }
-
-                        /*  Do the triple product tB*D*Bprime: */
-                        MatrixMultiply(&B[0][0],8,27,1,&D[0][0],8,8,0,&tBD[0][0],0);
-                        MatrixMultiply(&tBD[0][0],27,8,0,&B_prime_bubble[0][0],8,3,0,&Ke_gaussian[0][0],0);
-
-                        /*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */
-                        for(i=0;i<27;i++){
-                                for(j=0;j<3;j++){
-                                        Ke_temp[i][j]+=Ke_gaussian[i][j];
-                                }
-                        }
-                }
+        gauss=new GaussPenta(5,5);
+        for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                /*Compute strain rate and viscosity: */
+                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
+                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+
+                /* Get Jacobian determinant: */
+                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+
+                /* Get nodal functions */
+                GetNodalFunctionsMINI(&l1l7[0], gauss);
+
+                /* Build gaussian vector */
+                for(i=0;i<NUMVERTICES+1;i++){
+                        Pe_gaussian[i*NDOF4+2]=-rho_ice*gravity*Jdet*gauss->weight*l1l7[i];
+                }
+
+                /*Add Pe_gaussian to terms in Pe_temp. Watch out for column orientation from matlab: */
+                for(i=0;i<27;i++) Pe_temp[i]+=Pe_gaussian[i];
+
+                /*Get B and Bprime matrices: */
+                GetBStokes(&B[0][0],&xyz_list[0][0],gauss); 
+                GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss); 
+
+                /*Get bubble part of Bprime */
+                for(i=0;i<8;i++) for(j=0;j<3;j++) B_prime_bubble[i][j]=B_prime[i][j+24];
+
+                /* 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;
+                for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning;
+
+                /*  Do the triple product tB*D*Bprime: */
+                MatrixMultiply(&B[0][0],8,27,1,&D[0][0],8,8,0,&tBD[0][0],0);
+                MatrixMultiply(&tBD[0][0],27,8,0,&B_prime_bubble[0][0],8,3,0,&Ke_gaussian[0][0],0);
+
+                /*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */
+                for(i=0;i<27;i++) for(j=0;j<3;j++) Ke_temp[i][j]+=Ke_gaussian[i][j];
         }
 
@@ -3895 +3851 @@
                 }
 
-                xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); xfree((void**)&third_gauss_area_coord); xfree((void**)&area_gauss_weights);
-                GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
-
                 /* Start looping on the number of gauss 2d (nodes on the bedrock) */
-                for (igarea=0; igarea<num_area_gauss; igarea++){
-                        gauss_weight=*(area_gauss_weights+igarea);
-                        gauss_coord[0]=*(first_gauss_area_coord+igarea); 
-                        gauss_coord[1]=*(second_gauss_area_coord+igarea);
-                        gauss_coord[2]=*(third_gauss_area_coord+igarea);
-                        gauss_coord[3]=-1;
-
-                        gauss_coord_tria[0]=*(first_gauss_area_coord+igarea); 
-                        gauss_coord_tria[1]=*(second_gauss_area_coord+igarea);
-                        gauss_coord_tria[2]=*(third_gauss_area_coord+igarea);
+                delete gauss; //gauss of previous loop
+                gauss=new GaussPenta();
+                gauss->GaussFaceTria(0,1,2,2);
+                gauss_tria=new GaussTria();
+                for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                        gauss->GaussPoint(ig);
+                        gauss->SynchronizeGaussTria(gauss_tria);
 
                         /*Get the Jacobian determinant */
-                        tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria);
+                        tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_tria);
 
                         /* Get bed at gaussian point */
-                        bed_input->GetParameterValue(&bed, gauss_coord);
+                        bed_input->GetParameterValue(&bed, gauss);
 
                         /*Get L matrix : */
-                        GetNodalFunctionsP1(l1l6, gauss_coord);
+                        GetNodalFunctionsP1(l1l6, gauss);
 
                         /*Get water_pressure at gaussian point */
@@ -3927 +3878 @@
                         for(i=0;i<NUMVERTICES2D;i++){
                                 for(j=0;j<3;j++){
-                                        Pe_temp[i*NDOF4+j]+=water_pressure*gauss_weight*Jdet2d*l1l6[i]*bed_normal[j];
+                                        Pe_temp[i*NDOF4+j]+=water_pressure*gauss->weight*Jdet2d*l1l6[i]*bed_normal[j];
                                 }
                         }
@@ -3940 +3891 @@
 
         /*Free ressources:*/
-        xfree((void**)&first_gauss_area_coord);
-        xfree((void**)&second_gauss_area_coord);
-        xfree((void**)&third_gauss_area_coord);
-        xfree((void**)&area_gauss_weights);
-        xfree((void**)&vert_gauss_coord);
-        xfree((void**)&vert_gauss_weights); 
+        delete gauss;
+        delete gauss_tria;
         xfree((void**)&doflist);