-// Copyright (C) 2007-2014 CEA/DEN, EDF R&D
+// Copyright (C) 2007-2023 CEA, EDF
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
#include "INTERPKERNELDefines.hxx"
#include "InterpKernelException.hxx"
+#include "VolSurfUser.hxx"
#include "NormalizedUnstructuredMesh.hxx"
/* calcul la surface d'un triangle */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- inline double Surf_Tri(const double* P_1,const double* P_2,const double* P_3)
+ inline double Surf_Tri(const double *P_1,const double *P_2,const double *P_3)
{
double A=(P_3[1]-P_1[1])*(P_2[0]-P_1[0])-(P_2[1]-P_1[1])*(P_3[0]-P_1[0]);
double Surface = 0.5*fabs(A);
//fonction qui calcul le determinant des vecteurs: P3P1 et P3P2
//(cf doc CGAL).
- inline double mon_determinant(const double* P_1,
- const double* P_2,
- const double* P_3)
+ inline double mon_determinant(const double *P_1,
+ const double *P_2,
+ const double *P_3)
{
double mon_det=(P_1[0]-P_3[0])*(P_2[1]-P_3[1])-(P_2[0]-P_3[0])*(P_1[1]-P_3[1]);
return mon_det;
{
double X=P_1[0]-P_2[0];
double Y=P_1[1]-P_2[1];
- double norme=sqrt(X*X+Y*Y);
- return norme;
+ return sqrt(X*X+Y*Y);
+ }
+
+ inline void mid_of_seg2(const double P1[2], const double P2[2], double MID[2])
+ {
+ MID[0]=(P1[0]+P2[0])/2.;
+ MID[1]=(P1[1]+P1[1])/2.;
}
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
double tmp[SPACEDIM];
std::transform(triIn,triIn+SPACEDIM,triIn+SPACEDIM,tmp,std::plus<double>());
//2nd point
- std::transform(tmp,tmp+SPACEDIM,quadOut+SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
+ std::transform(tmp,tmp+SPACEDIM,quadOut+SPACEDIM,std::bind(std::multiplies<double>(),std::placeholders::_1,0.5));
std::transform(tmp,tmp+SPACEDIM,triIn+2*SPACEDIM,tmp,std::plus<double>());
//3rd point
- std::transform(tmp,tmp+SPACEDIM,quadOut+2*SPACEDIM,std::bind2nd(std::multiplies<double>(),1/3.));
+ std::transform(tmp,tmp+SPACEDIM,quadOut+2*SPACEDIM,std::bind(std::multiplies<double>(),std::placeholders::_1,1/3.));
//4th point
std::transform(triIn,triIn+SPACEDIM,triIn+2*SPACEDIM,tmp,std::plus<double>());
- std::transform(tmp,tmp+SPACEDIM,quadOut+3*SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
+ std::transform(tmp,tmp+SPACEDIM,quadOut+3*SPACEDIM,std::bind(std::multiplies<double>(),std::placeholders::_1,0.5));
}
/*!
* This method builds a potentially non-convex polygon cell built with the first point of 'triIn' the barycenter of two edges starting or ending with
* the first point of 'triIn' and the barycenter of 'triIn'.
*
- * @param triIn is a 6 doubles array in full interlace mode, that represents a triangle.
- * @param quadOut is a 8 doubles array filled after the following call.
+ * @param polygIn is a 6 doubles array in full interlace mode, that represents a triangle.
+ * @param polygOut is a 8 doubles array filled after the following call.
*/
template<int SPACEDIM>
- inline void fillDualCellOfPolyg(const double *polygIn, int nPtsPolygonIn, double *polygOut)
+ inline void fillDualCellOfPolyg(const double *polygIn, mcIdType nPtsPolygonIn, double *polygOut)
{
//1st point
std::copy(polygIn,polygIn+SPACEDIM,polygOut);
std::transform(polygIn,polygIn+SPACEDIM,polygIn+SPACEDIM,polygOut+SPACEDIM,std::plus<double>());
//2nd point
- std::transform(polygOut+SPACEDIM,polygOut+2*SPACEDIM,polygOut+SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
+ std::transform(polygOut+SPACEDIM,polygOut+2*SPACEDIM,polygOut+SPACEDIM,std::bind(std::multiplies<double>(),std::placeholders::_1,0.5));
double tmp[SPACEDIM];
//
- for(int i=0;i<nPtsPolygonIn-2;i++)
+ for(mcIdType i=0;i<nPtsPolygonIn-2;i++)
{
std::transform(polygIn,polygIn+SPACEDIM,polygIn+(i+2)*SPACEDIM,tmp,std::plus<double>());
- std::transform(tmp,tmp+SPACEDIM,polygOut+(2*i+3)*SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
+ std::transform(tmp,tmp+SPACEDIM,polygOut+(2*i+3)*SPACEDIM,std::bind(std::multiplies<double>(),std::placeholders::_1,0.5));
std::transform(polygIn+(i+1)*SPACEDIM,polygIn+(i+2)*SPACEDIM,tmp,tmp,std::plus<double>());
- std::transform(tmp,tmp+SPACEDIM,polygOut+(2*i+2)*SPACEDIM,std::bind2nd(std::multiplies<double>(),1./3.));
+ std::transform(tmp,tmp+SPACEDIM,polygOut+(2*i+2)*SPACEDIM,std::bind(std::multiplies<double>(),std::placeholders::_1,1./3.));
}
}
inline std::vector<double> bary_poly(const std::vector<double>& V)
{
std::vector<double> Bary;
- long taille=V.size();
+ std::size_t taille=V.size();
double x=0;
double y=0;
- for(long i=0;i<taille/2;i++)
+ for(std::size_t i=0;i<taille/2;i++)
{
x=x+V[2*i];
y=y+V[2*i+1];
}
- double A=2*x/((double)taille);
- double B=2*y/((double)taille);
+ double A=2*x/(static_cast<double>(taille));
+ double B=2*y/(static_cast<double>(taille));
Bary.push_back(A);//taille vecteur=2*nb de points.
Bary.push_back(B);
/*!
* \brief Solve system equation in matrix form using Gaussian elimination algorithm
- * \param M - N x N+NB_OF_VARS matrix
- * \param sol - vector of N solutions
+ * \param matrix - N x N+NB_OF_VARS matrix
+ * \param solutions - vector of N solutions
+ * \param eps - precision
* \retval bool - true if succeeded
*/
template<unsigned SZ, unsigned NB_OF_RES>
while (n<(int)SZ);
}
s=B[np];//s is the Pivot
- std::transform(B+k*nr,B+(k+1)*nr,B+k*nr,std::bind2nd(std::divides<double>(),s));
+ std::transform(B+k*nr,B+(k+1)*nr,B+k*nr,std::bind(std::divides<double>(),std::placeholders::_1,s));
for(j=0;j<SZ;j++)
{
if(j!=k)
double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11;
// vector
double r11 = p[0]-triaCoords[2*SPACEDIM], r12 = p[1]-triaCoords[2*SPACEDIM+1];
- // barycentric coordinates: mutiply matrix by vector
+ // 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_seg2(const std::vector<const double*>& 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<const double*>& 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<double>::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<const double*>& 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<const double*>& 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<double>::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<const double*>& 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<double>::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 verices.
+ * Calculate barycentric coordinates of a point p with respect to triangle or tetra vertices.
* This method makes 2 assumptions :
* - this is a simplex
* - spacedim == meshdim. For TRI3 and TRI6 spaceDim is expected to be equal to 2 and for TETRA4 spaceDim is expected to be equal to 3.
*/
inline void barycentric_coords(const std::vector<const double*>& n, const double *p, double *bc)
{
- enum { _X, _Y, _Z };
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][_X]-n[2][_X], T12 = n[1][_X]-n[2][_X],
- T21 = n[0][_Y]-n[2][_Y], T22 = n[1][_Y]-n[2][_Y];
- // matrix determinant
- double Tdet = T11*T22 - T12*T21;
- if ( (std::fabs( Tdet) ) < (std::numeric_limits<double>::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[_X]-n[2][_X], r12 = p[_Y]-n[2][_Y];
- // barycentric coordinates: mutiply 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][_X]-n[3][_X], n[1][_X]-n[3][_X], n[2][_X]-n[3][_X], p[_X]-n[3][_X] },
- { n[0][_Y]-n[3][_Y], n[1][_Y]-n[3][_Y], n[2][_Y]-n[3][_Y], p[_Y]-n[3][_Y] },
- { n[0][_Z]-n[3][_Z], n[1][_Z]-n[3][_Z], n[2][_Z]-n[3][_Z], p[_Z]-n[3][_Z] }};
-
- 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<double>::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<double>::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:
}
}
+ inline void barycentric_coords(INTERP_KERNEL::NormalizedCellType ct, const std::vector<const double*>& 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:
+ throw INTERP_KERNEL::Exception("INTERP_KERNEL::barycentric_coords : unrecognized simplex !");
+ }
+ }
+
+ /*!
+ * Compute barycentric coords of point \a point relative to segment S [segStart,segStop] in 3D space.
+ * If point \a point is further from S than eps false is returned.
+ * If point \a point projection on S is outside S false is also returned.
+ * If point \a point is closer from S than eps and its projection inside S true is returned and \a bc0 and \a bc1 output parameter set.
+ */
+ inline bool IsPointOn3DSeg(const double segStart[3], const double segStop[3], const double point[3], double eps, double& bc0, double& bc1)
+ {
+ double AB[3]={segStop[0]-segStart[0],segStop[1]-segStart[1],segStop[2]-segStart[2]},AP[3]={point[0]-segStart[0],point[1]-segStart[1],point[2]-segStart[2]};
+ double l_AB(sqrt(AB[0]*AB[0]+AB[1]*AB[1]+AB[2]*AB[2]));
+ double AP_dot_AB((AP[0]*AB[0]+AP[1]*AB[1]+AP[2]*AB[2])/(l_AB*l_AB));
+ double projOfPOnAB[3]={segStart[0]+AP_dot_AB*AB[0],segStart[1]+AP_dot_AB*AB[1],segStart[2]+AP_dot_AB*AB[2]};
+ double V_dist_P_AB[3]={point[0]-projOfPOnAB[0],point[1]-projOfPOnAB[1],point[2]-projOfPOnAB[2]};
+ double dist_P_AB(sqrt(V_dist_P_AB[0]*V_dist_P_AB[0]+V_dist_P_AB[1]*V_dist_P_AB[1]+V_dist_P_AB[2]*V_dist_P_AB[2]));
+ if(dist_P_AB>=eps)
+ return false;//to far from segment [segStart,segStop]
+ if(AP_dot_AB<-eps || AP_dot_AB>1.+eps)
+ return false;
+ AP_dot_AB=std::max(AP_dot_AB,0.); AP_dot_AB=std::min(AP_dot_AB,1.);
+ bc0=1.-AP_dot_AB; bc1=AP_dot_AB;
+ return true;
+ }
+
+ /*!
+ * Calculate pseudo barycentric coordinates of a point p with respect to the quadrangle vertices.
+ * This method makes the assumption that:
+ * - spacedim == meshdim (2 here).
+ * - the point is within the quad
+ * Quadratic elements are not supported yet.
+ *
+ * A quadrangle can be described as 3 vectors, one point being taken as the origin.
+ * Denoting A, B, C the three other points, any point P within the quad is written as
+ * P = xA+ yC + xy(B-A-C)
+ * This method solve those 2 equations (one per component) for x and y.
+ *
+
+ A------B
+ | |
+ | |
+ 0------C
+ */
+ inline void quad_mapped_coords(const std::vector<const double*>& n, const double *p, double *bc)
+ {
+ double prec = 1.0e-14;
+ enum { _XX=0, _YY, _ZZ };
+
+ if(n.size() != 4)
+ throw INTERP_KERNEL::Exception("INTERP_KERNEL::quad_mapped_coords : unrecognized geometric type! Only QUAD4 supported.");
+
+ double A[2] = {n[1][_XX] - n[0][_XX], n[1][_YY] - n[0][_YY]};
+ double B[2] = {n[2][_XX] - n[0][_XX], n[2][_YY] - n[0][_YY]};
+ double C[2] = {n[3][_XX] - n[0][_XX], n[3][_YY] - n[0][_YY]};
+ double N[2] = {B[_XX] - A[_XX] - C[_XX], B[_YY] - A[_YY] - C[_YY]};
+ double P[2] = {p[_XX] - n[0][_XX], p[_YY] - n[0][_YY]};
+
+ // degenerated case: a rectangle:
+ if (fabs(N[0]) < prec && fabs(N[1]) < prec)
+ {
+ double det = C[0]*A[1] -C[1]*A[0];
+ if (fabs(det) < prec)
+ throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords() has a degenerated 2x2 system!");
+ bc[0] = (P[0]*A[1]-P[1]*A[0])/det;
+ bc[1] = (P[1]*C[0]-P[0]*C[1])/det;
+ return;
+ }
+ double b,c ,a = A[1]*N[0]-A[0]*N[1];
+ bool cas1;
+ if (fabs(a) > 1.0e-14)
+ {
+ b = A[1]*C[0]+N[1]*P[0]-N[0]*P[1]-A[0]*C[1];
+ c = P[0]*C[1] - P[1]*C[0];
+ cas1 = true;
+ }
+ else
+ {
+ a = -C[1]*N[0]+C[0]*N[1];
+ b = A[1]*C[0]-N[1]*P[0]+N[0]*P[1]-A[0]*C[1];
+ c = -P[0]*A[1] + P[1]*A[0];
+ cas1 = false;
+ }
+ double delta = b*b - 4.0*a*c;
+ if (delta < 0.0)
+ throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords(): imaginary solutions!");
+ bc[1] = 0.5*(-b+sqrt(delta))/a;
+ if (bc[1] < -prec || bc[1] > (1.0+prec))
+ bc[1] = 0.5*(-b-sqrt(delta))/a;
+ if (bc[1] < -prec || bc[1] > (1.0+prec))
+ throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords(): point doesn't seem to be in quad4!");
+ if (cas1)
+ {
+ double denom = C[0]+bc[1]*N[0];
+ if (fabs(denom) < prec)
+ throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords(): point doesn't seem to be in quad4!");
+ bc[0] = (P[0]-bc[1]*A[0])/denom;
+ if (bc[0] < -prec || bc[0] > (1.0+prec))
+ throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords(): point doesn't seem to be in quad4!");
+ }
+ else
+ {
+ bc[0] = bc[1];
+ double denom = A[1]+bc[0]*N[1];
+ if (fabs(denom) < prec)
+ throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: cuboid_mapped_coord(): point doesn't seem to be in quad4!");
+ bc[1] = (P[1]-bc[0]*C[1])/denom;
+ if (bc[1] < -prec || bc[1] > (1.0+prec))
+ throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: cuboid_mapped_coord(): point doesn't seem to be in quad4!");
+ }
+ }
+
+ /*!
+ * Doing as in quad_mapped_coords() would lead to a 4th order equation ... So go simpler here:
+ * orthogonal distance to each pair of parallel faces is computed. The ratio gives a number in [0,1]
+ *
+ * Conventions:
+ * - for HEXA8, point F (5) is taken to be the origin (see med file ref connec):
+ * 0 ------ 3
+ /| /|
+ / | / |
+ 1 ------ 2 |
+ | | | |
+ | | | |
+ | 4-----|- 7
+ | / | /
+ 5 ------ 6
+
+ *
+ */
+
+ inline void cuboid_mapped_coords(const std::vector<const double*>& n, const double *p, double *bc)
+ {
+ double prec = 1.0e-14;
+ enum { _XX=0, _YY };
+ if (n.size() != 8)
+ throw INTERP_KERNEL::Exception("INTERP_KERNEL::cuboid_mapped_coords: unrecognized geometric type! Only HEXA8 supported.");
+
+ double dx1, dx2, dy1, dy2, dz1, dz2;
+ dx1 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[4],n[5],n[1]);
+ dx2 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[7],n[3],n[2]);
+
+ dy1 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[5],n[6],n[2]);
+ dy2 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[4],n[0],n[3]);
+
+ dz1 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[5],n[4],n[7]);
+ dz2 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[1],n[2],n[3]);
+
+ if (dx1 < -prec || dx2 < -prec || dy1 < -prec || dy2 < -prec || dz1 < -prec || dz2 < -prec)
+ throw INTERP_KERNEL::Exception("INTERP_KERNEL::cuboid_mapped_coords: point outside HEXA8");
+
+ bc[0] = dx1+dx2 < prec ? 0.5 : dx1/(dx1+dx2);
+ bc[1] = dy1+dy2 < prec ? 0.5 : dy1/(dy1+dy2);
+ bc[2] = dz1+dz2 < prec ? 0.5 : dz1/(dz1+dz2);
+ }
+
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
/* calcul la surface d'un polygone. */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
inline void verif_point_dans_vect(const double* P, std::vector<double>& V, double absolute_precision )
{
- long taille=V.size();
+ std::size_t taille=V.size();
bool isPresent=false;
- for(long i=0;i<taille/2;i++)
+ for(std::size_t i=0;i<taille/2;i++)
{
if (sqrt(((P[0]-V[2*i])*(P[0]-V[2*i])+(P[1]-V[2*i+1])*(P[1]-V[2*i+1])))<absolute_precision)
isPresent=true;
inline void verif_maill_dans_vect(int Num, std::vector<int>& V)
{
- long taille=V.size();
+ std::size_t taille=V.size();
int A=0;
- for(long i=0;i<taille;i++)
+ for(std::size_t i=0;i<taille;i++)
{
if(Num==V[i])
{
}
/*! Function that compares two angles from the values of the pairs (sin,cos)*/
- /*! Angles are considered in [0, 2Pi] bt are not computed explicitely */
+ /*! Angles are considered in [0, 2Pi] bt are not computed explicitly */
class AngleLess
{
public:
- bool operator()(std::pair<double,double>theta1, std::pair<double,double> theta2)
+ bool operator()(std::pair<double,double>theta1, std::pair<double,double> theta2) const
{
double norm1 = sqrt(theta1.first*theta1.first +theta1.second*theta1.second);
double norm2 = sqrt(theta2.first*theta2.first +theta2.second*theta2.second);
}
template<int DIM, NumberingPolicy numPol, class MyMeshType>
- inline void getElemBB(double* bb, const double *coordsOfMesh, int iP, int nb_nodes)
+ inline void getElemBB(double* bb, const double *coordsOfMesh, mcIdType iP, int nb_nodes)
{
bb[0]=std::numeric_limits<double>::max();
bb[1]=-std::numeric_limits<double>::max();
}
}
+ /*!
+ * Find a vector orthogonal to the input vector
+ */
+ inline void orthogonalVect3(const double inpVect[3], double outVect[3])
+ {
+ std::vector<bool> sw(3,false);
+ double inpVect2[3];
+ std::transform(inpVect,inpVect + 3,inpVect2,[](double c){return fabs(c);});
+ std::size_t posMin(std::distance(inpVect2,std::min_element(inpVect2,inpVect2+3)));
+ sw[posMin]=true;
+ std::size_t posMax(std::distance(inpVect2,std::max_element(inpVect2,inpVect2+3)));
+ if(posMax==posMin)
+ { posMax=(posMin+1)%3; }
+ sw[posMax]=true;
+ std::size_t posMid(std::distance(sw.begin(),std::find(sw.begin(),sw.end(),false)));
+ outVect[posMin]=0.; outVect[posMid]=1.; outVect[posMax]=-inpVect[posMid]/inpVect[posMax];
+ }
+
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
/* Computes the dot product of a and b */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
return result;
}
template<class T, int dim>
- inline double distance2( T * a, int inda, T * b, int indb)
+ inline double distance2( T * a, std::size_t inda, T * b, std::size_t indb)
{
double result =0;
for(int idim =0; idim<dim; idim++) result += ((*a)[inda+idim] - (*b)[indb+idim])* ((*a)[inda+idim] - (*b)[indb+idim]);
// just to be able to compile
}
- /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
- /* Checks wether point A is inside the quadrangle BCDE */
- /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
+ /* Checks whether point A is inside the quadrangle BCDE */
+ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
template<int dim> inline double check_inside(const double* A,const double* B,const double* C,const double* D,
const double* E,double* ABC, double* ADE)
}
return result;
}
+
+ /*!
+ * This method normalize input "tetrahedrized polyhedron" to put it around 0,0,0 point and with the normalization factor.
+ *
+ * \param [in,out] ptsOfTetrahedrizedPolyhedron nbfaces * 3 * 3 vector that store in full interlace all the 3 first points of each face of the input polyhedron
+ * \param [in] nbfaces number of faces in the tetrahedrized polyhedron to be normalized
+ * \param [out] centerPt the center of input tetrahedrized polyhedron
+ * \return the normalization factor corresponding to the max amplitude among all nbfaces*3 input points and among X, Y and Z axis.
+ */
+ inline double NormalizeTetrahedrizedPolyhedron(double *ptsOfTetrahedrizedPolyhedron, int nbfaces, double centerPt[3])
+ {
+ centerPt[0] = 0.0; centerPt[1] = 0.0; centerPt[2] = 0.0;
+ double minMax[6]={ std::numeric_limits<double>::max(), -std::numeric_limits<double>::max(),
+ std::numeric_limits<double>::max(), -std::numeric_limits<double>::max(),
+ std::numeric_limits<double>::max(), -std::numeric_limits<double>::max() };
+ for(int iPt = 0 ; iPt < 3 * nbfaces ; ++iPt)
+ {
+ for(int k = 0 ; k < 3 ; ++k)
+ {
+ minMax[2*k] = std::min(minMax[2*k],ptsOfTetrahedrizedPolyhedron[3*iPt+k]);
+ minMax[2*k+1] = std::max(minMax[2*k+1],ptsOfTetrahedrizedPolyhedron[3*iPt+k]);
+ }
+ }
+ double normFact = 0.0;
+ for(int k = 0 ; k < 3 ; ++k)
+ {
+ centerPt[k] = (minMax[2*k] + minMax[2*k+1]) / 2.0 ;
+ normFact = std::max(minMax[2*k+1] - minMax[2*k], normFact);
+ }
+ // renormalize into ptsOfTetrahedrizedPolyhedron
+ for(int iPt = 0 ; iPt < 3 * nbfaces ; ++iPt)
+ {
+ for(int k = 0 ; k < 3 ; ++k)
+ {
+ ptsOfTetrahedrizedPolyhedron[3*iPt+k] = ( ptsOfTetrahedrizedPolyhedron[3*iPt+k] - centerPt[k] ) / normFact;
+ }
+ }
+ return normFact;
+ }
- /*! Computes the triple product (XA^XB).XC (in 3D)*/
- inline double triple_product(const double* A, const double*B, const double*C, const double*X)
+ /*!
+ * Computes the triple product (XA^XB).XC/(norm(XA^XB)) (in 3D)
+ * Returned value is equal to the distance (positive or negative depending of the position of C relative to XAB plane) between XAB plane and C point.
+ */
+ inline double TripleProduct(const double *A, const double *B, const double *C, const double *X)
{
- double XA[3];
- XA[0]=A[0]-X[0];
- XA[1]=A[1]-X[1];
- XA[2]=A[2]-X[2];
- double XB[3];
- XB[0]=B[0]-X[0];
- XB[1]=B[1]-X[1];
- XB[2]=B[2]-X[2];
- double XC[3];
- XC[0]=C[0]-X[0];
- XC[1]=C[1]-X[1];
- XC[2]=C[2]-X[2];
+ double XA[3]={ A[0]-X[0], A[1]-X[1], A[2]-X[2] };
+ double XB[3]={ B[0]-X[0], B[1]-X[1], B[2]-X[2] };
+ double XC[3]={ C[0]-X[0], C[1]-X[1], C[2]-X[2] };
- return
- (XA[1]*XB[2]-XA[2]*XB[1])*XC[0]+
- (XA[2]*XB[0]-XA[0]*XB[2])*XC[1]+
- (XA[0]*XB[1]-XA[1]*XB[0])*XC[2];
+ double XA_cross_XB[3] = {XA[1]*XB[2]-XA[2]*XB[1], XA[2]*XB[0]-XA[0]*XB[2], XA[0]*XB[1]-XA[1]*XB[0]};
+ // norm is equal to double the area of the triangle
+ double norm = std::sqrt(XA_cross_XB[0]*XA_cross_XB[0]+XA_cross_XB[1]*XA_cross_XB[1]+XA_cross_XB[2]*XA_cross_XB[2]);
+
+ return ( XA_cross_XB[0]*XC[0]+ XA_cross_XB[1]*XC[1] + XA_cross_XB[2]*XC[2] ) / norm;
}
- /*! Subroutine of checkEqualPolygins that tests if two list of nodes (not necessarily distincts) describe the same polygon, assuming they share a comon point.*/
+ /*! Subroutine of checkEqualPolygins that tests if two list of nodes (not necessarily distincts) describe the same polygon, assuming they share a common point.*/
/*! Indexes istart1 and istart2 designate two points P1 in L1 and P2 in L2 that have identical coordinates. Generally called with istart1=0.*/
/*! Integer sign ( 1 or -1) indicate the direction used in going all over L2. */
template<class T, int dim>
- bool checkEqualPolygonsOneDirection(T * L1, T * L2, int size1, int size2, int istart1, int istart2, double epsilon, int sign)
+ bool checkEqualPolygonsOneDirection(T * L1, T * L2, std::size_t size1, std::size_t size2, std::size_t istart1, std::size_t istart2, double epsilon, int sign)
{
- int i1 = istart1;
- int i2 = istart2;
- int i1next = ( i1 + 1 ) % size1;
- int i2next = ( i2 + sign +size2) % size2;
-
+ std::size_t i1 = istart1;
+ std::size_t i2 = istart2;
+ std::size_t i1next = ( i1 + 1 ) % size1;
+ std::size_t i2next = ( i2 + sign +size2) % size2;
+
while(true)
{
while( i1next!=istart1 && distance2<T,dim>(L1,i1*dim, L1,i1next*dim) < epsilon ) i1next = ( i1next + 1 ) % size1;
while( i2next!=istart2 && distance2<T,dim>(L2,i2*dim, L2,i2next*dim) < epsilon ) i2next = ( i2next + sign +size2 ) % size2;
-
+
if(i1next == istart1)
{
if(i2next == istart2)
std::cout << "Warning InterpolationUtils.hxx:checkEqualPolygonsPointer: Null pointer " << std::endl;
throw(Exception("big error: not closed polygon..."));
}
-
- int size1 = (*L1).size()/dim;
- int size2 = (*L2).size()/dim;
- int istart1 = 0;
- int istart2 = 0;
-
+
+ std::size_t size1 = (*L1).size()/dim;
+ std::size_t size2 = (*L2).size()/dim;
+ std::size_t istart1 = 0;
+ std::size_t istart2 = 0;
+
while( istart2 < size2 && distance2<T,dim>(L1,istart1*dim, L2,istart2*dim) > epsilon ) istart2++;
-
+
if(istart2 == size2)
{
return (size1 == 0) && (size2 == 0);
