Changeset 5651

Timestamp:

09/01/10 16:04:42 (15 years ago)

Author:

seroussi

Message:

started to use new gauss points in penta

Location:

issm/trunk/src/c/objects

Files:

• Elements/Penta.cpp (modified) (21 diffs)
• Elements/Penta.h (modified) (1 diff)
• Gauss/GaussPenta.cpp (modified) (1 diff)
• Gauss/GaussPenta.h (modified) (2 diffs)

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

@@ -2384 +2384 @@
 
         /* 3d gaussian points: */
-        int     num_gauss,ig;
-        double* first_gauss_area_coord  =  NULL;
-        double* second_gauss_area_coord =  NULL;
-        double* third_gauss_area_coord  =  NULL;
-        double* fourth_gauss_vert_coord  =  NULL;
-        double* area_gauss_weights           =  NULL;
-        double* vert_gauss_weights           =  NULL;
-        int     ig1,ig2;
-        double  gauss_weight1,gauss_weight2;
-        double  gauss_coord[4];
-        int     order_area_gauss;
-        int     num_vert_gauss;
-        int     num_area_gauss;
-        double  gauss_weight;
-
-        /* 2d gaussian point: */
-        int     num_gauss2d;
-        double* first_gauss_area_coord2d  =  NULL;
-        double* second_gauss_area_coord2d =  NULL;
-        double* third_gauss_area_coord2d  =  NULL;
-        double* gauss_weights2d=NULL;
-        double  gauss_l1l2l3[3];
+        int     ig;
+        GaussPenta *gauss=NULL;
+        GaussTria  *gauss_tria=NULL;
 
         /* material data: */
@@ -2471 +2452 @@
         GetDofList(&doflistp,PattynApproximationEnum); //MacAyeal dof list
 
-        /*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.*/
-
-        order_area_gauss=5;
-        num_vert_gauss=5;
-
-        gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
-
         /*Retrieve all inputs we will be needing: */
         vx_input=inputs->GetInput(VxEnum);
@@ -2489 +2459 @@
 
         /* Start  looping on the number of gaussian points: */
-        for (ig1=0; ig1<num_area_gauss; ig1++){
-                for (ig2=0; ig2<num_vert_gauss; ig2++){
-
-                        /*Pick up the gaussian point: */
-                        gauss_weight1=*(area_gauss_weights+ig1);
-                        gauss_weight2=*(vert_gauss_weights+ig2);
-                        gauss_weight=gauss_weight1*gauss_weight2;
-
-                        gauss_coord[0]=*(first_gauss_area_coord+ig1); 
-                        gauss_coord[1]=*(second_gauss_area_coord+ig1);
-                        gauss_coord[2]=*(third_gauss_area_coord+ig1);
-                        gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
-
-                        /*Get strain rate from velocity: */
-                        this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input);
-                        this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,vxold_input,vyold_input);
-
-                        /*Get viscosity: */
-                        matice->GetViscosity3d(&viscosity, &epsilon[0]);
-                        matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
-
-                        /*Get B and Bprime matrices: */
-                        GetBMacAyealPattyn(&B[0][0], &xyz_list[0][0], gauss_coord);
-                        tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_coord);
-
-                        /* Get Jacobian determinant: */
-                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
-
-                        /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant 
-                          onto this scalar matrix, so that we win some computational time: */
-
-                        newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
-                        D_scalar=2*newviscosity*gauss_weight*Jdet;
-                        for (i=0;i<3;i++){
-                                D[i][i]=D_scalar;
-                        }
-
-                        /*  Do the triple product tB*D*Bprime: */
-                        TripleMultiply( &B[0][0],3,numdofp,1,
-                                                &D[0][0],3,3,0,
-                                                &Bprime[0][0],3,numdofm,0,
-                                                &Ke_gg_gaussian[0][0],0);
-
-                        /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
-                        for( i=0; i<numdofp; i++){
-                                for(j=0;j<numdofm; j++){
-                                        Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
-                                }
-                        }
-                } //for (ig2=0; ig2<num_vert_gauss; ig2++)
-        } //for (ig1=0; ig1<num_area_gauss; ig1++)
+        gauss=new GaussPenta(5,5);
+        gauss_tria=new GaussTria();
+        for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+                gauss->SynchronizeGaussTria(gauss_tria);
+
+                /*Get strain rate from velocity: */
+                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);
+
+                /*Get viscosity: */
+                matice->GetViscosity3d(&viscosity, &epsilon[0]);
+                matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
+
+                /*Get B and Bprime matrices: */
+                GetBMacAyealPattyn(&B[0][0], &xyz_list[0][0], gauss);
+                tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
+
+                /* Get Jacobian determinant: */
+                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+
+                /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant 
+                  onto this scalar matrix, so that we win some computational time: */
+
+                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
+                D_scalar=2*newviscosity*gauss->weight*Jdet;
+                for (i=0;i<3;i++) D[i][i]=D_scalar;
+
+                /*  Do the triple product tB*D*Bprime: */
+                TripleMultiply( &B[0][0],3,numdofp,1,
+                                        &D[0][0],3,3,0,
+                                        &Bprime[0][0],3,numdofm,0,
+                                        &Ke_gg_gaussian[0][0],0);
+
+                /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
+                for( i=0; i<numdofp; i++) for(j=0;j<numdofm; j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
+        } 
 
         /*Add Ke_gg and its transpose to global matrix Kgg: */
