Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5874) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5875) @@ -2862,25 +2862,22 @@ ElementMatrix* Tria::CreateKMatrixCouplingMacAyealPattynFriction(void){ - const int numdof = 2 *NUMVERTICES; - const int numdoftotal = 4 *NUMVERTICES; - int i,j; - int ig; - int analysis_type; + /*Constants*/ + const int numdof = NDOF2 *NUMVERTICES; + const int numdoftotal = NDOF4 *NUMVERTICES; + + /*Intermediaries */ + int i,j,ig,analysis_type,drag_type; + double Jdet,slope_magnitude,alpha2; double xyz_list[NUMVERTICES][3]; - int numberofdofspernode=2; + double slope[2]={0.0,0.0}; + double MAXSLOPE=.06; // 6 % + double MOUNTAINKEXPONENT=10; + double L[2][numdof]; + double DL[2][2] ={{ 0,0 },{0,0}}; //for basal drag + double DL_scalar; + double Ke_gg[numdof][numdof] ={0.0}; + double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag + Friction *friction = NULL; GaussTria *gauss=NULL; - double L[2][numdof]; - double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag - double DL_scalar; - double Ke_gg[numdof][numdof]={0.0}; - double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag - double Jdet; - double slope[2]={0.0,0.0}; - double slope_magnitude; - Friction *friction = NULL; - double alpha2; - double MAXSLOPE=.06; // 6 % - double MOUNTAINKEXPONENT=10; - int drag_type; /*Initialize Element matrix and return if necessary*/ @@ -2926,5 +2923,5 @@ /*Get L matrix: */ - GetL(&L[0][0], &xyz_list[0][0], gauss,numberofdofspernode); + GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2); @@ -2955,24 +2952,21 @@ ElementMatrix* Tria::CreateKMatrixDiagnosticPattynFriction(void){ + /*Constants*/ const int numdof = NDOF2*NUMVERTICES; - int i,j; - int ig; - int analysis_type; + + /*Intermediaries */ + int i,j,ig; + int analysis_type,drag_type; double xyz_list[NUMVERTICES][3]; - int numberofdofspernode=2; + double slope_magnitude,alpha2,Jdet; + double slope[2]={0.0,0.0}; + double MAXSLOPE=.06; // 6 % + double MOUNTAINKEXPONENT=10; + double L[2][numdof]; + double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag + double DL_scalar; + double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag + Friction *friction = NULL; GaussTria *gauss=NULL; - double L[2][numdof]; - double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag - double DL_scalar; - double Ke_gg[numdof][numdof]={0.0}; - double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag - double Jdet; - double slope[2]={0.0,0.0}; - double slope_magnitude; - Friction *friction = NULL; - double alpha2; - double MAXSLOPE=.06; // 6 % - double MOUNTAINKEXPONENT=10; - int drag_type; /*Initialize Element matrix and return if necessary*/ @@ -2999,17 +2993,16 @@ gauss->GaussPoint(ig); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2); + + surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss); friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT + slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); // If we have a slope > 6% for this element, it means we are on a mountain. In this particular case, //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */ - surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss); - slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); - if (slope_magnitude>MAXSLOPE){ alpha2=pow((double)10,MOUNTAINKEXPONENT); } - - GetL(&L[0][0], &xyz_list[0][0], gauss,numberofdofspernode); - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); DL_scalar=alpha2*gauss->weight*Jdet; @@ -3055,29 +3048,15 @@ void Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg){ - int i,j; - - /* node data: */ - const int numdof=NDOF1*NUMVERTICES; - double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; - - /* gaussian points: */ - int ig; + /*Constants*/ + const int numdof=NDOF1*NUMVERTICES; + + /*Intermediaries */ + int i,j,ig; + int* doflist=NULL; + double xyz_list[NUMVERTICES][3]; + double surface_normal[3]; + double Jdet,DL_scalar; + double L[3]; GaussTria *gauss=NULL; - - /* surface normal: */ - double x4,y4,z4; - double x5,y5,z5; - double x6,y6,z6; - double v46[3]; - double v56[3]; - double normal[3]; - double norm_normal; - double nz; - - /*Matrices: */ - double DL_scalar; - double L[3]; - double Jdet; /* local element matrices: */ @@ -3088,32 +3067,5 @@ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - /*Build normal vector to the surface:*/ - x4=xyz_list[0][0]; - y4=xyz_list[0][1]; - z4=xyz_list[0][2]; - - x5=xyz_list[1][0]; - y5=xyz_list[1][1]; - z5=xyz_list[1][2]; - - x6=xyz_list[2][0]; - y6=xyz_list[2][1]; - z6=xyz_list[2][2]; - - v46[0]=x4-x6; - v46[1]=y4-y6; - v46[2]=z4-z6; - - v56[0]=x5-x6; - v56[1]=y5-y6; - v56[2]=z5-z6; - - normal[0]=(y4-y6)*(z5-z6)-(z4-z6)*(y5-y6); - normal[1]=(z4-z6)*(x5-x6)-(x4-x6)*(z5-z6); - normal[2]=(x4-x6)*(y5-y6)-(y4-y6)*(x5-x6); - - norm_normal=sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2)); - nz=1.0/norm_normal*normal[2]; + SurfaceNormal(&surface_normal[0],xyz_list); /* Start looping on the number of gaussian points: */ @@ -3123,15 +3075,11 @@ gauss->GaussPoint(ig); - /* Get Jacobian determinant: */ GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss); - - //Get L matrix if viscous basal drag present: GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); /**********************Do not forget the sign**********************************/ - DL_scalar=- gauss->weight*Jdet*nz; + DL_scalar=- gauss->weight*Jdet*surface_normal[2]; /******************************************************************************/ - /* Do the triple producte tL*D*L: */ TripleMultiply( L,1,3,1, &DL_scalar,1,1,0, @@ -3139,11 +3087,7 @@ &Ke_gg_gaussian[0][0],0); - /* Add the Ke_gg_gaussian, onto Ke_gg: */ for( i=0; i