//
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
-#ifndef __BBTREE_TXX__
-#define __BBTREE_TXX__
+
+#pragma once
#include <vector>
#include <algorithm>
#include <limits>
#include <cmath>
+constexpr double BBTREE_DFT_EPSILON = 1e-12;
+
template <int dim, class ConnType = int>
class BBTree
{
\endcode
*/
- BBTree(const double* bbs, ConnType* elems, int level, ConnType nbelems, double epsilon=1e-12):
+ BBTree(const double* bbs, ConnType* elems, int level, ConnType nbelems, double epsilon=BBTREE_DFT_EPSILON):
_left(0), _right(0), _level(level), _bb(bbs), _terminal(false),_nbelems(nbelems),_epsilon(epsilon)
{
if (nbelems < MIN_NB_ELEMS || level> MAX_LEVEL)
return _left->size()+_right->size();
}
};
-#endif
void createBoundingBoxes(const MyMeshType& mesh, std::vector<double>& bbox);
void adjustBoundingBoxes(std::vector<double>& bbox, double adjustmentEpsAbs);
static void getElemBB(double* bb, const MyMeshType& mesh, ConnType iP, ConnType nb_nodes);
- static bool ComputeBaryCoordsOf(double startOfSeg, double endOfSeg, double pt, double& startPos, double& endPos);
+ static void ComputeBaryCoordsOf(double startOfSeg, double endOfSeg, double pt, double& startPos, double& endPos);
protected :
+ bool projectionThis(const double *coordsT, const double *coordsS, double& xs0, double& xs1, double& xt) const;
bool projectionThis(const double *coordsT, const double *coordsS, double& xs0, double& xs1, double& xt0, double& xt1) const;
bool getRealTargetCoordinates(ConnType icellT, std::vector<double>& coordsT) const;
typename MyMeshType::MyConnType getNodeIdOfTargetCellAt(ConnType icellT, ConnType nodeIdInCellT) const;
typename MyMeshType::MyConnType getNodeIdOfSourceCellAt(ConnType icellT, ConnType nodeIdInCellT) const;
double intersectSegments(const double *coordsT, const double *coordsS) const;
double intersectSegmentsInternal(const double *coordsT, const double *coordsS, double& xs0, double& xs1, double& xt0, double& xt1) const;
+ bool isPtIncludedInSeg(const double *coordsT, const double *coordsS, double& xs0, double& xs1, double& xt) const;
struct TDualSegment
{
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
// Author : Anthony Geay (CEA/DEN)
-#ifndef __CURVEINTERSECTOR_TXX__
-#define __CURVEINTERSECTOR_TXX__
+
+#pragma once
#include "CurveIntersector.hxx"
#include "InterpolationUtils.hxx"
+#include "PointLocatorAlgos.txx"
#include <limits>
/*!
* \param [in] startOfSeg - input coming from intersectSegments or intersectSegmentsInternal
- * \param [in] endOfSeg - input coming from intersectSegments or intersectSegmentsInternal. Assume that endOfSeg>startOfSeg.
+ * \param [in] endOfSeg - input coming from intersectSegments or intersectSegmentsInternal. NO Assume about sort
* \param [in] pt - position of point that the method computes the bary coords for.
*/
template<class MyMeshType, class MyMatrix>
- bool CurveIntersector<MyMeshType,MyMatrix>::ComputeBaryCoordsOf(double startOfSeg, double endOfSeg, double pt, double& startPos, double& endPos)
+ void CurveIntersector<MyMeshType,MyMatrix>::ComputeBaryCoordsOf(double startOfSeg, double endOfSeg, double pt, double& startPos, double& endPos)
{
double deno(endOfSeg-startOfSeg);
- startPos=(endOfSeg-pt)/deno;
- endPos=1.-startPos;
- return startPos>=0. && endPos>=0.;
+ startPos = (endOfSeg-pt)/deno;
+ startPos = std::max(startPos,0.); startPos = std::min(startPos,1.);
+ endPos=1.-startPos;
}
-
+
/*! Readjusts a set of bounding boxes so that they are extended
in all dimensions for avoiding missing interesting intersections
}
}
}
+
+ template<class MyMeshType, class MyMatrix>
+ bool CurveIntersector<MyMeshType,MyMatrix>::projectionThis(const double *coordsT, const double *coordsS, double& xs0, double& xs1, double& xt) const
+ {
+ enum { X=0, Y };
+ switch(SPACEDIM)
+ {
+ case 1:
+ {
+ xt = coordsT[0];
+ xs0 = coordsS[0]; xs1 = coordsS[1];
+ return true;
+ }
+ case 2:
+ {
+ const double *s0(coordsS),*s1(coordsS + 2);
+ double s01[2] = { s1[X]-s0[X], s1[Y]-s0[Y] }; // src segment direction
+ double sSize = sqrt( s01[X]*s01[X] + s01[Y]*s01[Y] ); // src segment size
+ if( sSize < this->_precision )
+ return false;
+ s01[X] /= sSize; s01[Y] /= sSize; // normalize s01
+ double t[2] = { coordsT[X]-s0[X], coordsT[Y]-s0[Y] };
+ xs0 = 0. ; xs1 = sSize; xt = s01[X]*t[X] + s01[Y]*t[Y];
+ double proj_t_on_s[2] = { s0[X]+xt*s01[X], s0[Y]+xt*s01[Y] };
+ double dist_t_s_vect[2] = { coordsT[X]-proj_t_on_s[X], coordsT[Y]-proj_t_on_s[Y] };
+ double dist_t_s = sqrt( dist_t_s_vect[X]*dist_t_s_vect[X]+dist_t_s_vect[Y]*dist_t_s_vect[Y] );
+ return dist_t_s < this->_precision;
+ }
+ default:
+ throw Exception("CurveIntersector::projectionThis : space dimension of mesh must be 1 or 2");
+ }
+ }
template<class MyMeshType, class MyMatrix>
bool CurveIntersector<MyMeshType,MyMatrix>::projectionThis(const double *coordsT, const double *coordsS,
double x1 = std::min( xt1, xs1 );
return ( x0 < x1 ) ? ( x1 - x0 ) : 0.;
}
+
+ template<class MyMeshType>
+ class DummyMyMeshType1D
+ {
+ public:
+ static const int MY_SPACEDIM=1;
+ static const int MY_MESHDIM=8;
+ typedef mcIdType MyConnType;
+ static const INTERP_KERNEL::NumberingPolicy My_numPol=MyMeshType::My_numPol;
+ // begin
+ // useless, but for windows compilation ...
+ const double *getCoordinatesPtr() const { return nullptr; }
+ const MyConnType *getConnectivityPtr() const { return nullptr; }
+ const MyConnType *getConnectivityIndexPtr() const { return nullptr; }
+ INTERP_KERNEL::NormalizedCellType getTypeOfElement(MyConnType) const { return (INTERP_KERNEL::NormalizedCellType)0; }
+ // end
+ };
+
+ /*!
+ * This method determines if a target point ( \a coordsT ) is in source seg2 contained in \a coordsS. To do so _precision attribute is used.
+ * If target point is in, \a xs0, \a xs1 and \a xt are set to 1D referential for a further barycentric computation.
+ * This method deals with SPACEDIM == 2 (see projectionThis).
+ */
+ template<class MyMeshType, class MyMatrix>
+ bool CurveIntersector<MyMeshType,MyMatrix>::isPtIncludedInSeg(const double *coordsT, const double *coordsS, double& xs0, double& xs1, double& xt) const
+ {
+ if(!projectionThis(coordsT,coordsS,xs0,xs1,xt))
+ return false;
+ constexpr ConnType TAB[2]={0,1};
+ const double coordsS_1D[2]={xs0,xs1};
+ const double *coordsT_1D(&xt);
+ return PointLocatorAlgos<DummyMyMeshType1D<MyMeshType>>::isElementContainsPoint(coordsT_1D,NORM_SEG2,coordsS_1D,TAB,2,this->_precision);
+ }
/*!
* \brief Return length of intersection of two segments
double xs0,xs1,xt0,xt1;
return intersectSegmentsInternal(coordsT,coordsS,xs0,xs1,xt0,xt1);
}
-
}
-
-#endif
//
// Author : Anthony Geay (EDF R&D)
-#ifndef __CURVEINTERSECTORP1P1PL_TXX__
-#define __CURVEINTERSECTORP1P1PL_TXX__
+#pragma once
#include "CurveIntersectorP1P1PL.hxx"
#include "CurveIntersector.txx"
template<class MyMeshType, class MyMatrix>
typename MyMeshType::MyConnType CurveIntersectorP1P1PL<MyMeshType,MyMatrix>::getNumberOfRowsOfResMatrix() const
{
- return CurveIntersector<MyMeshType,MyMatrix>::_meshT.getNumberOfNodes();
+ return this->_meshT.getNumberOfNodes();
}
template<class MyMeshType, class MyMatrix>
typename MyMeshType::MyConnType CurveIntersectorP1P1PL<MyMeshType,MyMatrix>::getNumberOfColsOfResMatrix() const
{
- return CurveIntersector<MyMeshType,MyMatrix>::_meshS.getNumberOfNodes();
+ return this->_meshS.getNumberOfNodes();
}
template<class MyMeshType, class MyMatrix>
void CurveIntersectorP1P1PL<MyMeshType,MyMatrix>::AppendValueInMatrix(MyMatrix& res, ConnType nodeIdT, ConnType nodeIdS0, double val0, ConnType nodeIdS1, double val1)
{
typename MyMatrix::value_type& resRow(res[nodeIdT]);
- AppendValueInMatrix2(resRow,nodeIdS0,val0);
- AppendValueInMatrix2(resRow,nodeIdS1,val1);
+ if(std::abs(val0)>0.)
+ AppendValueInMatrix2(resRow,nodeIdS0,val0);
+ if(std::abs(val1)>0.)
+ AppendValueInMatrix2(resRow,nodeIdS1,val1);
}
template<class MyMeshType, class MyMatrix>
void CurveIntersectorP1P1PL<MyMeshType,MyMatrix>::intersectCells(ConnType icellT, const std::vector<ConnType>& icellsS, MyMatrix& res)
{
std::vector<double> coordsT;
- if(CurveIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates(icellT,coordsT))
+ if(this->getRealTargetCoordinates(icellT,coordsT))
throw INTERP_KERNEL::Exception("Invalid target cell detected for meshdim==1. Only SEG2 supported !");
- assert(coordsT.size()/SPACEDIM==2);
+ std::size_t nbOfNodesT(coordsT.size()/SPACEDIM);
for(typename std::vector<ConnType>::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++)
{
std::vector<double> coordsS;
- if(CurveIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinates(*iter,coordsS))
+ if(this->getRealSourceCoordinates(*iter,coordsS))
throw INTERP_KERNEL::Exception("Invalid source cell detected for meshdim==1. Only SEG2 supported !");
assert(coordsS.size()/SPACEDIM==2);
- double xs0,xs1,xt0,xt1;
- double lgth(CurveIntersector<MyMeshType,MyMatrix>::intersectSegmentsInternal(&coordsT[0],&coordsS[0],xs0,xs1,xt0,xt1));
- ConnType nodeIdS0(CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfSourceCellAt(*iter,0));
- ConnType nodeIdS1(CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfSourceCellAt(*iter,1));
- if(lgth>0.)
+ ConnType nodeIdS0(this->getNodeIdOfSourceCellAt(*iter,0));
+ ConnType nodeIdS1(this->getNodeIdOfSourceCellAt(*iter,1));
+ for(std::size_t nodeIdT = 0 ; nodeIdT<nbOfNodesT ; ++nodeIdT)
{
- double a,b;
- // for first
- ConnType nodeIdT0(CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfTargetCellAt(icellT,0));
- if(CurveIntersector<MyMeshType,MyMatrix>::ComputeBaryCoordsOf(xs0,xs1,xt0,a,b))
+ ConnType nodeIdT0(this->getNodeIdOfTargetCellAt(icellT,nodeIdT));
+ double xs0,xs1,xt;
+ if( this->isPtIncludedInSeg(coordsT.data()+nodeIdT*SPACEDIM,coordsS.data(),xs0,xs1,xt) )
{
- a*=lgth; b*=lgth;
+ double a,b;
+ this->ComputeBaryCoordsOf(xs0,xs1,xt,a,b);
AppendValueInMatrix(res,nodeIdT0,nodeIdS0,a,nodeIdS1,b);
}
- //
- ConnType nodeIdT1(CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfTargetCellAt(icellT,1));
- typename MyMatrix::value_type& resRow1=res[nodeIdT1];
- if(CurveIntersector<MyMeshType,MyMatrix>::ComputeBaryCoordsOf(xs0,xs1,xt1,a,b))
- {
- a*=lgth; b*=lgth;
- AppendValueInMatrix(res,nodeIdT1,nodeIdS0,a,nodeIdS1,b);
- }
}
}
}
}
-
-#endif
#include "Interpolation.hxx"\r
#include "InterpolationOptions.hxx"\r
\r
+#include "BBTree.txx"\r
+\r
+#include <functional>\r
+\r
namespace INTERP_KERNEL\r
{\r
template<class RealCurve>\r
\r
// Main function to interpolate\r
template<class MyMeshType, class MatrixType>\r
+ typename MyMeshType::MyConnType interpolateMeshesInternal(const MyMeshType& meshS, const MyMeshType& meshT,\r
+ MatrixType& result, const std::string& method,\r
+ std::function< void(const BBTree< MyMeshType::MY_SPACEDIM , typename MyMeshType::MyConnType>&, const double*, std::vector<typename MyMeshType::MyConnType>&) > bbtreeMethod);\r
+ template<class MyMeshType, class MatrixType>\r
typename MyMeshType::MyConnType interpolateMeshes(const MyMeshType& meshS, const MyMeshType& meshT,\r
- MatrixType& result, const std::string& method);\r
- \r
+ MatrixType& result, const std::string& method) { return this->interpolateMeshesInternal(meshS,meshT,result,method,[](const BBTree< MyMeshType::MY_SPACEDIM ,\r
+ typename MyMeshType::MyConnType>& bbtree, const double *bb,\r
+ std::vector<typename MyMeshType::MyConnType>& intersecting_elems)\r
+ { bbtree.getIntersectingElems(bb, intersecting_elems); }); }\r
+ template<class MyMeshType, class MatrixType>\r
+ typename MyMeshType::MyConnType interpolateMeshes0D(const MyMeshType& meshS, const MyMeshType& meshT,\r
+ MatrixType& result, const std::string& method) { return this->interpolateMeshesInternal(meshS,meshT,result,method,[](const BBTree< MyMeshType::MY_SPACEDIM ,\r
+ typename MyMeshType::MyConnType>& bbtree, const double *bb,\r
+ std::vector<typename MyMeshType::MyConnType>& intersecting_elems)\r
+ {\r
+ double TMP[MyMeshType::MY_SPACEDIM];\r
+ for(int i=0;i<MyMeshType::MY_SPACEDIM;++i)\r
+ TMP[i] = bb[2*i];\r
+ bbtree.getElementsAroundPoint(TMP, intersecting_elems);\r
+ }\r
+ ); }\r
};\r
}\r
\r
#include "CurveIntersectorP0P1.txx"
#include "CurveIntersectorP1P1.txx"
#include "CurveIntersectorP1P1PL.txx"
-#include "BBTree.txx"
#include <time.h>
+#include <memory>
namespace INTERP_KERNEL
{
*/
template<class RealCurve>
template<class MyMeshType, class MatrixType>
- typename MyMeshType::MyConnType InterpolationCurve<RealCurve>::interpolateMeshes (const MyMeshType& myMeshS,
- const MyMeshType& myMeshT,
- MatrixType& result,
- const std::string& method)
+ typename MyMeshType::MyConnType InterpolationCurve<RealCurve>::interpolateMeshesInternal (const MyMeshType& myMeshS,
+ const MyMeshType& myMeshT,
+ MatrixType& result,
+ const std::string& method,
+ std::function< void(const BBTree< MyMeshType::MY_SPACEDIM ,
+ typename MyMeshType::MyConnType>&, const double*,
+ std::vector<typename MyMeshType::MyConnType>&) > bbtreeMethod
+ )
{
static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
typedef typename MyMeshType::MyConnType ConnType;
ConnType nbMailleS = myMeshS.getNumberOfElements();
ConnType nbMailleT = myMeshT.getNumberOfElements();
- CurveIntersector<MyMeshType,MatrixType>* intersector=0;
+ std::unique_ptr< CurveIntersector<MyMeshType,MatrixType> > intersector;
if(method=="P0P0")
{
switch (InterpolationOptions::getIntersectionType())
{
case Triangulation:
{
- intersector = new CurveIntersectorP0P0<MyMeshType,MatrixType>(myMeshT, myMeshS,
- InterpolationOptions::getPrecision(),
- InterpolationOptions::getBoundingBoxAdjustmentAbs(),
- InterpolationOptions::getMedianPlane(),
- InterpolationOptions::getPrintLevel());
+ intersector.reset( new CurveIntersectorP0P0<MyMeshType,MatrixType>(myMeshT, myMeshS,
+ InterpolationOptions::getPrecision(),
+ InterpolationOptions::getBoundingBoxAdjustmentAbs(),
+ InterpolationOptions::getMedianPlane(),
+ InterpolationOptions::getPrintLevel()) );
break;
}
default:
{
case Triangulation:
{
- intersector = new CurveIntersectorP0P1<MyMeshType,MatrixType>(myMeshT, myMeshS,
- InterpolationOptions::getPrecision(),
- InterpolationOptions::getBoundingBoxAdjustmentAbs(),
- InterpolationOptions::getMedianPlane(),
- InterpolationOptions::getPrintLevel());
+ intersector.reset( new CurveIntersectorP0P1<MyMeshType,MatrixType>(myMeshT, myMeshS,
+ InterpolationOptions::getPrecision(),
+ InterpolationOptions::getBoundingBoxAdjustmentAbs(),
+ InterpolationOptions::getMedianPlane(),
+ InterpolationOptions::getPrintLevel()) );
break;
}
default:
{
case Triangulation:
{
- intersector = new CurveIntersectorP1P0<MyMeshType,MatrixType>(myMeshT, myMeshS,
- InterpolationOptions::getPrecision(),
- InterpolationOptions::getBoundingBoxAdjustmentAbs(),
- InterpolationOptions::getMedianPlane(),
- InterpolationOptions::getPrintLevel());
+ intersector.reset( new CurveIntersectorP1P0<MyMeshType,MatrixType>(myMeshT, myMeshS,
+ InterpolationOptions::getPrecision(),
+ InterpolationOptions::getBoundingBoxAdjustmentAbs(),
+ InterpolationOptions::getMedianPlane(),
+ InterpolationOptions::getPrintLevel()) );
break;
}
- default:
- throw INTERP_KERNEL::Exception("For P1P0 in 1D or 2D curve only Triangulation supported for the moment !");
+ default:
+ throw INTERP_KERNEL::Exception("For P1P0 in 1D or 2D curve only Triangulation supported for the moment !");
}
}
else if(method=="P1P1")
switch (InterpolationOptions::getIntersectionType())
{
case Triangulation:
- intersector = new CurveIntersectorP1P1<MyMeshType,MatrixType>
- (myMeshT, myMeshS,
- InterpolationOptions::getPrecision(),
- InterpolationOptions::getBoundingBoxAdjustmentAbs(),
- InterpolationOptions::getMedianPlane(),
- InterpolationOptions::getPrintLevel());
+ intersector.reset( new CurveIntersectorP1P1<MyMeshType,MatrixType>
+ (myMeshT, myMeshS,
+ InterpolationOptions::getPrecision(),
+ InterpolationOptions::getBoundingBoxAdjustmentAbs(),
+ InterpolationOptions::getMedianPlane(),
+ InterpolationOptions::getPrintLevel()) );
break;
case PointLocator:
- intersector = new CurveIntersectorP1P1PL<MyMeshType,MatrixType>
- (myMeshT, myMeshS,
- InterpolationOptions::getPrecision(),
- InterpolationOptions::getBoundingBoxAdjustmentAbs(),
- InterpolationOptions::getMedianPlane(),
- InterpolationOptions::getPrintLevel());
+ intersector.reset( new CurveIntersectorP1P1PL<MyMeshType,MatrixType>
+ (myMeshT, myMeshS,
+ InterpolationOptions::getPrecision(),
+ InterpolationOptions::getBoundingBoxAdjustmentAbs(),
+ InterpolationOptions::getMedianPlane(),
+ InterpolationOptions::getPrintLevel()) );
break;
default:
throw INTERP_KERNEL::Exception("For P1P1 in 1D or 2D curve only Triangulation and PointLocator supported !");
std::vector<ConnType> intersecting_elems;
double bb[2*SPACEDIM];
intersector->getElemBB(bb,myMeshT,OTT<ConnType,numPol>::indFC(iT),nb_nodesT);
- my_tree.getIntersectingElems(bb, intersecting_elems);
+ bbtreeMethod(my_tree,bb,intersecting_elems);
intersector->intersectCells(iT,intersecting_elems,result);
counter += intersecting_elems.size();
}
- ConnType ret = intersector->getNumberOfColsOfResMatrix();
- delete intersector;
if (InterpolationOptions::getPrintLevel() >= 1)
{
std::cout << "Number of computed intersections = " << counter << std::endl;
std::cout << "Global time= " << global_end - global_start << std::endl;
}
- return ret;
+ return intersector->getNumberOfColsOfResMatrix();
}
}
void PlanarIntersectorP1P1PL<MyMeshType,MyMatrix>::intersectCells(ConnType icellT, const std::vector<ConnType>& icellsS, MyMatrix& res)
{
std::vector<double> CoordsT;
- PlanarIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates(OTT<ConnType,numPol>::indFC(icellT),CoordsT);
+ this->getRealTargetCoordinates(OTT<ConnType,numPol>::indFC(icellT),CoordsT);
ConnType nbOfNodesT=ToConnType(CoordsT.size())/SPACEDIM;
for(typename std::vector<ConnType>::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++)
{
- NormalizedCellType tS=PlanarIntersector<MyMeshType,MyMatrix>::_meshS.getTypeOfElement(OTT<ConnType,numPol>::indFC(*iter));
+ NormalizedCellType tS=this->_meshS.getTypeOfElement(OTT<ConnType,numPol>::indFC(*iter));
if(tS!=NORM_TRI3)
throw INTERP_KERNEL::Exception("Invalid source cell detected for meshdim==2. Only TRI3 supported !");
std::vector<double> CoordsS;
- PlanarIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinates(OTT<ConnType,numPol>::indFC(*iter),CoordsS);
+ this->getRealSourceCoordinates(OTT<ConnType,numPol>::indFC(*iter),CoordsS);
std::vector<double> CoordsTTmp(CoordsT);
if(SPACEDIM==3)
- PlanarIntersector<MyMeshType,MyMatrix>::projectionThis(&CoordsS[0],&CoordsTTmp[0],ToConnType(CoordsS.size())/SPACEDIM,nbOfNodesT);
- const ConnType *startOfCellNodeConnT=PlanarIntersector<MyMeshType,MyMatrix>::_connectT+OTT<ConnType,numPol>::conn2C(PlanarIntersector<MyMeshType,MyMatrix>::_connIndexT[icellT]);
+ this->projectionThis(CoordsS.data(),CoordsTTmp.data(),ToConnType(CoordsS.size())/SPACEDIM,nbOfNodesT);
+ const ConnType *startOfCellNodeConnT=this->_connectT+OTT<ConnType,numPol>::conn2C(this->_connIndexT[icellT]);
for(int nodeIdT=0;nodeIdT<nbOfNodesT;nodeIdT++)
{
typename MyMatrix::value_type& resRow=res[OTT<ConnType,numPol>::ind2C(startOfCellNodeConnT[nodeIdT])];
- if( PointLocatorAlgos<MyMeshType>::isElementContainsPointAlg2D(&CoordsTTmp[nodeIdT*SPACEDIM],&CoordsS[0],3,PlanarIntersector<MyMeshType,MyMatrix>::_precision) )
+ if( PointLocatorAlgos<MyMeshType>::isElementContainsPointAlg2D(CoordsTTmp.data()+nodeIdT*SPACEDIM,CoordsS.data(),3,this->_precision) )
{
double resLoc[3];
barycentric_coords<SPACEDIM>(&CoordsS[0],&CoordsTTmp[nodeIdT*SPACEDIM],resLoc);
- const ConnType *startOfCellNodeConnS=PlanarIntersector<MyMeshType,MyMatrix>::_connectS+OTT<ConnType,numPol>::conn2C(PlanarIntersector<MyMeshType,MyMatrix>::_connIndexS[*iter]);
+ const ConnType *startOfCellNodeConnS=this->_connectS+OTT<ConnType,numPol>::conn2C(this->_connIndexS[*iter]);
for(int nodeIdS=0;nodeIdS<3;nodeIdS++)
{
- if(fabs(resLoc[nodeIdS])>PlanarIntersector<MyMeshType,MyMatrix>::_precision)
+ if(fabs(resLoc[nodeIdS])>this->_precision)
{
ConnType curNodeSInCmode=OTT<ConnType,numPol>::coo2C(startOfCellNodeConnS[nodeIdS]);
typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT<ConnType,numPol>::indFC(curNodeSInCmode));
template<class MyMeshType, class MyMatrix>
typename MyMeshType::MyConnType PlanarIntersectorP1P1PL<MyMeshType,MyMatrix>::getNumberOfRowsOfResMatrix() const
{
- return PlanarIntersector<MyMeshType,MyMatrix>::_meshT.getNumberOfNodes();
+ return this->_meshT.getNumberOfNodes();
}
template<class MyMeshType, class MyMatrix>
typename MyMeshType::MyConnType PlanarIntersectorP1P1PL<MyMeshType,MyMatrix>::getNumberOfColsOfResMatrix() const
{
- return PlanarIntersector<MyMeshType,MyMatrix>::_meshS.getNumberOfNodes();
+ return this->_meshS.getNumberOfNodes();
}
}
INTERP_KERNEL::Interpolation1D interpolation(*this);
nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
}
+ else if(srcMeshDim==1 && trgMeshDim==0 && srcSpaceDim==1 )
+ {
+ if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
+ throw INTERP_KERNEL::Exception("Invalid interpolation requested between 1D and 0D into 1D space ! Select PointLocator as intersection type !");
+ MEDCouplingNormalizedUnstructuredMesh<1,1> source_mesh_wrapper(src_mesh);
+ MEDCouplingNormalizedUnstructuredMesh<1,1> target_mesh_wrapper(target_mesh);
+ INTERP_KERNEL::Interpolation1D interpolation(*this);
+ nbCols=interpolation.interpolateMeshes0D(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
+ }
+ else if(srcMeshDim==1 && trgMeshDim==0 && srcSpaceDim==2 )
+ {
+ if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
+ throw INTERP_KERNEL::Exception("Invalid interpolation requested between 1D and 0D into 2D space ! Select PointLocator as intersection type !");
+ MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
+ MEDCouplingNormalizedUnstructuredMesh<2,1> target_mesh_wrapper(target_mesh);
+ INTERP_KERNEL::Interpolation1D interpolation(*this);
+ nbCols=interpolation.interpolateMeshes0D(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
+ }
else if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==2)
{
MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
INTERP_KERNEL::Interpolation1D interpolation1D(*this);
if(interpolation1D.getIntersectionType()==INTERP_KERNEL::Geometric2D)// For intersection type of 1D, Geometric2D do not deal with it ! -> make interpolation1D not inherite from this
interpolation1D.setIntersectionType(INTERP_KERNEL::Triangulation);//
- mcIdType nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC);
+ mcIdType nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC);
s1D->decrRef();
t1D->decrRef();
buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
m=rem.getCrudeCSRMatrix()
row=array([1,1,2,2,3,3])
col=array([0,1,1,2,5,6])
- data=array([0.9,0.1,0.3,0.7,0.5,0.5])
+ data=array([1.8,0.2,0.6,1.4,1.0,1.0])
mExp2=csr_matrix((data,(row,col)),shape=(5,11))
diff=abs(m-mExp2)
self.assertAlmostEqual(diff.sum(),0.,14)
self.checkMatrix(rem.getCrudeMatrix(),[{0:1.0}],src.getNumberOfNodes(),1e-12)
pass
+ def test1D0DPointLocator(self):
+ """
+ For pointlocator fans, Remapper support following intersection
+ IntersectionType == PointLocator
+ - source == 1D
+ - target == 0D
+ """
+ # P1P1 - 0
+ src = MEDCouplingUMesh("src",1)
+ src.allocateCells()
+ src.setCoords( DataArrayDouble([0,1]) )
+ src.insertNextCell(NORM_SEG2,[0,1])
+ trg = MEDCouplingUMesh.Build0DMeshFromCoords( DataArrayDouble([0.4]) )
+ rem=MEDCouplingRemapper()
+ rem.setIntersectionType(PointLocator)
+ rem.prepare(src,trg,"P1P1")
+ self.checkMatrix(rem.getCrudeMatrix(),[{0:0.6,1:0.4}],src.getNumberOfNodes(),1e-12)
+ # P1P1 - 1
+ src = MEDCouplingUMesh("src",1)
+ src.allocateCells()
+ src.setCoords( DataArrayDouble([0,1]) )
+ src.insertNextCell(NORM_SEG2,[1,0]) # permutation
+ trg = MEDCouplingUMesh.Build0DMeshFromCoords( DataArrayDouble([0.4]) )
+ rem=MEDCouplingRemapper()
+ rem.setIntersectionType(PointLocator)
+ rem.prepare(src,trg,"P1P1")
+ self.checkMatrix(rem.getCrudeMatrix(),[{0:0.6,1:0.4}],src.getNumberOfNodes(),1e-12)
+ # P1P1 - 2
+ src = MEDCouplingUMesh("src",1)
+ src.allocateCells()
+ src.setCoords( DataArrayDouble([1,0]) )
+ src.insertNextCell(NORM_SEG2,[0,1])
+ trg = MEDCouplingUMesh.Build0DMeshFromCoords( DataArrayDouble([0.4]) )
+ rem=MEDCouplingRemapper()
+ rem.setIntersectionType(PointLocator)
+ rem.prepare(src,trg,"P1P1")
+ self.checkMatrix(rem.getCrudeMatrix(),[{0:0.4,1:0.6}],src.getNumberOfNodes(),1e-12)
+ # P1P1 - 3 - 2DCurve
+ src = MEDCouplingUMesh("src",1)
+ src.allocateCells()
+ src.setCoords( DataArrayDouble([0,1]) )
+ src.insertNextCell(NORM_SEG2,[0,1])
+ trg = MEDCouplingUMesh.Build0DMeshFromCoords( DataArrayDouble([0.4]) )
+ src.changeSpaceDimension(2) ; trg.changeSpaceDimension(2)
+ src.rotate([-1.,-1.],1.2)
+ trg.rotate([-1.,-1.],1.2)
+ rem=MEDCouplingRemapper()
+ rem.setIntersectionType(PointLocator)
+ rem.prepare(src,trg,"P1P1")
+ self.checkMatrix(rem.getCrudeMatrix(),[{0:0.6,1:0.4}],src.getNumberOfNodes(),1e-12)
+ # P1P1 - 4
+ src = MEDCouplingUMesh("src",1)
+ src.allocateCells()
+ src.setCoords( DataArrayDouble([1.1,7.6,2.3,5.4]) )
+ src.insertNextCell(NORM_SEG2,[0,2])
+ src.insertNextCell(NORM_SEG2,[2,3])
+ src.insertNextCell(NORM_SEG2,[3,1])
+ for eps in [0,1e-13,-1e-13]:
+ trg = MEDCouplingUMesh.Build0DMeshFromCoords( DataArrayDouble([0.4,2.3+eps,4.,7.]) )
+ rem=MEDCouplingRemapper()
+ rem.setIntersectionType(PointLocator)
+ rem.prepare(src,trg,"P1P1")
+ rem.nullifiedTinyCoeffInCrudeMatrixAbs(1e-12)
+ self.checkMatrix(rem.getCrudeMatrix(),[{}, {2: 2.0}, {2: 0.4516129032258065, 3: 0.5483870967741935}, {1: 0.7272727272727273, 3: 0.27272727272727265}],src.getNumberOfNodes(),1e-12)
+ # P1P1 - 5 - descending order of coords in source mesh
+ src = MEDCouplingUMesh("src",1)
+ src.allocateCells()
+ src.setCoords( DataArrayDouble([3.,1.]) )
+ src.insertNextCell(NORM_SEG2,[0,1])
+ trg = MEDCouplingUMesh.Build0DMeshFromCoords( DataArrayDouble([2.3]) )
+ rem=MEDCouplingRemapper()
+ rem.setIntersectionType(PointLocator)
+ rem.prepare(src,trg,"P1P1")
+ self.checkMatrix(rem.getCrudeMatrix(),[{0:0.65,1:0.35}],src.getNumberOfNodes(),1e-12)
+ pass
+
def checkMatrix(self,mat1,mat2,nbCols,eps):
self.assertEqual(len(mat1),len(mat2))
for i in range(len(mat1)):
- self.assertTrue(max(mat2[i].keys())<nbCols)
- self.assertTrue(max(mat1[i].keys())<nbCols)
- self.assertTrue(min(mat2[i].keys())>=0)
- self.assertTrue(min(mat1[i].keys())>=0)
+ if len(mat2[i].keys())>0:
+ self.assertTrue(max(mat2[i].keys())<nbCols)
+ if len(mat1[i].keys())>0:
+ self.assertTrue(max(mat1[i].keys())<nbCols)
+ if len(mat2[i].keys())>0:
+ self.assertTrue(min(mat2[i].keys())>=0)
+ if len(mat1[i].keys())>0:
+ self.assertTrue(min(mat1[i].keys())>=0)
s1=set(mat1[i].keys()) ; s2=set(mat2[i].keys())
for elt in s1.intersection(s2):
self.assertTrue(abs(mat1[i][elt]-mat2[i][elt])<eps)