@@ -2558 +2515 @@
         delete tria->matice; delete tria;
 
-        xfree((void**)&first_gauss_area_coord);
-        xfree((void**)&second_gauss_area_coord);
-        xfree((void**)&third_gauss_area_coord);
-        xfree((void**)&fourth_gauss_vert_coord);
-        xfree((void**)&area_gauss_weights);
-        xfree((void**)&vert_gauss_weights);
-        xfree((void**)&first_gauss_area_coord2d);
-        xfree((void**)&second_gauss_area_coord2d);
-        xfree((void**)&third_gauss_area_coord2d);
-        xfree((void**)&gauss_weights2d);
         xfree((void**)&doflistm);
         xfree((void**)&doflistp);
+        delete gauss;
+        delete gauss_tria;
 }
 /*}}}*/
@@ -3113 +3062 @@
 
         /* 3d gaussian points: */
-        int     num_gauss,ig;
-        double* first_gauss_area_coord  =  NULL;
-        double* second_gauss_area_coord =  NULL;
-        double* third_gauss_area_coord  =  NULL;
-        double* fourth_gauss_vert_coord  =  NULL;
-        double* area_gauss_weights           =  NULL;
-        double* vert_gauss_weights           =  NULL;
-        int     ig1,ig2;
-        double  gauss_weight1,gauss_weight2;
-        double  gauss_coord[4];
-        int     order_area_gauss;
-        int     num_vert_gauss;
-        int     num_area_gauss;
-        double  gauss_weight;
+        int     ig;
+        GaussPenta *gauss=NULL;
 
         /* matrices: */
@@ -3164 +3101 @@
         GetDofList(&doflist);
 
-        /*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.*/
-
-        order_area_gauss=2;
-        num_vert_gauss=2;
-
-        gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
-
         /* Start  looping on the number of gaussian points: */
