-// Copyright (C) 2007-2019 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
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)
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]);
}
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 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);
