From 8ae58c830281271b65d8afcce3258ce333b843d2 Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Mon, 14 Jan 2019 13:49:12 +0100 Subject: [PATCH] Protect MEDCouplingFieldDouble::getValueOnMulti from SIGEGV on PENTA6 --- src/INTERP_KERNEL/InterpolationUtils.hxx | 233 +++++++++++------- .../MEDCouplingFieldDiscretization.cxx | 3 +- 2 files changed, 149 insertions(+), 87 deletions(-) diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx index d037a1fa7..1ee1175d5 100644 --- a/src/INTERP_KERNEL/InterpolationUtils.hxx +++ b/src/INTERP_KERNEL/InterpolationUtils.hxx @@ -418,6 +418,113 @@ namespace INTERP_KERNEL bc[2] = 1. - bc[0] - bc[1]; } + inline void barycentric_coords_seg2(const std::vector& n, const double *p, double *bc) + { + double delta=n[0][0]-n[1][0]; + bc[0]=fabs((*p-n[1][0])/delta); + bc[1]=fabs((*p-n[0][0])/delta); + } + + inline void barycentric_coords_tri3(const std::vector& n, const double *p, double *bc) + { + enum { _XX=0, _YY, _ZZ }; + // matrix 2x2 + double + T11 = n[0][_XX]-n[2][_XX], T12 = n[1][_XX]-n[2][_XX], + T21 = n[0][_YY]-n[2][_YY], T22 = n[1][_YY]-n[2][_YY]; + // matrix determinant + double Tdet = T11*T22 - T12*T21; + if ( (std::fabs( Tdet) ) < (std::numeric_limits::min()) ) + { + bc[0]=1; bc[1]=bc[2]=0; // no solution + return; + } + // matrix inverse + double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11; + // vector + double r11 = p[_XX]-n[2][_XX], r12 = p[_YY]-n[2][_YY]; + // barycentric coordinates: multiply matrix by vector + bc[0] = (t11 * r11 + t12 * r12)/Tdet; + bc[1] = (t21 * r11 + t22 * r12)/Tdet; + bc[2] = 1. - bc[0] - bc[1]; + } + + inline void barycentric_coords_quad4(const std::vector& n, const double *p, double *bc) + { + enum { _XX=0, _YY, _ZZ }; + // Find bc by solving system of 3 equations using Gaussian elimination algorithm + // bc1*( x1 - x4 ) + bc2*( x2 - x4 ) + bc3*( x3 - x4 ) = px - x4 + // bc1*( y1 - y4 ) + bc2*( y2 - y4 ) + bc3*( y3 - y4 ) = px - y4 + // bc1*( z1 - z4 ) + bc2*( z2 - z4 ) + bc3*( z3 - z4 ) = px - z4 + + double T[3][4]= + {{ n[0][_XX]-n[3][_XX], n[1][_XX]-n[3][_XX], n[2][_XX]-n[3][_XX], p[_XX]-n[3][_XX] }, + { n[0][_YY]-n[3][_YY], n[1][_YY]-n[3][_YY], n[2][_YY]-n[3][_YY], p[_YY]-n[3][_YY] }, + { n[0][_ZZ]-n[3][_ZZ], n[1][_ZZ]-n[3][_ZZ], n[2][_ZZ]-n[3][_ZZ], p[_ZZ]-n[3][_ZZ] }}; + + if ( !solveSystemOfEquations<3>( T, bc ) ) + bc[0]=1., bc[1] = bc[2] = bc[3] = 0; + else + bc[ 3 ] = 1. - bc[0] - bc[1] - bc[2]; + } + + inline void barycentric_coords_tri6(const std::vector& n, const double *p, double *bc) + { + double matrix2[48]={1., 0., 0., 0., 0., 0., 0., 0., + 1., 0., 0., 0., 0., 0., 1., 0., + 1., 0., 0., 0., 0., 0., 0., 1., + 1., 0., 0., 0., 0., 0., 0.5, 0., + 1., 0., 0., 0., 0., 0., 0.5, 0.5, + 1., 0., 0., 0., 0., 0., 0.,0.5}; + for(int i=0;i<6;i++) + { + matrix2[8*i+1]=n[i][0]; + matrix2[8*i+2]=n[i][1]; + matrix2[8*i+3]=n[i][0]*n[i][0]; + matrix2[8*i+4]=n[i][0]*n[i][1]; + matrix2[8*i+5]=n[i][1]*n[i][1]; + } + double res[12]; + solveSystemOfEquations2<6,2>(matrix2,res,std::numeric_limits::min()); + double refCoo[2]; + refCoo[0]=computeTria6RefBase(res,p); + refCoo[1]=computeTria6RefBase(res+6,p); + computeWeightedCoeffsInTria6FromRefBase(refCoo,bc); + } + + inline void barycentric_coords_tetra10(const std::vector& n, const double *p, double *bc) + { + double matrix2[130]={1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., + 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., + 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., + 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0., 0., + 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0.5, 0., + 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,0.5, 0., + 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, + 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0., 0.5, + 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0.5}; + for(int i=0;i<10;i++) + { + matrix2[13*i+1]=n[i][0]; + matrix2[13*i+2]=n[i][1]; + matrix2[13*i+3]=n[i][2]; + matrix2[13*i+4]=n[i][0]*n[i][0]; + matrix2[13*i+5]=n[i][0]*n[i][1]; + matrix2[13*i+6]=n[i][0]*n[i][2]; + matrix2[13*i+7]=n[i][1]*n[i][1]; + matrix2[13*i+8]=n[i][1]*n[i][2]; + matrix2[13*i+9]=n[i][2]*n[i][2]; + } + double res[30]; + solveSystemOfEquations2<10,3>(matrix2,res,std::numeric_limits::min()); + double refCoo[3]; + refCoo[0]=computeTetra10RefBase(res,p); + refCoo[1]=computeTetra10RefBase(res+10,p); + refCoo[2]=computeTetra10RefBase(res+20,p); + computeWeightedCoeffsInTetra10FromRefBase(refCoo,bc); + } + /*! * Calculate barycentric coordinates of a point p with respect to triangle or tetra vertices. * This method makes 2 assumptions : @@ -427,113 +534,67 @@ namespace INTERP_KERNEL */ inline void barycentric_coords(const std::vector& n, const double *p, double *bc) { - enum { _XX=0, _YY, _ZZ }; switch(n.size()) { case 2: {// SEG 2 - double delta=n[0][0]-n[1][0]; - bc[0]=fabs((*p-n[1][0])/delta); - bc[1]=fabs((*p-n[0][0])/delta); + barycentric_coords_seg2(n,p,bc); break; } case 3: { // TRIA3 - // matrix 2x2 - double - T11 = n[0][_XX]-n[2][_XX], T12 = n[1][_XX]-n[2][_XX], - T21 = n[0][_YY]-n[2][_YY], T22 = n[1][_YY]-n[2][_YY]; - // matrix determinant - double Tdet = T11*T22 - T12*T21; - if ( (std::fabs( Tdet) ) < (std::numeric_limits::min()) ) - { - bc[0]=1; bc[1]=bc[2]=0; // no solution - return; - } - // matrix inverse - double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11; - // vector - double r11 = p[_XX]-n[2][_XX], r12 = p[_YY]-n[2][_YY]; - // barycentric coordinates: multiply matrix by vector - bc[0] = (t11 * r11 + t12 * r12)/Tdet; - bc[1] = (t21 * r11 + t22 * r12)/Tdet; - bc[2] = 1. - bc[0] - bc[1]; + barycentric_coords_tri3(n,p,bc); break; } case 4: { // TETRA4 - // Find bc by solving system of 3 equations using Gaussian elimination algorithm - // bc1*( x1 - x4 ) + bc2*( x2 - x4 ) + bc3*( x3 - x4 ) = px - x4 - // bc1*( y1 - y4 ) + bc2*( y2 - y4 ) + bc3*( y3 - y4 ) = px - y4 - // bc1*( z1 - z4 ) + bc2*( z2 - z4 ) + bc3*( z3 - z4 ) = px - z4 - - double T[3][4]= - {{ n[0][_XX]-n[3][_XX], n[1][_XX]-n[3][_XX], n[2][_XX]-n[3][_XX], p[_XX]-n[3][_XX] }, - { n[0][_YY]-n[3][_YY], n[1][_YY]-n[3][_YY], n[2][_YY]-n[3][_YY], p[_YY]-n[3][_YY] }, - { n[0][_ZZ]-n[3][_ZZ], n[1][_ZZ]-n[3][_ZZ], n[2][_ZZ]-n[3][_ZZ], p[_ZZ]-n[3][_ZZ] }}; - - if ( !solveSystemOfEquations<3>( T, bc ) ) - bc[0]=1., bc[1] = bc[2] = bc[3] = 0; - else - bc[ 3 ] = 1. - bc[0] - bc[1] - bc[2]; + barycentric_coords_quad4(n,p,bc); break; } case 6: { // TRIA6 - double matrix2[48]={1., 0., 0., 0., 0., 0., 0., 0., - 1., 0., 0., 0., 0., 0., 1., 0., - 1., 0., 0., 0., 0., 0., 0., 1., - 1., 0., 0., 0., 0., 0., 0.5, 0., - 1., 0., 0., 0., 0., 0., 0.5, 0.5, - 1., 0., 0., 0., 0., 0., 0.,0.5}; - for(int i=0;i<6;i++) - { - matrix2[8*i+1]=n[i][0]; - matrix2[8*i+2]=n[i][1]; - matrix2[8*i+3]=n[i][0]*n[i][0]; - matrix2[8*i+4]=n[i][0]*n[i][1]; - matrix2[8*i+5]=n[i][1]*n[i][1]; - } - double res[12]; - solveSystemOfEquations2<6,2>(matrix2,res,std::numeric_limits::min()); - double refCoo[2]; - refCoo[0]=computeTria6RefBase(res,p); - refCoo[1]=computeTria6RefBase(res+6,p); - computeWeightedCoeffsInTria6FromRefBase(refCoo,bc); + barycentric_coords_tri6(n,p,bc); break; } case 10: {//TETRA10 - double matrix2[130]={1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., - 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., - 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., - 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0., 0., - 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0.5, 0., - 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,0.5, 0., - 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, - 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0., 0.5, - 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0.5}; - for(int i=0;i<10;i++) - { - matrix2[13*i+1]=n[i][0]; - matrix2[13*i+2]=n[i][1]; - matrix2[13*i+3]=n[i][2]; - matrix2[13*i+4]=n[i][0]*n[i][0]; - matrix2[13*i+5]=n[i][0]*n[i][1]; - matrix2[13*i+6]=n[i][0]*n[i][2]; - matrix2[13*i+7]=n[i][1]*n[i][1]; - matrix2[13*i+8]=n[i][1]*n[i][2]; - matrix2[13*i+9]=n[i][2]*n[i][2]; - } - double res[30]; - solveSystemOfEquations2<10,3>(matrix2,res,std::numeric_limits::min()); - double refCoo[3]; - refCoo[0]=computeTetra10RefBase(res,p); - refCoo[1]=computeTetra10RefBase(res+10,p); - refCoo[2]=computeTetra10RefBase(res+20,p); - computeWeightedCoeffsInTetra10FromRefBase(refCoo,bc); + barycentric_coords_tetra10(n,p,bc); + break; + } + default: + throw INTERP_KERNEL::Exception("INTERP_KERNEL::barycentric_coords : unrecognized simplex !"); + } + } + + inline void barycentric_coords(INTERP_KERNEL::NormalizedCellType ct, const std::vector& n, const double *p, double *bc) + { + switch(ct) + { + case NORM_SEG2: + {// SEG 2 + barycentric_coords_seg2(n,p,bc); + break; + } + case NORM_TRI3: + { // TRIA3 + barycentric_coords_tri3(n,p,bc); + break; + } + case NORM_TETRA4: + { // TETRA4 + barycentric_coords_quad4(n,p,bc); + break; + } + case NORM_TRI6: + { + // TRIA6 + barycentric_coords_tri6(n,p,bc); + break; + } + case NORM_TETRA10: + {//TETRA10 + barycentric_coords_tetra10(n,p,bc); break; } default: diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx index 0c93445c5..acce38a74 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx @@ -1059,7 +1059,8 @@ void MEDCouplingFieldDiscretizationP1::getValueInCell(const MEDCouplingMesh *mes for(std::size_t i=0;i tmp=new double[nbOfNodes]; - INTERP_KERNEL::barycentric_coords(vec,loc,tmp); + INTERP_KERNEL::NormalizedCellType ct(mesh->getTypeOfCell(cellId)); + INTERP_KERNEL::barycentric_coords(ct,vec,loc,tmp); int sz=arr->getNumberOfComponents(); INTERP_KERNEL::AutoPtr tmp2=new double[sz]; std::fill(res,res+sz,0.); -- 2.39.2