-        for (ig1=0; ig1<num_area_gauss; ig1++){
-                for (ig2=0; ig2<num_vert_gauss; ig2++){
-
-                        /*Pick up the gaussian point: */
-                        gauss_weight1=*(area_gauss_weights+ig1);
-                        gauss_weight2=*(vert_gauss_weights+ig2);
-                        gauss_weight=gauss_weight1*gauss_weight2;
-
-                        gauss_coord[0]=*(first_gauss_area_coord+ig1); 
-                        gauss_coord[1]=*(second_gauss_area_coord+ig1);
-                        gauss_coord[2]=*(third_gauss_area_coord+ig1);
-                        gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
-
-                        /*Get B and Bprime matrices: */
-                        GetBVert(&B[0][0], &xyz_list[0][0], gauss_coord);
-                        GetBprimeVert(&Bprime[0][0], &xyz_list[0][0], gauss_coord);
-
-                        /* Get Jacobian determinant: */
-                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
-                        DL_scalar=gauss_weight*Jdet;
-
-                        /*  Do the triple product tB*D*Bprime: */
-                        TripleMultiply( &B[0][0],1,numgrids,1,
-                                                &DL_scalar,1,1,0,
-                                                &Bprime[0][0],1,numgrids,0,
-                                                &Ke_gg_gaussian[0][0],0);
-
-                        /* Add the Ke_gg_gaussian, and optionally 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];
-                                }
-                        }       
-                } //for (ig2=0; ig2<num_vert_gauss; ig2++)
-        } //for (ig1=0; ig1<num_area_gauss; ig1++)
+        gauss=new GaussPenta(2,2);
+        for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                /*Get B and Bprime matrices: */
+                GetBVert(&B[0][0], &xyz_list[0][0], gauss);
+                GetBprimeVert(&Bprime[0][0], &xyz_list[0][0], gauss);
+
+                /* Get Jacobian determinant: */
+                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+                DL_scalar=gauss->weight*Jdet;
+
+                /*  Do the triple product tB*D*Bprime: */
+                TripleMultiply( &B[0][0],1,numgrids,1,
+                                        &DL_scalar,1,1,0,
+                                        &Bprime[0][0],1,numgrids,0,
+                                        &Ke_gg_gaussian[0][0],0);
+
+                /* Add the Ke_gg_gaussian, and optionally 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);
 
-        xfree((void**)&first_gauss_area_coord);
-        xfree((void**)&second_gauss_area_coord);
-        xfree((void**)&third_gauss_area_coord);
-        xfree((void**)&fourth_gauss_vert_coord);
-        xfree((void**)&area_gauss_weights);
-        xfree((void**)&vert_gauss_weights);
         xfree((void**)&doflist);
+        delete gauss;
 }
 /*}}}*/
@@ -3333 +3241 @@
 
         /* gaussian points: */
-        int     num_area_gauss,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;
-        double  gauss_weight,area_gauss_weight,vert_gauss_weight;
-        double  gauss_coord[4];
-        double  gauss_l1l2l3[3];
-
-        int area_order=5;
-        int num_vert_gauss=5;
+        int     ig;
+        GaussPenta *gauss=NULL;
 
         double  K[2][2]={0.0};
@@ -3392 +3289 @@
         bool onbed;
         bool shelf;
+        Input* vx_input=NULL;
+        Input* vy_input=NULL;
+        Input* vz_input=NULL;
 
         /*retrieve inputs :*/
@@ -3417 +3317 @@
         this->parameters->FindParam(&epsvel,EpsVelEnum);
 
-        /* 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: */
-        area_order=2;
-        num_vert_gauss=2;
-
-        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);
+        vx_input=inputs->GetInput(VxEnum);
+        vy_input=inputs->GetInput(VyEnum);
+        vz_input=inputs->GetInput(VzEnum);
 
         /* 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);
-
-                        /* Get Jacobian determinant: */
-                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
-
-                        /*Conduction: */
-
-                        /*Get B_conduct matrix: */
-                        GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss_coord); 
-
-                        /*Build D: */
-                        D_scalar=gauss_weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity));
-
-                        if(dt){
-                                D_scalar=D_scalar*dt;
-                        }
-
-                        D[0][0]=D_scalar; D[0][1]=0; D[0][2]=0;
-                        D[1][0]=0; D[1][1]=D_scalar; D[1][2]=0;
-                        D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar;
-
-                        /*  Do the triple product B'*D*B: */
-                        MatrixMultiply(&B_conduct[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_conduct[0][0],0);
-                        MatrixMultiply(&tBD_conduct[0][0],numdof,3,0,&B_conduct[0][0],3,numdof,0,&Ke_gaussian_conduct[0][0],0);
-
-                        /*Advection: */
-
-                        /*Get B_advec and Bprime_advec matrices: */
-                        GetBAdvec(&B_advec[0][0],&xyz_list[0][0],gauss_coord); 
-                        GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss_coord); 
-
-                        //Build the D matrix
-                        inputs->GetParameterValue(&u, gauss_coord,VxEnum);
-                        inputs->GetParameterValue(&v, gauss_coord,VyEnum);
-                        inputs->GetParameterValue(&w, gauss_coord,VzEnum);
-
-                        D_scalar=gauss_weight*Jdet;
-
-                        if(dt){
-                                D_scalar=D_scalar*dt;
-                        }
-
-                        D[0][0]=D_scalar*u;D[0][1]=0;         D[0][2]=0;
-                        D[1][0]=0;         D[1][1]=D_scalar*v;D[1][2]=0;
-                        D[2][0]=0;         D[2][1]=0;         D[2][2]=D_scalar*w;
-
-                        /*  Do the triple product B'*D*Bprime: */
-                        MatrixMultiply(&B_advec[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_advec[0][0],0);
-                        MatrixMultiply(&tBD_advec[0][0],numdof,3,0,&Bprime_advec[0][0],3,numdof,0,&Ke_gaussian_advec[0][0],0);
-
-                        /*Transient: */
-                        if(dt){
-                                GetNodalFunctionsP1(&L[0], gauss_coord);
-                                D_scalar=gauss_weight*Jdet;
-                                D_scalar=D_scalar;
-
-                                /*  Do the triple product L'*D*L: */
-                                MatrixMultiply(&L[0],numdof,1,0,&D_scalar,1,1,0,&tLD[0],0);
-                                MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian_transient[0][0],0);
-                        }
-                        else{
-                                for(i=0;i<numdof;i++){
-                                        for(j=0;j<numdof;j++){
-                                                Ke_gaussian_transient[i][j]=0;
-                                        }
-                                }
-                        }
-
-                        /*Artifficial diffusivity*/
-                        if(artdiff){
-                                /*Build K: */
-                                D_scalar=gauss_weight*Jdet/(pow(u,2)+pow(v,2)+epsvel);
-                                if(dt){
-                                        D_scalar=D_scalar*dt;
-                                }
-                                K[0][0]=D_scalar*pow(u,2);       K[0][1]=D_scalar*fabs(u)*fabs(v);
-                                K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2);
-
-                                /*Get B_artdiff: */
-                                GetBArtdiff(&B_artdiff[0][0],&xyz_list[0][0],gauss_coord); 
-
-                                /*  Do the triple product B'*K*B: */
-                                MatrixMultiply(&B_artdiff[0][0],2,numdof,1,&K[0][0],2,2,0,&tBD_artdiff[0][0],0);
-                                MatrixMultiply(&tBD_artdiff[0][0],numdof,2,0,&B_artdiff[0][0],2,numdof,0,&Ke_gaussian_artdiff[0][0],0);
-                        }
-                        else{
-                                for(i=0;i<numdof;i++){
-                                        for(j=0;j<numdof;j++){
-                                                Ke_gaussian_artdiff[i][j]=0;
-                                        }
-                                }
-                        }
-
-                        /*Add Ke_gaussian to pKe: */
-                        for(i=0;i<numdof;i++){
-                                for(j=0;j<numdof;j++){
-                                        K_terms[i][j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j]+Ke_gaussian_artdiff[i][j];
-                                }
-                        }
-                }
+        gauss=new GaussPenta(2,2);
+        for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+
+                /*Conduction: */
+
+                /*Get B_conduct matrix: */
+                GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); 
+
+                /*Build D: */
+                D_scalar=gauss->weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity));
+
+                if(dt) D_scalar=D_scalar*dt;
+
+                D[0][0]=D_scalar; D[0][1]=0; D[0][2]=0;
+                D[1][0]=0; D[1][1]=D_scalar; D[1][2]=0;
+                D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar;
+
+                /*  Do the triple product B'*D*B: */
+                MatrixMultiply(&B_conduct[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_conduct[0][0],0);
+                MatrixMultiply(&tBD_conduct[0][0],numdof,3,0,&B_conduct[0][0],3,numdof,0,&Ke_gaussian_conduct[0][0],0);
+
+                /*Advection: */
+
+                /*Get B_advec and Bprime_advec matrices: */
+                GetBAdvec(&B_advec[0][0],&xyz_list[0][0],gauss); 
+                GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); 
+
+                //Build the D matrix
+                vx_input->GetParameterValue(&u, gauss);
+                vy_input->GetParameterValue(&v, gauss);
+                vz_input->GetParameterValue(&w, gauss);
+
+                D_scalar=gauss->weight*Jdet;
+
+                if(dt) D_scalar=D_scalar*dt;
+
+                D[0][0]=D_scalar*u;D[0][1]=0;         D[0][2]=0;
+                D[1][0]=0;         D[1][1]=D_scalar*v;D[1][2]=0;
+                D[2][0]=0;         D[2][1]=0;         D[2][2]=D_scalar*w;
+
+                /*  Do the triple product B'*D*Bprime: */
+                MatrixMultiply(&B_advec[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_advec[0][0],0);
+                MatrixMultiply(&tBD_advec[0][0],numdof,3,0,&Bprime_advec[0][0],3,numdof,0,&Ke_gaussian_advec[0][0],0);
+
+                /*Transient: */
+                if(dt){
+                        GetNodalFunctionsP1(&L[0], gauss);
+                        D_scalar=gauss->weight*Jdet;
+                        D_scalar=D_scalar;
+
+                        /*  Do the triple product L'*D*L: */
+                        MatrixMultiply(&L[0],numdof,1,0,&D_scalar,1,1,0,&tLD[0],0);
+                        MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian_transient[0][0],0);
+                }
+                else{
+                        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gaussian_transient[i][j]=0;
+                }
+
+                /*Artifficial diffusivity*/
+                if(artdiff){
+                        /*Build K: */
+                        D_scalar=gauss->weight*Jdet/(pow(u,2)+pow(v,2)+epsvel);
+                        if(dt) D_scalar=D_scalar*dt;
+                        K[0][0]=D_scalar*pow(u,2);       K[0][1]=D_scalar*fabs(u)*fabs(v);
+                        K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2);
+
+                        /*Get B_artdiff: */
+                        GetBArtdiff(&B_artdiff[0][0],&xyz_list[0][0],gauss); 
+
+                        /*  Do the triple product B'*K*B: */
+                        MatrixMultiply(&B_artdiff[0][0],2,numdof,1,&K[0][0],2,2,0,&tBD_artdiff[0][0],0);
+                        MatrixMultiply(&tBD_artdiff[0][0],numdof,2,0,&B_artdiff[0][0],2,numdof,0,&Ke_gaussian_artdiff[0][0],0);
+                }
+                else{
+                        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gaussian_artdiff[i][j]=0;
+                }
+
+                /*Add Ke_gaussian to pKe: */
+                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) K_terms[i][j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j]+Ke_gaussian_artdiff[i][j];
         }
 
