From: Anthony Geay Date: Tue, 31 Dec 2019 14:54:57 +0000 (+0100) Subject: Deal with 1D/0D pointlocator into Remapper X-Git-Tag: V9_5_0a1~17 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=f92459c4cde6b0f4a9d12e9b0361591fc514ab2c;p=tools%2Fmedcoupling.git Deal with 1D/0D pointlocator into Remapper --- 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..847915e79 100644 --- a/src/INTERP_KERNEL/CurveIntersector.hxx +++ b/src/INTERP_KERNEL/CurveIntersector.hxx @@ -41,8 +41,9 @@ namespace INTERP_KERNEL void createBoundingBoxes(const MyMeshType& mesh, std::vector& bbox); void adjustBoundingBoxes(std::vector& 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& 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..c7860b292 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 @@ -138,18 +139,18 @@ namespace INTERP_KERNEL /*! * \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 - bool CurveIntersector::ComputeBaryCoordsOf(double startOfSeg, double endOfSeg, double pt, double& startPos, double& endPos) + void CurveIntersector::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 @@ -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,39 @@ namespace INTERP_KERNEL double x1 = std::min( xt1, xs1 ); return ( x0 < x1 ) ? ( x1 - x0 ) : 0.; } + + template + 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 + 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}; + const double coordsS_1D[2]={xs0,xs1}; + const double *coordsT_1D(&xt); + return PointLocatorAlgos>::isElementContainsPoint(coordsT_1D,NORM_SEG2,coordsS_1D,TAB,2,this->_precision); + } /*! * \brief Return length of intersection of two segments @@ -424,7 +490,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..820fbd525 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 @@ -65,48 +64,38 @@ namespace INTERP_KERNEL void CurveIntersectorP1P1PL::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 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; + double a,b; + this->ComputeBaryCoordsOf(xs0,xs1,xt,a,b); 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); - } } } } } - -#endif diff --git a/src/INTERP_KERNEL/InterpolationCurve.hxx b/src/INTERP_KERNEL/InterpolationCurve.hxx index af9332319..40eecfc03 100755 --- a/src/INTERP_KERNEL/InterpolationCurve.hxx +++ b/src/INTERP_KERNEL/InterpolationCurve.hxx @@ -24,6 +24,10 @@ #include "Interpolation.hxx" #include "InterpolationOptions.hxx" +#include "BBTree.txx" + +#include + namespace INTERP_KERNEL { template @@ -35,9 +39,27 @@ namespace INTERP_KERNEL // Main function to interpolate template + typename MyMeshType::MyConnType interpolateMeshesInternal(const MyMeshType& meshS, const MyMeshType& meshT, + MatrixType& result, const std::string& method, + std::function< void(const BBTree< MyMeshType::MY_SPACEDIM , typename MyMeshType::MyConnType>&, const double*, std::vector&) > bbtreeMethod); + template typename MyMeshType::MyConnType interpolateMeshes(const MyMeshType& meshS, const MyMeshType& meshT, - MatrixType& result, const std::string& method); - + MatrixType& result, const std::string& method) { return this->interpolateMeshesInternal(meshS,meshT,result,method,[](const BBTree< MyMeshType::MY_SPACEDIM , + typename MyMeshType::MyConnType>& bbtree, const double *bb, + std::vector& intersecting_elems) + { bbtree.getIntersectingElems(bb, intersecting_elems); }); } + template + typename MyMeshType::MyConnType interpolateMeshes0D(const MyMeshType& meshS, const MyMeshType& meshT, + MatrixType& result, const std::string& method) { return this->interpolateMeshesInternal(meshS,meshT,result,method,[](const BBTree< MyMeshType::MY_SPACEDIM , + typename MyMeshType::MyConnType>& bbtree, const double *bb, + std::vector& intersecting_elems) + { + double TMP[MyMeshType::MY_SPACEDIM]; + for(int i=0;i +#include namespace INTERP_KERNEL { @@ -73,10 +73,14 @@ namespace INTERP_KERNEL */ template template - typename MyMeshType::MyConnType InterpolationCurve::interpolateMeshes (const MyMeshType& myMeshS, - const MyMeshType& myMeshT, - MatrixType& result, - const std::string& method) + typename MyMeshType::MyConnType InterpolationCurve::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&) > bbtreeMethod + ) { static const int SPACEDIM=MyMeshType::MY_SPACEDIM; typedef typename MyMeshType::MyConnType ConnType; @@ -88,18 +92,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 +116,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 +133,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 +149,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 !"); @@ -193,12 +197,10 @@ namespace INTERP_KERNEL std::vector intersecting_elems; double bb[2*SPACEDIM]; intersector->getElemBB(bb,myMeshT,OTT::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) { @@ -209,7 +211,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..13128690b 100755 --- a/src/MEDCoupling/MEDCouplingRemapper.cxx +++ b/src/MEDCoupling/MEDCouplingRemapper.cxx @@ -427,6 +427,24 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() 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); @@ -703,7 +721,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyEE() 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, diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py index 773e0a147..189887c21 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) @@ -1370,13 +1370,93 @@ class MEDCouplingBasicsTest(unittest.TestCase): 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())=0) - self.assertTrue(min(mat1[i].keys())>=0) + if len(mat2[i].keys())>0: + self.assertTrue(max(mat2[i].keys())0: + self.assertTrue(max(mat1[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])