Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5807) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5808) @@ -4300,4 +4300,21 @@ } /*}}}*/ +/*FUNCTION Penta::GetZcoord {{{1*/ +double Penta::GetZcoord(GaussPenta* gauss){ + + int i; + double xyz_list[NUMVERTICES][3]; + double z_list[NUMVERTICES]; + double z; + + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + for(i=0;i=0 && index1<6); + ISSMASSERT(index2>=0 && index2<6); + ISSMASSERT(index3>=0 && index3<6); + ISSMASSERT(index4>=0 && index4<6); + + l1l4[0]=BasisFunctions[index1]; + l1l4[1]=BasisFunctions[index2]; + l1l4[2]=BasisFunctions[index3]; + l1l4[3]=BasisFunctions[index4]; + +} +/*}}}*/ +/*FUNCTION PentaRef::GetQuadJacobianDeterminant{{{1*/ +void PentaRef::GetQuadJacobianDeterminant(double* Jdet,double xyz_list[4][3],GaussPenta* gauss){ + /*This routine returns the values of the nodal functions at the gaussian point.*/ + + double x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4; + + x1=xyz_list[0][0]; + y1=xyz_list[0][1]; + z1=xyz_list[0][2]; + x2=xyz_list[1][0]; + y2=xyz_list[1][1]; + z2=xyz_list[1][2]; + x3=xyz_list[2][0]; + y3=xyz_list[2][1]; + z3=xyz_list[2][2]; + x4=xyz_list[3][0]; + y4=xyz_list[3][1]; + z4=xyz_list[3][2]; + + /*Jdet = (Area of the trapezoid)/(Area trapezoid ref) with AreaRef = 4*/ + /*Area of a trabezoid = altitude * (base1 + base2)/2 */ + *Jdet= pow(pow(x2-x1,2.) + pow(y2-y1,2.),0.5) * (z4-z1 + z3-z2)/8; + if(*Jdet<0) ISSMERROR("negative jacobian determinant!"); + } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/PentaRef.h =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaRef.h (revision 5807) +++ /issm/trunk/src/c/objects/Elements/PentaRef.h (revision 5808) @@ -29,4 +29,6 @@ void GetNodalFunctionsP1DerivativesReference(double* dl1dl6,GaussPenta* gauss); void GetNodalFunctionsMINIDerivativesReference(double* dl1dl7,GaussPenta* gauss); + void GetQuadNodalFunctions(double* l1l4,GaussPenta* gauss,int index1,int index2,int index3,int index4); + void GetQuadJacobianDeterminant(double* Jdet, double xyz_list[4][3],GaussPenta* gauss); void GetJacobian(double* J, double* xyz_list,GaussPenta* gauss); void GetJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss); Index: /issm/trunk/src/c/objects/Elements/TriaRef.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/TriaRef.cpp (revision 5807) +++ /issm/trunk/src/c/objects/Elements/TriaRef.cpp (revision 5808) @@ -265,9 +265,8 @@ } /*}}}*/ -/*FUNCTION TriaRef::GetSegmentJacobian{{{1*/ -void TriaRef::GetSegmentJacobian(double* J, double* xyz_list,GaussTria* gauss){ - /*The Jacobian is constant over the element, discard the gaussian points. - * J is assumed to have been allocated.*/ - +/*FUNCTION TriaRef::GetSegmentJacobianDeterminant{{{1*/ +void TriaRef::GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss){ + /*The Jacobian determinant is constant over the element, discard the gaussian points. + * J is assumed to have been allocated*/ double x1,y1,x2,y2; @@ -277,13 +276,5 @@ y2=*(xyz_list+3*1+1); - *J=1.0/2.0*sqrt(pow(x2-x1,2.) + pow(y2-y1,2.)); -} -/*}}}*/ -/*FUNCTION TriaRef::GetSegmentJacobianDeterminant{{{1*/ -void TriaRef::GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss){ - /*The Jacobian determinant is constant over the element, discard the gaussian points. - * J is assumed to have been allocated*/ - - GetSegmentJacobian(Jdet,xyz_list,gauss); + *Jdet=1.0/2.0*sqrt(pow(x2-x1,2.) + pow(y2-y1,2.)); if(*Jdet<0) ISSMERROR("negative jacobian determinant!"); Index: /issm/trunk/src/c/objects/Elements/TriaRef.h =================================================================== --- /issm/trunk/src/c/objects/Elements/TriaRef.h (revision 5807) +++ /issm/trunk/src/c/objects/Elements/TriaRef.h (revision 5808) @@ -31,5 +31,4 @@ void GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof); void GetJacobian(double* J, double* xyz_list,GaussTria* gauss); - void GetSegmentJacobian(double* J, double* xyz_list,GaussTria* gauss); void GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss); void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss); Index: /issm/trunk/src/c/objects/Gauss/GaussPenta.cpp =================================================================== --- /issm/trunk/src/c/objects/Gauss/GaussPenta.cpp (revision 5807) +++ /issm/trunk/src/c/objects/Gauss/GaussPenta.cpp (revision 5808) @@ -179,34 +179,34 @@ /*Quads: get the gauss points using the product of two line rules */ - if(index1==0 && index2==1 && index3==3 && index4==4){ + if(index1==0 && index2==1 && index3==4 && index4==3){ for(i=0;iGetG(); - GetNormal(&normal[0],xyz_list); + GetSegmentNormal(&normal[0],xyz_list); index1=tria->GetNodeIndex(nodes[0]); @@ -542,235 +542,158 @@ void Icefront::CreatePVectorDiagnosticPattyn( Vec pg){ - int i,j; - - const int numnodes=NUMVERTICESQUA; - const int numdofs=numnodes*NDOF2; - int* doflist=NULL; - - double xyz_list[numnodes][3]; - double xyz_list_quad[numnodes+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]={0.0}; - int element_type; - - /*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]; - - - element_type=element->Enum(); - - /*check icefront is associated to a pentaelem: */ - if(element_type!=PentaEnum){ - ISSMERROR(" quad icefront is associated to a TriaElem!"); - } + /*Constants*/ + const int numdofs = NUMVERTICESQUA *NDOF2; + + /*Intermediaries*/ + int i,j,ig,index1,index2,index3,index4; + int fill; + double surface,pressure,ice_pressure,rho_water,rho_ice,gravity; + double water_pressure,air_pressure; + double Jdet,z_g; + double xyz_list[NUMVERTICESQUA][3]; + double normal[3]; + double pe_g[numdofs] = {0.0}; + double l1l4[4]; + int *doflist = NULL; + GaussPenta *gauss = NULL; + Penta *penta = NULL; /* Get dof list and node coordinates: */ + penta=(Penta*)element; GetDofList(&doflist,PattynApproximationEnum,GsetEnum); - GetVerticesCoordinates(&xyz_list[0][0], nodes, numnodes); - - //Identify which grids are comprised in the quad: - grid1=element->GetNodeIndex(nodes[0]); - grid2=element->GetNodeIndex(nodes[1]); - grid3=element->GetNodeIndex(nodes[2]); - grid4=element->GetNodeIndex(nodes[3]); - - /*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->GetParameterListOnVertices(&thickness_list[0],ThicknessEnum); - element->GetParameterListOnVertices(&bed_list[0],BedEnum); - - 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]); - - /*Plug pe_g into vector: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICESQUA); + + /*Recover inputs and parameters */ + Input* surface_input =penta->inputs->GetInput(SurfaceEnum); ISSMASSERT(surface_input); + inputs->GetParameterValue(&fill,FillEnum); + rho_water=matpar->GetRhoWater(); + rho_ice =matpar->GetRhoIce(); + gravity =matpar->GetG(); + GetQuadNormal(&normal[0],xyz_list); + + /*Identify which nodes are in the quad: */ + index1=element->GetNodeIndex(nodes[0]); + index2=element->GetNodeIndex(nodes[1]); + index3=element->GetNodeIndex(nodes[2]); + index4=element->GetNodeIndex(nodes[3]); + + /* Start looping on the number of gaussian points: */ + double zmax=xyz_list[0][2]; for(i=1;izmax) zmax=xyz_list[i][2]; + double zmin=xyz_list[0][2]; for(i=1;i0 && zmin<0) gauss=new GaussPenta(index1,index2,index3,index4,3,10); //refined in vertical because of the sea level discontinuity + else gauss=new GaussPenta(index1,index2,index3,index4,3,3); + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + penta->GetQuadNodalFunctions(l1l4,gauss,index1,index2,index3,index4); + penta->GetQuadJacobianDeterminant(&Jdet,xyz_list,gauss); + z_g=penta->GetZcoord(gauss); + surface_input->GetParameterValue(&surface,gauss); + + switch(fill){ + case WaterEnum: + water_pressure=rho_water*gravity*min(0,z_g);//0 if the gaussian point is above water level + break; + case AirEnum: + water_pressure=0; + break; + default: + ISSMERROR("fill type %s not supported yet",EnumToString(fill)); + } + ice_pressure=rho_ice*gravity*(surface-z_g); + air_pressure=0; + pressure = ice_pressure + water_pressure + air_pressure; + + for(i=0;iweight*pressure*l1l4[i]*normal[j]; + } + } + } + VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES); /*Free ressources:*/ + delete gauss; xfree((void**)&doflist); - -} -/*}}}*/ -/*FUNCTION Icefront::CreatePVectorDiagnosticStokes {{{1*/ +} +/*}}}*/ +/*FUNCTION Icefront::CreatePVectorDiagnosticStokes{{{1*/ void Icefront::CreatePVectorDiagnosticStokes( Vec pg){ - int i,j; - - const int numnodes=4; - const int numdofs=numnodes*NDOF4; - int* doflist=NULL; - - double xyz_list[numnodes][3]; - double xyz_list_quad[numnodes+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]={0.0}; - Penta* penta=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]; - int approximation; - - /*check icefront is associated to a pentaelem: */ - if(element->Enum()!=PentaEnum) ISSMERROR("Only Penta supported yet"); + /*Constants*/ + const int numdofs = NUMVERTICESQUA *NDOF4; + + /*Intermediaries*/ + int i,j,ig,index1,index2,index3,index4; + int fill; + double pressure,rho_water,gravity; + double water_pressure,air_pressure; + double Jdet,z_g; + double xyz_list[NUMVERTICESQUA][3]; + double normal[3]; + double pe_g[numdofs] = {0.0}; + double l1l4[4]; + int *doflist = NULL; + GaussPenta *gauss = NULL; + Penta *penta = NULL; + + /* Get dof list and node coordinates: */ penta=(Penta*)element; - - /*If not Stokes, return*/ - penta->inputs->GetParameterValue(&approximation,ApproximationEnum); - if(approximation!=StokesApproximationEnum) return; - - /* Get dof list and node coordinates: */ GetDofList(&doflist,StokesApproximationEnum,GsetEnum); - GetVerticesCoordinates(&xyz_list[0][0], nodes, numnodes); - - //Identify which grids are comprised in the quad: - grid1=element->GetNodeIndex(nodes[0]); - grid2=element->GetNodeIndex(nodes[1]); - grid3=element->GetNodeIndex(nodes[2]); - grid4=element->GetNodeIndex(nodes[3]); - - /*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->GetParameterListOnVertices(&thickness_list[0],ThicknessEnum); - element->GetParameterListOnVertices(&bed_list[0],BedEnum); - - 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: - QuadPressureLoadStokes(pe_g,matpar->GetRhoWater(),matpar->GetRhoIce(),matpar->GetG(),thickness_list_quad,bed_list_quad,normal1,normal2,normal3,normal4,&xyz_list_quad[0][0]); - - /*Plug pe_g into vector: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICESQUA); + + /*Recover inputs and parameters */ + inputs->GetParameterValue(&fill,FillEnum); + rho_water=matpar->GetRhoWater(); + gravity =matpar->GetG(); + GetQuadNormal(&normal[0],xyz_list); + + /*Identify which nodes are in the quad: */ + index1=element->GetNodeIndex(nodes[0]); + index2=element->GetNodeIndex(nodes[1]); + index3=element->GetNodeIndex(nodes[2]); + index4=element->GetNodeIndex(nodes[3]); + + /* Start looping on the number of gaussian points: */ + double zmax=xyz_list[0][2]; for(i=1;izmax) zmax=xyz_list[i][2]; + double zmin=xyz_list[0][2]; for(i=1;i0 && zmin<0) gauss=new GaussPenta(index1,index2,index3,index4,3,30); //refined in vertical because of the sea level discontinuity + else gauss=new GaussPenta(index1,index2,index3,index4,3,3); + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + penta->GetQuadNodalFunctions(l1l4,gauss,index1,index2,index3,index4); + penta->GetQuadJacobianDeterminant(&Jdet,xyz_list,gauss); + z_g=penta->GetZcoord(gauss); + + switch(fill){ + case WaterEnum: + water_pressure=rho_water*gravity*min(0,z_g);//0 if the gaussian point is above water level + break; + case AirEnum: + water_pressure=0; + break; + default: + ISSMERROR("fill type %s not supported yet",EnumToString(fill)); + } + air_pressure=0; + pressure = water_pressure + air_pressure; //no ice pressure fore Stokes + + for(i=0;iweight*pressure*l1l4[i]*normal[j]; + else pe_g[i*NDOF4+j]+=0; //pressure term + } + } + } + VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES); /*Free ressources:*/ + delete gauss; xfree((void**)&doflist); - } /*}}}*/ @@ -949,489 +872,6 @@ } /*}}}*/ -/*FUNCTION Icefront::QuadPressureLoad {{{1*/ -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){ - - - //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 ig; - GaussTria *gauss=NULL; - - 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]; - int fill; - - - /*To use tria functionalities: */ - Tria* tria=NULL; - - /*Recover inputs: */ - inputs->GetParameterValue(&fill,FillEnum); - - /*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]; - - //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 number of gaussian points: */ - gauss=new GaussTria(2); - for(ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - tria->GetNodalFunctions(l1l2l3_tria,gauss); - - //Get the coordinates of gauss point for each triangle in penta/quad - r_tria=gauss->coord2-gauss->coord1; - s_tria=-3/sqrt((double)3)*(gauss->coord1+gauss->coord2-2/3); - - //Coordinates of gauss points in the reference triangle - r_quad[0]=r_tria; - s_quad[0]=1/sqrt((double)3)*s_tria-2/3; - r_quad[1]=-1/sqrt((double)3)*s_tria+2/3; - s_quad[1]=r_tria; - r_quad[2]=-r_tria; - s_quad[2]=-1/sqrt((double)3)*s_tria+2/3; - r_quad[3]=1/sqrt((double)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->GetJacobianDeterminant2d(&J[i],&complete_list[0][0],gauss); - } - - //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==WaterEnum){ //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==AirEnum){ - water_pressure_tria=0; - } - else{ - ISSMERROR("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++) - } - - delete tria; - delete gauss; -} -/*}}}*/ -/*FUNCTION Icefront::QuadPressureLoadStokes {{{1*/ -void Icefront::QuadPressureLoadStokes(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){ - - - //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]; - double nz[4]; - - /* gaussian points: */ - int ig; - GaussTria* gauss=NULL; - - 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 air_pressure_tria; - double water_level_above_g_tria; - double water_pressure_tria; - double pressure_tria; - int fill; - - /*To use tria functionalities: */ - Tria* tria=NULL; - - /*Recover inputs: */ - inputs->GetParameterValue(&fill,FillEnum); - - /*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]; - nz[0]=normal1[2]; nz[1]=normal2[2]; nz[2]=normal3[2]; nz[3]=normal4[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 number of gaussian points: */ - gauss=new GaussTria(2); - for(ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - tria->GetNodalFunctions(l1l2l3_tria,gauss); - - //Get the coordinates of gauss point for each triangle in penta/quad - r_tria=gauss->coord2-gauss->coord1; - s_tria=-3/sqrt((double)3)*(gauss->coord1+gauss->coord2-2/3); - - //Coordinates of gauss points in the reference triangle - r_quad[0]=r_tria; - s_quad[0]=1/sqrt((double)3)*s_tria-2/3; - r_quad[1]=-1/sqrt((double)3)*s_tria+2/3; - s_quad[1]=r_tria; - r_quad[2]=-r_tria; - s_quad[2]=-1/sqrt((double)3)*s_tria+2/3; - r_quad[3]=1/sqrt((double)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->GetJacobianDeterminant2d(&J[i],&complete_list[0][0],gauss); - } - - //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 - air_pressure_tria=0; - - //Now deal with water pressure: - if(fill==WaterEnum){ //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==AirEnum){ - water_pressure_tria=0; - } - else{ - ISSMERROR("QuadPressureLoadStokes error message: unknow fill type for icefront boundary condition"); - } - - //Add pressure - pressure_tria = water_pressure_tria + air_pressure_tria; - - pe_g[0]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][0]*nx[i]; - pe_g[1]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][0]*ny[i]; - pe_g[2]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][0]*nz[i]; - pe_g[3]+= 0; - pe_g[4]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][1]*nx[i]; - pe_g[5]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][1]*ny[i]; - pe_g[6]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][1]*nz[i]; - pe_g[7]+= 0; - pe_g[8]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][2]*nx[i]; - pe_g[9]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][2]*ny[i]; - pe_g[10]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][2]*nz[i]; - pe_g[11]+= 0; - pe_g[12]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][3]*nx[i]; - pe_g[13]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][3]*ny[i]; - pe_g[14]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][3]*nz[i]; - pe_g[15]+= 0; - - } //for(i=0;i<4;i++) - } - - delete tria; - delete gauss; -} -/*}}}*/ -/*FUNCTION Icefront::GetNormal {{{1*/ -void Icefront:: GetNormal(double* normal,double xyz_list[2][3]){ +/*FUNCTION Icefront::GetSegmentNormal {{{1*/ +void Icefront:: GetSegmentNormal(double* normal,double xyz_list[4][3]){ /*Build unit outward pointing vector*/ @@ -1447,4 +887,25 @@ normal[0]= + vector[1]/norm; normal[1]= - vector[0]/norm; +} +/*}}}*/ +/*FUNCTION Icefront::GetQuadNormal {{{1*/ +void Icefront:: GetQuadNormal(double* normal,double xyz_list[4][3]){ + + /*Build unit outward pointing vector*/ + double AB[3]; + double AC[3]; + double norm; + + AB[0]=xyz_list[1][0] - xyz_list[0][0]; + AB[1]=xyz_list[1][1] - xyz_list[0][1]; + AB[2]=xyz_list[1][2] - xyz_list[0][2]; + AC[0]=xyz_list[2][0] - xyz_list[0][0]; + AC[1]=xyz_list[2][1] - xyz_list[0][1]; + AC[2]=xyz_list[2][2] - xyz_list[0][2]; + + cross(normal,AB,AC); + norm=sqrt(pow(normal[0],2.0)+pow(normal[1],2.0)+pow(normal[2],2.0)); + + for(int i=0;i<3;i++) normal[i]=normal[i]/norm; } /*}}}*/ Index: /issm/trunk/src/c/objects/Loads/Icefront.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Icefront.h (revision 5807) +++ /issm/trunk/src/c/objects/Loads/Icefront.h (revision 5808) @@ -85,9 +85,6 @@ int* GetExternalDofList(int approximation_enum,int setenum); int GetNumberOfDofs(int approximation_enum,int setenum); - void 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); - void QuadPressureLoadStokes(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); - void GetNormal(double* normal,double xyz_list[2][3]); + void GetSegmentNormal(double* normal,double xyz_list[2][3]); + void GetQuadNormal(double* normal,double xyz_list[4][3]); ElementVector* NewElementVector(int approximation); /*}}}*/ Index: /issm/trunk/src/c/shared/Numerics/GaussPoints.cpp =================================================================== --- /issm/trunk/src/c/shared/Numerics/GaussPoints.cpp (revision 5807) +++ /issm/trunk/src/c/shared/Numerics/GaussPoints.cpp (revision 5808) @@ -1720,11 +1720,2 @@ GaussLegendreLinear(pzgaus, pzwgt, nkgaus); }/*}}}1*/ -/*FUNCTION gaussPenta{{{1*/ -void gaussPenta( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, double** pzgaus, double** pzwgt, int iord, int nkgaus ) { - /*Gauss quadrature points for the pentahedron.*/ - - /* get the gauss points using the product of triangular and line rules */ - GaussLegendreTria(pngaus, pl1, pl2, pl3, pwgt, iord); - GaussLegendreLinear(pzgaus, pzwgt, nkgaus); - -}/*}}}1*/ Index: /issm/trunk/src/c/shared/Numerics/GaussPoints.h =================================================================== --- /issm/trunk/src/c/shared/Numerics/GaussPoints.h (revision 5807) +++ /issm/trunk/src/c/shared/Numerics/GaussPoints.h (revision 5808) @@ -19,5 +19,4 @@ void gaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus ); void gaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus ); -void gaussPenta( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, double** pzgaus, double** pzwgt, int iord, int nkgaus ); #endif