@@ -3553 +3419 @@
         
         /*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_weights);
-        xfree((void**)&vert_gauss_coord);
+        delete gauss;
         xfree((void**)&doflist);
 
@@ -3906 +3767 @@
 
         /* gaussian points: */
-        int     num_gauss,ig;
-        double* first_gauss_area_coord  =  NULL;
-        double* second_gauss_area_coord =  NULL;
-        double* third_gauss_area_coord  =  NULL;
-        double* fourth_gauss_vert_coord  =  NULL;
-        double* area_gauss_weights      =  NULL;
-        double* vert_gauss_weights      =  NULL;
-        double  gauss_coord[4];
-        int     order_area_gauss;
-        int     num_vert_gauss;
-        int     num_area_gauss;
-        int     ig1,ig2;
-        double  gauss_weight1,gauss_weight2;
-        double  gauss_weight;
+        int     ig;
+        GaussPenta *gauss=NULL;
 
         /* Jacobian: */
@@ -3948 +3797 @@
         GetDofList(&doflist,PattynApproximationEnum);
 
-        /*Get gaussian points and weights :*/
-        order_area_gauss=2;
-        num_vert_gauss=3;
-
-        gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
-
         /*Retrieve all inputs we will be needing: */
         thickness_input=inputs->GetInput(ThicknessEnum);
