Changeset 5726

Timestamp:

09/09/10 15:02:52 (15 years ago)

Author:

Mathieu Morlighem

Message:

Use GaussPenta for Hutter 3d

Location:

issm/trunk/src/c/objects/Elements

Files:

• Penta.cpp (modified) (5 diffs)
• PentaRef.cpp (modified) (1 diff)
• PentaRef.h (modified) (1 diff)

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

@@ -2740 +2740 @@
         int     ig;
         GaussPenta *gauss=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;
-
         double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
         double  viscosity;
@@ -3416 +3408 @@
         double    pe_g[numdofs]={0.0};
         double    xyz_list[NUMVERTICES][3];
+        double    xyz_list_segment[2][3];
         double    z_list[NUMVERTICES];
         double    z_segment[2];
         double    slope2,constant_part;
-        double    gauss[numdofs][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
         int       node0,node1;
-        BeamRef*  beam=NULL;
-
-        /*gaussian points: */
-        int      num_gauss;
-        double*  segment_gauss_coord=NULL;
-        double   gauss_coord;
-        double*  gauss_weights=NULL;
-        double   gauss_weight;
+        GaussPenta* gauss=NULL;
         int      ig;
         double   slope[2];
@@ -3444 +3429 @@
         int  connectivity[2];
 
-        /*Inputs*/
-        Input* thickness_input;
-        Input* surface_input;
-        Input* slopex_input;
-        Input* slopey_input;
-
         /*recover some inputs: */
         inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
@@ -3469 +3448 @@
         B=matice->GetB();
 
-        //Initilaize beam
-        beam=new BeamRef();
-
         //recover extra inputs
-        thickness_input=inputs->GetInput(ThicknessEnum);
-        surface_input=inputs->GetInput(SurfaceEnum);
-        slopex_input=inputs->GetInput(SurfaceSlopeXEnum);
-        slopey_input=inputs->GetInput(SurfaceSlopeYEnum);
+        Input* thickness_input=inputs->GetInput(ThicknessEnum);
+        Input* surface_input=inputs->GetInput(SurfaceEnum);
+        Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum);
+        Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum);
 
         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
         for(i=0;i<NUMVERTICES;i++)z_list[i]=xyz_list[i][2];
 
-        //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
-        num_gauss=3;
-        gaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
-
         /*Loop on the three segments*/
         for(i=0;i<3;i++){
-                slopex_input->GetParameterValue(&slope[0], &gauss[i][0]);
-                slopey_input->GetParameterValue(&slope[1], &gauss[i][0]);
-                surface_input->GetParameterValue(&surface, &gauss[i][0]);
-                thickness_input->GetParameterValue(&thickness, &gauss[i][0]);
-
-                //compute slope2 slopex and slopey
-                slope2=pow(slope[0],2)+pow(slope[1],2);
-
-                //%compute constant_part
-                constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));
-
-                z_segment[0]=z_list[i];
-                z_segment[1]=z_list[i+3];
 
                 node0=i;
                 node1=i+3;
 
+                for(j=0;j<3;j++){
+                        xyz_list_segment[0][j]=xyz_list[node0][j];
+                        xyz_list_segment[1][j]=xyz_list[node1][j];
+                }
+
                 connectivity[0]=nodes[node0]->GetConnectivity();
                 connectivity[1]=nodes[node1]->GetConnectivity();
 
                 /*Loop on the Gauss points: */
-                for(ig=0;ig<num_gauss;ig++){
-
-                        //Pick up the gaussian point and its weight:
-                        gauss_weight=gauss_weights[ig];
-                        gauss_coord=segment_gauss_coord[ig];
-
-                        beam->GetParameterValue(&z_g, &z_segment[0],gauss_coord);
-
-                        //Get Jacobian determinant:
-                        beam->GetJacobianDeterminant(&Jdet, &z_segment[0],gauss_coord);
+                gauss=new GaussPenta(node0,node1,3);
+                for(ig=gauss->begin();ig<gauss->end();ig++){
+                        gauss->GaussPoint(ig);
+
+                        slopex_input->GetParameterValue(&slope[0],gauss);
+                        slopey_input->GetParameterValue(&slope[1],gauss);
+                        surface_input->GetParameterValue(&surface,gauss);
+                        thickness_input->GetParameterValue(&thickness,gauss);
+
+                        slope2=pow(slope[0],2)+pow(slope[1],2);
+                        constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));
+
+                        PentaRef::GetParameterValue(&z_g,&z_list[0],gauss);
+                        GetSegmentJacobianDeterminant(&Jdet,&xyz_list_segment[0][0],gauss);
 
                         if (onsurface){
-                                for(j=0;j<NDOF2;j++){
-                                        pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight/(double)connectivity[1];
-                                }
+                                for(j=0;j<NDOF2;j++) pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight/(double)connectivity[1];
                         }
                         else{//connectivity is too large, should take only half on it
-                                for(j=0;j<NDOF2;j++){
-                                        pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight*2/(double)connectivity[1];
-                                }
+                                for(j=0;j<NDOF2;j++) pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight*2/(double)connectivity[1];
                         }
                 }
+                delete gauss;
 
                 //Deal with lower surface
                 if (onbed){
-
-                        //compute ub
                         constant_part=-1.58*pow((double)10.0,-(double)10.0)*rho_ice*gravity*thickness;
                         ub=constant_part*slope[0];
                         vb=constant_part*slope[1];
 
-                        //Add to pe: 
                         pe_g[2*node0]+=ub/(double)connectivity[0];
                         pe_g[2*node0+1]+=vb/(double)connectivity[0];
@@ -3549 +3511 @@
 
         /*Clean up */
-        delete beam;
-        xfree((void**)&gauss_weights);
-        xfree((void**)&segment_gauss_coord);
         xfree((void**)&doflist);
 }

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

@@ -1699 +1699 @@
 }
 /*}}}*/
+/*FUNCTION PentaRef::GetSegmentJacobianDeterminant{{{1*/
+void PentaRef::GetSegmentJacobianDeterminant(double*  Jdet, double* xyz_list,GaussPenta* gauss){
+        /*The Jacobian determinant is constant over the element, discard the gaussian points. 
+         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+
+        double x1,x2,y1,y2,z1,z2;
+
+        x1=*(xyz_list+3*0+0);
+        y1=*(xyz_list+3*0+1);
+        z1=*(xyz_list+3*0+2);
+        x2=*(xyz_list+3*1+0);
+        y2=*(xyz_list+3*1+1);
+        z2=*(xyz_list+3*1+2);
+
+        *Jdet=1.0/2.0*sqrt(pow(x2-x1,2.) + pow(y2-y1,2.) + pow(z2-z1,2.));
+        if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
+
+}
+/*}}}*/
 /*FUNCTION PentaRef::GetJacobianInvert {{{1*/
 void PentaRef::GetJacobianInvert(double* Jinv, double* xyz_list,GaussPenta* gauss){

issm/trunk/src/c/objects/Elements/PentaRef.h

@@ -59 +59 @@
                 void GetJacobianDeterminant(double*  Jdet, double* xyz_list,GaussPenta* gauss);
                 void GetTriaJacobianDeterminant(double*  Jdet, double* xyz_list,GaussPenta* gauss);
+                void GetSegmentJacobianDeterminant(double*  Jdet, double* xyz_list,GaussPenta* gauss);
                 void GetJacobianInvert(double*  Jinv, double* xyz_list,GaussPenta* gauss);
                 void GetBMacAyealPattyn(double* B, double* xyz_list, GaussPenta* gauss);