Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5636) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5637) @@ -3247,11 +3247,6 @@ /* 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]; + int ig; + GaussTria *gauss=NULL; /* matrices: */ @@ -3274,5 +3269,5 @@ double alpha2; double alpha2_list[numvertices]; //TO BE DELETED - double gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED + double gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED double MAXSLOPE=.06; // 6 % @@ -3314,25 +3309,20 @@ /*COMPUT alpha2_list (TO BE DELETED)*/ for(i=0;iGetAlpha2(&alpha2_list[i],&gauss[i][0],VxEnum,VyEnum,VzEnum); - } - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + friction->GetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum); + } /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); /*Friction: */ - // friction->GetAlpha2(&alpha2, gauss_l1l2l3,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT - TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss_l1l2l3); // TO BE DELETED + // friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT + TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss); // TO BE DELETED // 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_l1l2l3[0]); + surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss); slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); @@ -3342,11 +3332,11 @@ /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); /*Get L matrix: */ - GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); + GetL(&L[0][0], &xyz_list[0][0], gauss,numberofdofspernode); - DL_scalar=alpha2*gauss_weight*Jdet; + DL_scalar=alpha2*gauss->weight*Jdet; for (i=0;i<2;i++){ DL[i][i]=DL_scalar; @@ -3367,11 +3357,9 @@ MatSetValues(Kgg,numdof,doflistp,numdof,doflistm,(const double*)Ke_gg,ADD_VALUES); - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); + /*Clean up and return*/ + delete gauss; + delete friction; xfree((void**)&doflistm); xfree((void**)&doflistp); - delete friction; } /*}}}*/ @@ -3391,11 +3379,6 @@ /* 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]; + int ig; + GaussTria *gauss=NULL; /* matrices: */ @@ -3418,5 +3401,5 @@ double alpha2; double alpha2_list[numvertices]; //TO BE DELETED - double gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED + double gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED double MAXSLOPE=.06; // 6 % @@ -3457,25 +3440,20 @@ /*COMPUT alpha2_list (TO BE DELETED)*/ for(i=0;iGetAlpha2(&alpha2_list[i],&gauss[i][0],VxEnum,VyEnum,VzEnum); - } - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + friction->GetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum); + } /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); /*Friction: */ - // friction->GetAlpha2(&alpha2, gauss_l1l2l3,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT - TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss_l1l2l3); // TO BE DELETED + // friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT + TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss); // TO BE DELETED // 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_l1l2l3[0]); + surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss); slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); @@ -3485,11 +3463,11 @@ /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); /*Get L matrix: */ - GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); + GetL(&L[0][0], &xyz_list[0][0], gauss,numberofdofspernode); - DL_scalar=alpha2*gauss_weight*Jdet; + DL_scalar=alpha2*gauss->weight*Jdet; for (i=0;i<2;i++){ DL[i][i]=DL_scalar; @@ -3509,10 +3487,8 @@ MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); + /*Clean up and return*/ + delete gauss; + delete friction; xfree((void**)&doflist); - delete friction; } /*}}}*/ @@ -3532,11 +3508,6 @@ /* 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]; + int ig; + GaussTria *gauss=NULL; /* matrices: */ @@ -3559,5 +3530,5 @@ double alpha2; double alpha2_list[numvertices]; //TO BE DELETED - double gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED + double gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED double MAXSLOPE=.06; // 6 % @@ -3598,25 +3569,20 @@ /*COMPUT alpha2_list (TO BE DELETED)*/ for(i=0;iGetAlpha2(&alpha2_list[i],&gauss[i][0],VxEnum,VyEnum,VzEnum); - } - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + friction->GetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum); + } /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); /*Friction: */ - // friction->GetAlpha2(&alpha2, gauss_l1l2l3,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT - TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss_l1l2l3); // TO BE DELETED + // friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT + TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss); // TO BE DELETED // 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_l1l2l3[0]); + surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss); slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); @@ -3626,11 +3592,11 @@ /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); /*Get L matrix: */ - GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); + GetL(&L[0][0], &xyz_list[0][0], gauss,numberofdofspernode); - DL_scalar=alpha2*gauss_weight*Jdet; + DL_scalar=alpha2*gauss->weight*Jdet; for (i=0;i<2;i++){ DL[i][i]=DL_scalar; @@ -3650,10 +3616,8 @@ MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); + /*Clean up and return*/ + delete gauss; + delete friction; xfree((void**)&doflist); - delete friction; } /*}}}*/ Index: /issm/trunk/src/c/objects/Gauss/GaussTria.cpp =================================================================== --- /issm/trunk/src/c/objects/Gauss/GaussTria.cpp (revision 5636) +++ /issm/trunk/src/c/objects/Gauss/GaussTria.cpp (revision 5637) @@ -102,4 +102,14 @@ } /*}}}*/ +/*FUNCTION GaussTria::GaussCenter{{{1*/ +void GaussTria::GaussCenter(void){ + + /*update static arrays*/ + coord1=ONETHIRD; + coord2=ONETHIRD; + coord3=ONETHIRD; + +} +/*}}}*/ /*FUNCTION GaussTria::GaussPoint{{{1*/ void GaussTria::GaussPoint(int ig){ Index: /issm/trunk/src/c/objects/Gauss/GaussTria.h =================================================================== --- /issm/trunk/src/c/objects/Gauss/GaussTria.h (revision 5636) +++ /issm/trunk/src/c/objects/Gauss/GaussTria.h (revision 5637) @@ -40,4 +40,5 @@ void GaussPoint(int ig); void GaussVertex(int iv); + void GaussCenter(void); }; #endif /* _GAUSSTRIA_H_ */