@@ -3959 +3802 @@
 
         /* Start  looping on the number of gaussian points: */
-        for (ig1=0; ig1<num_area_gauss; ig1++){
-                for (ig2=0; ig2<num_vert_gauss; ig2++){
-
-                        /*Pick up the gaussian point: */
-                        gauss_weight1=*(area_gauss_weights+ig1);
-                        gauss_weight2=*(vert_gauss_weights+ig2);
-                        gauss_weight=gauss_weight1*gauss_weight2;
-
-                        gauss_coord[0]=*(first_gauss_area_coord+ig1); 
-                        gauss_coord[1]=*(second_gauss_area_coord+ig1);
-                        gauss_coord[2]=*(third_gauss_area_coord+ig1);
-                        gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
-
-                        /*Compute thickness at gaussian point: */
-                        thickness_input->GetParameterValue(&thickness, gauss_coord);
-
-                        /*Compute slope at gaussian point: */
-                        surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord);
-
-                        /* Get Jacobian determinant: */
-                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
-
-                        /*Get nodal functions: */
-                        GetNodalFunctionsP1(l1l6, gauss_coord);
-
-                        /*Compute driving stress: */
-                        driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG();
-
-                        /*Build pe_g_gaussian vector: */
-                        for (i=0;i<numgrids;i++){
-                                for (j=0;j<NDOF2;j++){
-                                        pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss_weight*l1l6[i];
-                                }
+        gauss=new GaussPenta(2,3);
+        for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                /*Compute thickness at gaussian point: */
+                thickness_input->GetParameterValue(&thickness, gauss);
+
+                /*Compute slope at gaussian point: */
+                surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
+
+                /* Get Jacobian determinant: */
+                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+
+                /*Get nodal functions: */
+                GetNodalFunctionsP1(l1l6, gauss);
+
+                /*Compute driving stress: */
+                driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG();
+
+                /*Build pe_g_gaussian vector: */
+                for (i=0;i<numgrids;i++){
+                        for (j=0;j<NDOF2;j++){
+                                pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*l1l6[i];
                         }
-
-                        /*Add pe_g_gaussian vector to pe_g: */
-                        for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
-
-                } //for (ig2=0; ig2<num_vert_gauss; ig2++)
-        } //for (ig1=0; ig1<num_area_gauss; ig1++)
+                }
+
+                /*Add pe_g_gaussian vector to pe_g: */
+                for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
+        }
 
         /*Add pe_g to global vector pg: */
         VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
 
