Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5721) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5722) @@ -508,14 +508,8 @@ int dofp[1]={3}; double Jdet2d; - Tria* tria=NULL; /*Gauss*/ - 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[4]; + int ig; + GaussPenta* gauss; /*Output*/ @@ -529,8 +523,4 @@ bool onbed; int approximation; - Input* pressure_input=NULL; - Input* vx_input=NULL; - Input* vy_input=NULL; - Input* vz_input=NULL; /*parameters: */ @@ -569,48 +559,41 @@ } - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussLegendreTria(&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights,2); - /*Retrieve all inputs we will be needing: */ - pressure_input=inputs->GetInput(PressureEnum); - vx_input=inputs->GetInput(VxEnum); - vy_input=inputs->GetInput(VyEnum); - vz_input=inputs->GetInput(VzEnum); + Input* pressure_input=inputs->GetInput(PressureEnum); + Input* vx_input=inputs->GetInput(VxEnum); + Input* vy_input=inputs->GetInput(VyEnum); + Input* vz_input=inputs->GetInput(VzEnum); /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(0,1,2,2); for (ig=0; igGetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - pressure_input->GetParameterValue(&pressure, &gauss_coord[0]); - - /*Compute Stress*/ - sigma_xx=2*viscosity*epsilon[0]-pressure*stokesreconditioning; // sigma = nu eps - pressure - sigma_yy=2*viscosity*epsilon[1]-pressure*stokesreconditioning; - sigma_zz=2*viscosity*epsilon[2]-pressure*stokesreconditioning; - sigma_xy=2*viscosity*epsilon[3]; - sigma_xz=2*viscosity*epsilon[4]; - sigma_yz=2*viscosity*epsilon[5]; - - /*Get normal vector to the bed */ - BedNormal(&bed_normal[0],xyz_list_tria); - - /*basalforce*/ - basalforce[0] += sigma_xx*bed_normal[0] + sigma_xy*bed_normal[1] + sigma_xz*bed_normal[2]; - basalforce[1] += sigma_xy*bed_normal[0] + sigma_yy*bed_normal[1] + sigma_yz*bed_normal[2]; - basalforce[2] += sigma_xz*bed_normal[0] + sigma_yz*bed_normal[1] + sigma_zz*bed_normal[2]; - - /*Get the Jacobian determinant */ - tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0],gauss_coord); - value+=sigma_zz*Jdet2d*gauss_weight; - surface+=Jdet2d*gauss_weight; + gauss->GaussPoint(ig); + + /*Compute strain rate viscosity and pressure: */ + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + pressure_input->GetParameterValue(&pressure, &gauss_coord[0]); + + /*Compute Stress*/ + sigma_xx=2*viscosity*epsilon[0]-pressure*stokesreconditioning; // sigma = nu eps - pressure + sigma_yy=2*viscosity*epsilon[1]-pressure*stokesreconditioning; + sigma_zz=2*viscosity*epsilon[2]-pressure*stokesreconditioning; + sigma_xy=2*viscosity*epsilon[3]; + sigma_xz=2*viscosity*epsilon[4]; + sigma_yz=2*viscosity*epsilon[5]; + + /*Get normal vector to the bed */ + BedNormal(&bed_normal[0],xyz_list_tria); + + /*basalforce*/ + basalforce[0] += sigma_xx*bed_normal[0] + sigma_xy*bed_normal[1] + sigma_xz*bed_normal[2]; + basalforce[1] += sigma_xy*bed_normal[0] + sigma_yy*bed_normal[1] + sigma_yz*bed_normal[2]; + basalforce[2] += sigma_xz*bed_normal[0] + sigma_yz*bed_normal[1] + sigma_zz*bed_normal[2]; + + /*Get the Jacobian determinant */ + GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss); + value+=sigma_zz*Jdet2d*gauss_weight; + surface+=Jdet2d*gauss_weight; } value=value/surface; @@ -2735,7 +2718,4 @@ const int numdof2d=NUMVERTICES2D*NDOF4; - /*Collapsed formulation: */ - Tria* tria=NULL; - /*Grid data: */ double xyz_list[NUMVERTICES][3]; @@ -2761,6 +2741,4 @@ int ig; GaussPenta *gauss=NULL; - GaussTria *gauss_tria=NULL; - /* specific gaussian point: */ @@ -2786,9 +2764,4 @@ bool shelf; - /*inputs: */ - Input* vx_input=NULL; - Input* vy_input=NULL; - Input* vz_input=NULL; - /*retrive parameters: */ parameters->FindParam(&analysis_type,AnalysisTypeEnum); @@ -2808,14 +2781,8 @@ GetDofList(&doflist,StokesApproximationEnum); - /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore - get tria gaussian points as well as segment gaussian points. For tria gaussian - points, order of integration is 2, because we need to integrate the product tB*D*B' - which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian - points, same deal, which yields 3 gaussian points.*/ - /*Retrieve all inputs we will be needing: */ - vx_input=inputs->GetInput(VxEnum); - vy_input=inputs->GetInput(VyEnum); - vz_input=inputs->GetInput(VzEnum); + Input* vx_input=inputs->GetInput(VxEnum); + Input* vy_input=inputs->GetInput(VyEnum); + Input* vz_input=inputs->GetInput(VzEnum); /* Start looping on the number of gaussian points: */ @@ -2853,4 +2820,5 @@ for(i=0;i<27;i++) for(j=0;j<27;j++) Ke_temp[i][j]+=Ke_gaussian[i][j]; } + delete gauss; //gauss of previous loop if((onbed==1) && (shelf==0)){ @@ -2866,15 +2834,11 @@ /* Start looping on the number of gaussian points: */ - delete gauss; //gauss of previous loop - gauss=new GaussPenta(); - gauss->GaussFaceTria(0,1,2,2); - gauss_tria=new GaussTria(); + gauss=new GaussPenta(0,1,2,2); for (ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - gauss->SynchronizeGaussTria(gauss_tria); /*Get the Jacobian determinant */ - tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_tria); + GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss); /*Get L matrix if viscous basal drag present: */ @@ -2920,4 +2884,5 @@ /*Free ressources:*/ delete friction; + delete gauss; } //if ( (onbed==1) && (shelf==0)) @@ -2930,9 +2895,5 @@ /*Free ressources:*/ - delete gauss; - delete gauss_tria; xfree((void**)&doflist); - - return; } /*}}}*/ @@ -3742,5 +3703,4 @@ int ig; GaussPenta *gauss=NULL; - GaussTria *gauss_tria=NULL; double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ @@ -3854,4 +3814,5 @@ for(i=0;i<27;i++) for(j=0;j<3;j++) Ke_temp[i][j]+=Ke_gaussian[i][j]; } + delete gauss; /*Deal with 2d friction at the bedrock interface: */ @@ -3865,15 +3826,11 @@ /* Start looping on the number of gauss 2d (nodes on the bedrock) */ - delete gauss; //gauss of previous loop - gauss=new GaussPenta(); - gauss->GaussFaceTria(0,1,2,2); - gauss_tria=new GaussTria(); + gauss=new GaussPenta(0,1,2,2) for (ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - gauss->SynchronizeGaussTria(gauss_tria); /*Get the Jacobian determinant */ - tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_tria); + GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); /* Get bed at gaussian point */ @@ -3895,4 +3852,6 @@ } } + /*Free ressources:*/ + delete gauss; } //if ( (onbed==1) && (shelf==1)) @@ -3905,5 +3864,4 @@ /*Free ressources:*/ delete gauss; - delete gauss_tria; xfree((void**)&doflist); Index: /issm/trunk/src/c/objects/Elements/PentaRef.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 5721) +++ /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 5722) @@ -1677,4 +1677,26 @@ } /*}}}*/ +/*FUNCTION PentaRef::GetTriaJacobianDeterminant{{{1*/ +void PentaRef::GetTriaJacobianDeterminant(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,x3,y1,y2,y3,z1,z2,z3; + + 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); + x3=*(xyz_list+3*2+0); + y3=*(xyz_list+3*2+1); + z3=*(xyz_list+3*2+2); + + *Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5); + 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 5721) +++ /issm/trunk/src/c/objects/Elements/PentaRef.h (revision 5722) @@ -58,4 +58,5 @@ void GetJacobian(double* J, double* xyz_list,GaussPenta* gauss); void GetJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss); + void GetTriaJacobianDeterminant(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); Index: /issm/trunk/src/c/objects/Gauss/GaussPenta.cpp =================================================================== --- /issm/trunk/src/c/objects/Gauss/GaussPenta.cpp (revision 5721) +++ /issm/trunk/src/c/objects/Gauss/GaussPenta.cpp (revision 5722) @@ -80,4 +80,23 @@ } /*}}}*/ +/*FUNCTION GaussPenta::GaussPenta{{{1*/ +GaussPenta::GaussPenta(int index1, int index2, int index3, int order){ + + /*Basal Tria*/ + if(index1==0 && index2==1 && index3==2){ + + /*Get GaussTria*/ + GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order); + + /*compute z coordinate*/ + coords4=(double*)xmalloc(numgauss*sizeof(double)); + for(int i=0;iGetParameterValue(&value,nodes[0],nodes[1],gauss_coord,input_in); - - /*Assign output pointers:*/ - *pvalue=value; -} -/*}}}*/ /*FUNCTION Icefront::GetNormal {{{1*/ void Icefront:: GetNormal(double* normal,double xyz_list[2][3]){ Index: /issm/trunk/src/c/objects/Loads/Icefront.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Icefront.h (revision 5721) +++ /issm/trunk/src/c/objects/Loads/Icefront.h (revision 5722) @@ -85,5 +85,4 @@ 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 GetParameterValue(double* pvalue, double gauss_coord,Input* input_in); void GetNormal(double* normal,double xyz_list[2][3]); /*}}}*/