Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5646) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5647) @@ -2713,26 +2713,6 @@ /* 3d 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* fourth_gauss_vert_coord = NULL; - double* area_gauss_weights = NULL; - double* vert_gauss_weights = NULL; - int ig1,ig2; - double gauss_weight1,gauss_weight2; - double gauss_coord[4]; - int order_area_gauss; - int num_vert_gauss; - int num_area_gauss; - double gauss_weight; - - /* 2d gaussian point: */ - int num_gauss2d; - double* first_gauss_area_coord2d = NULL; - double* second_gauss_area_coord2d = NULL; - double* third_gauss_area_coord2d = NULL; - double* gauss_weights2d=NULL; - double gauss_l1l2l3[3]; + int ig; + GaussPenta *gauss=NULL; /* material data: */ @@ -2802,15 +2782,4 @@ GetDofList(&doflist,PattynApproximationEnum); - /*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.*/ - - order_area_gauss=5; - num_vert_gauss=5; - - gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss); - /*Retrieve all inputs we will be needing: */ vx_input=inputs->GetInput(VxEnum); @@ -2820,57 +2789,40 @@ /* Start looping on the number of gaussian points: */ - for (ig1=0; ig1GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input); - this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,vxold_input,vyold_input); - - /*Get viscosity: */ - matice->GetViscosity3d(&viscosity, &epsilon[0]); - matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); - - /*Get B and Bprime matrices: */ - GetBPattyn(&B[0][0], &xyz_list[0][0], gauss_coord); - GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss_coord); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); - - /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant - onto this scalar matrix, so that we win some computational time: */ - - newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); - D_scalar=2*newviscosity*gauss_weight*Jdet; - for (i=0;i<5;i++){ - D[i][i]=D_scalar; - } - - /* Do the triple product tB*D*Bprime: */ - TripleMultiply( &B[0][0],5,numdof,1, - &D[0][0],5,5,0, - &Bprime[0][0],5,numdof,0, - &Ke_gg_gaussian[0][0],0); - - /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ - for( i=0; ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + /*Get strain rate from velocity: */ + this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); + + /*Get viscosity: */ + matice->GetViscosity3d(&viscosity, &epsilon[0]); + matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); + + /*Get B and Bprime matrices: */ + GetBPattyn(&B[0][0], &xyz_list[0][0], gauss); + GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + + /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant + onto this scalar matrix, so that we win some computational time: */ + + newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); + D_scalar=2*newviscosity*gauss->weight*Jdet; + for (i=0;i<5;i++) D[i][i]=D_scalar; + + /* Do the triple product tB*D*Bprime: */ + TripleMultiply( &B[0][0],5,numdof,1, + &D[0][0],5,5,0, + &Bprime[0][0],5,numdof,0, + &Ke_gg_gaussian[0][0],0); + + /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ + for(i=0;iGetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss); + vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss); + + /*Sum all contributions*/ + for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; + +} +/*}}}*/ +/*FUNCTION Penta::GetStrainRate3dPattyn{{{1*/ +void Penta::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input){ /*Compute the 3d Blatter/PattynStrain Rate (5 components): * Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 5646) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 5647) @@ -168,4 +168,5 @@ void GetSolutionFromInputsThermal(Vec solutiong); void GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input); + void GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input); void GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input); Penta* GetUpperElement(void); Index: /issm/trunk/src/c/objects/Elements/PentaRef.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 5646) +++ /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 5647) @@ -981,2 +981,835 @@ } /*}}}*/ + +/*FUNCTION PentaRef::GetBMacAyealPattyn {{{1*/ +void PentaRef::GetBMacAyealPattyn(double* B, double* xyz_list, GaussPenta* gauss){ + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. + * For grid i, Bi can be expressed in the actual coordinate system + * by: + * Bi=[ dh/dx 0 ] + * [ 0 dh/dy ] + * [ 1/2*dh/dy 1/2*dh/dx ] + * where h is the interpolation function for grid i. + * + * We assume B has been allocated already, of size: 5x(NDOF2*numgrids) + */ + + int i; + const int numgrids=6; + const int NDOF3=3; + const int NDOF2=2; + + double dh1dh6[NDOF3][numgrids]; + + /*Get dh1dh6 in actual coordinate system: */ + GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss); + + /*Build B: */ + for (i=0;icoord1; + A2=gauss->coord2; + A3=gauss->coord3; + + xi=A2-A1; + eta=SQRT3*A3; + zi=gauss->coord4; + + x1=*(xyz_list+3*0+0); + x2=*(xyz_list+3*1+0); + x3=*(xyz_list+3*2+0); + x4=*(xyz_list+3*3+0); + x5=*(xyz_list+3*4+0); + x6=*(xyz_list+3*5+0); + + y1=*(xyz_list+3*0+1); + y2=*(xyz_list+3*1+1); + y3=*(xyz_list+3*2+1); + y4=*(xyz_list+3*3+1); + y5=*(xyz_list+3*4+1); + y6=*(xyz_list+3*5+1); + + z1=*(xyz_list+3*0+2); + z2=*(xyz_list+3*1+2); + z3=*(xyz_list+3*2+2); + z4=*(xyz_list+3*3+2); + z5=*(xyz_list+3*4+2); + z6=*(xyz_list+3*5+2); + + *(J+NDOF3*0+0)=0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5); + *(J+NDOF3*1+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+SQRT3/12.0*(-x1-x2+2*x3-x4-x5+2*x6); + *(J+NDOF3*2+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4); + + *(J+NDOF3*0+1)=0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5); + *(J+NDOF3*1+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+SQRT3/12.0*(-y1-y2+2*y3-y4-y5+2*y6); + *(J+NDOF3*2+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2); + + *(J+NDOF3*0+2)=0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5); + *(J+NDOF3*1+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+SQRT3/12.0*(-z1-z2+2*z3-z4-z5+2*z6); + *(J+NDOF3*2+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4); + +} +/*}}}*/ +/*FUNCTION PentaRef::GetJacobianDeterminant {{{1*/ +void PentaRef::GetJacobianDeterminant(double* Jdet, double* xyz_list,GaussPenta* gauss){ + /*On a penta, Jacobian varies according to coordinates. We need to get the Jacobian, and take + * the determinant of it: */ + double J[3][3]; + + /*Get Jacobian*/ + GetJacobian(&J[0][0],xyz_list,gauss); + + /*Get Determinant*/ + Matrix3x3Determinant(Jdet,&J[0][0]); + if(*Jdet<0) ISSMERROR("negative jacobian determinant!"); + +} +/*}}}*/ +/*FUNCTION PentaRef::GetJacobianInvert {{{1*/ +void PentaRef::GetJacobianInvert(double* Jinv, double* xyz_list,GaussPenta* gauss){ + + /*Jacobian*/ + double J[3][3]; + + /*Call Jacobian routine to get the jacobian:*/ + GetJacobian(&J[0][0], xyz_list, gauss); + + /*Invert Jacobian matrix: */ + Matrix3x3Invert(Jinv,&J[0][0]); +} +/*}}}*/ +/*FUNCTION PentaRef::GetNodalFunctionsMINI{{{1*/ +void PentaRef::GetNodalFunctionsMINI(double* l1l7, GaussPenta* gauss){ + /*This routine returns the values of the nodal functions at the gaussian point.*/ + + l1l7[0]=gauss->coord1*(1.0-gauss->coord4)/2.0; + l1l7[1]=gauss->coord2*(1.0-gauss->coord4)/2.0; + l1l7[2]=gauss->coord3*(1.0-gauss->coord4)/2.0; + l1l7[3]=gauss->coord1*(1.0+gauss->coord4)/2.0; + l1l7[4]=gauss->coord2*(1.0+gauss->coord4)/2.0; + l1l7[5]=gauss->coord3*(1.0+gauss->coord4)/2.0; + l1l7[6]=27*gauss->coord1*gauss->coord2*gauss->coord3*(1.0+gauss->coord4)*(1.0-gauss->coord4); + +} +/*}}}*/ +/*FUNCTION PentaRef::GetNodalFunctionsMINIDerivatives{{{1*/ +void PentaRef::GetNodalFunctionsMINIDerivatives(double* dh1dh7,double* xyz_list, GaussPenta* gauss){ + + /*This routine returns the values of the nodal functions derivatives (with respect to the + * actual coordinate system): */ + + int i; + const int numgrids = 7; + double dh1dh7_ref[3][numgrids]; + double Jinv[3][3]; + + /*Get derivative values with respect to parametric coordinate system: */ + GetNodalFunctionsMINIDerivativesReference(&dh1dh7_ref[0][0], gauss); + + /*Get Jacobian invert: */ + GetJacobianInvert(&Jinv[0][0], xyz_list, gauss); + + /*Build dh1dh6: + * + * [dhi/dx]= Jinv'*[dhi/dr] + * [dhi/dy] [dhi/ds] + * [dhi/dz] [dhi/dzeta] + */ + + for (i=0;icoord2-gauss->coord1; + double s=-3.0/SQRT3*(gauss->coord1+gauss->coord2-2.0/3.0); + double zeta=gauss->coord4; + + /*First nodal function: */ + *(dl1dl7+numgrids*0+0)=-0.5*(1.0-zeta)/2.0; + *(dl1dl7+numgrids*1+0)=-SQRT3/6.0*(1.0-zeta)/2.0; + *(dl1dl7+numgrids*2+0)=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD); + + /*Second nodal function: */ + *(dl1dl7+numgrids*0+1)=0.5*(1.0-zeta)/2.0; + *(dl1dl7+numgrids*1+1)=-SQRT3/6.0*(1.0-zeta)/2.0; + *(dl1dl7+numgrids*2+1)=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD); + + /*Third nodal function: */ + *(dl1dl7+numgrids*0+2)=0; + *(dl1dl7+numgrids*1+2)=SQRT3/3.0*(1.0-zeta)/2.0; + *(dl1dl7+numgrids*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD); + + /*Fourth nodal function: */ + *(dl1dl7+numgrids*0+3)=-0.5*(1.0+zeta)/2.0; + *(dl1dl7+numgrids*1+3)=-SQRT3/6.0*(1.0+zeta)/2.0; + *(dl1dl7+numgrids*2+3)=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD); + + /*Fith nodal function: */ + *(dl1dl7+numgrids*0+4)=0.5*(1.0+zeta)/2.0; + *(dl1dl7+numgrids*1+4)=-SQRT3/6.0*(1.0+zeta)/2.0; + *(dl1dl7+numgrids*2+4)=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD); + + /*Sixth nodal function: */ + *(dl1dl7+numgrids*0+5)=0; + *(dl1dl7+numgrids*1+5)=SQRT3/3.0*(1.0+zeta)/2.0; + *(dl1dl7+numgrids*2+5)=0.5*(SQRT3/3.0*s+ONETHIRD); + + /*Seventh nodal function: */ + *(dl1dl7+numgrids*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0); + *(dl1dl7+numgrids*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0)); + *(dl1dl7+numgrids*2+6)=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta); + +} +/*}}}*/ +/*FUNCTION PentaRef::GetNodalFunctionsP1 {{{1*/ +void PentaRef::GetNodalFunctionsP1(double* l1l6, GaussPenta* gauss){ + /*This routine returns the values of the nodal functions at the gaussian point.*/ + + l1l6[0]=gauss->coord1*(1-gauss->coord4)/2.0; + l1l6[1]=gauss->coord2*(1-gauss->coord4)/2.0; + l1l6[2]=gauss->coord3*(1-gauss->coord4)/2.0; + l1l6[3]=gauss->coord1*(1+gauss->coord4)/2.0; + l1l6[4]=gauss->coord2*(1+gauss->coord4)/2.0; + l1l6[5]=gauss->coord3*(1+gauss->coord4)/2.0; + +} +/*}}}*/ +/*FUNCTION PentaRef::GetNodalFunctionsP1Derivatives {{{1*/ +void PentaRef::GetNodalFunctionsP1Derivatives(double* dh1dh6,double* xyz_list, GaussPenta* gauss){ + + /*This routine returns the values of the nodal functions derivatives (with respect to the + * actual coordinate system): */ + int i; + const int NDOF3 = 3; + const int numgrids = 6; + double dh1dh6_ref[NDOF3][numgrids]; + double Jinv[NDOF3][NDOF3]; + + /*Get derivative values with respect to parametric coordinate system: */ + GetNodalFunctionsP1DerivativesReference(&dh1dh6_ref[0][0], gauss); + + /*Get Jacobian invert: */ + GetJacobianInvert(&Jinv[0][0], xyz_list, gauss); + + /*Build dh1dh3: + * + * [dhi/dx]= Jinv*[dhi/dr] + * [dhi/dy] [dhi/ds] + * [dhi/dz] [dhi/dn] + */ + + for (i=0;icoord1; ISSMASSERT(A1>=0 && A1<=1);//first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*SQRT3); + A2=gauss->coord2; ISSMASSERT(A2>=0 && A2<=1);//second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*SQRT3); + A3=gauss->coord3; ISSMASSERT(A3>=0 && A3<=1);//third area coordinate value In term of xi and eta: A3=y/SQRT3; + z =gauss->coord4; ISSMASSERT(z>=-1 && z<=1);//fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2 + + /*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/ + *(dl1dl6+numgrids*0+0)=-0.5*(1.0-z)/2.0; + *(dl1dl6+numgrids*1+0)=-0.5/SQRT3*(1.0-z)/2.0; + *(dl1dl6+numgrids*2+0)=-0.5*A1; + + /*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/ + *(dl1dl6+numgrids*0+1)=0.5*(1.0-z)/2.0; + *(dl1dl6+numgrids*1+1)=-0.5/SQRT3*(1.0-z)/2.0; + *(dl1dl6+numgrids*2+1)=-0.5*A2; + + /*Third nodal function: The corresponding nodal function is N=A3*(1-z)/2. Its derivatives follow*/ + *(dl1dl6+numgrids*0+2)=0.0; + *(dl1dl6+numgrids*1+2)=1.0/SQRT3*(1.0-z)/2.0; + *(dl1dl6+numgrids*2+2)=-0.5*A3; + + /*Fourth nodal function: The corresponding nodal function is N=A1*(1+z)/2. Its derivatives follow*/ + *(dl1dl6+numgrids*0+3)=-0.5*(1.0+z)/2.0; + *(dl1dl6+numgrids*1+3)=-0.5/SQRT3*(1.0+z)/2.0; + *(dl1dl6+numgrids*2+3)=0.5*A1; + + /*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/ + *(dl1dl6+numgrids*0+4)=0.5*(1.0+z)/2.0; + *(dl1dl6+numgrids*1+4)=-0.5/SQRT3*(1.0+z)/2.0; + *(dl1dl6+numgrids*2+4)=0.5*A2; + + /*Sixth nodal function: The corresponding nodal function is N=A3*(1+z)/2. Its derivatives follow*/ + *(dl1dl6+numgrids*0+5)=0.0; + *(dl1dl6+numgrids*1+5)=1.0/SQRT3*(1.0+z)/2.0; + *(dl1dl6+numgrids*2+5)=0.5*A3; +} +/*}}}*/ +/*FUNCTION PentaRef::GetParameterValue{{{1*/ +void PentaRef::GetParameterValue(double* pvalue,double* plist,GaussPenta* gauss){ + /*P1 interpolation on Gauss point*/ + + /*intermediary*/ + double l1l6[6]; + + /*nodal functions: */ + GetNodalFunctionsP1(&l1l6[0],gauss); + + /*Assign output pointers:*/ + *pvalue=l1l6[0]*plist[0]+l1l6[1]*plist[1]+l1l6[2]*plist[2]+l1l6[3]*plist[3]+l1l6[4]*plist[4]+l1l6[5]*plist[5]; + +} +/*}}}*/ +/*FUNCTION PentaRef::GetParameterDerivativeValue{{{1*/ +void PentaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, GaussPenta* gauss){ + /*From grid values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_coord: + * dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx; + * dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy; + * dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz; + * + * p is a vector of size 3x1 already allocated. + */ + double dh1dh6[3][6]; + + /*Get nodal funnctions derivatives in actual coordinate system: */ + GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss); + + /*Assign output*/ + p[0]=plist[0]*dh1dh6[0][0]+plist[1]*dh1dh6[0][1]+plist[2]*dh1dh6[0][2]+plist[3]*dh1dh6[0][3]+plist[4]*dh1dh6[0][4]+plist[5]*dh1dh6[0][5]; + p[1]=plist[0]*dh1dh6[1][0]+plist[1]*dh1dh6[1][1]+plist[2]*dh1dh6[1][2]+plist[3]*dh1dh6[1][3]+plist[4]*dh1dh6[1][4]+plist[5]*dh1dh6[1][5]; + p[2]=plist[0]*dh1dh6[2][0]+plist[1]*dh1dh6[2][1]+plist[2]*dh1dh6[2][2]+plist[3]*dh1dh6[2][3]+plist[4]*dh1dh6[2][4]+plist[5]*dh1dh6[2][5]; + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/PentaRef.h =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaRef.h (revision 5646) +++ /issm/trunk/src/c/objects/Elements/PentaRef.h (revision 5647) @@ -50,4 +50,31 @@ void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussTria* gauss){ISSMERROR("only PentaGauss are supported");}; + void GetNodalFunctionsP1(double* l1l6, GaussPenta* gauss); + void GetNodalFunctionsMINI(double* l1l7, GaussPenta* gauss); + void GetNodalFunctionsP1Derivatives(double* dh1dh6,double* xyz_list, GaussPenta* gauss); + void GetNodalFunctionsMINIDerivatives(double* dh1dh7,double* xyz_list, GaussPenta* gauss); + void GetNodalFunctionsP1DerivativesReference(double* dl1dl6,GaussPenta* gauss); + void GetNodalFunctionsMINIDerivativesReference(double* dl1dl7,GaussPenta* gauss); + void GetJacobian(double* J, double* xyz_list,GaussPenta* gauss); + void GetJacobianDeterminant(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); + void GetBPattyn(double* B, double* xyz_list, GaussPenta* gauss); + void GetBprimePattyn(double* B, double* xyz_list, GaussPenta* gauss); + void GetBStokes(double* B, double* xyz_list, GaussPenta* gauss); + void GetBprimeStokes(double* B_prime, double* xyz_list, GaussPenta* gauss); + void GetBprimeVert(double* B, double* xyz_list, GaussPenta* gauss); + void GetBAdvec(double* B_advec, double* xyz_list, GaussPenta* gauss); + void GetBArtdiff(double* B_artdiff, double* xyz_list, GaussPenta* gauss); + void GetBConduct(double* B_conduct, double* xyz_list, GaussPenta* gauss); + void GetBVert(double* B, double* xyz_list, GaussPenta* gauss); + void GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss); + //void GetLStokes(double* LStokes, double* gauss_tria); + void GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss_tria, GaussPenta* gauss); + void GetParameterValue(double* pvalue,double* plist,GaussPenta* gauss); + //void GetParameterValue(double* pvalue,double* plist,GaussTria* gauss){ISSMERROR("only PentaGauss are supported");}; + void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussPenta* gauss); + //void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussTria* gauss){ISSMERROR("only PentaGauss are supported");}; + }; #endif Index: /issm/trunk/src/c/objects/Inputs/BoolInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/BoolInput.cpp (revision 5646) +++ /issm/trunk/src/c/objects/Inputs/BoolInput.cpp (revision 5647) @@ -172,4 +172,7 @@ void BoolInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");} /*}}}*/ +/*FUNCTION BoolInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/ +void BoolInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR(" not supported yet!");} +/*}}}*/ /*FUNCTION BoolInput::GetParameterValues{{{1*/ void BoolInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){ISSMERROR(" not supported yet!");} @@ -180,4 +183,7 @@ /*FUNCTION BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/ void BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){ISSMERROR(" not supported yet!");} +/*}}}*/ +/*FUNCTION BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){{{1*/ +void BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){ISSMERROR(" not supported yet!");} /*}}}*/ /*FUNCTION BoolInput::ChangeEnum{{{1*/ Index: /issm/trunk/src/c/objects/Inputs/BoolInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/BoolInput.h (revision 5646) +++ /issm/trunk/src/c/objects/Inputs/BoolInput.h (revision 5647) @@ -49,7 +49,9 @@ void GetParameterValue(double* pvalue,double* gauss); void GetParameterValue(double* pvalue,GaussTria* gauss); + void GetParameterValue(double* pvalue,GaussPenta* gauss); void GetParameterValues(double* values,double* gauss_pointers, int numgauss); void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss); void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss); + void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss); void GetParameterAverage(double* pvalue){ISSMERROR("not implemented yet");}; void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss){ISSMERROR("not implemented yet");}; @@ -62,9 +64,9 @@ void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; + void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; + void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; + void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; + void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; + void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; void ChangeEnum(int newenumtype); void SquareMin(double* psquaremin, bool process_units,Parameters* parameters); Index: /issm/trunk/src/c/objects/Inputs/DoubleInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/DoubleInput.cpp (revision 5646) +++ /issm/trunk/src/c/objects/Inputs/DoubleInput.cpp (revision 5647) @@ -185,4 +185,7 @@ void DoubleInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");} /*}}}*/ +/*FUNCTION DoubleInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/ +void DoubleInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR(" not supported yet!");} +/*}}}*/ /*FUNCTION DoubleInput::GetParameterValues{{{1*/ void DoubleInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){ISSMERROR(" not supported yet!");} @@ -193,4 +196,7 @@ /*FUNCTION DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/ void DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){ISSMERROR(" not supported yet!");} +/*}}}*/ +/*FUNCTION DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){{{1*/ +void DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){ISSMERROR(" not supported yet!");} /*}}}*/ /*FUNCTION DoubleInput::ChangeEnum{{{1*/ Index: /issm/trunk/src/c/objects/Inputs/DoubleInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/DoubleInput.h (revision 5646) +++ /issm/trunk/src/c/objects/Inputs/DoubleInput.h (revision 5647) @@ -48,7 +48,9 @@ void GetParameterValue(double* pvalue,double* gauss); void GetParameterValue(double* pvalue,GaussTria* gauss); + void GetParameterValue(double* pvalue,GaussPenta* gauss); void GetParameterValues(double* values,double* gauss_pointers, int numgauss); void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss); void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss); + void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss); void GetParameterAverage(double* pvalue); void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss){ISSMERROR("not implemented yet");}; @@ -61,9 +63,9 @@ void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; + void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; + void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; + void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; + void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; + void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; void ChangeEnum(int newenumtype); void SquareMin(double* psquaremin, bool process_units,Parameters* parameters); Index: /issm/trunk/src/c/objects/Inputs/Input.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/Input.h (revision 5646) +++ /issm/trunk/src/c/objects/Inputs/Input.h (revision 5647) @@ -27,7 +27,9 @@ virtual void GetParameterValue(double* pvalue,double* gauss)=0; virtual void GetParameterValue(double* pvalue,GaussTria* gauss)=0; + virtual void GetParameterValue(double* pvalue,GaussPenta* gauss)=0; virtual void GetParameterValues(double* values,double* gauss_pointers, int numgauss)=0; virtual void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss)=0; virtual void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss)=0; + virtual void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss)=0; virtual void GetParameterAverage(double* pvalue)=0; virtual void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss)=0; @@ -40,9 +42,9 @@ virtual void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss)=0; virtual void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss)=0; - virtual void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussTria* gauss)=0; - virtual void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussTria* gauss)=0; - virtual void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussTria* gauss)=0; - virtual void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussTria* gauss)=0; - virtual void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussTria* gauss)=0; + virtual void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss)=0; + virtual void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss)=0; + virtual void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss)=0; + virtual void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss)=0; + virtual void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss)=0; virtual void ChangeEnum(int newenumtype)=0; Index: /issm/trunk/src/c/objects/Inputs/IntInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/IntInput.cpp (revision 5646) +++ /issm/trunk/src/c/objects/Inputs/IntInput.cpp (revision 5647) @@ -173,4 +173,7 @@ void IntInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");} /*}}}*/ +/*FUNCTION IntInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/ +void IntInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR(" not supported yet!");} +/*}}}*/ /*FUNCTION IntInput::GetParameterValues{{{1*/ void IntInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){ISSMERROR(" not supported yet!");} @@ -181,4 +184,7 @@ /*FUNCTION IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/ void IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){ISSMERROR(" not supported yet!");} +/*}}}*/ +/*FUNCTION IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){{{1*/ +void IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){ISSMERROR(" not supported yet!");} /*}}}*/ /*FUNCTION IntInput::ChangeEnum{{{1*/ Index: /issm/trunk/src/c/objects/Inputs/IntInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/IntInput.h (revision 5646) +++ /issm/trunk/src/c/objects/Inputs/IntInput.h (revision 5647) @@ -49,7 +49,9 @@ void GetParameterValue(double* pvalue,double* gauss); void GetParameterValue(double* pvalue,GaussTria* gauss); + void GetParameterValue(double* pvalue,GaussPenta* gauss); void GetParameterValues(double* values,double* gauss_pointers, int numgauss); void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss); void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss); + void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss); void GetParameterAverage(double* pvalue){ISSMERROR("not implemented yet");}; void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss){ISSMERROR("not implemented yet");}; @@ -62,9 +64,9 @@ void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; + void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; + void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; + void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; + void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; + void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; void ChangeEnum(int newenumtype); void SquareMin(double* psquaremin, bool process_units,Parameters* parameters); Index: /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp (revision 5646) +++ /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp (revision 5647) @@ -184,6 +184,6 @@ } /*}}}*/ -/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/ -void PentaVertexInput::GetParameterValue(double* pvalue,GaussTria* gauss){ +/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/ +void PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ /*Call PentaRef function*/ @@ -216,6 +216,6 @@ } /*}}}*/ -/*FUNCTION PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){{{1*/ -void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){ +/*FUNCTION PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussPenta* gauss){{{1*/ +void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussPenta* gauss){ /*Call PentaRef function*/ @@ -386,4 +386,190 @@ /*FUNCTION PentaVertexInput::GetVyStrainRate3dPattyn{{{1*/ void PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, double* gauss){ + + int i; + const int numgrids=6; + const int NDOF2=2; + double B[5][NDOF2*numgrids]; + double velocity[numgrids][NDOF2]; + + /*Get B matrix: */ + GetBPattyn(&B[0][0], xyz_list, gauss); + + /*Here, we are computing the strain rate of (0,vy)*/ + for(i=0;ivalues[i]; + } + + /*Multiply B by velocity, to get strain rate: */ + MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0, + &velocity[0][0],NDOF2*numgrids,1,0, + epsilonvy,0); + +} +/*}}}*/ +/*FUNCTION PentaVertexInput::GetVxStrainRate3d{{{1*/ +void PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ + int i,j; + + const int numgrids=6; + const int DOFVELOCITY=3; + double B[8][27]; + double B_reduced[6][DOFVELOCITY*numgrids]; + double velocity[numgrids][DOFVELOCITY]; + + /*Get B matrix: */ + GetBStokes(&B[0][0], xyz_list, gauss); + /*Create a reduced matrix of B to get rid of pressure */ + for (i=0;i<6;i++){ + for (j=0;j<3;j++){ + B_reduced[i][j]=B[i][j]; + } + for (j=4;j<7;j++){ + B_reduced[i][j-1]=B[i][j]; + } + for (j=8;j<11;j++){ + B_reduced[i][j-2]=B[i][j]; + } + for (j=12;j<15;j++){ + B_reduced[i][j-3]=B[i][j]; + } + for (j=16;j<19;j++){ + B_reduced[i][j-4]=B[i][j]; + } + for (j=20;j<23;j++){ + B_reduced[i][j-5]=B[i][j]; + } + } + + /*Here, we are computing the strain rate of (vx,0,0)*/ + for(i=0;ivalues[i]; + velocity[i][1]=0.0; + velocity[i][2]=0.0; + } + /*Multiply B by velocity, to get strain rate: */ + MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvx,0); + +} +/*}}}*/ +/*FUNCTION PentaVertexInput::GetVyStrainRate3d{{{1*/ +void PentaVertexInput::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ + int i,j; + + const int numgrids=6; + const int DOFVELOCITY=3; + double B[8][27]; + double B_reduced[6][DOFVELOCITY*numgrids]; + double velocity[numgrids][DOFVELOCITY]; + + /*Get B matrix: */ + GetBStokes(&B[0][0], xyz_list, gauss); + /*Create a reduced matrix of B to get rid of pressure */ + for (i=0;i<6;i++){ + for (j=0;j<3;j++){ + B_reduced[i][j]=B[i][j]; + } + for (j=4;j<7;j++){ + B_reduced[i][j-1]=B[i][j]; + } + for (j=8;j<11;j++){ + B_reduced[i][j-2]=B[i][j]; + } + for (j=12;j<15;j++){ + B_reduced[i][j-3]=B[i][j]; + } + for (j=16;j<19;j++){ + B_reduced[i][j-4]=B[i][j]; + } + for (j=20;j<23;j++){ + B_reduced[i][j-5]=B[i][j]; + } + } + + /*Here, we are computing the strain rate of (0,vy,0)*/ + for(i=0;ivalues[i]; + velocity[i][2]=0.0; + } + /*Multiply B by velocity, to get strain rate: */ + MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvy,0); + +} +/*}}}*/ +/*FUNCTION PentaVertexInput::GetVzStrainRate3d{{{1*/ +void PentaVertexInput::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){ + int i,j; + + const int numgrids=6; + const int DOFVELOCITY=3; + double B[8][27]; + double B_reduced[6][DOFVELOCITY*numgrids]; + double velocity[numgrids][DOFVELOCITY]; + + /*Get B matrix: */ + GetBStokes(&B[0][0], xyz_list, gauss); + /*Create a reduced matrix of B to get rid of pressure */ + for (i=0;i<6;i++){ + for (j=0;j<3;j++){ + B_reduced[i][j]=B[i][j]; + } + for (j=4;j<7;j++){ + B_reduced[i][j-1]=B[i][j]; + } + for (j=8;j<11;j++){ + B_reduced[i][j-2]=B[i][j]; + } + for (j=12;j<15;j++){ + B_reduced[i][j-3]=B[i][j]; + } + for (j=16;j<19;j++){ + B_reduced[i][j-4]=B[i][j]; + } + for (j=20;j<23;j++){ + B_reduced[i][j-5]=B[i][j]; + } + } + + /*Here, we are computing the strain rate of (0,0,vz)*/ + for(i=0;ivalues[i]; + } + + /*Multiply B by velocity, to get strain rate: */ + MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvz,0); + +} +/*}}}*/ +/*FUNCTION PentaVertexInput::GetVxStrainRate3dPattyn{{{1*/ +void PentaVertexInput::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ + + int i; + const int numgrids=6; + const int NDOF2=2; + double B[5][NDOF2*numgrids]; + double velocity[numgrids][NDOF2]; + + /*Get B matrix: */ + GetBPattyn(&B[0][0], xyz_list, gauss); + + /*Here, we are computing the strain rate of (vx,0)*/ + for(i=0;ivalues[i]; + velocity[i][1]=0.0; + } + + /*Multiply B by velocity, to get strain rate: */ + MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0, + &velocity[0][0],NDOF2*numgrids,1,0, + epsilonvx,0); + +} +/*}}}*/ +/*FUNCTION PentaVertexInput::GetVyStrainRate3dPattyn{{{1*/ +void PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ int i; Index: /issm/trunk/src/c/objects/Inputs/PentaVertexInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/PentaVertexInput.h (revision 5646) +++ /issm/trunk/src/c/objects/Inputs/PentaVertexInput.h (revision 5647) @@ -48,8 +48,10 @@ void GetParameterValue(double* pvalue){ISSMERROR("not implemented yet");}; void GetParameterValue(double* pvalue,double* gauss); - void GetParameterValue(double* pvalue,GaussTria* gauss); + void GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR("not implemented yet");}; + void GetParameterValue(double* pvalue,GaussPenta* gauss); void GetParameterValues(double* values,double* gauss_pointers, int numgauss); void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss); - void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss); + void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; + void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss); void GetParameterAverage(double* pvalue); void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss){ISSMERROR("not implemented yet");}; @@ -62,9 +64,9 @@ void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; + void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss); + void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss); + void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss); + void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss); + void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss); void ChangeEnum(int newenumtype); Index: /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h (revision 5646) +++ /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h (revision 5647) @@ -49,7 +49,9 @@ void GetParameterValue(double* pvalue,double* gauss); void GetParameterValue(double* pvalue,GaussTria* gauss); + void GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR("not implemented yet");}; void GetParameterValues(double* values,double* gauss_pointers, int numgauss); void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss); void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss); + void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; void GetParameterAverage(double* pvalue); void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss); @@ -62,9 +64,9 @@ void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss); void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss); - void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; - void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; + void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; + void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; + void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; + void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; + void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; void ChangeEnum(int newenumtype);