-        xfree((void**)&first_gauss_area_coord);
-        xfree((void**)&second_gauss_area_coord);
-        xfree((void**)&third_gauss_area_coord);
-        xfree((void**)&fourth_gauss_vert_coord);
-        xfree((void**)&area_gauss_weights);
-        xfree((void**)&vert_gauss_weights);
         xfree((void**)&doflist);
+        delete gauss;
 }
 /*}}}*/
@@ -4294 +4122 @@
 
         /* gaussian points: */
-        int     num_gauss,ig;
-        double* first_gauss_area_coord  =  NULL;
-        double* second_gauss_area_coord =  NULL;
-        double* third_gauss_area_coord  =  NULL;
-        double* fourth_gauss_vert_coord  =  NULL;
-        double* area_gauss_weights      =  NULL;
-        double* vert_gauss_weights      =  NULL;
-        double  gauss_coord[4];
-        int     order_area_gauss;
-        int     num_vert_gauss;
-        int     num_area_gauss;
-        int     ig1,ig2;
-        double  gauss_weight1,gauss_weight2;
-        double  gauss_weight;
+        int     ig;
+        GaussPenta *gauss=NULL;
 
         /* Jacobian: */
@@ -4351 +4167 @@
         GetDofList(&doflist);
 
-        /*Get gaussian points and weights :*/
-        order_area_gauss=2;
-        num_vert_gauss=2;
-
-        gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
-
         /*Retrieve all inputs we will be needing: */
         vx_input=inputs->GetInput(VxEnum);
