Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5725) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5726) @@ -2740,12 +2740,4 @@ 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,17 +3408,10 @@ 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,10 +3429,4 @@ int connectivity[2]; - /*Inputs*/ - Input* thickness_input; - Input* surface_input; - Input* slopex_input; - Input* slopey_input; - /*recover some inputs: */ inputs->GetParameterValue(&onwater,ElementOnWaterEnum); @@ -3469,75 +3448,58 @@ 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;iGetParameterValue(&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;igGetParameterValue(&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();igend();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;jweight/(double)connectivity[1]; } else{//connectivity is too large, should take only half on it - for(j=0;jweight*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,7 +3511,4 @@ /*Clean up */ - delete beam; - xfree((void**)&gauss_weights); - xfree((void**)&segment_gauss_coord); xfree((void**)&doflist); } Index: /issm/trunk/src/c/objects/Elements/PentaRef.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 5725) +++ /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 5726) @@ -1699,4 +1699,23 @@ } /*}}}*/ +/*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){ Index: /issm/trunk/src/c/objects/Elements/PentaRef.h =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaRef.h (revision 5725) +++ /issm/trunk/src/c/objects/Elements/PentaRef.h (revision 5726) @@ -59,4 +59,5 @@ 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);