Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5637) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5638) @@ -1007,11 +1007,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: */ @@ -1035,5 +1030,5 @@ double l1l2l3[3]; double alpha2complement_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 /* strain rate: */ @@ -1075,9 +1070,6 @@ /*COMPUT alpha2_list (TO BE DELETED)*/ for(i=0;iGetAlphaComplement(&alpha2complement_list[i],&gauss[i][0],VxEnum,VyEnum); - } - - /* 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, 4); + friction->GetAlphaComplement(&alpha2complement_list[i],&gaussgrids[i][0],VxEnum,VyEnum); + } /*Retrieve all inputs we will be needing: */ @@ -1089,38 +1081,36 @@ /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); /*Build alpha_complement_list: */ - //if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxEnum,VyEnum); // TO BE UNCOMMENTED + //if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum); // TO BE UNCOMMENTED //else alpha_complement=0; - TriaRef::GetParameterValue(&alpha_complement,&alpha2complement_list[0],gauss_l1l2l3); // TO BE DELETED + TriaRef::GetParameterValue(&alpha_complement,&alpha2complement_list[0],gauss); // TO BE DELETED /*Recover alpha_complement and k: */ - dragcoefficient_input->GetParameterValue(&drag, gauss_l1l2l3); + dragcoefficient_input->GetParameterValue(&drag, gauss); /*recover lambda and mu: */ - adjointx_input->GetParameterValue(&lambda, gauss_l1l2l3); - adjointy_input->GetParameterValue(&mu, gauss_l1l2l3); + adjointx_input->GetParameterValue(&lambda, gauss); + adjointy_input->GetParameterValue(&mu, gauss); /*recover vx and vy: */ - inputs->GetParameterValue(&vx, gauss_l1l2l3,VxEnum); - inputs->GetParameterValue(&vy, gauss_l1l2l3,VyEnum); + vx_input->GetParameterValue(&vx,gauss); + vy_input->GetParameterValue(&vy,gauss); /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); /* Get nodal functions value at gaussian point:*/ - GetNodalFunctions(l1l2l3, gauss_l1l2l3); + GetNodalFunctions(l1l2l3, gauss); /*Get nodal functions derivatives*/ - GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3); + GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss); /*Get k derivative: dk/dx */ - dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0]); + dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ @@ -1128,8 +1118,8 @@ //standard term dJ/dki - grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss_weight*l1l2l3[i]; + grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss->weight*l1l2l3[i]; //noise dampening d/dki(1/2*(dk/dx)^2) - grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); + grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); } @@ -1141,8 +1131,6 @@ VecSetValues(gradient,numvertices,doflist1,(const double*)grade_g,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; @@ -3672,12 +3660,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; /* surface normal: */ @@ -3704,9 +3686,5 @@ GetDofList(&doflist); - /* 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); - /*Build normal vector to the surface:*/ - x4=xyz_list[0][0]; y4=xyz_list[0][1]; @@ -3737,19 +3715,17 @@ /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); /* Get Jacobian determinant: */ - GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss); //Get L matrix if viscous basal drag present: - GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,NDOF1); + 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*nz; /******************************************************************************/ @@ -3769,8 +3745,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); } @@ -3794,11 +3768,6 @@ /* gaussian points: */ - int num_area_gauss,ig; - double* gauss_weights = NULL; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double gauss_weight; - double gauss_coord[3]; + int ig; + GaussTria *gauss=NULL; /*matrices: */ @@ -3818,22 +3787,18 @@ GetDofList(&doflist); - /* Get gaussian points and weights: */ - GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /* Start looping on the number of gauss (nodes on the bedrock) */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss); /*Get L matrix : */ - GetL(&L[0], &xyz_list[0][0], gauss_coord,NDOF1); + GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); /*Calculate DL on gauss point */ - D_scalar=latentheat/heatcapacity*gauss_weight*Jdet; + D_scalar=latentheat/heatcapacity*gauss->weight*Jdet; /* Do the triple product tL*D*L: */ @@ -3851,8 +3816,6 @@ MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,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); } @@ -3873,11 +3836,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: */ @@ -3933,6 +3891,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) @@ -3944,23 +3904,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); /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); /*Get L matrix: */ - GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); - - DL_scalar=gauss_weight*Jdettria; + GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); + + DL_scalar=gauss->weight*Jdettria; /* Do the triple product tL*D*L: */ @@ -3971,18 +3925,18 @@ /*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=dt*gauss_weight*Jdettria; + DL_scalar=dt*gauss->weight*Jdettria; //Create DL and DLprime matrix @@ -4032,8 +3986,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); } @@ -4054,11 +4006,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: */ @@ -4092,7 +4039,4 @@ GetDofList(&doflist); - /* 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){ @@ -4106,19 +4050,16 @@ /* 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 L matrix: */ - GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); - - DL_scalar=gauss_weight*Jdettria; + GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); + + DL_scalar=gauss->weight*Jdettria; /* Do the triple product tL*D*L: */ @@ -4130,12 +4071,12 @@ /*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 - vxaverage_input->GetParameterValue(&vx,&gauss_l1l2l3[0]); - vyaverage_input->GetParameterValue(&vy,&gauss_l1l2l3[0]); - - DL_scalar=-dt*gauss_weight*Jdettria; + vxaverage_input->GetParameterValue(&vx,gauss); + vyaverage_input->GetParameterValue(&vy,gauss); + + DL_scalar=-dt*gauss->weight*Jdettria; DLprime[0][0]=DL_scalar*vx; @@ -4157,8 +4098,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); } @@ -4231,11 +4170,6 @@ double heatcapacity; - 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_coord[3]; + int ig; + GaussTria *gauss=NULL; /*matrices: */ @@ -4264,21 +4198,18 @@ heatcapacity=matpar->GetHeatCapacity(); - GaussLegendreTria (&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - /* Start looping on the number of gauss (nodes on the bedrock) */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); //Get the Jacobian determinant - GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss_coord); + GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss); /*Get nodal functions values: */ - GetNodalFunctions(&l1l2l3[0], gauss_coord); + GetNodalFunctions(&l1l2l3[0], gauss); /*Calculate DL on gauss point */ - D_scalar=gauss_weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice); + D_scalar=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice); if(dt){ D_scalar=dt*D_scalar; @@ -4299,8 +4230,7 @@ MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,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); } @@ -4321,11 +4251,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; /* matrix */ @@ -4347,7 +4272,4 @@ GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices); GetDofList(&doflist); - - /* 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 inputs :*/ @@ -4357,24 +4279,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 L matrix: */ - GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); + GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); /* Get accumulation, melting at gauss point */ - accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]); - melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]); - dhdt_input->GetParameterValue(&dhdt_g, &gauss_l1l2l3[0]); + accumulation_input->GetParameterValue(&accumulation_g,gauss); + melting_input->GetParameterValue(&melting_g,gauss); + dhdt_input->GetParameterValue(&dhdt_g,gauss); /* Add value into pe_g: */ - for( i=0; iweight*(accumulation_g-melting_g-dhdt_g)*L[i]; } // for (ig=0; igGetInput(AccumulationRateEnum); @@ -4441,24 +4351,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 L matrix: */ - GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); + GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); /* Get accumulation, melting and thickness at gauss point */ - accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]); - melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]); - dhdt_input->GetParameterValue(&dhdt_g, &gauss_l1l2l3[0]); + accumulation_input->GetParameterValue(&accumulation_g,gauss); + melting_input->GetParameterValue(&melting_g,gauss); + dhdt_input->GetParameterValue(&dhdt_g,gauss); /* Add value into pe_g: */ - for( i=0; iweight*(accumulation_g-melting_g-dhdt_g)*L[i]; } // for (ig=0; igGetInput(AccumulationRateEnum); @@ -4522,23 +4420,21 @@ /* 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 L matrix: */ - GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); + GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); /* Get accumulation, melting at gauss point */ - accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]); - melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]); + accumulation_input->GetParameterValue(&accumulation_g,gauss); + melting_input->GetParameterValue(&melting_g,gauss); /* Add value into pe_g: */ - for( i=0; iweight*(accumulation_g-melting_g)*L[i]; } // for (ig=0; igGetInput(MeltingRateEnum); @@ -4616,33 +4501,30 @@ /*For icesheets: */ /* Start looping on the number of gaussian points: */ - for (ig=0; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); /*Get melting at gaussian point: */ - melting_input->GetParameterValue(&meltingvalue, &gauss_l1l2l3[0]); + melting_input->GetParameterValue(&meltingvalue, gauss); /*Get velocity at gaussian point: */ - vx_input->GetParameterValue(&vx, &gauss_l1l2l3[0]); - vy_input->GetParameterValue(&vy, &gauss_l1l2l3[0]); + vx_input->GetParameterValue(&vx, gauss); + vy_input->GetParameterValue(&vy, gauss); /*Get bed slope: */ - bed_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0]); + bed_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss); dbdx=slope[0]; dbdy=slope[1]; /* Get Jacobian determinant: */ - GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss); //Get L matrix if viscous basal drag present: - GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,NDOF1); + GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); /*Build gaussian vector: */ for(i=0;iweight*(vx*dbdx+vy*dbdy-meltingvalue)*L[i]; } @@ -4655,9 +4537,6 @@ VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); - /*Free ressources:*/ - 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); } Index: /issm/trunk/src/c/objects/Elements/TriaRef.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/TriaRef.cpp (revision 5637) +++ /issm/trunk/src/c/objects/Elements/TriaRef.cpp (revision 5638) @@ -487,4 +487,26 @@ } /*}}}*/ +/*FUNCTION TriaRef::GetJacobianDeterminant3d {{{1*/ +void TriaRef::GetJacobianDeterminant3d(double* Jdet, double* xyz_list,GaussTria* gauss){ + /*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,z1,z2,z3; + + x1=*(xyz_list+3*0+0); + y1=*(xyz_list+3*0+1); + z1=*(xyz_list+3*0+2); + x2=*(xyz_list+3*1+0); + y2=*(xyz_list+3*1+1); + z2=*(xyz_list+3*1+2); + x3=*(xyz_list+3*2+0); + y3=*(xyz_list+3*2+1); + z3=*(xyz_list+3*2+2); + + *Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5); + if(*Jdet<0) ISSMERROR("negative jacobian determinant!"); + +} +/*}}}*/ /*FUNCTION TriaRef::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss) {{{1*/ void TriaRef::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss){ Index: /issm/trunk/src/c/objects/Elements/TriaRef.h =================================================================== --- /issm/trunk/src/c/objects/Elements/TriaRef.h (revision 5637) +++ /issm/trunk/src/c/objects/Elements/TriaRef.h (revision 5638) @@ -40,4 +40,5 @@ void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss); void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,double* gauss); + void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,GaussTria* gauss); void GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss); void GetJacobianInvert(double* Jinv, double* xyz_list,GaussTria* gauss);