Index: /issm/trunk/src/c/objects/Icefront.cpp =================================================================== --- /issm/trunk/src/c/objects/Icefront.cpp (revision 382) +++ /issm/trunk/src/c/objects/Icefront.cpp (revision 383) @@ -899,5 +899,5 @@ complete_list[j][2]=0; } - tria->GetJacobianDeterminant(&J[i],&complete_list[0][0],l1l2l3_tria); + tria->GetJacobianDeterminant3d(&J[i],&complete_list[0][0],l1l2l3_tria); } Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 382) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 383) @@ -410,5 +410,5 @@ /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant @@ -640,5 +640,5 @@ /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); /*Get L matrix: */ @@ -732,5 +732,5 @@ /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); DL_scalar=gauss_weight*Jdet; @@ -882,5 +882,5 @@ /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); /*Get nodal functions: */ @@ -1014,5 +1014,5 @@ /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); /*Get nodal functions: */ @@ -1194,6 +1194,32 @@ #undef __FUNCT__ -#define __FUNCT__ "Tria::GetJacobianDeterminant" -void Tria::GetJacobianDeterminant(double* Jdet, double* xyz_list,double* gauss_l1l2l3){ +#define __FUNCT__ "Tria::GetJacobianDeterminant2d" +void Tria::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss_l1l2l3){ + + /*The Jacobian determinant is constant over the element, discard the gaussian points. + * J is assumed to have been allocated of size NDOF2xNDOF2.*/ + + double x1,x2,x3,y1,y2,y3; + + x1=*(xyz_list+3*0+0); + y1=*(xyz_list+3*0+1); + x2=*(xyz_list+3*1+0); + y2=*(xyz_list+3*1+1); + x3=*(xyz_list+3*2+0); + y3=*(xyz_list+3*2+1); + + + *Jdet=sqrt(3.0)/6.0*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)); + + + if(Jdet<0){ + printf("%s%s\n",__FUNCT__," error message: negative jacobian determinant!"); + } + +} + +#undef __FUNCT__ +#define __FUNCT__ "Tria::GetJacobianDeterminant3d" +void Tria::GetJacobianDeterminant3d(double* Jdet, double* xyz_list,double* gauss_l1l2l3){ /*The Jacobian determinant is constant over the element, discard the gaussian points. @@ -1664,5 +1690,5 @@ /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); #ifdef _DEBUG_ printf("Element id %i Jacobian determinant: %lf\n",TriaElementGetID(this),Jdet); @@ -1912,5 +1938,5 @@ /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); #ifdef _DEBUG_ printf("Element id %i Jacobian determinant: %lf\n",TriaElementGetID(this),Jdet); @@ -2053,5 +2079,5 @@ /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); /* Get nodal functions value at gaussian point:*/ @@ -2179,5 +2205,5 @@ /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); #ifdef _DEBUG_ printf("Element id %i Jacobian determinant: %lf\n",TriaElementGetID(this),Jdet); @@ -2245,89 +2271,4 @@ - -#undef __FUNCT__ -#define __FUNCT__ "Tria::CreateKMatrixDiagnosticBaseVert" -void Tria::CreateKMatrixDiagnosticBaseVert(Mat Kgg,void* vinputs,int analysis_type){ - - int i,j; - - /* node data: */ - const int numgrids=3; - const int NDOF1=1; - const int numdofs=NDOF1*numgrids; - double xyz_list[numgrids][3]; - int doflist[numdofs]; - int numberofdofspernode; - - /* 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* gauss_weights = NULL; - double gauss_weight; - double gauss_l1l2l3[3]; - - - /* matrices: */ - double L[1][3]; - double DL_scalar; - - double Ke_gg[3][3]; //stiffness matrix - double Ke_gg_gaussian[3][3]; //stiffness matrix evaluated at the gaussian point. - double Jdet; - - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - - /* Get node coordinates and dof list: */ - GetElementNodeData( &xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* Set Ke_gg to 0: */ - for(i=0;i