-// Copyright (C) 2007-2016 CEA/DEN, EDF R&D
+// Copyright (C) 2007-2019 CEA/DEN, EDF R&D
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
/* 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.;
}
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
* @param quadOut 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(polygOut+SPACEDIM,polygOut+2*SPACEDIM,polygOut+SPACEDIM,std::bind2nd(std::multiplies<double>(),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));
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);
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 vertices.
* This method makes 2 assumptions :
*/
inline void barycentric_coords(const std::vector<const double*>& 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<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: 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][_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<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:
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:
}
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();
// 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)
(XA[0]*XB[1]-XA[1]*XB[0])*XC[2];
}
- /*! 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>
