//
// 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
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);
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>
}
}
}
+
+ 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.;
}
+
+ /*!
+ * 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};
+ return PointLocatorAlgos<MyMeshType>::isElementContainsPoint(coordsT,NORM_SEG2,coordsS,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>::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;
- 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);
+ double a,b;
+ if(this->ComputeBaryCoordsOf(xs0,xs1,xt,a,b))
+ {
+ AppendValueInMatrix(res,nodeIdT0,nodeIdS0,a,nodeIdS1,b);
+ }
}
}
}
}
}
-
-#endif
#include "BBTree.txx"
#include <time.h>
+#include <memory>
namespace INTERP_KERNEL
{
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 !");
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();
}
}
{
MEDCouplingNormalizedUnstructuredMesh<1,1> source_mesh_wrapper(src_mesh);
MEDCouplingNormalizedUnstructuredMesh<1,1> target_mesh_wrapper(target_mesh);
+ //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.interpolateMeshes(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.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
+ }
else if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==2)
{
MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
throw INTERP_KERNEL::Exception("Invalid interpolation requested between 1D and 0D into 3D space ! Select PointLocator as intersection type !");
MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
- INTERP_KERNEL::Interpolation1D0D interpolation(*this);
- nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
+ //INTERP_KERNEL::Interpolation1D0D interpolation(*this);
+ //nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
}
else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
{
MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
std::vector<std::map<mcIdType,double> > matrix1D;
- INTERP_KERNEL::Interpolation1D interpolation1D(*this);
+ /*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);*/
+ mcIdType nbCols1D=-1;
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)
ref=float(m.getMeasureField(True).getArray())
self.assertTrue(abs(res-ref)/ref<1e-12)
pass
-
+
def test3D0DPointLocator(self):
"""
For pointlocator fans, Remapper support following intersection
self.checkMatrix(rem.getCrudeMatrix(),[{0:1.0}],src.getNumberOfNodes(),1e-12)
pass
+class MEDCouplingBasicsTest(unittest.TestCase):
+ 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])
+ trg = MEDCouplingUMesh.Build0DMeshFromCoords( DataArrayDouble([0.4,2.3,4.,7.]) )
+ rem=MEDCouplingRemapper()
+ rem.setIntersectionType(PointLocator)
+ rem.prepare(src,trg,"P1P1")
+ print(rem.getCrudeMatrix())
+ #self.checkMatrix(rem.getCrudeMatrix(),[{0:0.4,1:0.6}],src.getNumberOfNodes(),1e-12)
+ pass
+
def checkMatrix(self,mat1,mat2,nbCols,eps):
self.assertEqual(len(mat1),len(mat2))
for i in range(len(mat1)):
