Changeset 93

Timestamp:

04/28/09 15:18:46 (16 years ago)

Author:

Eric.Larour

Message:

New quad pressure

File:

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

issm/trunk/src/c/objects/Icefront.cpp

@@ -67 +67 @@
         printf("   matpar_offset: %i\n",matpar_offset);
         printf("   matpar: \n");
-        if(matpar)matpar->Echo();
         
         printf("   element_type: %i\n",element_type);
@@ -73 +72 @@
         printf("   element_offset: %i\n",eid);
         printf("   element\n");
-        if(element)element->Echo();
                 
         if (strcmp(type,"segment")==0){
@@ -319 +317 @@
                 }
         }
-
         /* Set pe_g to 0: */
         for(i=0;i<numdofs;i++) pe_g[i]=0.0;
@@ -327 +324 @@
         if(element_type==PentaEnum())element_nodes=(Node**)xmalloc(6*sizeof(Node*));
         element->GetNodes(element_nodes);
-                
+
         /*go through first 3 grids (all grids for tria, bed grids for penta) and identify 1st and 2nd grid: */
         for(i=0;i<3;i++){
@@ -341 +338 @@
         element->GetThicknessList(&thickness_list_element[0]);
         element->GetBedList(&bed_list_element[0]);
-        
+                
         /*Build thickness_list and bed_list: */
         thickness_list[0]=thickness_list_element[grid1];
@@ -394 +391 @@
 void Icefront::CreatePVectorDiagnosticHorizQuad( Vec pg,ParameterInputs* inputs, int analysis_type){
 
-        throw ErrorException(__FUNCT__," not implemented yet!");
-}
+        int i,j;
+
+        const int numgrids=4;
+        const int NDOF2=2;
+        const int numdofs=numgrids*NDOF2;
+        int   numberofdofspernode;
+        int   doflist[numdofs];
+
+        int    fill;
+        double xyz_list[numgrids][3];
+        double xyz_list_quad[numgrids+1][3]; //center node
+
+        double thickness_list[6]; //from penta element
+        double thickness_list_quad[6]; //selection of 4 thickness from penta elements
+
+        double bed_list[6]; //from penta element
+        double bed_list_quad[6]; //selection of 4 beds from penta elements
+
+        /*Objects: */
+        double  pe_g[numdofs];
+        Matpar* matpar=NULL;
+        Node**  element_nodes=NULL;
+
+        /*quad grids: */
+        int grid1,grid2,grid3,grid4;
+
+        /*normals: */
+        double normal1[3];
+        double normal2[3];
+        double normal3[3];
+        double normal4[3];
+        double normal_norm;
+        double v15[3];
+        double v25[3];
+        double v35[3];
+        double v45[3];
+                
+
+        /*check icefront is associated to a pentaelem: */
+        if(element_type!=PentaEnum()){
+                throw ErrorException(__FUNCT__," quad icefront is associated to a TriaElem!");
+        }
+        /*Recover material and fill parameters: */
+        matpar=element->GetMatPar();
+        fill=element->GetShelf();
+
+        /* Set pe_g to 0: */
+        for(i=0;i<numdofs;i++) pe_g[i]=0.0;
+
+        /* Get dof list and node coordinates: */
+        GetDofList(&doflist[0],&numberofdofspernode);
+        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+        
+        //Identify which grids are comprised in the quad: 
+        if(element_type==PentaEnum())element_nodes=(Node**)xmalloc(6*sizeof(Node*));
+        element->GetNodes(element_nodes);
+
+        grid1=-1; grid2=-1; grid3=-1; grid4=-1;
+        for(i=0;i<6;i++){
+                if (nodes[0]==element_nodes[i])grid1=i;
+                if (nodes[1]==element_nodes[i])grid2=i;
+                if (nodes[2]==element_nodes[i])grid3=i;
+                if (nodes[3]==element_nodes[i])grid4=i;
+        }
+        
+        if((grid1==-1) || (grid2==-1)|| (grid3==-1)||(grid4==-1)){
+                throw ErrorException(__FUNCT__,"could not find element grids corresponding to quad icefront!");
+        }
+
+        /*Build new xyz, bed and thickness lists for grids 1 to 4: */
+        for(i=0;i<4;i++){
+                for(j=0;j<3;j++){
+                        xyz_list_quad[i][j]=xyz_list[i][j];
+                }
+        }
+
+        element->GetThicknessList(&thickness_list[0]);
+        element->GetBedList(&bed_list[0]);
+
+        thickness_list_quad[0]=thickness_list[grid1];
+        thickness_list_quad[1]=thickness_list[grid2];
+        thickness_list_quad[2]=thickness_list[grid3];
+        thickness_list_quad[3]=thickness_list[grid4];
+
+        bed_list_quad[0]=bed_list[grid1];
+        bed_list_quad[1]=bed_list[grid2];
+        bed_list_quad[2]=bed_list[grid3];
+        bed_list_quad[3]=bed_list[grid4];
+
+        //Create a new grid in the midle of the quad and add it at the end of the list
+        xyz_list_quad[4][0] = (xyz_list_quad[0][0]+xyz_list_quad[1][0]+xyz_list_quad[2][0]+xyz_list_quad[3][0])/4.0;
+        xyz_list_quad[4][1] = (xyz_list_quad[0][1]+xyz_list_quad[1][1]+xyz_list_quad[2][1]+xyz_list_quad[3][1])/4.0;
+        xyz_list_quad[4][2] = (xyz_list_quad[0][2]+xyz_list_quad[1][2]+xyz_list_quad[2][2]+xyz_list_quad[3][2])/4.0;
+
+        //Compute four normals (the quad is divided into four triangles)
+        v15[0]=xyz_list_quad[0][0]-xyz_list_quad[4][0];
+        v15[1]=xyz_list_quad[0][1]-xyz_list_quad[4][1];
+        v15[2]=xyz_list_quad[0][2]-xyz_list_quad[4][2];
+        
+        v25[0]=xyz_list_quad[1][0]-xyz_list_quad[4][0];
+        v25[1]=xyz_list_quad[1][1]-xyz_list_quad[4][1];
+        v25[2]=xyz_list_quad[1][2]-xyz_list_quad[4][2];
+        
+        v35[0]=xyz_list_quad[2][0]-xyz_list_quad[4][0];
+        v35[1]=xyz_list_quad[2][1]-xyz_list_quad[4][1];
+        v35[2]=xyz_list_quad[2][2]-xyz_list_quad[4][2];
+        
+        v45[0]=xyz_list_quad[3][0]-xyz_list_quad[4][0];
+        v45[1]=xyz_list_quad[3][1]-xyz_list_quad[4][1];
+        v45[2]=xyz_list_quad[3][2]-xyz_list_quad[4][2];
+
+        cross(normal1,v15,v25); 
+        cross(normal2,v25,v35);
+        cross(normal3,v35,v45);
+        cross(normal4,v45,v15);
+
+        normal_norm=norm(normal1);normal1[0]=normal1[0]/normal_norm;normal1[1]=normal1[1]/normal_norm;normal1[2]=normal1[2]/normal_norm;
+        normal_norm=norm(normal2);normal2[0]=normal2[0]/normal_norm;normal2[1]=normal2[1]/normal_norm;normal2[2]=normal2[2]/normal_norm;
+        normal_norm=norm(normal3);normal3[0]=normal3[0]/normal_norm;normal3[1]=normal3[1]/normal_norm;normal3[2]=normal3[2]/normal_norm;
+        normal_norm=norm(normal4);normal4[0]=normal4[0]/normal_norm;normal4[1]=normal4[1]/normal_norm;normal4[2]=normal4[2]/normal_norm;
+
+
+        //Compute load contribution for this quad:
+        QuadPressureLoad(pe_g,matpar->GetRhoWater(),matpar->GetRhoIce(),matpar->GetG(),thickness_list_quad,bed_list_quad,normal1,normal2,normal3,normal4,&xyz_list_quad[0][0],fill);
+
+
+        #ifdef _DEBUGELEMENTS_
+        if(my_rank==RANK && eid==ELID){ 
+                printf("\nicefront load\n");
+                printf("grids %i %i %i %i\n",grid1,grid2,grid3,grid4);
+                printf("rho_water %g\n",rho_water);
+                printf("rho_ice %g\n",rho_ice);
+                printf("gravity %g\n",gravity);
+                printf("thickness_list (%g,%g,%g,%g)\n",thickness_list[0],thickness_list[1],thickness_list[2],thickness_list[3]);
+                printf("bed_list (%g,%g,%g,%g)\n",bed_list[0],bed_list[1],bed_list[2],bed_list[3]);
+                printf("normal1 (%g,%g,%g)\n",normal1[0],normal1[1],normal1[2]);
+                printf("normal2 (%g,%g,%g)\n",normal2[0],normal2[1],normal2[2]);
+                printf("normal3 (%g,%g,%g)\n",normal3[0],normal3[1],normal3[2]);
+                printf("normal4 (%g,%g,%g)\n",normal4[0],normal4[1],normal4[2]);
+                printf("fill %i\n",fill);
+                printf("pe_g->terms\n");
+                for(i=0;i<numgridsload*NDOF2;i++){
+                        printf("%g ",*(pe_g->terms+i));
+                }
+        }
+        #endif
+
+        /*Plug pe_g into vector: */
+        VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
+
+        /*Free ressources:*/
+        xfree((void**)&element_nodes);
+
+}
+
 
 #undef __FUNCT__ 
