Index: /issm/trunk/src/c/Container/Inputs.cpp =================================================================== --- /issm/trunk/src/c/Container/Inputs.cpp (revision 5742) +++ /issm/trunk/src/c/Container/Inputs.cpp (revision 5743) @@ -43,6 +43,6 @@ /*Object management*/ -/*FUNCTION Inputs::GetParameterValue(double* pvalue,double* gauss,int enum_type){{{1*/ -void Inputs::GetParameterValue(double* pvalue,double* gauss, int enum_type){ +/*FUNCTION Inputs::GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type){{{1*/ +void Inputs::GetParameterValue(double* pvalue,GaussTria* gauss, int enum_type){ vector::iterator object; @@ -72,5 +72,5 @@ /*}}}*/ /*FUNCTION Inputs::GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type){{{1*/ -void Inputs::GetParameterValue(double* pvalue,GaussTria* gauss, int enum_type){ +void Inputs::GetParameterValue(double* pvalue,GaussPenta* gauss, int enum_type){ vector::iterator object; @@ -99,58 +99,4 @@ } /*}}}*/ -/*FUNCTION Inputs::GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type){{{1*/ -void Inputs::GetParameterValue(double* pvalue,GaussPenta* gauss, int enum_type){ - - vector::iterator object; - Input* input=NULL; - bool found=false; - - /*Go through inputs and check whether any input with the same name is already in: */ - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - input=(Input*)(*object); - if (input->EnumType()==enum_type){ - found=true; - break; - } - } - - if (!found){ - /*we could not find an input with the correct enum type. No defaults values were provided, - * error out: */ - ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumToString(enum_type)); - } - - /*Ok, we have an input if we made it here, request the input to return the values: */ - input->GetParameterValue(pvalue,gauss); - -} -/*}}}*/ -/*FUNCTION Inputs::GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue){{{1*/ -void Inputs::GetParameterValue(double* pvalue,double* gauss, int enum_type,double defaultvalue){ - - vector::iterator object; - Input* input=NULL; - bool found=false; - - /*Go through inputs and check whether any input with the same name is already in: */ - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - input=(Input*)(*object); - if (input->EnumType()==enum_type){ - found=true; - break; - } - } - - if (!found){ - /*we could not find an input with the correct enum type. Return the default value: */ - *pvalue=defaultvalue; - } - else{ - input->GetParameterValue(pvalue,gauss); - } -} -/*}}}*/ /*FUNCTION Inputs::GetParameterValue(bool* pvalue,int enum-type){{{1*/ void Inputs::GetParameterValue(bool* pvalue,int enum_type){ @@ -263,31 +209,4 @@ input->GetParameterAverage(pvalue); -} -/*}}}*/ -/*FUNCTION Inputs::GetParameterDerivativeValue{{{1*/ -void Inputs::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type){ - - vector::iterator object; - Input* input=NULL; - bool found=false; - - /*Go through inputs and check whether any input with the same name is already in: */ - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - input=(Input*)(*object); - if (input->EnumType()==enum_type){ - found=true; - break; - } - } - - if (!found){ - /*we could not find an input with the correct enum type. No defaults values were provided, - * error out: */ - ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumToString(enum_type)); - } - - /*Ok, we have an input if we made it here, request the input to return the value: */ - input->GetParameterDerivativeValue(derivativevalues,xyz_list,gauss); } /*}}}*/ Index: /issm/trunk/src/c/Container/Inputs.h =================================================================== --- /issm/trunk/src/c/Container/Inputs.h (revision 5742) +++ /issm/trunk/src/c/Container/Inputs.h (revision 5743) @@ -48,9 +48,6 @@ void GetParameterValue(int* pvalue,int enum_type); void GetParameterValue(double* pvalue,int enum_type); - void GetParameterValue(double* pvalue,double* gauss,int enum_type); void GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type); void GetParameterValue(double* pvalue,GaussPenta* gauss,int enum_type); - void GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue); - void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type); /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 5742) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 5743) @@ -31,5 +31,4 @@ virtual void CreatePVector(Vec pg)=0; virtual void GetSolutionFromInputs(Vec solution)=0; - virtual void GetNodes(void** nodes)=0; virtual int GetNodeIndex(Node* node)=0; virtual bool GetShelf()=0; Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5742) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5743) @@ -884,18 +884,4 @@ this->results=new Results(); -} -/*}}}*/ -/*FUNCTION Penta::GetNodes {{{1*/ -void Penta::GetNodes(void** vpnodes){ - - int i; - Node** pnodes=NULL; - - /*virtual object: */ - pnodes=(Node**)vpnodes; - - for(i=0;i<6;i++){ - pnodes[i]=nodes[i]; - } } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 5742) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 5743) @@ -78,5 +78,4 @@ void CreatePVector(Vec pg); void DeleteResults(void); - void GetNodes(void** nodes); int GetNodeIndex(Node* node); bool GetOnBed(); @@ -149,11 +148,5 @@ void GetDofList(int** pdoflist,int approximation_enum=0); void GetDofList1(int* doflist); - int GetElementType(void); - void GetNodalFunctions(double* l1l6, double* gauss_coord); - void GetNodalFunctionsDerivatives(double* dh1dh6,double* xyz_list, double* gauss_coord); - void GetNodalFunctionsDerivativesReference(double* dl1dl6,double* gauss_coord); - void GetNodalFunctionsDerivativesReferenceStokes(double* dl1dl7,double* gauss_coord); - void GetNodalFunctionsDerivativesStokes(double* dh1dh7,double* xyz_list, double* gauss_coord); - void GetNodalFunctionsStokes(double* l1l7, double* gauss_coord); + int GetElementType(void); void GetParameterListOnVertices(double* pvalue,int enumtype); void GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue); Index: /issm/trunk/src/c/objects/Elements/PentaRef.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 5742) +++ /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 5743) @@ -53,5 +53,5 @@ /*Reference Element numerics*/ /*FUNCTION PentaRef::GetBMacAyealPattyn {{{1*/ -void PentaRef::GetBMacAyealPattyn(double* B, double* xyz_list, double* gauss){ +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 @@ -91,5 +91,5 @@ /*}}}*/ /*FUNCTION PentaRef::GetBPattyn {{{1*/ -void PentaRef::GetBPattyn(double* B, double* xyz_list, double* gauss){ +void PentaRef::GetBPattyn(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 @@ -136,5 +136,5 @@ /*}}}*/ /*FUNCTION PentaRef::GetBprimePattyn {{{1*/ -void PentaRef::GetBprimePattyn(double* B, double* xyz_list, double* gauss_coord){ +void PentaRef::GetBprimePattyn(double* B, double* xyz_list, GaussPenta* gauss_coord){ /*Compute B prime 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 @@ -180,5 +180,5 @@ /*}}}*/ /*FUNCTION PentaRef::GetBStokes {{{1*/ -void PentaRef::GetBStokes(double* B, double* xyz_list, double* gauss){ +void PentaRef::GetBStokes(double* B, double* xyz_list, GaussPenta* gauss){ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*DOFPERGRID. @@ -251,5 +251,5 @@ /*}}}*/ /*FUNCTION PentaRef::GetBprimeStokes {{{1*/ -void PentaRef::GetBprimeStokes(double* B_prime, double* xyz_list, double* gauss){ +void PentaRef::GetBprimeStokes(double* B_prime, double* xyz_list, GaussPenta* gauss){ /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. * For grid i, Bi' can be expressed in the actual coordinate system @@ -323,5 +323,5 @@ /*}}}*/ /*FUNCTION PentaRef::GetBArtdiff {{{1*/ -void PentaRef::GetBArtdiff(double* B_artdiff, double* xyz_list, double* gauss){ +void PentaRef::GetBArtdiff(double* B_artdiff, double* xyz_list, GaussPenta* gauss){ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. * For grid i, Bi' can be expressed in the actual coordinate system @@ -353,5 +353,5 @@ /*}}}*/ /*FUNCTION PentaRef::GetBAdvec{{{1*/ -void PentaRef::GetBAdvec(double* B_advec, double* xyz_list, double* gauss){ +void PentaRef::GetBAdvec(double* B_advec, double* xyz_list, GaussPenta* gauss){ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. * For grid i, Bi' can be expressed in the actual coordinate system @@ -385,5 +385,5 @@ /*}}}*/ /*FUNCTION PentaRef::GetBConduct{{{1*/ -void PentaRef::GetBConduct(double* B_conduct, double* xyz_list, double* gauss){ +void PentaRef::GetBConduct(double* B_conduct, double* xyz_list, GaussPenta* gauss){ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. * For grid i, Bi' can be expressed in the actual coordinate system @@ -417,5 +417,5 @@ /*}}}*/ /*FUNCTION PentaRef::GetBVert{{{1*/ -void PentaRef::GetBVert(double* B, double* xyz_list, double* gauss){ +void PentaRef::GetBVert(double* B, double* xyz_list, GaussPenta* gauss){ /* Compute B matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz]; where hi is the interpolation function for grid i.*/ @@ -437,5 +437,5 @@ /*}}}*/ /*FUNCTION PentaRef::GetBprimeAdvec{{{1*/ -void PentaRef::GetBprimeAdvec(double* Bprime_advec, double* xyz_list, double* gauss){ +void PentaRef::GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss){ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. * For grid i, Bi' can be expressed in the actual coordinate system @@ -469,5 +469,5 @@ /*}}}*/ /*FUNCTION PentaRef::GetBprimeVert{{{1*/ -void PentaRef::GetBprimeVert(double* B, double* xyz_list, double* gauss){ +void PentaRef::GetBprimeVert(double* B, double* xyz_list, GaussPenta* gauss){ /* Compute Bprime matrix. Bprime=[L1 L2 L3 L4 L5 L6] where Li is the nodal function for grid i*/ @@ -477,5 +477,5 @@ /*}}}*/ /*FUNCTION PentaRef::GetLStokes {{{1*/ -void PentaRef::GetLStokes(double* LStokes, double* gauss){ +void PentaRef::GetLStokes(double* LStokes, GaussPenta* gauss){ /* * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. @@ -507,7 +507,7 @@ /*Get l1l2l3 in actual coordinate system: */ - l1l2l3[0]=gauss[0]; - l1l2l3[1]=gauss[1]; - l1l2l3[2]=gauss[2]; + l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; + l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0; + l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0; /*Build LStokes: */ @@ -574,5 +574,5 @@ /*}}}*/ /*FUNCTION PentaRef::GetLprimeStokes {{{1*/ -void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss){ +void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){ /* @@ -604,7 +604,7 @@ /*Get l1l2l3 in actual coordinate system: */ - l1l2l3[0]=gauss[0]; - l1l2l3[1]=gauss[1]; - l1l2l3[2]=gauss[2]; + l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; + l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0; + l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0; GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss); @@ -673,934 +673,4 @@ /*}}}*/ /*FUNCTION PentaRef::GetJacobian {{{1*/ -void PentaRef::GetJacobian(double* J, double* xyz_list,double* gauss){ - - const int NDOF3=3; - int i,j; - - /*The Jacobian is constant over the element, discard the gaussian points. - * J is assumed to have been allocated of size NDOF2xNDOF2.*/ - - double A1,A2,A3; //area coordinates - double xi,eta,zi; //parametric coordinates - - double x1,x2,x3,x4,x5,x6; - double y1,y2,y3,y4,y5,y6; - double z1,z2,z3,z4,z5,z6; - - /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */ - A1=gauss[0]; - A2=gauss[1]; - A3=gauss[2]; - - xi=A2-A1; - eta=SQRT3*A3; - zi=gauss[3]; - - 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,double* 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,double* 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, double* gauss){ - /*This routine returns the values of the nodal functions at the gaussian point.*/ - - l1l7[0]=gauss[0]*(1.0-gauss[3])/2.0; - l1l7[1]=gauss[1]*(1.0-gauss[3])/2.0; - l1l7[2]=gauss[2]*(1.0-gauss[3])/2.0; - l1l7[3]=gauss[0]*(1.0+gauss[3])/2.0; - l1l7[4]=gauss[1]*(1.0+gauss[3])/2.0; - l1l7[5]=gauss[2]*(1.0+gauss[3])/2.0; - l1l7[6]=27*gauss[0]*gauss[1]*gauss[2]*(1.0+gauss[3])*(1.0-gauss[3]); - -} -/*}}}*/ -/*FUNCTION PentaRef::GetNodalFunctionsMINIDerivatives{{{1*/ -void PentaRef::GetNodalFunctionsMINIDerivatives(double* dh1dh7,double* xyz_list, double* 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;i=0 && A1<=1);//first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*SQRT3); - A2=gauss[1]; ISSMASSERT(A2>=0 && A2<=1);//second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*SQRT3); - A3=gauss[2]; ISSMASSERT(A3>=0 && A3<=1);//third area coordinate value In term of xi and eta: A3=y/SQRT3; - z=gauss[3]; 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,double* 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, double* 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]; - -} -/*}}}*/ - -/*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*(1-gauss->coord4)/2.0; - l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0; - l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0; - - /*Build LStokes: */ - for (i=0;i<3;i++){ - *(LStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i]; - *(LStokes+num_dof*numgrids2d*0+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*0+num_dof*i+2)=0; - *(LStokes+num_dof*numgrids2d*0+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*1+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*1+num_dof*i+2)=0; - *(LStokes+num_dof*numgrids2d*1+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*2+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*2+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*2+num_dof*i+2)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*2+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*3+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*3+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*3+num_dof*i+2)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*3+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*4+num_dof*i)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*4+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*4+num_dof*i+2)=0; - *(LStokes+num_dof*numgrids2d*4+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*5+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*5+num_dof*i+1)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*5+num_dof*i+2)=0; - *(LStokes+num_dof*numgrids2d*5+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*6+num_dof*i)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*6+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*6+num_dof*i+2)=0; - *(LStokes+num_dof*numgrids2d*6+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*7+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*7+num_dof*i+1)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*7+num_dof*i+2)=0; - *(LStokes+num_dof*numgrids2d*7+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*8+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*8+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*8+num_dof*i+2)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*8+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*9+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*9+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*9+num_dof*i+2)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*9+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*10+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*10+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*10+num_dof*i+2)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*10+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*11+num_dof*i)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*11+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*11+num_dof*i+2)=0; - *(LStokes+num_dof*numgrids2d*11+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*12+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*12+num_dof*i+1)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*12+num_dof*i+2)=0; - *(LStokes+num_dof*numgrids2d*12+num_dof*i+3)=0; - *(LStokes+num_dof*numgrids2d*13+num_dof*i)=0; - *(LStokes+num_dof*numgrids2d*13+num_dof*i+1)=0; - *(LStokes+num_dof*numgrids2d*13+num_dof*i+2)=l1l2l3[i]; - *(LStokes+num_dof*numgrids2d*13+num_dof*i+3)=0; - - } -} -/*}}}*/ -/*FUNCTION PentaRef::GetLprimeStokes {{{1*/ -void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){ - - /* - * Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. - * For grid i, Lpi can be expressed in the actual coordinate system - * by: - * Lpi=[ h 0 0 0] - * [ 0 h 0 0] - * [ h 0 0 0] - * [ 0 h 0 0] - * [ 0 0 h 0] - * [ 0 0 h 0] - * [ 0 0 dh/dz 0] - * [ 0 0 dh/dz 0] - * [ 0 0 dh/dz 0] - * [dh/dz 0 dh/dx 0] - * [ 0 dh/dz dh/dy 0] - * [ 0 0 0 h] - * [ 0 0 0 h] - * [ 0 0 0 h] - * where h is the interpolation function for grid i. - */ - int i; - const int numgrids2d=3; - int num_dof=4; - - double l1l2l3[numgrids2d]; - double dh1dh6[3][6]; - - /*Get l1l2l3 in actual coordinate system: */ - l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; - l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0; - l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0; - - GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss); - - /*Build LprimeStokes: */ - for (i=0;i<3;i++){ - *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i]; - *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+2)=0; - *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i]; - *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+2)=0; - *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i)=l1l2l3[i]; - *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+2)=0; - *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+1)=l1l2l3[i]; - *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+2)=0; - *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+2)=l1l2l3[i]; - *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+2)=l1l2l3[i]; - *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+2)=dh1dh6[2][i]; - *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+2)=dh1dh6[2][i]; - *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+2)=dh1dh6[2][i]; - *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i)=dh1dh6[2][i]; - *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+2)=dh1dh6[0][i]; - *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+1)=dh1dh6[2][i]; - *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+2)=dh1dh6[1][i]; - *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+3)=0; - *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+2)=0; - *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+3)=l1l2l3[i]; - *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+2)=0; - *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+3)=l1l2l3[i]; - *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i)=0; - *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+1)=0; - *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+2)=0; - *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+3)=l1l2l3[i]; - - } -} -/*}}}*/ -/*FUNCTION PentaRef::GetJacobian {{{1*/ void PentaRef::GetJacobian(double* J, double* xyz_list,GaussPenta* gauss){ Index: /issm/trunk/src/c/objects/Elements/PentaRef.h =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaRef.h (revision 5742) +++ /issm/trunk/src/c/objects/Elements/PentaRef.h (revision 5743) @@ -23,31 +23,4 @@ /*Numerics*/ - void GetNodalFunctionsP1(double* l1l6, double* gauss); - void GetNodalFunctionsMINI(double* l1l7, double* gauss); - void GetNodalFunctionsP1Derivatives(double* dh1dh6,double* xyz_list, double* gauss); - void GetNodalFunctionsMINIDerivatives(double* dh1dh7,double* xyz_list, double* gauss); - void GetNodalFunctionsP1DerivativesReference(double* dl1dl6,double* gauss); - void GetNodalFunctionsMINIDerivativesReference(double* dl1dl7,double* gauss); - void GetJacobian(double* J, double* xyz_list,double* gauss); - void GetJacobianDeterminant(double* Jdet, double* xyz_list,double* gauss); - void GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss); - void GetBMacAyealPattyn(double* B, double* xyz_list, double* gauss); - void GetBPattyn(double* B, double* xyz_list, double* gauss); - void GetBprimePattyn(double* B, double* xyz_list, double* gauss); - void GetBStokes(double* B, double* xyz_list, double* gauss); - void GetBprimeStokes(double* B_prime, double* xyz_list, double* gauss); - void GetBprimeVert(double* B, double* xyz_list, double* gauss); - void GetBAdvec(double* B_advec, double* xyz_list, double* gauss); - void GetBArtdiff(double* B_artdiff, double* xyz_list, double* gauss); - void GetBConduct(double* B_conduct, double* xyz_list, double* gauss); - void GetBVert(double* B, double* xyz_list, double* gauss); - void GetBprimeAdvec(double* Bprime_advec, double* xyz_list, double* gauss); - void GetLStokes(double* LStokes, double* gauss); - void GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss); - void GetParameterValue(double* pvalue,double* plist,double* gauss); - void GetParameterValue(double* pvalue,double* plist,GaussTria* gauss){ISSMERROR("only PentaGauss are supported");}; - void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, double* gauss); - 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); @@ -75,7 +48,7 @@ void GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss); void GetParameterValue(double* pvalue,double* plist, GaussPenta* gauss); - //void GetParameterValue(double* pvalue,double* plist,GaussTria* gauss){ISSMERROR("only PentaGauss are supported");}; + 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");}; + void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussTria* gauss){ISSMERROR("only PentaGauss are supported");}; }; Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5742) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5743) @@ -807,17 +807,4 @@ } /*}}}*/ -/*FUNCTION Tria::GetNodes {{{1*/ -void Tria::GetNodes(void** vpnodes){ - int i; - Node** pnodes=NULL; - - /*recover nodes: */ - pnodes=(Node**)vpnodes; - - for(i=0;i<3;i++){ - pnodes[i]=nodes[i]; - } -} -/*}}}*/ /*FUNCTION Tria::GetNodeIndex {{{1*/ int Tria::GetNodeIndex(Node* node){ @@ -5502,124 +5489,4 @@ } /*}}}*/ -/*FUNCTION Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype) {{{1*/ -void Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype){ - - /*Output*/ - double value; - - /*Intermediaries*/ - const int numnodes=3; - int grid1=-1,grid2=-1; - int grid3; - int i; - double gauss_tria[numnodes]; - - /*go through 3 nodes (all nodes for tria) and identify 1st and 2nd nodes: */ - ISSMASSERT(nodes); - for(i=0;igrid2){ - - /*Reverse grid1 and grid2*/ - grid3=grid1; grid1=grid2; grid2=grid3; - - /*Change segment gauss point (between -1 and +1)*/ - gauss_seg=-gauss_seg; - } - - /*Build Triangle Gauss point*/ - if (grid1==0 && grid2==1){ - /*gauss_seg is between 0 and 1*/ - gauss_tria[0]=0.5*(1.-gauss_seg); - gauss_tria[1]=1.-0.5*(1.-gauss_seg); - gauss_tria[2]=0.; - } - else if (grid1==0 && grid2==2){ - /*gauss_seg is between 0 and 2*/ - gauss_tria[0]=0.5*(1.-gauss_seg); - gauss_tria[1]=0.; - gauss_tria[2]=1.-0.5*(1.-gauss_seg); - } - else if (grid1==1 && grid2==2){ - /*gauss_seg is between 1 and 2*/ - gauss_tria[0]=0.; - gauss_tria[1]=0.5*(1.-gauss_seg); - gauss_tria[2]=1.-0.5*(1.-gauss_seg); - } - else{ - ISSMERROR("The 2 nodes provided are either the same or did not match current Tria nodes"); - } - - /*OK, now we can call input method*/ - this->inputs->GetParameterValue(&value, &gauss_tria[0],enumtype); - - /*Assign output pointers:*/ - *pvalue=value; -} -/*}}}*/ -/*FUNCTION Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in) {{{1*/ -void Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in){ - - /*Output*/ - double value; - - /*Intermediaries*/ - const int numnodes=3; - int grid1=-1,grid2=-1; - int grid3; - int i; - double gauss_tria[numnodes]; - - /*go through 3 nodes (all nodes for tria) and identify 1st and 2nd nodes: */ - ISSMASSERT(nodes); - for(i=0;igrid2){ - - /*Reverse grid1 and grid2*/ - grid3=grid1; grid1=grid2; grid2=grid3; - - /*Change segment gauss point (between -1 and +1)*/ - gauss_seg=-gauss_seg; - } - - /*Build Triangle Gauss point*/ - if (grid1==0 && grid2==1){ - /*gauss_seg is between 0 and 1*/ - gauss_tria[0]=0.5*(1.-gauss_seg); - gauss_tria[1]=1.-0.5*(1.-gauss_seg); - gauss_tria[2]=0.; - } - else if (grid1==0 && grid2==2){ - /*gauss_seg is between 0 and 2*/ - gauss_tria[0]=0.5*(1.-gauss_seg); - gauss_tria[1]=0.; - gauss_tria[2]=1.-0.5*(1.-gauss_seg); - } - else if (grid1==1 && grid2==2){ - /*gauss_seg is between 1 and 2*/ - gauss_tria[0]=0.; - gauss_tria[1]=0.5*(1.-gauss_seg); - gauss_tria[2]=1.-0.5*(1.-gauss_seg); - } - else{ - ISSMERROR("The 2 nodes provided are either the same or did not match current Tria nodes"); - } - - /*OK, now we can call input method*/ - input_in->GetParameterValue(&value, &gauss_tria[0]); - - /*Assign output pointers:*/ - *pvalue=value; -} -/*}}}*/ /*FUNCTION Tria::GetSolutionFromInputsDiagnosticHoriz{{{1*/ void Tria::GetSolutionFromInputsDiagnosticHoriz(Vec solution){ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 5742) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 5743) @@ -74,5 +74,4 @@ void CreateKMatrix(Mat Kgg); void CreatePVector(Vec pg); - void GetNodes(void** nodes); int GetNodeIndex(Node* node); bool GetOnBed(); @@ -153,7 +152,5 @@ void GetParameterListOnVertices(double* pvalue,int enumtype); void GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue); - void GetParameterValue(double* pvalue,Node* node,int enumtype); - void GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype); - void GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in); + void GetParameterValue(double* pvalue,Node* node,int enumtype); void GetSolutionFromInputsDiagnosticHoriz(Vec solution); void GetSolutionFromInputsDiagnosticHutter(Vec solution); Index: /issm/trunk/src/c/objects/Elements/TriaRef.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/TriaRef.cpp (revision 5742) +++ /issm/trunk/src/c/objects/Elements/TriaRef.cpp (revision 5743) @@ -52,38 +52,4 @@ /*Reference Element numerics*/ -/*FUNCTION TriaRef::GetBMacAyeal {{{1*/ -void TriaRef::GetBMacAyeal(double* B, double* xyz_list, double* gauss){ - /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*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: 3x(NDOF2*numgrids) - */ - - int i; - const int NDOF2=2; - const int numgrids=3; - - double dh1dh3[NDOF2][numgrids]; - - - /*Get dh1dh2dh3 in actual coordinate system: */ - GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss); - - /*Build B: */ - for (i=0;iGetParameterValue(&drag_p,DragPEnum); - inputs->GetParameterValue(&drag_q,DragQEnum); - inputs->GetParameterValue(&thickness, gauss,ThicknessEnum); - inputs->GetParameterValue(&bed, gauss,BedEnum); - inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum); - - /*Get material parameters: */ - gravity=matpar->GetG(); - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); - - //compute r and q coefficients: */ - r=drag_q/drag_p; - s=1./drag_p; - - //From bed and thickness, compute effective pressure when drag is viscous: - Neff=gravity*(rho_ice*thickness+rho_water*bed); - - /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because - the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour - for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should - replace it by Neff=0 (ie, equival it to an ice shelf)*/ - if (Neff<0)Neff=0; - - if(strcmp(element_type,"2d")==0){ - inputs->GetParameterValue(&vx, gauss,vxenum); - inputs->GetParameterValue(&vy, gauss,vyenum); - vmag=sqrt(pow(vx,2)+pow(vy,2)); - } - else if (strcmp(element_type,"3d")==0){ - inputs->GetParameterValue(&vx, gauss,vxenum); - inputs->GetParameterValue(&vy, gauss,vyenum); - inputs->GetParameterValue(&vz, gauss,vzenum); - vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2)); - } - else ISSMERROR("element_type %s not supported yet",element_type); - - alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1)); - - /*Assign output pointers:*/ - *palpha2=alpha2; -} -/*}}}*/ /*FUNCTION Friction::GetAlphaComplement {{{1*/ -void Friction::GetAlphaComplement(double* palpha_complement, double* gauss,int vxenum,int vyenum){ - +void Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum){ + /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. * FrictionGetAlphaComplement is used in control methods on drag, and it computes: @@ -274,5 +213,5 @@ r=drag_q/drag_p; s=1./drag_p; - + //From bed and thickness, compute effective pressure when drag is viscous: Neff=gravity*(rho_ice*thickness+rho_water*bed); @@ -288,5 +227,5 @@ inputs->GetParameterValue(&vy, gauss,vyenum); vmag=sqrt(pow(vx,2)+pow(vy,2)); - + alpha_complement=pow(Neff,r)*pow(vmag,(s-1)); @@ -296,60 +235,2 @@ } /*}}}*/ -/*FUNCTION Friction::GetAlphaComplement {{{1*/ -void Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum){ - - /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. - * FrictionGetAlphaComplement is used in control methods on drag, and it computes: - * alpha_complement= Neff ^r * vel ^s*/ - - /*diverse: */ - int i; - double Neff; - double r,s; - double vx; - double vy; - double vmag; - double drag_p,drag_q; - double drag_coefficient; - double bed,thickness; - double gravity,rho_ice,rho_water; - double alpha_complement; - - /*Recover parameters: */ - inputs->GetParameterValue(&drag_p,DragPEnum); - inputs->GetParameterValue(&drag_q,DragQEnum); - inputs->GetParameterValue(&thickness, gauss,ThicknessEnum); - inputs->GetParameterValue(&bed, gauss,BedEnum); - inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum); - - /*Get material parameters: */ - gravity=matpar->GetG(); - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); - - - //compute r and q coefficients: */ - r=drag_q/drag_p; - s=1./drag_p; - - //From bed and thickness, compute effective pressure when drag is viscous: - Neff=gravity*(rho_ice*thickness+rho_water*bed); - - /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because - the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour - for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should - replace it by Neff=0 (ie, equival it to an ice shelf)*/ - if (Neff<0)Neff=0; - - //We need the velocity magnitude to evaluate the basal stress: - inputs->GetParameterValue(&vx, gauss,vxenum); - inputs->GetParameterValue(&vy, gauss,vyenum); - vmag=sqrt(pow(vx,2)+pow(vy,2)); - - alpha_complement=pow(Neff,r)*pow(vmag,(s-1)); - - /*Assign output pointers:*/ - *palpha_complement=alpha_complement; - -} -/*}}}*/ Index: /issm/trunk/src/c/objects/Loads/Friction.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Friction.h (revision 5742) +++ /issm/trunk/src/c/objects/Loads/Friction.h (revision 5743) @@ -27,8 +27,6 @@ void Echo(void); - void GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum); void GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum); void GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum); - void GetAlphaComplement(double* alpha_complement, double* gauss,int vxenum,int vyenum); void GetAlphaComplement(double* alpha_complement, GaussTria* gauss,int vxenum,int vyenum); Index: /issm/trunk/src/c/objects/Loads/Icefront.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Icefront.cpp (revision 5742) +++ /issm/trunk/src/c/objects/Loads/Icefront.cpp (revision 5743) @@ -549,5 +549,4 @@ /*Objects: */ double pe_g[numdofs]={0.0}; - Node** element_nodes=NULL; int element_type; @@ -579,18 +578,8 @@ //Identify which grids are comprised in the quad: - if(element_type==PentaEnum)element_nodes=(Node**)xmalloc(6*sizeof(Node*)); - element->GetNodes((void**)element_nodes); - - grid1=-1; grid2=-1; grid3=-1; grid4=-1; - for(i=0;i<6;i++){ - if (nodes[0]==element_nodes[i])grid1=i; - if (nodes[1]==element_nodes[i])grid2=i; - if (nodes[2]==element_nodes[i])grid3=i; - if (nodes[3]==element_nodes[i])grid4=i; - } - - if((grid1==-1) || (grid2==-1)|| (grid3==-1)||(grid4==-1)){ - ISSMERROR("could not find element grids corresponding to quad icefront!"); - } + 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: */ @@ -653,5 +642,4 @@ /*Free ressources:*/ - xfree((void**)&element_nodes); xfree((void**)&doflist); @@ -678,5 +666,4 @@ /*Objects: */ double pe_g[numdofs]={0.0}; - Node** element_nodes=NULL; Penta* penta=NULL; @@ -709,18 +696,8 @@ //Identify which grids are comprised in the quad: - element_nodes=(Node**)xmalloc(6*sizeof(Node*)); - element->GetNodes((void**)element_nodes); - - grid1=-1; grid2=-1; grid3=-1; grid4=-1; - for(i=0;i<6;i++){ - if (nodes[0]==element_nodes[i])grid1=i; - if (nodes[1]==element_nodes[i])grid2=i; - if (nodes[2]==element_nodes[i])grid3=i; - if (nodes[3]==element_nodes[i])grid4=i; - } - - if((grid1==-1) || (grid2==-1)|| (grid3==-1)||(grid4==-1)){ - ISSMERROR("could not find element grids corresponding to quad icefront!"); - } + 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: */ @@ -784,5 +761,4 @@ /*Free ressources:*/ - xfree((void**)&element_nodes); xfree((void**)&doflist); @@ -1019,5 +995,5 @@ complete_list[j][2]=0; } - tria->GetJacobianDeterminant2d(&J[i],&complete_list[0][0],l1l2l3_tria); + tria->GetJacobianDeterminant2d(&J[i],&complete_list[0][0],gauss); } @@ -1261,5 +1237,5 @@ complete_list[j][2]=0; } - tria->GetJacobianDeterminant2d(&J[i],&complete_list[0][0],l1l2l3_tria); + tria->GetJacobianDeterminant2d(&J[i],&complete_list[0][0],gauss); } @@ -1310,5 +1286,5 @@ } //for(i=0;i<4;i++) - } //for(ig=0;ig