Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 3839) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 3840) @@ -545,7 +545,52 @@ } /*}}}*/ -/*FUNCTION Penta::GetSolutionFromInputs(Vec solution, int analysis_type,int sub_analysis_type);{{{1*/ -void Penta::GetSolutionFromInputs(Vec solution, int analysis_type,int sub_analysis_type){ - ISSMERROR(" not supported yet!"); +/*FUNCTION Penta::GetSolutionFromInputs(Vec solution,int analysis_type,int sub_analysis_type){{{1*/ +void Penta::GetSolutionFromInputs(Vec solution,int analysis_type,int sub_analysis_type){ + /*Just branch to the correct UpdateInputsFromSolution generator, according to the type of analysis we are carrying out: */ + if (analysis_type==DiagnosticAnalysisEnum){ + if (sub_analysis_type==HorizAnalysisEnum){ + GetSolutionFromInputsDiagnosticHoriz(solution,analysis_type,sub_analysis_type); + } + else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); + } + else{ + ISSMERROR("%s%i%s\n","analysis: ",analysis_type," not supported yet"); + } +} +/*}}}*/ +/*FUNCTION Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution,int analysis_type,int sub_analysis_type){{{1*/ +void Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution,int analysis_type,int sub_analysis_type){ + + int i; + + const int numvertices=6; + const int numdofpervertex=2; + const int numdof=numdofpervertex*numvertices; + double gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + + int doflist[numdof]; + double values[numdof]; + double vx; + double vy; + + int dummy; + + /*Get dof list: */ + GetDofList(&doflist[0],&dummy); + + /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ + /*P1 element only for now*/ + for(i=0;iGetParameterValue(&vx,&gauss[i][0],VxEnum); + inputs->GetParameterValue(&vy,&gauss[i][0],VyEnum); + values[i*numdofpervertex+0]=vx; + values[i*numdofpervertex+1]=vy; + } + + /*Add value to global vector*/ + VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + } /*}}}*/ @@ -1169,5 +1214,5 @@ double* area_gauss_weights = NULL; double* vert_gauss_weights = NULL; - double gaussgrids[numgrids][numgrids]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; /* specific gaussian point: */ @@ -3064,5 +3109,5 @@ const int numgrids=6; - double gaussgrids[numgrids][numgrids]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; inputs->GetParameterValues(bedlist,&gaussgrids[0][0],6,BedEnum); @@ -3989,5 +4034,5 @@ const int numgrids=6; - double gaussgrids[numgrids][numgrids]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; inputs->GetParameterValues(thickness_list,&gaussgrids[0][0],6,ThicknessEnum); Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 3839) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 3840) @@ -63,4 +63,5 @@ void CreatePVector(Vec pg, int analysis_type,int sub_analysis_type); void GetSolutionFromInputs(Vec solution, int analysis_type,int sub_analysis_type); + void GetSolutionFromInputsDiagnosticHoriz(Vec solution,int analysis_type,int sub_analysis_type); void GetDofList(int* doflist,int* pnumberofdofs); void GetDofList1(int* doflist); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 3839) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 3840) @@ -483,6 +483,6 @@ } /*}}}*/ -/*FUNCTION Tria::GetSolutionFromInputs(Vec solution, int analysis_type,int sub_analysis_type){{{1*/ -void Tria::GetSolutionFromInputs(Vec solution, int analysis_type,int sub_analysis_type){ +/*FUNCTION Tria::GetSolutionFromInputs(Vec solution,int analysis_type,int sub_analysis_type){{{1*/ +void Tria::GetSolutionFromInputs(Vec solution,int analysis_type,int sub_analysis_type){ /*Just branch to the correct UpdateInputsFromSolution generator, according to the type of analysis we are carrying out: */ if (analysis_type==DiagnosticAnalysisEnum){ @@ -497,5 +497,5 @@ } /*}}}*/ -/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHoriz(Vec solution, int analysis_type,int sub_analysis_type){{{1*/ +/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHoriz(Vec solution,int analysis_type,int sub_analysis_type){{{1*/ void Tria::GetSolutionFromInputsDiagnosticHoriz(Vec solution,int analysis_type,int sub_analysis_type){ Index: /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp (revision 3839) +++ /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp (revision 3840) @@ -152,5 +152,17 @@ /*}}}*/ /*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){{{1*/ -void PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){ISSMERROR(" not supported yet!");} +void PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){ + /*P1 interpolation on Gauss point*/ + + /*intermediary*/ + double l1l6[6]; + + /*nodal functions: */ + GetNodalFunctionsP1(&l1l6[0],gauss); + + /*Assign output pointers:*/ + *pvalue=l1l6[0]*values[0]+l1l6[1]*values[1]+l1l6[2]*values[2]+l1l6[3]*values[3]+l1l6[4]*values[4]+l1l6[5]*values[5]; + +} /*}}}*/ /*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,double* gauss,double defaultvalue){{{1*/ @@ -158,8 +170,280 @@ /*}}}*/ /*FUNCTION PentaVertexInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){{{1*/ -void PentaVertexInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){ISSMERROR(" not supported yet!");} +void PentaVertexInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){ + /*It is assumed that output values has been correctly allocated*/ + + int i,j; + double gauss[4]; + + for (i=0;ivalues[0]*dh1dh6[0][0]+this->values[1]*dh1dh6[0][1]+this->values[2]*dh1dh6[0][2]+this->values[3]*dh1dh6[0][3]+this->values[4]*dh1dh6[0][4]+this->values[5]*dh1dh6[0][5]; + p[1]=this->values[0]*dh1dh6[1][0]+this->values[1]*dh1dh6[1][1]+this->values[2]*dh1dh6[1][2]+this->values[3]*dh1dh6[1][3]+this->values[4]*dh1dh6[1][4]+this->values[5]*dh1dh6[1][5]; + p[2]=this->values[0]*dh1dh6[2][0]+this->values[1]*dh1dh6[2][1]+this->values[2]*dh1dh6[2][2]+this->values[3]*dh1dh6[2][3]+this->values[4]*dh1dh6[2][4]+this->values[5]*dh1dh6[2][5]; + +} +/*}}}*/ +/*FUNCTION PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, double* gauss,int formulation_enum) {{{1*/ +void PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, double* gauss,int formulation_enum){ + + if(formulation_enum==StokesAnalysisEnum){ + GetVxStrainRate3dStokes(epsilonvx,xyz_list,gauss); + } + else if(formulation_enum==PattynFormulationEnum){ + GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss); + } + else{ + ISSMERROR("Formulation enum %i (%s) not supported yet",formulation_enum,EnumAsString(formulation_enum)); + } +} +/*}}}*/ +/*FUNCTION PentaVertexInput::GetVxStrainRate3dStokes(double* epsilonvx,double* xyz_list, double* gauss) {{{1*/ +void PentaVertexInput::GetVxStrainRate3dStokes(double* epsilonvx,double* xyz_list, double* gauss){ + int i,j; + + const int numgrids=6; + const int DOFVELOCITY=3; + double B[8][27]; + double B_reduced[6][DOFVELOCITY*numgrids]; + double velocity[3][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::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, double* gauss) {{{1*/ +void PentaVertexInput::GetVxStrainRate3dPattyn(double* epsilonvx,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 (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::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, double* gauss,int formulation_enum) {{{1*/ +void PentaVertexInput::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, double* gauss,int formulation_enum){ + + if(formulation_enum==StokesAnalysisEnum){ + GetVyStrainRate3dStokes(epsilonvy,xyz_list,gauss); + } + else if(formulation_enum==PattynFormulationEnum){ + GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss); + } + else{ + ISSMERROR("Formulation enum %i (%s) not supported yet",formulation_enum,EnumAsString(formulation_enum)); + } +} +/*}}}*/ +/*FUNCTION PentaVertexInput::GetVyStrainRate3dStokes(double* epsilonvy,double* xyz_list, double* gauss) {{{1*/ +void PentaVertexInput::GetVyStrainRate3dStokes(double* epsilonvy,double* xyz_list, double* gauss){ + int i,j; + + const int numgrids=6; + const int DOFVELOCITY=3; + double B[8][27]; + double B_reduced[6][DOFVELOCITY*numgrids]; + double velocity[3][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::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, double* gauss) {{{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::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, double* gauss,int formulation_enum) {{{1*/ +void PentaVertexInput::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, double* gauss,int formulation_enum){ + + if(formulation_enum==StokesAnalysisEnum){ + GetVzStrainRate3dStokes(epsilonvz,xyz_list,gauss); + } + else if(formulation_enum==PattynFormulationEnum){ + GetVzStrainRate3dPattyn(epsilonvz,xyz_list,gauss); + } + else{ + ISSMERROR("Formulation enum %i (%s) not supported yet",formulation_enum,EnumAsString(formulation_enum)); + } +} +/*}}}*/ +/*FUNCTION PentaVertexInput::GetVzStrainRate3dStokes(double* epsilonvz,double* xyz_list, double* gauss) {{{1*/ +void PentaVertexInput::GetVzStrainRate3dStokes(double* epsilonvz,double* xyz_list, double* gauss){ + int i,j; + + const int numgrids=6; + const int DOFVELOCITY=3; + double B[8][27]; + double B_reduced[6][DOFVELOCITY*numgrids]; + double velocity[3][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::GetVzStrainRate3dPattyn(double* epsilonvz,double* xyz_list, double* gauss) {{{1*/ +void PentaVertexInput::GetVzStrainRate3dPattyn(double* epsilonvz,double* xyz_list, double* gauss){ + + /*vz does not contribute to the strain rate in Pattyn/Blatter's model*/ + for (int i=0;i<5;i++){ + epsilonvz[i]=0.0; + } + +} /*}}}*/ /*FUNCTION PentaVertexInput::ChangeEnum(int newenumtype){{{1*/ @@ -173,2 +457,401 @@ } /*}}}*/ + +/*Intermediary*/ +/*FUNCTION PentaVertexInput::GetNodalFunctionsP1 {{{1*/ +void PentaVertexInput::GetNodalFunctionsP1(double* l1l6, double* gauss_coord){ + + /*This routine returns the values of the nodal functions at the gaussian point.*/ + + l1l6[0]=gauss_coord[0]*(1-gauss_coord[3])/2.0; + + l1l6[1]=gauss_coord[1]*(1-gauss_coord[3])/2.0; + + l1l6[2]=gauss_coord[2]*(1-gauss_coord[3])/2.0; + + l1l6[3]=gauss_coord[0]*(1+gauss_coord[3])/2.0; + + l1l6[4]=gauss_coord[1]*(1+gauss_coord[3])/2.0; + + l1l6[5]=gauss_coord[2]*(1+gauss_coord[3])/2.0; + +} +/*}}}*/ +/*FUNCTION PentaVertexInput::GetNodalFunctionsMINI{{{1*/ +void PentaVertexInput::GetNodalFunctionsMINI(double* l1l7, double* gauss_coord){ + + /*This routine returns the values of the nodal functions at the gaussian point.*/ + + /*First nodal function: */ + l1l7[0]=gauss_coord[0]*(1.0-gauss_coord[3])/2.0; + + /*Second nodal function: */ + l1l7[1]=gauss_coord[1]*(1.0-gauss_coord[3])/2.0; + + /*Third nodal function: */ + l1l7[2]=gauss_coord[2]*(1.0-gauss_coord[3])/2.0; + + /*Fourth nodal function: */ + l1l7[3]=gauss_coord[0]*(1.0+gauss_coord[3])/2.0; + + /*Fifth nodal function: */ + l1l7[4]=gauss_coord[1]*(1.0+gauss_coord[3])/2.0; + + /*Sixth nodal function: */ + l1l7[5]=gauss_coord[2]*(1.0+gauss_coord[3])/2.0; + + /*Seventh nodal function: */ + l1l7[6]=27*gauss_coord[0]*gauss_coord[1]*gauss_coord[2]*(1.0+gauss_coord[3])*(1.0-gauss_coord[3]); + +} +/*}}}*/ +/*FUNCTION PentaVertexInput::GetNodalFunctionsP1Derivatives {{{1*/ +void PentaVertexInput::GetNodalFunctionsP1Derivatives(double* dh1dh6,double* xyz_list, double* gauss_coord){ + + /*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_coord); + + /*Get Jacobian invert: */ + GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_coord); + + /*Build dh1dh3: + * + * [dhi/dx]= Jinv*[dhi/dr] + * [dhi/dy] [dhi/ds] + * [dhi/dz] [dhi/dn] + */ + + for (i=0;i