@@ -428 +578 @@
                 }
         }
+
 
         /*Assign output pointers:*/
@@ -521 +672 @@
 }
 
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Icefront::QuadPressureLoad"
+void Icefront::QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, 
+                                              double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list, int fill){
+        
+        
+        //The quad is divided into four triangles tria1 tria2 tria3 and tria4 as follows
+        //
+        //   grid 4 +-----------------+ grid 3
+        //          |\2            1 /|
+        //          |1\    tria3    /2|
+        //          |  \           /  |
+        //          |   \         /   |
+        //          |    \       /    |
+        //          |     \     /     |
+        //          |      \ 3 /      |
+        //          |tria4  \ / 3     |
+        //          |      3 \grid5   |
+        //          |       / \       | 
+        //          |      / 3 \ tria2|
+        //          |     /     \     |
+        //          |    /       \    |
+        //          |   /         \   |
+        //          |  /   tria1   \  |
+        //          |2/1           2\1|
+        //   grid1  +-----------------+ grid 2
+        //
+        //
+
+        int      i,j;
+
+        double nx[4];
+        double ny[4];
+
+        /* 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* gauss_weights           =  NULL;
+        double  gauss_weight;
+        double  gauss_coord[3];
+
+        double  surface_list[5];
+        double  x[5];
+        double  y[5];
+        double  z[5];
+        double l12,l14,l15,l23,l25,l34,l35,l45;
+        double cos_theta_triangle1,cos_theta_triangle2,cos_theta_triangle3,cos_theta_triangle4;
+
+        double xyz_tria[12][2];
+        double l1l2l3_tria[3];
+        double r_tria,s_tria;
+        double r_quad[4];
+        double s_quad[4];
+        double l1l4_tria[4][4];
+
+        double J[4];
+        double z_g[4];
+        double ice_pressure_tria[4];
+        double air_pressure_tria;
+        double water_level_above_g_tria;
+        double water_pressure_tria;
+        double pressure_tria[4];
+
+        /*To use tria functionalities: */
+        Tria* tria=NULL;
+
+        /*Initialize fake tria: */
+        tria=new Tria();
+
+        //Build the four normal vectors 
+        nx[0]=normal1[0]; nx[1]=normal2[0]; nx[2]=normal3[0]; nx[3]=normal4[0];
+        ny[0]=normal1[1]; ny[1]=normal2[1]; ny[2]=normal3[1]; ny[3]=normal4[1];
+
+        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+        
+        //Recover the surface of the four nodes
+        for(i=0;i<4;i++){
+                surface_list[i]=thickness_list[i]+bed_list[i];
+        }
+        //Add surface sor the fifth point (average of the surfaces)
+        surface_list[4]=(surface_list[0]+surface_list[1]+surface_list[2]+surface_list[3])/4.0;
+
+        //Recover grid coordinates
+        for(i=0;i<5;i++){
+                x[i]=*(xyz_list+3*i+0);
+                y[i]=*(xyz_list+3*i+1);
+                z[i]=*(xyz_list+3*i+2);
+        }
+        
+        //Build triangles in a 2D plan before using reference elements
+        l12=pow((x[1]-x[0]),2)+pow((y[1]-y[0]),2)+pow((z[1]-z[0]),2);
+        l14=pow((x[3]-x[0]),2)+pow((y[3]-y[0]),2)+pow((z[3]-z[0]),2);
+        l15=pow((x[4]-x[0]),2)+pow((y[4]-y[0]),2)+pow((z[4]-z[0]),2);
+        l23=pow((x[2]-x[1]),2)+pow((y[2]-y[1]),2)+pow((z[2]-z[1]),2);
+        l25=pow((x[4]-x[1]),2)+pow((y[4]-y[1]),2)+pow((z[4]-z[1]),2);
+        l34=pow((x[3]-x[2]),2)+pow((y[3]-y[2]),2)+pow((z[3]-z[2]),2);
+        l35=pow((x[4]-x[2]),2)+pow((y[4]-y[2]),2)+pow((z[4]-z[2]),2);
+        l45=pow((x[4]-x[3]),2)+pow((y[4]-y[3]),2)+pow((z[4]-z[3]),2);
+
+        cos_theta_triangle1=(l15+l12-l25)/(2*sqrt(l12*l15));
+        cos_theta_triangle2=(l25+l23-l35)/(2*sqrt(l23*l25));
+        cos_theta_triangle3=(l35+l34-l45)/(2*sqrt(l34*l35));
+        cos_theta_triangle4=(l45+l14-l15)/(2*sqrt(l14*l45));
+
+        //First triangle : nodes 1, 2 and 5
+        xyz_tria[0][0]=0;
+        xyz_tria[0][1]=0;
+        xyz_tria[1][0]=sqrt(l12);
+        xyz_tria[1][1]=0;
+        xyz_tria[2][0]=cos_theta_triangle1*sqrt(l15);
+        xyz_tria[2][1]=sqrt(1-pow(cos_theta_triangle1,2))*sqrt(l15);
+
+        //Second triangle : nodes 2, 3 and 5
+        xyz_tria[3][0]=0;
+        xyz_tria[3][1]=0;
+        xyz_tria[4][0]=sqrt(l23);
+        xyz_tria[4][1]=0;
+        xyz_tria[5][0]=cos_theta_triangle2*sqrt(l25);
+        xyz_tria[5][1]=sqrt(1-pow(cos_theta_triangle2,2))*sqrt(l25);
+        
+        //Third triangle : nodes 3, 4 and 5
+        xyz_tria[6][0]=0;
+        xyz_tria[6][1]=0;
+        xyz_tria[7][0]=sqrt(l34);
+        xyz_tria[7][1]=0;
+        xyz_tria[8][0]=cos_theta_triangle3*sqrt(l35);
+        xyz_tria[8][1]=sqrt(1-pow(cos_theta_triangle3,2))*sqrt(l35);
+
+        //Fourth triangle : nodes 4, 1 and 5
+        xyz_tria[9][0]=0;
+        xyz_tria[9][1]=0;
+        xyz_tria[10][0]=sqrt(l14);
+        xyz_tria[10][1]=0;
+        xyz_tria[11][0]=cos_theta_triangle4*sqrt(l45);
+        xyz_tria[11][1]=sqrt(1-pow(cos_theta_triangle4,2))*sqrt(l45);
+
+        
+
+        //Start  looping on the triangle's gaussian points: 
+        for(ig=0;ig<num_gauss;ig++){
+
+                /*Pick up the gaussian point: */
+                gauss_weight=*(gauss_weights+ig);
+                gauss_coord[0]=*(first_gauss_area_coord+ig);
+                gauss_coord[1]=*(second_gauss_area_coord+ig);
+                gauss_coord[2]=*(third_gauss_area_coord+ig);
+        
+                tria->GetNodalFunctions(l1l2l3_tria,gauss_coord);
+
+                //Get the coordinates of gauss point for each triangle in penta/quad
+                r_tria=gauss_coord[1]-gauss_coord[0];
+                s_tria=-3/sqrt(3)*(gauss_coord[0]+gauss_coord[1]-2/3);
+
+                //Coordinates of gauss points in the reference triangle
+                r_quad[0]=r_tria;
+                s_quad[0]=1/sqrt(3)*s_tria-2/3;
+                r_quad[1]=-1/sqrt(3)*s_tria+2/3;
+                s_quad[1]=r_tria;
+                r_quad[2]=-r_tria;
+                s_quad[2]=-1/sqrt(3)*s_tria+2/3;
+                r_quad[3]=1/sqrt(3)*s_tria-2/3;
+                s_quad[3]=-r_tria;
+
+                //Get the nodal function of the quad for the gauss points of each triangle
+                for(i=0;i<4;i++){
+                        l1l4_tria[i][0]=1.0/4.0*(
+                                        (r_quad[i]-1.0)*(s_quad[i]-1.0) 
+                                        ); 
+                        
+                        l1l4_tria[i][1]=1.0/4.0*(
+                                         -(r_quad[i]+1.0)*(s_quad[i]-1.0)
+                                        ); 
+
+                        l1l4_tria[i][2]=1.0/4.0*(
+                                        (r_quad[i]+1.0)*(s_quad[i]+1.0) 
+                                        );
+                        
+                        l1l4_tria[i][3]=1.0/4.0*(
+                                         -(r_quad[i]-1.0)*(s_quad[i]+1.0)
+                                        );
+
+                } //for(i=0;i<4;i++)
+                
+                
+                //Compute jacobian of each triangle
+                for (i=0;i<4;i++){
+                        double complete_list[3][3]; //a third coordinate is needed for the jacobian determinant calculation, here it is zero
+                        for(j=0;j<3;j++){
+                                complete_list[j][0]=xyz_tria[3*i+j][0];
+                                complete_list[j][1]=xyz_tria[3*i+j][1];
+                                complete_list[j][2]=0;
+                        }
+                        tria->GetJacobianDeterminant(&J[i],&complete_list[0][0],l1l2l3_tria);
+                }
+
+                //Calculation of the z coordinate for the gaussian point ig for each triangle
+                z_g[0]=z[0]*l1l2l3_tria[0]+z[1]*l1l2l3_tria[1]+z[4]*l1l2l3_tria[2];
+                z_g[1]=z[1]*l1l2l3_tria[0]+z[2]*l1l2l3_tria[1]+z[4]*l1l2l3_tria[2];
+                z_g[2]=z[2]*l1l2l3_tria[0]+z[3]*l1l2l3_tria[1]+z[4]*l1l2l3_tria[2];
+                z_g[3]=z[3]*l1l2l3_tria[0]+z[0]*l1l2l3_tria[1]+z[4]*l1l2l3_tria[2];
+                
+                //Loop on the triangles
+                for(i=0;i<4;i++){
+
+                        //Loop on the grids of the quad
+                        //Calculate the ice pressure
+                        for(j=0;j<4;j++){
+                                ice_pressure_tria[j]=gravity*rho_ice*(surface_list[j]-z_g[i]);
+                        }
+                        air_pressure_tria=0;
+
+                        //Now deal with water pressure: 
+                        if(fill==1){ //icefront ends in water
+                                water_level_above_g_tria=min(0,z_g[i]);//0 if the gaussian point is above water level
+                                water_pressure_tria=rho_water*gravity*water_level_above_g_tria;
+                        }
+                        else if(fill==0){
+                                water_pressure_tria=0;
+                        }
+                        else{
+                                throw ErrorException(__FUNCT__,"QuadPressureLoad error message: unknow fill type for icefront boundary condition");
+                        }
+
+                        //Add pressure from triangle i
+                        //Loop on the grids of the quad
+                        for(j=0;j<4;j++){
+                                pressure_tria[j] = ice_pressure_tria[j] + water_pressure_tria + air_pressure_tria;
+                        }
+
+
+                        pe_g[0]+= J[i]*gauss_weight* pressure_tria[0]*l1l4_tria[i][0]*nx[i];
+                        pe_g[1]+= J[i]*gauss_weight* pressure_tria[0]*l1l4_tria[i][0]*ny[i];
+                        pe_g[2]+= J[i]*gauss_weight* pressure_tria[1]*l1l4_tria[i][1]*nx[i];
+                        pe_g[3]+= J[i]*gauss_weight* pressure_tria[1]*l1l4_tria[i][1]*ny[i];
+                        pe_g[4]+= J[i]*gauss_weight* pressure_tria[2]*l1l4_tria[i][2]*nx[i];
+                        pe_g[5]+= J[i]*gauss_weight* pressure_tria[2]*l1l4_tria[i][2]*ny[i];
+                        pe_g[6]+= J[i]*gauss_weight* pressure_tria[3]*l1l4_tria[i][3]*nx[i];
+                        pe_g[7]+= J[i]*gauss_weight* pressure_tria[3]*l1l4_tria[i][3]*ny[i];
+
+                        
+
+                } //for(i=0;i<4;i++)
+        } //for(ig=0;ig<num_gauss;ig++)
+
+        /*Delete fake tria: */
+        delete tria;
+}