From: Anthony Geay Date: Tue, 21 Jul 2015 07:28:25 +0000 (+0200) Subject: PointLocator implementation for P1P1 for 1D and 2D curve. X-Git-Tag: V7_7_0a1~7 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;ds=sidebyside;h=2f5cf6aec964845422b9a03c838dd26009c1f386;p=modules%2Fmed.git PointLocator implementation for P1P1 for 1D and 2D curve. --- diff --git a/src/INTERP_KERNEL/CurveIntersector.hxx b/src/INTERP_KERNEL/CurveIntersector.hxx index f2f846e1a..e1c3e69ce 100644 --- a/src/INTERP_KERNEL/CurveIntersector.hxx +++ b/src/INTERP_KERNEL/CurveIntersector.hxx @@ -41,12 +41,16 @@ 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); protected : - bool getRealTargetCoordinates(ConnType icellT, std::vector& coordsT); - bool getRealSourceCoordinates(ConnType icellS, std::vector& coordsS); - double intersectSegments(double *Coords_T, double *Coords_S); - + 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; + bool getRealSourceCoordinates(ConnType icellS, std::vector& coordsS) 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; + struct TDualSegment { std::vector _coords; diff --git a/src/INTERP_KERNEL/CurveIntersector.txx b/src/INTERP_KERNEL/CurveIntersector.txx index 3dc7ff60a..c311605f7 100644 --- a/src/INTERP_KERNEL/CurveIntersector.txx +++ b/src/INTERP_KERNEL/CurveIntersector.txx @@ -103,12 +103,9 @@ namespace INTERP_KERNEL } } - //================================================================================ /*! Computes the bouding box of a given element. iP in numPol mode. */ - //================================================================================ - template void CurveIntersector::getElemBB (double* bb, const MyMeshType& mesh, @@ -138,15 +135,26 @@ 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] 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) + { + double deno(endOfSeg-startOfSeg); + startPos=(endOfSeg-pt)/deno; + endPos=1.-startPos; + return startPos>=0. && endPos>=0.; + } - //================================================================================ /*! Readjusts a set of bounding boxes so that they are extended in all dimensions for avoiding missing interesting intersections \param bbox vector containing the bounding boxes */ - //================================================================================ - template void CurveIntersector::adjustBoundingBoxes (std::vector& bbox, double adjustmentEpsAbs) @@ -162,7 +170,6 @@ namespace INTERP_KERNEL } } - //================================================================================ /*! * @param icellT id in target mesh in format of MyMeshType. * @param coordsT output val that stores coordinates of the target cell @@ -170,14 +177,10 @@ namespace INTERP_KERNEL * @return true if segment is quadratic and in this case coordinates of medium node * are placed in the middle of coordsT */ - //================================================================================ - template - bool CurveIntersector::getRealTargetCoordinates - (ConnType icellT, std::vector& coordsT) + bool CurveIntersector::getRealTargetCoordinates(ConnType icellT, std::vector& coordsT) const { - int nbNodesT = _connIndexT[OTT::ind2C(icellT)+1] - - _connIndexT[OTT::ind2C(icellT)]; + int nbNodesT(_connIndexT[OTT::ind2C(icellT)+1] - _connIndexT[OTT::ind2C(icellT)]); coordsT.resize(SPACEDIM*nbNodesT); for (ConnType iT=0; iT + typename MyMeshType::MyConnType CurveIntersector::getNodeIdOfTargetCellAt(ConnType icellT, ConnType nodeIdInCellT) const + { + int nbNodesT(_connIndexT[OTT::ind2C(icellT)+1] - _connIndexT[OTT::ind2C(icellT)]); + if(nodeIdInCellT>=0 && nodeIdInCellT::coo2C(_connectT[OTT::conn2C(_connIndexT[OTT::ind2C(icellT)]+nodeIdInCellT)]); + else + throw Exception("getNodeIdOfTargetCellAt : error in nodeId in cell"); + } - //================================================================================ /*! * @param icellS id in source mesh in format of MyMeshType. * @param coordsS output val that stores coordinates of the source cell automatically resized to the right length. * @return true if segment is quadratic and in this case coordinates of medium node * are placed in the middle of coordsS */ - //================================================================================ - template - bool CurveIntersector::getRealSourceCoordinates - (ConnType icellS, std::vector& coordsS) + bool CurveIntersector::getRealSourceCoordinates(ConnType icellS, std::vector& coordsS) const { - int nbNodesS = _connIndexS[OTT::ind2C(icellS)+1] - - _connIndexS[OTT::ind2C(icellS)]; + int nbNodesS = _connIndexS[OTT::ind2C(icellS)+1] - _connIndexS[OTT::ind2C(icellS)]; coordsS.resize(SPACEDIM*nbNodesS); for(ConnType iS=0; iS + typename MyMeshType::MyConnType CurveIntersector::getNodeIdOfSourceCellAt(ConnType icellS, ConnType nodeIdInCellS) const + { + int nbNodesS(_connIndexS[OTT::ind2C(icellS)+1] - _connIndexS[OTT::ind2C(icellS)]); + if(nodeIdInCellS>=0 && nodeIdInCellS::coo2C(_connectS[OTT::conn2C(_connIndexS[OTT::ind2C(icellS)]+nodeIdInCellS)]); + else + throw Exception("getNodeIdOfSourceCellAt : error in nodeId in cell"); + } + /*! * \brief Return dual segments of given segment * \param icell - given segment in C mode * \param mesh - mesh * \param segments - dual segments */ - //================================================================================ - template void CurveIntersector::getDualSegments(ConnType icell, const MyMeshType& mesh, @@ -300,18 +315,12 @@ namespace INTERP_KERNEL } } - //================================================================================ - /*! - * \brief Return length of intersection of two segments - */ - //================================================================================ - template - double CurveIntersector::intersectSegments(double *Coords_T, - double *Coords_S) + bool CurveIntersector::projectionThis(const double *coordsT, const double *coordsS, + double& xs0, double& xs1, double& xt0, double& xt1) const { - double xt0 = Coords_T[0], xt1 = Coords_T[1]; - double xs0 = Coords_S[0], xs1 = Coords_S[1]; + xt0 = coordsT[0]; xt1 = coordsT[1]; + xs0 = coordsS[0]; xs1 = coordsS[1]; if ( SPACEDIM == 2 ) { // Pass 2D->1D @@ -320,15 +329,16 @@ namespace INTERP_KERNEL // check if two segments overlap in 2D within tolerance - double* t0 = Coords_T; - double* t1 = Coords_T + 2; + const double* t0 = coordsT; + const double* t1 = coordsT + 2; double t01[2] = { t1[X]-t0[X], t1[Y]-t0[Y] }; // tgt segment direction double tSize = sqrt( t01[X]*t01[X] + t01[Y]*t01[Y] ); // tgt segment size - if ( tSize < _precision ) return 0; // degenerated segment + if ( tSize < _precision ) + return false; // degenerated segment t01[X] /= tSize, t01[Y] /= tSize; // normalize t01 - double* s0 = Coords_S; - double* s1 = Coords_S + 2; + const double* s0 = coordsS; + const double* s1 = coordsS + 2; double t0s0[2] = { s0[X]-t0[X], s0[Y]-t0[Y] }; double t0s1[2] = { s1[X]-t0[X], s1[Y]-t0[Y] }; double nt01_x_t0s0 = t0s0[X] * t01[Y] - t0s0[Y] * t01[X]; // t0s0 dot norm of t01 @@ -338,7 +348,7 @@ namespace INTERP_KERNEL bool s0_out_of_tol = ( dist_ts0 > _tolerance ); bool s1_out_of_tol = ( dist_ts1 > _tolerance ); if ( nt01_x_t0s0 * nt01_x_t0s1 > 0 && ( s0_out_of_tol || s1_out_of_tol )) - return 0; // tgt segment is to far from src segment + return false; // tgt segment is to far from src segment double S0[2] = { s0[X], s0[Y] }; double S1[2] = { s1[X], s1[Y] }; @@ -361,7 +371,8 @@ namespace INTERP_KERNEL 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 < _precision ) return 0; // degenerated segment + if ( sSize < _precision ) + return false; // degenerated segment s01[X] /= sSize, s01[Y] /= sSize; // normalize s01 // make t01 and s01 codirected @@ -375,7 +386,7 @@ namespace INTERP_KERNEL }; double medianSize = sqrt( medianDir[X]*medianDir[X] + medianDir[Y]*medianDir[Y] ); if ( medianSize < std::numeric_limits::min() ) - return 0; // strange... + return false; // strange... medianDir[X] /= medianSize, medianDir[Y] /= medianSize; xt0 = t0[X] * medianDir[X] + t0[Y] * medianDir[Y]; @@ -384,7 +395,18 @@ namespace INTERP_KERNEL xs1 = S1[X] * medianDir[X] + S1[Y] * medianDir[Y]; } // if ( SPACEDIM == 2 ) - + return true; + } + + /*! + * \brief Return length of intersection of two segments + */ + template + double CurveIntersector::intersectSegmentsInternal(const double *coordsT, const double *coordsS, double& xs0, double& xs1, double& xt0, double& xt1) const + { + if(!projectionThis(coordsT,coordsS,xs0,xs1,xt0,xt1)) + return 0.; + if ( xt0 > xt1 ) std::swap( xt0, xt1 ); if ( xs0 > xs1 ) std::swap( xs0, xs1 ); @@ -392,6 +414,16 @@ namespace INTERP_KERNEL double x1 = std::min( xt1, xs1 ); return ( x0 < x1 ) ? ( x1 - x0 ) : 0.; } + + /*! + * \brief Return length of intersection of two segments + */ + template + double CurveIntersector::intersectSegments(const double *coordsT, const double *coordsS) const + { + double xs0,xs1,xt0,xt1; + return intersectSegmentsInternal(coordsT,coordsS,xs0,xs1,xt0,xt1); + } } diff --git a/src/INTERP_KERNEL/CurveIntersectorP1P1PL.hxx b/src/INTERP_KERNEL/CurveIntersectorP1P1PL.hxx new file mode 100644 index 000000000..d9bab059b --- /dev/null +++ b/src/INTERP_KERNEL/CurveIntersectorP1P1PL.hxx @@ -0,0 +1,50 @@ +// Copyright (C) 2007-2015 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author : Anthony Geay (EDF R&D) + +#ifndef __CURVEINTERSECTORP1P1PL_HXX__ +#define __CURVEINTERSECTORP1P1PL_HXX__ + +#include "CurveIntersector.hxx" + +namespace INTERP_KERNEL +{ + template + class CurveIntersectorP1P1PL : public CurveIntersector + { + public: + static const int SPACEDIM=MyMeshType::MY_SPACEDIM; + static const int MESHDIM=MyMeshType::MY_MESHDIM; + typedef typename MyMeshType::MyConnType ConnType; + static const NumberingPolicy numPol=MyMeshType::My_numPol; + + CurveIntersectorP1P1PL(const MyMeshType& meshT, const MyMeshType& meshS, + double precision, double tolerance, + double medianLine, int printLevel); + public: + void intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res); + int getNumberOfRowsOfResMatrix() const; + int getNumberOfColsOfResMatrix() const; + private: + static void AppendValueInMatrix(MyMatrix& res, ConnType nodeIdT, ConnType nodeIdS0, double val0, ConnType nodeIdS1, double val1); + static void AppendValueInMatrix2(typename MyMatrix::value_type& resRow, ConnType nodeIdS0, double val0); + }; +} + +#endif diff --git a/src/INTERP_KERNEL/CurveIntersectorP1P1PL.txx b/src/INTERP_KERNEL/CurveIntersectorP1P1PL.txx new file mode 100644 index 000000000..02b09eb91 --- /dev/null +++ b/src/INTERP_KERNEL/CurveIntersectorP1P1PL.txx @@ -0,0 +1,112 @@ +// Copyright (C) 2007-2015 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author : Anthony Geay (EDF R&D) + +#ifndef __CURVEINTERSECTORP1P1PL_TXX__ +#define __CURVEINTERSECTORP1P1PL_TXX__ + +#include "CurveIntersectorP1P1PL.hxx" +#include "CurveIntersector.txx" + +#include + +namespace INTERP_KERNEL +{ + template + CurveIntersectorP1P1PL::CurveIntersectorP1P1PL(const MyMeshType& meshT, const MyMeshType& meshS, double precision, double tolerance, double medianLine, int printLevel):CurveIntersector(meshT, meshS, precision, tolerance, medianLine, printLevel) + { + } + + template + int CurveIntersectorP1P1PL::getNumberOfRowsOfResMatrix() const + { + return CurveIntersector::_meshT.getNumberOfNodes(); + } + + template + int CurveIntersectorP1P1PL::getNumberOfColsOfResMatrix() const + { + return CurveIntersector::_meshS.getNumberOfNodes(); + } + + template + void CurveIntersectorP1P1PL::AppendValueInMatrix2(typename MyMatrix::value_type& resRow, ConnType nodeIdS0, double val0) + { + typename MyMatrix::value_type::const_iterator iterRes(resRow.find(OTT::indFC(nodeIdS0))); + if(iterRes==resRow.end()) + { + resRow.insert(std::make_pair(OTT::indFC(nodeIdS0),val0)); + } + else + { + double val((*iterRes).second+val0); + resRow.erase(OTT::indFC(nodeIdS0)); + resRow.insert(std::make_pair(OTT::indFC(nodeIdS0),val)); + } + } + + template + 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); + } + + template + void CurveIntersectorP1P1PL::intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res) + { + std::vector coordsT; + if(CurveIntersector::getRealTargetCoordinates(icellT,coordsT)) + throw INTERP_KERNEL::Exception("Invalid target cell detected for meshdim==1. Only SEG2 supported !"); + assert(coordsT.size()/SPACEDIM==2); + for(typename std::vector::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++) + { + std::vector coordsS; + if(CurveIntersector::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.) + { + double a,b; + // for first + ConnType nodeIdT0(CurveIntersector::getNodeIdOfTargetCellAt(icellT,0)); + if(CurveIntersector::ComputeBaryCoordsOf(xs0,xs1,xt0,a,b)) + { + 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); + } + } + } + } +} + +#endif diff --git a/src/INTERP_KERNEL/InterpolationCurve.txx b/src/INTERP_KERNEL/InterpolationCurve.txx index 35e789172..239370a61 100644 --- a/src/INTERP_KERNEL/InterpolationCurve.txx +++ b/src/INTERP_KERNEL/InterpolationCurve.txx @@ -26,6 +26,7 @@ #include "CurveIntersectorP1P0.txx" #include "CurveIntersectorP0P1.txx" #include "CurveIntersectorP1P1.txx" +#include "CurveIntersectorP1P1PL.txx" #include "BBTree.txx" #include @@ -90,39 +91,78 @@ namespace INTERP_KERNEL CurveIntersector* intersector=0; if(method=="P0P0") { - intersector = new CurveIntersectorP0P0 - (myMeshT, myMeshS, - InterpolationOptions::getPrecision(), - InterpolationOptions::getBoundingBoxAdjustmentAbs(), - InterpolationOptions::getMedianPlane(), - InterpolationOptions::getPrintLevel()); + switch (InterpolationOptions::getIntersectionType()) + { + case Triangulation: + { + intersector = new CurveIntersectorP0P0(myMeshT, myMeshS, + InterpolationOptions::getPrecision(), + InterpolationOptions::getBoundingBoxAdjustmentAbs(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getPrintLevel()); + break; + } + default: + throw INTERP_KERNEL::Exception("For P0P0 in 1D or 2D curve only Triangulation supported for the moment !"); + } } else if(method=="P0P1") { - intersector = new CurveIntersectorP0P1 - (myMeshT, myMeshS, - InterpolationOptions::getPrecision(), - InterpolationOptions::getBoundingBoxAdjustmentAbs(), - InterpolationOptions::getMedianPlane(), - InterpolationOptions::getPrintLevel()); + switch (InterpolationOptions::getIntersectionType()) + { + case Triangulation: + { + intersector = new CurveIntersectorP0P1(myMeshT, myMeshS, + InterpolationOptions::getPrecision(), + InterpolationOptions::getBoundingBoxAdjustmentAbs(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getPrintLevel()); + break; + } + default: + throw INTERP_KERNEL::Exception("For P0P1 in 1D or 2D curve only Triangulation supported for the moment !"); + } } else if(method=="P1P0") { - intersector = new CurveIntersectorP1P0 - (myMeshT, myMeshS, - InterpolationOptions::getPrecision(), - InterpolationOptions::getBoundingBoxAdjustmentAbs(), - InterpolationOptions::getMedianPlane(), - InterpolationOptions::getPrintLevel()); + switch (InterpolationOptions::getIntersectionType()) + { + case Triangulation: + { + intersector = 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 !"); + } } else if(method=="P1P1") { - intersector = new CurveIntersectorP1P1 - (myMeshT, myMeshS, - InterpolationOptions::getPrecision(), - InterpolationOptions::getBoundingBoxAdjustmentAbs(), - InterpolationOptions::getMedianPlane(), - InterpolationOptions::getPrintLevel()); + switch (InterpolationOptions::getIntersectionType()) + { + case Triangulation: + intersector = 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()); + break; + default: + throw INTERP_KERNEL::Exception("For P1P1 in 1D or 2D curve only Triangulation and PointLocator supported !"); + } } else throw INTERP_KERNEL::Exception("Invalid method specified ! Must be in : \"P0P0\" \"P0P1\" \"P1P0\" or \"P1P1\""); diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py index 469f92463..ba5ba6428 100644 --- a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py @@ -786,6 +786,74 @@ class MEDCouplingBasicsTest(unittest.TestCase): self.assertTrue(coarse.isEqual(trgField.getArray(),1e-12)) pass + @unittest.skipUnless(MEDCouplingHasNumPyBindings() and MEDCouplingHasSciPyBindings(),"requires numpy AND scipy") + def test1DPointLocator1(self): + """This test focuses on PointLocator for P1P1 in 1D and 2DCurve.""" + from numpy import array + from scipy.sparse import diags,csr_matrix,identity + ## basic case 1D + arrS=DataArrayInt.Range(0,11,1).convertToDblArr() + arrT=DataArrayDouble([0.1,1.7,5.5,9.6]) + mS=MEDCouplingCMesh() ; mS.setCoords(arrS) + mT=MEDCouplingCMesh() ; mT.setCoords(arrT) + rem=MEDCouplingRemapper() + rem.setIntersectionType(PointLocator) + self.assertEqual(rem.prepare(mS.buildUnstructured(),mT.buildUnstructured(),"P1P1"),1) + m=rem.getCrudeCSRMatrix() + rowSum=m.sum(axis=1) + m=diags(array(1/rowSum.transpose()),[0])*m + # expected matrix + row=array([0,0,1,1,2,2,3,3]) + col=array([0,1,1,2,5,6,9,10]) + data=array([0.9,0.1,0.3,0.7,0.5,0.5,0.4,0.6]) + mExp0=csr_matrix((data,(row,col)),shape=(4,11)) + # compute diff and check + diff=abs(m-mExp0) + self.assertAlmostEqual(diff.max(),0.,14) + ## full specific case 1D where target=source + rem=MEDCouplingRemapper() + rem.setIntersectionType(PointLocator) + self.assertEqual(rem.prepare(mS.buildUnstructured(),mS.buildUnstructured(),"P1P1"),1) + m=rem.getCrudeCSRMatrix() + rowSum=m.sum(axis=1) + m=diags(array(1/rowSum.transpose()),[0])*m + # expected matrix + mExp1=identity(11) + diff=abs(m-mExp1) + self.assertAlmostEqual(diff.max(),0.,14) + ## case where some points in target are not in source + arrT=DataArrayDouble([-0.2,0.1,1.7,5.5,10.3]) + mT=MEDCouplingCMesh() ; mT.setCoords(arrT) + mT=mT.buildUnstructured() + rem=MEDCouplingRemapper() + rem.setIntersectionType(PointLocator) + self.assertEqual(rem.prepare(mS.buildUnstructured(),mT,"P1P1"),1) + 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]) + mExp2=csr_matrix((data,(row,col)),shape=(5,11)) + diff=abs(m-mExp2) + self.assertAlmostEqual(diff.max(),0.,14) + ## basic case 2D Curve + arrS=DataArrayInt.Range(0,11,1).convertToDblArr() + arrT=DataArrayDouble([0.1,1.7,5.5,9.6]) + mS=MEDCouplingCMesh() ; mS.setCoords(arrS) + mT=MEDCouplingCMesh() ; mT.setCoords(arrT) + mS=mS.buildUnstructured() ; mS.changeSpaceDimension(2) + mT=mT.buildUnstructured() ; mT.changeSpaceDimension(2) + mS.rotate([-1.,-1.],1.2) + mT.rotate([-1.,-1.],1.2) + rem=MEDCouplingRemapper() + rem.setIntersectionType(PointLocator) + self.assertEqual(rem.prepare(mS,mT,"P1P1"),1) + m=rem.getCrudeCSRMatrix() + rowSum=m.sum(axis=1) + m=diags(array(1/rowSum.transpose()),[0])*m + diff=abs(m-mExp0) + self.assertAlmostEqual(diff.max(),0.,14) + pass + def build2DSourceMesh_1(self): sourceCoords=[-0.3,-0.3, 0.7,-0.3, -0.3,0.7, 0.7,0.7] sourceConn=[0,3,1,0,2,3]