@@ -4362 +4172 @@
 
         /* Start  looping on the number of gaussian points: */
-        for (ig1=0; ig1<num_area_gauss; ig1++){
-                for (ig2=0; ig2<num_vert_gauss; ig2++){
-
-                        /*Pick up the gaussian point: */
-                        gauss_weight1=*(area_gauss_weights+ig1);
-                        gauss_weight2=*(vert_gauss_weights+ig2);
-                        gauss_weight=gauss_weight1*gauss_weight2;
-
-                        gauss_coord[0]=*(first_gauss_area_coord+ig1); 
-                        gauss_coord[1]=*(second_gauss_area_coord+ig1);
-                        gauss_coord[2]=*(third_gauss_area_coord+ig1);
-                        gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
-
-                        /*Get velocity derivative, with respect to x and y: */
-
-                        vx_input->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord);
-                        vy_input->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss_coord);
-                        dudx=du[0];
-                        dvdy=dv[1];
-
-                        /* Get Jacobian determinant: */
-                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
-
-                        /*Get nodal functions: */
-                        GetNodalFunctionsP1(l1l6, gauss_coord);
-
-                        /*Build pe_g_gaussian vector: */
-                        for (i=0;i<numgrids;i++){
-                                pe_g_gaussian[i]=(dudx+dvdy)*Jdet*gauss_weight*l1l6[i];
-                        }
-
-                        /*Add pe_g_gaussian vector to pe_g: */
-                        for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
-
-                } //for (ig2=0; ig2<num_vert_gauss; ig2++)
-        } //for (ig1=0; ig1<num_area_gauss; ig1++)
+        gauss=new GaussPenta(2,2);
+        for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                /*Get velocity derivative, with respect to x and y: */
+                vx_input->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss);
+                vy_input->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss);
+                dudx=du[0];
+                dvdy=dv[1];
+
+                /* Get Jacobian determinant: */
+                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+
+                /*Get nodal functions: */
+                GetNodalFunctionsP1(l1l6, gauss);
+
+                /*Build pe_g_gaussian vector: */
+                for (i=0;i<numgrids;i++){
+                        pe_g_gaussian[i]=(dudx+dvdy)*Jdet*gauss->weight*l1l6[i];
+                }
+
+                /*Add pe_g_gaussian vector to pe_g: */
+                for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
+        }
 
         /*Add pe_g to global vector pg: */
         VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
 
-        xfree((void**)&first_gauss_area_coord);
-        xfree((void**)&second_gauss_area_coord);
-        xfree((void**)&third_gauss_area_coord);
-        xfree((void**)&fourth_gauss_vert_coord);
-        xfree((void**)&area_gauss_weights);
-        xfree((void**)&vert_gauss_weights);
         xfree((void**)&doflist);
+        delete gauss;
 }
 /*}}}*/
@@ -4495 +4289 @@
 
         /* gaussian points: */
-        int     num_area_gauss,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;
-        double  gauss_weight,area_gauss_weight,vert_gauss_weight;
-        double  gauss_coord[4];
-        int     area_order=2;
-        int       num_vert_gauss=3;
+        int     ig;
+        GaussPenta *gauss=NULL;
 
         double temperature_list[numgrids];
@@ -4550 +4335 @@
         Input* vy_input=NULL;
         Input* vz_input=NULL;
+        Input* temperature_input=NULL;
 
         /*retrieve inputs :*/
@@ -4574 +4360 @@
         meltingpoint=matpar->GetMeltingPoint();
 
-        /* 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: */
-        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);
         vy_input=inputs->GetInput(VyEnum);
         vz_input=inputs->GetInput(VzEnum);
+        if (dt) temperature_input=inputs->GetInput(TemperatureEnum);
 
         /* 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 */
