Index: /issm/trunk/src/c/objects/Icefront.cpp =================================================================== --- /issm/trunk/src/c/objects/Icefront.cpp (revision 92) +++ /issm/trunk/src/c/objects/Icefront.cpp (revision 93) @@ -67,5 +67,4 @@ printf(" matpar_offset: %i\n",matpar_offset); printf(" matpar: \n"); - if(matpar)matpar->Echo(); printf(" element_type: %i\n",element_type); @@ -73,5 +72,4 @@ printf(" element_offset: %i\n",eid); printf(" element\n"); - if(element)element->Echo(); if (strcmp(type,"segment")==0){ @@ -319,5 +317,4 @@ } } - /* Set pe_g to 0: */ for(i=0;iGetNodes(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,5 +338,5 @@ 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,6 +391,159 @@ 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;iGetNodes(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;iterms+i)); + } + } + #endif + + /*Plug pe_g into vector: */ + VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES); + + /*Free ressources:*/ + xfree((void**)&element_nodes); + +} + #undef __FUNCT__ @@ -428,4 +578,5 @@ } } + /*Assign output pointers:*/ @@ -521,2 +672,253 @@ } + +#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;igGetNodalFunctions(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