Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5635) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5636) @@ -1931,11 +1931,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; /* parameters: */ @@ -2008,31 +2003,17 @@ } - /* 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); - /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); + Jelem+=misfit*Jdet*gauss->weight; + } + + /*clean up and Return: */ + delete gauss; return Jelem; } @@ -2062,11 +2043,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; /* parameters: */ @@ -2083,5 +2059,4 @@ /*inputs: */ bool onwater; - double meanvel, epsvel; @@ -2137,31 +2112,17 @@ } - /* 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); - /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); + Jelem+=misfit*Jdet*gauss->weight; + } + + /*clean-up and Return: */ + delete gauss; return Jelem; } @@ -2191,11 +2152,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; /* parameters: */ @@ -2268,31 +2224,17 @@ } - /* 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); - /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); + Jelem+=misfit*Jdet*gauss->weight; + } + + /*clean-up and Return: */ + delete gauss; return Jelem; } @@ -2322,11 +2264,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; /* parameters: */ @@ -2398,31 +2335,17 @@ } - /* 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); - /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); + Jelem+=misfit*Jdet*gauss->weight; + } + + /*clean-up and Return: */ + delete gauss; return Jelem; } @@ -2434,5 +2357,4 @@ int row; int vertices_ids[3]; - /*recover pointer: */ @@ -2779,11 +2701,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: */ @@ -2839,6 +2756,8 @@ if(artdiff){ //Get the Jacobian determinant - gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; - GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); + gauss=new GaussTria(); + gauss->GaussCenter(); + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); + delete gauss; //Build K matrix (artificial diffusivity matrix) @@ -2850,33 +2769,28 @@ } - /* 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); - /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); /*Get B and B prime matrix: */ - GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss_l1l2l3); - GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3); + GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss); + GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss); //Get vx, vy and their derivatives at gauss point - vxaverage_input->GetParameterValue(&vx, &gauss_l1l2l3[0]); - vyaverage_input->GetParameterValue(&vy, &gauss_l1l2l3[0]); - - vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0]); - vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0]); + vxaverage_input->GetParameterValue(&vx,gauss); + vyaverage_input->GetParameterValue(&vy,gauss); + + vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); + vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); dvxdx=dvx[0]; dvydy=dvy[1]; - DL_scalar=gauss_weight*Jdettria; + DL_scalar=gauss->weight*Jdettria; //Create DL and DLprime matrix @@ -2924,8 +2838,6 @@ 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*/ + delete gauss; xfree((void**)&doflist); } @@ -2945,11 +2857,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: */ @@ -2977,7 +2884,4 @@ this->parameters->FindParam(&dim,DimEnum); - /* 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); - /*Retrieve all inputs we will be needed*/ if(dim==2){ @@ -2987,24 +2891,22 @@ /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); /*Get B and B prime matrix: */ /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/ - GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3); - GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss_l1l2l3); + GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss); + GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss); //Get vx, vy and their derivatives at gauss point - vx_input->GetParameterValue(&vx, &gauss_l1l2l3[0]); - vy_input->GetParameterValue(&vy, &gauss_l1l2l3[0]); - - DL_scalar=-gauss_weight*Jdettria; + vx_input->GetParameterValue(&vx,gauss); + vy_input->GetParameterValue(&vy,gauss); + + DL_scalar=-gauss->weight*Jdettria; DLprime[0][0]=DL_scalar*vx; @@ -3025,8 +2927,6 @@ 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; xfree((void**)&doflist); } @@ -3046,11 +2946,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: */ @@ -3129,6 +3024,8 @@ if(artdiff){ //Get the Jacobian determinant - gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; - GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); + gauss=new GaussTria(); + gauss->GaussCenter(); + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); + delete gauss; //Build K matrix (artificial diffusivity matrix) @@ -3140,33 +3037,28 @@ } - /* 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); - /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); /*Get B and B prime matrix: */ - GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss_l1l2l3); - GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3); + GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss); + GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss); //Get vx, vy and their derivatives at gauss point - vxaverage_input->GetParameterValue(&vx,&gauss_l1l2l3[0]); - vyaverage_input->GetParameterValue(&vy,&gauss_l1l2l3[0]); - - vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0]); - vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0]); + vxaverage_input->GetParameterValue(&vx,gauss); + vyaverage_input->GetParameterValue(&vy,gauss); + + vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); + vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); dvxdx=dvx[0]; dvydy=dvy[1]; - DL_scalar=gauss_weight*Jdettria; + DL_scalar=gauss->weight*Jdettria; DLprime[0][0]=DL_scalar*nx; @@ -3204,8 +3096,6 @@ 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; xfree((void**)&doflist); } @@ -3224,11 +3114,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; /* material data: */ @@ -3281,7 +3166,4 @@ GetDofList(&doflist,MacAyealApproximationEnum); - /* 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); - /*Retrieve all inputs we will be needing: */ thickness_input=inputs->GetInput(ThicknessEnum); @@ -3292,18 +3174,15 @@ /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); /*Compute thickness at gaussian point: */ - thickness_input->GetParameterValue(&thickness, gauss_l1l2l3); + thickness_input->GetParameterValue(&thickness, gauss); /*Get strain rate from velocity: */ - this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss_l1l2l3,vx_input,vy_input); - this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss_l1l2l3,vxold_input,vyold_input); + this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); /*Get viscosity: */ @@ -3312,10 +3191,10 @@ /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant onto this scalar matrix, so that we win some computational time: */ newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); - D_scalar=2*newviscosity*thickness*gauss_weight*Jdet; + D_scalar=2*newviscosity*thickness*gauss->weight*Jdet; for (i=0;i<3;i++){ @@ -3324,6 +3203,6 @@ /*Get B and Bprime matrices: */ - GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss_l1l2l3); - GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3); + GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss); + GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss); /* Do the triple product tB*D*Bprime: */ @@ -3346,8 +3225,7 @@ } - 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; xfree((void**)&doflist); } @@ -6949,5 +6827,4 @@ /*}}}1*/ /*FUNCTION Tria::SurfaceNormal{{{1*/ - void Tria::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){ @@ -6972,5 +6849,4 @@ *(surface_normal+1)=normal[1]/normal_norm; *(surface_normal+2)=normal[2]/normal_norm; - -} -/*}}}*/ +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/TriaRef.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/TriaRef.cpp (revision 5635) +++ /issm/trunk/src/c/objects/Elements/TriaRef.cpp (revision 5636) @@ -149,4 +149,33 @@ } /*}}}*/ +/*FUNCTION TriaRef::GetBPrognostic{{{1*/ +void TriaRef::GetBPrognostic(double* B_prog, double* xyz_list, GaussTria* 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=[ h ] + * [ h ] + * where h is the interpolation function for grid i. + * + * We assume B_prog has been allocated already, of size: 2x(NDOF1*numgrids) + */ + + int i; + const int NDOF1=1; + const int numgrids=3; + + double l1l2l3[numgrids]; + + + /*Get dh1dh2dh3 in actual coordinate system: */ + GetNodalFunctions(&l1l2l3[0],gauss); + + /*Build B_prog: */ + for (i=0;i