From: Anthony Geay Date: Tue, 31 Dec 2019 14:54:57 +0000 (+0100) Subject: En route X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=941de44d86c9d3f7fc87763d72d94f224ff6fa5d;p=tools%2Fmedcoupling.git En route --- diff --git a/src/INTERP_KERNEL/BBTree.txx b/src/INTERP_KERNEL/BBTree.txx index c312a646b..582ec5e50 100644 --- a/src/INTERP_KERNEL/BBTree.txx +++ b/src/INTERP_KERNEL/BBTree.txx @@ -16,8 +16,8 @@ // // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -#ifndef __BBTREE_TXX__ -#define __BBTREE_TXX__ + +#pragma once #include #include @@ -26,6 +26,8 @@ #include #include +constexpr double BBTREE_DFT_EPSILON = 1e-12; + template class BBTree { @@ -66,7 +68,7 @@ public: \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) @@ -278,4 +280,3 @@ public: return _left->size()+_right->size(); } }; -#endif diff --git a/src/INTERP_KERNEL/CurveIntersector.hxx b/src/INTERP_KERNEL/CurveIntersector.hxx index 339343820..9fe10e7ce 100644 --- a/src/INTERP_KERNEL/CurveIntersector.hxx +++ b/src/INTERP_KERNEL/CurveIntersector.hxx @@ -43,6 +43,7 @@ namespace INTERP_KERNEL 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& coordsT) const; typename MyMeshType::MyConnType getNodeIdOfTargetCellAt(ConnType icellT, ConnType nodeIdInCellT) const; @@ -50,6 +51,7 @@ namespace INTERP_KERNEL 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 { diff --git a/src/INTERP_KERNEL/CurveIntersector.txx b/src/INTERP_KERNEL/CurveIntersector.txx index 50df7a994..329c6030c 100644 --- a/src/INTERP_KERNEL/CurveIntersector.txx +++ b/src/INTERP_KERNEL/CurveIntersector.txx @@ -17,11 +17,12 @@ // 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 @@ -314,6 +315,38 @@ namespace INTERP_KERNEL } } } + + template + bool CurveIntersector::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 bool CurveIntersector::projectionThis(const double *coordsT, const double *coordsS, @@ -414,6 +447,20 @@ namespace INTERP_KERNEL 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 + bool CurveIntersector::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::isElementContainsPoint(coordsT,NORM_SEG2,coordsS,TAB,2,this->_precision); + } /*! * \brief Return length of intersection of two segments @@ -424,7 +471,4 @@ namespace INTERP_KERNEL double xs0,xs1,xt0,xt1; return intersectSegmentsInternal(coordsT,coordsS,xs0,xs1,xt0,xt1); } - } - -#endif diff --git a/src/INTERP_KERNEL/CurveIntersectorP1P1PL.txx b/src/INTERP_KERNEL/CurveIntersectorP1P1PL.txx index e9bd3e5c2..6b6e3d2df 100644 --- a/src/INTERP_KERNEL/CurveIntersectorP1P1PL.txx +++ b/src/INTERP_KERNEL/CurveIntersectorP1P1PL.txx @@ -18,8 +18,7 @@ // // Author : Anthony Geay (EDF R&D) -#ifndef __CURVEINTERSECTORP1P1PL_TXX__ -#define __CURVEINTERSECTORP1P1PL_TXX__ +#pragma once #include "CurveIntersectorP1P1PL.hxx" #include "CurveIntersector.txx" @@ -36,13 +35,13 @@ namespace INTERP_KERNEL template typename MyMeshType::MyConnType CurveIntersectorP1P1PL::getNumberOfRowsOfResMatrix() const { - return CurveIntersector::_meshT.getNumberOfNodes(); + return this->_meshT.getNumberOfNodes(); } template typename MyMeshType::MyConnType CurveIntersectorP1P1PL::getNumberOfColsOfResMatrix() const { - return CurveIntersector::_meshS.getNumberOfNodes(); + return this->_meshS.getNumberOfNodes(); } template @@ -73,40 +72,30 @@ namespace INTERP_KERNEL void CurveIntersectorP1P1PL::intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res) { std::vector coordsT; - if(CurveIntersector::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::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++) { std::vector coordsS; - if(CurveIntersector::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::intersectSegmentsInternal(&coordsT[0],&coordsS[0],xs0,xs1,xt0,xt1)); - ConnType nodeIdS0(CurveIntersector::getNodeIdOfSourceCellAt(*iter,0)); - ConnType nodeIdS1(CurveIntersector::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::getNodeIdOfTargetCellAt(icellT,0)); - if(CurveIntersector::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::getNodeIdOfTargetCellAt(icellT,1)); - typename MyMatrix::value_type& resRow1=res[nodeIdT1]; - if(CurveIntersector::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 diff --git a/src/INTERP_KERNEL/InterpolationCurve.txx b/src/INTERP_KERNEL/InterpolationCurve.txx index 972aa5f91..08dfaefe3 100755 --- a/src/INTERP_KERNEL/InterpolationCurve.txx +++ b/src/INTERP_KERNEL/InterpolationCurve.txx @@ -30,6 +30,7 @@ #include "BBTree.txx" #include +#include namespace INTERP_KERNEL { @@ -88,18 +89,18 @@ namespace INTERP_KERNEL ConnType nbMailleS = myMeshS.getNumberOfElements(); ConnType nbMailleT = myMeshT.getNumberOfElements(); - CurveIntersector* intersector=0; + std::unique_ptr< CurveIntersector > intersector; if(method=="P0P0") { switch (InterpolationOptions::getIntersectionType()) { case Triangulation: { - intersector = new CurveIntersectorP0P0(myMeshT, myMeshS, - InterpolationOptions::getPrecision(), - InterpolationOptions::getBoundingBoxAdjustmentAbs(), - InterpolationOptions::getMedianPlane(), - InterpolationOptions::getPrintLevel()); + intersector.reset( new CurveIntersectorP0P0(myMeshT, myMeshS, + InterpolationOptions::getPrecision(), + InterpolationOptions::getBoundingBoxAdjustmentAbs(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getPrintLevel()) ); break; } default: @@ -112,11 +113,11 @@ namespace INTERP_KERNEL { case Triangulation: { - intersector = new CurveIntersectorP0P1(myMeshT, myMeshS, - InterpolationOptions::getPrecision(), - InterpolationOptions::getBoundingBoxAdjustmentAbs(), - InterpolationOptions::getMedianPlane(), - InterpolationOptions::getPrintLevel()); + intersector.reset( new CurveIntersectorP0P1(myMeshT, myMeshS, + InterpolationOptions::getPrecision(), + InterpolationOptions::getBoundingBoxAdjustmentAbs(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getPrintLevel()) ); break; } default: @@ -129,15 +130,15 @@ namespace INTERP_KERNEL { case Triangulation: { - intersector = new CurveIntersectorP1P0(myMeshT, myMeshS, - InterpolationOptions::getPrecision(), - InterpolationOptions::getBoundingBoxAdjustmentAbs(), - InterpolationOptions::getMedianPlane(), - InterpolationOptions::getPrintLevel()); + intersector.reset( new CurveIntersectorP1P0(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") @@ -145,20 +146,20 @@ namespace INTERP_KERNEL switch (InterpolationOptions::getIntersectionType()) { case Triangulation: - intersector = new CurveIntersectorP1P1 - (myMeshT, myMeshS, - InterpolationOptions::getPrecision(), - InterpolationOptions::getBoundingBoxAdjustmentAbs(), - InterpolationOptions::getMedianPlane(), - InterpolationOptions::getPrintLevel()); + intersector.reset( new CurveIntersectorP1P1 + (myMeshT, myMeshS, + InterpolationOptions::getPrecision(), + InterpolationOptions::getBoundingBoxAdjustmentAbs(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getPrintLevel()) ); break; case PointLocator: - intersector = new CurveIntersectorP1P1PL - (myMeshT, myMeshS, - InterpolationOptions::getPrecision(), - InterpolationOptions::getBoundingBoxAdjustmentAbs(), - InterpolationOptions::getMedianPlane(), - InterpolationOptions::getPrintLevel()); + intersector.reset( new CurveIntersectorP1P1PL + (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 !"); @@ -197,8 +198,6 @@ namespace INTERP_KERNEL intersector->intersectCells(iT,intersecting_elems,result); counter += intersecting_elems.size(); } - ConnType ret = intersector->getNumberOfColsOfResMatrix(); - delete intersector; if (InterpolationOptions::getPrintLevel() >= 1) { @@ -209,7 +208,7 @@ namespace INTERP_KERNEL std::cout << "Number of computed intersections = " << counter << std::endl; std::cout << "Global time= " << global_end - global_start << std::endl; } - return ret; + return intersector->getNumberOfColsOfResMatrix(); } } diff --git a/src/INTERP_KERNEL/PlanarIntersectorP1P1PL.txx b/src/INTERP_KERNEL/PlanarIntersectorP1P1PL.txx index ed82199f6..e62c6b9e8 100644 --- a/src/INTERP_KERNEL/PlanarIntersectorP1P1PL.txx +++ b/src/INTERP_KERNEL/PlanarIntersectorP1P1PL.txx @@ -41,30 +41,30 @@ namespace INTERP_KERNEL void PlanarIntersectorP1P1PL::intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res) { std::vector CoordsT; - PlanarIntersector::getRealTargetCoordinates(OTT::indFC(icellT),CoordsT); + this->getRealTargetCoordinates(OTT::indFC(icellT),CoordsT); ConnType nbOfNodesT=ToConnType(CoordsT.size())/SPACEDIM; for(typename std::vector::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++) { - NormalizedCellType tS=PlanarIntersector::_meshS.getTypeOfElement(OTT::indFC(*iter)); + NormalizedCellType tS=this->_meshS.getTypeOfElement(OTT::indFC(*iter)); if(tS!=NORM_TRI3) throw INTERP_KERNEL::Exception("Invalid source cell detected for meshdim==2. Only TRI3 supported !"); std::vector CoordsS; - PlanarIntersector::getRealSourceCoordinates(OTT::indFC(*iter),CoordsS); + this->getRealSourceCoordinates(OTT::indFC(*iter),CoordsS); std::vector CoordsTTmp(CoordsT); if(SPACEDIM==3) - PlanarIntersector::projectionThis(&CoordsS[0],&CoordsTTmp[0],ToConnType(CoordsS.size())/SPACEDIM,nbOfNodesT); - const ConnType *startOfCellNodeConnT=PlanarIntersector::_connectT+OTT::conn2C(PlanarIntersector::_connIndexT[icellT]); + this->projectionThis(CoordsS.data(),CoordsTTmp.data(),ToConnType(CoordsS.size())/SPACEDIM,nbOfNodesT); + const ConnType *startOfCellNodeConnT=this->_connectT+OTT::conn2C(this->_connIndexT[icellT]); for(int nodeIdT=0;nodeIdT::ind2C(startOfCellNodeConnT[nodeIdT])]; - if( PointLocatorAlgos::isElementContainsPointAlg2D(&CoordsTTmp[nodeIdT*SPACEDIM],&CoordsS[0],3,PlanarIntersector::_precision) ) + if( PointLocatorAlgos::isElementContainsPointAlg2D(CoordsTTmp.data()+nodeIdT*SPACEDIM,CoordsS.data(),3,this->_precision) ) { double resLoc[3]; barycentric_coords(&CoordsS[0],&CoordsTTmp[nodeIdT*SPACEDIM],resLoc); - const ConnType *startOfCellNodeConnS=PlanarIntersector::_connectS+OTT::conn2C(PlanarIntersector::_connIndexS[*iter]); + const ConnType *startOfCellNodeConnS=this->_connectS+OTT::conn2C(this->_connIndexS[*iter]); for(int nodeIdS=0;nodeIdS<3;nodeIdS++) { - if(fabs(resLoc[nodeIdS])>PlanarIntersector::_precision) + if(fabs(resLoc[nodeIdS])>this->_precision) { ConnType curNodeSInCmode=OTT::coo2C(startOfCellNodeConnS[nodeIdS]); typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT::indFC(curNodeSInCmode)); @@ -86,13 +86,13 @@ namespace INTERP_KERNEL template typename MyMeshType::MyConnType PlanarIntersectorP1P1PL::getNumberOfRowsOfResMatrix() const { - return PlanarIntersector::_meshT.getNumberOfNodes(); + return this->_meshT.getNumberOfNodes(); } template typename MyMeshType::MyConnType PlanarIntersectorP1P1PL::getNumberOfColsOfResMatrix() const { - return PlanarIntersector::_meshS.getNumberOfNodes(); + return this->_meshS.getNumberOfNodes(); } } diff --git a/src/MEDCoupling/MEDCouplingRemapper.cxx b/src/MEDCoupling/MEDCouplingRemapper.cxx index 884132a50..f31c9036e 100755 --- a/src/MEDCoupling/MEDCouplingRemapper.cxx +++ b/src/MEDCoupling/MEDCouplingRemapper.cxx @@ -424,9 +424,27 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() { 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); @@ -479,8 +497,8 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() 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) { @@ -700,10 +718,11 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyEE() MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D); MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D); std::vector > 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, diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py index 773e0a147..bdb243515 100644 --- a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py @@ -955,7 +955,7 @@ class MEDCouplingBasicsTest(unittest.TestCase): 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) @@ -1301,7 +1301,7 @@ class MEDCouplingBasicsTest(unittest.TestCase): 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 @@ -1370,6 +1370,72 @@ class MEDCouplingBasicsTest(unittest.TestCase): 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)):