-                        GetNodalFunctionsP1(&L[0], gauss_coord);
-
-                        /*Build deformational heating: */
-                        GetPhi(&phi, &epsilon[0], viscosity);
-
-                        /*Build pe_gaussian */
-                        scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss_weight;
-                        if(dt){
-                                scalar_def=scalar_def*dt;
-                        }
-
-                        for(i=0;i<numgrids;i++){
-                                P_terms[i]+=scalar_def*L[i];
-                        }
-
-                        /* Build transient now */
-                        if(dt){
-                                inputs->GetParameterValue(&temperature, gauss_coord,TemperatureEnum);
-                                scalar_transient=temperature*Jdet*gauss_weight;
-                                for(i=0;i<numgrids;i++){
-                                        P_terms[i]+=scalar_transient*L[i];
-                                }
-                        }
+        gauss=new GaussPenta(2,3);
+        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 */
+                GetNodalFunctionsP1(&L[0], gauss);
+
+                /*Build deformational heating: */
+                GetPhi(&phi, &epsilon[0], viscosity);
+
+                /*Build pe_gaussian */
+                scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss->weight;
+                if(dt) scalar_def=scalar_def*dt;
+
+                for(i=0;i<numgrids;i++) P_terms[i]+=scalar_def*L[i];
+
+                /* Build transient now */
+                if(dt){
+                        temperature_input->GetParameterValue(&temperature, gauss);
+                        scalar_transient=temperature*Jdet*gauss->weight;
+                        for(i=0;i<numgrids;i++) P_terms[i]+=scalar_transient*L[i];
                 }
         }
@@ -4653 +4418 @@
         }
 
-        xfree((void**)&first_gauss_area_coord);
-        xfree((void**)&second_gauss_area_coord);
-        xfree((void**)&third_gauss_area_coord);
-        xfree((void**)&vert_gauss_coord);
-        xfree((void**)&area_gauss_weights);
-        xfree((void**)&vert_gauss_weights);
         xfree((void**)&doflist);
+        delete gauss;
 
 }
@@ -5213 +4973 @@
 /*FUNCTION Penta::GetStrainRate3d{{{1*/
 void Penta::GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input){
+        /*Compute the 3d Strain Rate (6 components):
+         *
+         * epsilon=[exx eyy ezz exy exz eyz]
+         */
+
+        int i;
+
+        double epsilonvx[6];
+        double epsilonvy[6];
+        double epsilonvz[6];
+
+        /*Check that both inputs have been found*/
+        if (!vx_input || !vy_input || !vz_input){
+                ISSMERROR("Input missing. Here are the input pointers we have for vx: %p, vy: %p, vz: %p\n",vx_input,vy_input,vz_input);
+        }
+
+        /*Get strain rate assuming that epsilon has been allocated*/
+        vx_input->GetVxStrainRate3d(epsilonvx,xyz_list,gauss);
+        vy_input->GetVyStrainRate3d(epsilonvy,xyz_list,gauss);
+        vz_input->GetVzStrainRate3d(epsilonvz,xyz_list,gauss);
+
+        /*Sum all contributions*/
+        for(i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i];
+
+}
+/*}}}*/
+/*FUNCTION Penta::GetStrainRate3d{{{1*/
+void Penta::GetStrainRate3d(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input){
         /*Compute the 3d Strain Rate (6 components):
          *

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

@@ -170 +170 @@
                 void    GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);
                 void    GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
+                void    GetStrainRate3d(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
                 Penta*  GetUpperElement(void);
                 Penta*  GetLowerElement(void);

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

@@ -225 +225 @@
 }
 /*}}}*/
+/*FUNCTION GaussPenta::SynchronizeGaussTria{{{1*/
+void GaussPenta::SynchronizeGaussTria(GaussTria* gauss_tria){
+
+        gauss_tria->coord1=this->coord1;
+        gauss_tria->coord2=this->coord2;
+        gauss_tria->coord3=this->coord3;
+        gauss_tria->weight=UNDEF;
+}
+/*}}}*/

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

@@ -9 +9 @@
 /*{{{1*/
 #include "./../../shared/shared.h"
+class GaussTria;
 /*}}}*/
 
@@ -42 +43 @@
                 void GaussVertex(int iv);
                 void GaussCenter(void);
+                void SynchronizeGaussTria(GaussTria* gauss_tria);
 };
 #endif