-// Copyright (C) 2007-2016 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 "CellModel.hxx"
#include "BBTree.txx"
+#include "InterpKernelGeo2DNode.hxx"
+#include "InterpKernelGeo2DQuadraticPolygon.hxx"
+
#include <list>
+#include <vector>
#include <set>
#include <limits>
+#include <memory>
namespace INTERP_KERNEL
{
{
public:
virtual ~GenericPointLocatorAlgos() { }
- virtual std::list<int> locates(const double* x, double eps) = 0;
+ virtual std::list<mcIdType> locates(const double* x, double eps) = 0;
};
template<class MyMeshType>
typedef typename MyMeshType::MyConnType ConnType;
const int SPACEDIM=MyMeshType::MY_SPACEDIM;
const NumberingPolicy numPol=MyMeshType::My_numPol;
- int nelem = _mesh.getNumberOfElements();
+ ConnType nelem = _mesh.getNumberOfElements();
_bb = new double[SPACEDIM*2*nelem];
const ConnType* conn = _mesh.getConnectivityPtr();
const ConnType* conn_index = _mesh.getConnectivityIndexPtr();
const double* coords=_mesh.getCoordinatesPtr();
- for (int i=0; i<nelem; i++)
+ for (ConnType i=0; i<nelem; i++)
{
for (int idim=0; idim<SPACEDIM; idim++)
{
_bb[2*(i*SPACEDIM+idim)]=std::numeric_limits<double>::max();
_bb[2*(i*SPACEDIM+idim)+1]=-std::numeric_limits<double>::max();
}
- for (int index= conn_index[i]; index < conn_index[i+1];index++)
+ for (ConnType index= conn_index[i]; index < conn_index[i+1];index++)
{
//coordelem points to the coordinates of the current node of the i-th element
const double* coordelem = coords+OTT<ConnType,numPol>::ind2C(conn[OTT<ConnType,numPol>::ind2C(index)])*SPACEDIM;
std::list<ConnType> retlist;
for(unsigned int i=0; i< candidates.size(); i++)
{
- int ielem=candidates[i];
+ ConnType ielem=candidates[i];
if (elementContainsPoint(ielem,x,eps))
retlist.push_back(OTT<ConnType,numPol>::indFC(ielem));
}
return retlist;
}
- static bool isElementContainsPointAlg2D(const double *ptToTest, const double *cellPts, int nbEdges, double eps)
+ static bool isElementContainsPointAlg2DSimple(const double *ptToTest, const double *cellPts, mcIdType nbEdges, double eps)
{
/* with dimension 2, it suffices to check all the edges
and see if the sign of double products from the point
here XA^XC and XC^XB have different signs*/
const int SPACEDIM=MyMeshType::MY_SPACEDIM;
- int* sign = new int[nbEdges];
- for (int iedge=0; iedge<nbEdges; iedge++)
+ std::unique_ptr<char[]> sign( new char[nbEdges] );
+ for (mcIdType iedge=0; iedge<nbEdges; iedge++)
{
const double* A=cellPts+SPACEDIM*iedge;
const double* B=cellPts+SPACEDIM*((iedge+1)%nbEdges);
else
sign[iedge]=0;
}
- bool ret=decideFromSign(sign, nbEdges);
- delete [] sign;
+ bool ret=decideFromSign(sign.get(), nbEdges);
+ return ret;
+ }
+
+ // Same as isElementContainsPointAlg2DSimple() with a different input format ...
+ static bool isElementContainsPointAlgo2DSimple2(const double *ptToTest, NormalizedCellType type,
+ const double *coords, const typename MyMeshType::MyConnType *conn_elem,
+ typename MyMeshType::MyConnType conn_elem_sz, double eps)
+ {
+ const int SPACEDIM=MyMeshType::MY_SPACEDIM;
+ typedef typename MyMeshType::MyConnType ConnType;
+ const NumberingPolicy numPol=MyMeshType::My_numPol;
+
+ const CellModel& cmType=CellModel::GetCellModel(type);
+ bool ret=false;
+
+ int nbEdges=cmType.getNumberOfSons();
+ double *pts = new double[nbEdges*SPACEDIM];
+ for (int iedge=0; iedge<nbEdges; iedge++)
+ {
+ const double* a=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[iedge]));
+ std::copy(a,a+SPACEDIM,pts+iedge*SPACEDIM);
+ }
+ ret=isElementContainsPointAlg2DSimple(ptToTest,pts,nbEdges,eps);
+ delete [] pts;
return ret;
}
- static bool isElementContainsPointAlg3D(const double *ptToTest, const typename MyMeshType::MyConnType *conn_elem, int conn_elem_sz, const double *coords, const CellModel& cmType, double eps)
+ static bool isElementContainsPointAlgo2DPolygon(const double *ptToTest, NormalizedCellType type,
+ const double *coords, const typename MyMeshType::MyConnType *conn_elem,
+ typename MyMeshType::MyConnType conn_elem_sz, double eps)
+ {
+ // Override precision for this method only:
+ INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
+
+ const int SPACEDIM=MyMeshType::MY_SPACEDIM;
+ typedef typename MyMeshType::MyConnType ConnType;
+ const NumberingPolicy numPol=MyMeshType::My_numPol;
+
+ std::vector<INTERP_KERNEL::Node *> nodes(conn_elem_sz);
+ INTERP_KERNEL::QuadraticPolygon *pol(0);
+ for(mcIdType j=0;j<conn_elem_sz;j++)
+ {
+ mcIdType nodeId(OTT<ConnType,numPol>::ind2C(conn_elem[j]));
+ nodes[j]=new INTERP_KERNEL::Node(coords[nodeId*SPACEDIM],coords[nodeId*SPACEDIM+1]);
+ }
+ if(!INTERP_KERNEL::CellModel::GetCellModel(type).isQuadratic())
+ pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes);
+ else
+ pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes);
+ INTERP_KERNEL::Node *n(new INTERP_KERNEL::Node(ptToTest[0],ptToTest[1]));
+ double a(0.),b(0.),c(0.);
+ a=pol->normalizeMe(b,c); n->applySimilarity(b,c,a);
+ bool ret=pol->isInOrOut2(n);
+ delete pol; n->decrRef();
+ return ret;
+ }
+
+ static bool isElementContainsPointAlg3D(const double *ptToTest, const typename MyMeshType::MyConnType *conn_elem, typename MyMeshType::MyConnType conn_elem_sz, const double *coords, const CellModel& cmType, double eps)
{
const int SPACEDIM=MyMeshType::MY_SPACEDIM;
typedef typename MyMeshType::MyConnType ConnType;
const NumberingPolicy numPol=MyMeshType::My_numPol;
int nbfaces = cmType.getNumberOfSons2(conn_elem,conn_elem_sz);
- int *sign = new int[nbfaces];
- int *connOfSon = new int[conn_elem_sz];
+ std::unique_ptr<char[]> sign( new char[nbfaces] );
+ std::unique_ptr<ConnType[]> connOfSon( new ConnType[conn_elem_sz] );
+ std::unique_ptr<double[]> ptsOfTetrahedrizedPolyhedron( new double[3*3*nbfaces] );
for (int iface=0; iface<nbfaces; iface++)
{
NormalizedCellType typeOfSon;
- cmType.fillSonCellNodalConnectivity2(iface,conn_elem,conn_elem_sz,connOfSon,typeOfSon);
+ cmType.fillSonCellNodalConnectivity2(iface,conn_elem,conn_elem_sz,connOfSon.get(),typeOfSon);
const double* AA=coords+SPACEDIM*(OTT<ConnType,numPol>::coo2C(connOfSon[0]));
+ std::copy(AA,AA+3,ptsOfTetrahedrizedPolyhedron.get()+9*iface);
const double* BB=coords+SPACEDIM*(OTT<ConnType,numPol>::coo2C(connOfSon[1]));
+ std::copy(BB,BB+3,ptsOfTetrahedrizedPolyhedron.get()+9*iface+3);
const double* CC=coords+SPACEDIM*(OTT<ConnType,numPol>::coo2C(connOfSon[2]));
- double Vol=triple_product(AA,BB,CC,ptToTest);
- if (Vol<-eps)
+ std::copy(CC,CC+3,ptsOfTetrahedrizedPolyhedron.get()+9*iface+6);
+ }
+ double centerPt[3];
+ double normalizeFact = NormalizeTetrahedrizedPolyhedron(ptsOfTetrahedrizedPolyhedron.get(),nbfaces,centerPt);
+ double ptToTestNorm[3] = {(ptToTest[0]-centerPt[0])/normalizeFact,
+ (ptToTest[1]-centerPt[1])/normalizeFact,
+ (ptToTest[2]-centerPt[2])/normalizeFact};
+ for (int iface=0; iface<nbfaces; iface++)
+ {
+ double lengthNorm = ( TripleProduct(
+ ptsOfTetrahedrizedPolyhedron.get() + 9*iface + 3,
+ ptsOfTetrahedrizedPolyhedron.get() + 9*iface + 6,
+ ptToTestNorm,
+ ptsOfTetrahedrizedPolyhedron.get() + 9*iface) );
+ // assigne sign[iface] : 0 means on the face within eps, 1 or -1 means far from tetra representing face
+ if( lengthNorm<-eps )
sign[iface]=-1;
- else if (Vol>eps)
+ else if( lengthNorm>eps )
sign[iface]=1;
- else
- sign[iface]=0;
+ else
+ sign[iface]=0;
}
- bool ret=decideFromSign(sign, nbfaces);
- delete [] sign;
- delete [] connOfSon;
+ bool ret=decideFromSign(sign.get(), nbfaces);
return ret;
}
- static bool isElementContainsPoint(const double *ptToTest, NormalizedCellType type, const double *coords, const typename MyMeshType::MyConnType *conn_elem, int conn_elem_sz, double eps)
+ /*!
+ * Precondition : spacedim==meshdim. To be checked upstream to this call.
+ */
+ static bool isElementContainsPoint(const double *ptToTest, NormalizedCellType type, const double *coords,
+ const typename MyMeshType::MyConnType *conn_elem,
+ typename MyMeshType::MyConnType conn_elem_sz, double eps)
{
const int SPACEDIM=MyMeshType::MY_SPACEDIM;
typedef typename MyMeshType::MyConnType ConnType;
//
if (SPACEDIM==2)
{
- int nbEdges=cmType.getNumberOfSons();
- double *pts = new double[nbEdges*SPACEDIM];
- for (int iedge=0; iedge<nbEdges; iedge++)
- {
- const double* a=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[iedge]));
- std::copy(a,a+SPACEDIM,pts+iedge*SPACEDIM);
- }
- bool ret=isElementContainsPointAlg2D(ptToTest,pts,nbEdges,eps);
- delete [] pts;
- return ret;
+ if(type != INTERP_KERNEL::NORM_POLYGON && !cmType.isQuadratic())
+ return isElementContainsPointAlgo2DSimple2(ptToTest, type, coords, conn_elem, conn_elem_sz, eps);
+ else
+ return isElementContainsPointAlgo2DPolygon(ptToTest, type, coords, conn_elem, conn_elem_sz, eps);
}
-
if (SPACEDIM==3)
{
return isElementContainsPointAlg3D(ptToTest,conn_elem,conn_elem_sz,coords,cmType,eps);
return isElementContainsPoint(x,type,coords,conn_elem,conn_elem_sz,eps);
}
- static bool decideFromSign(const int* sign, int nbelem)
+ static bool decideFromSign(const char *sign, mcIdType nbelem)
{
- int min_sign = 1;
- int max_sign = -1;
+ char min_sign = 1;
+ char max_sign = -1;
for (int i=0; i<nbelem;i++)
{
min_sign=(sign[i]<min_sign)?sign[i]:min_sign;
*/
//================================================================================
- virtual std::list<int> locates(const double* x, double eps)
+ virtual std::list<typename MyMeshType::MyConnType> locates(const double* x, double eps)
{
typedef typename MyMeshType::MyConnType ConnType;
const NumberingPolicy numPol=MyMeshType::My_numPol;
- std::list<int> simplexNodes;
- std::list<int> candidates = PointLocatorAlgos<MyMeshType>::locates(x,eps);
- std::list<int>::iterator eIt = candidates.begin();
+ std::list<ConnType> simplexNodes;
+ std::list<ConnType> candidates = PointLocatorAlgos<MyMeshType>::locates(x,eps);
+ typename std::list<ConnType>::iterator eIt = candidates.begin();
for ( ; eIt != candidates.end(); ++eIt )
{
- const int i = OTT<ConnType,numPol>::ind2C( *eIt );
+ const ConnType i = OTT<ConnType,numPol>::ind2C( *eIt );
const double* coords= _mesh.getCoordinatesPtr();
const ConnType* conn=_mesh.getConnectivityPtr();
const ConnType* conn_index= _mesh.getConnectivityIndexPtr();
else
{
NormalizedCellType simlexType = cell.getDimension()==3 ? NORM_TETRA4 : NORM_TRI3;
- std::vector<int> sonNodes;
+ std::vector<mcIdType> sonNodes;
NormalizedCellType sonType;
const unsigned nbSons = cell.getNumberOfSons2( conn_elem, conn_elem_sz );
for ( unsigned s = 0; s < nbSons; ++s )
{
sonNodes.resize( cell.getNumberOfNodesConstituentTheSon2( s, conn_elem, conn_elem_sz ));
cell.fillSonCellNodalConnectivity2( s, conn_elem, conn_elem_sz, &sonNodes[0], sonType );
- std::set<int> sonNodesSet( sonNodes.begin(), sonNodes.end() );
+ std::set<mcIdType> sonNodesSet( sonNodes.begin(), sonNodes.end() );
std::set< std::set< ConnType > > checkedSonSimplex;
for ( unsigned sn = 0; sn < sonNodes.size(); ++sn )
