-// Copyright (C) 2007-2016 CEA/DEN, EDF R&D
+// Copyright (C) 2007-2021 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
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)
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: 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<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:
+ throw INTERP_KERNEL::Exception("INTERP_KERNEL::barycentric_coords : unrecognized simplex !");
+ }
+ }
+
+ 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:
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])
{
}
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();
{
std::vector<bool> sw(3,false);
double inpVect2[3];
- std::transform(inpVect,inpVect+3,inpVect2,std::ptr_fun<double,double>(fabs));
+ 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)));
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]);
/*! 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);
