void createBoundingBoxes(const MyMeshType& mesh, std::vector<double>& bbox);
void adjustBoundingBoxes(std::vector<double>& bbox, double adjustmentEpsAbs);
static void getElemBB(double* bb, const MyMeshType& mesh, ConnType iP, ConnType nb_nodes);
-
+ static bool ComputeBaryCoordsOf(double startOfSeg, double endOfSeg, double pt, double& startPos, double& endPos);
protected :
- bool getRealTargetCoordinates(ConnType icellT, std::vector<double>& coordsT);
- bool getRealSourceCoordinates(ConnType icellS, std::vector<double>& 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<double>& coordsT) const;
+ typename MyMeshType::MyConnType getNodeIdOfTargetCellAt(ConnType icellT, ConnType nodeIdInCellT) const;
+ bool getRealSourceCoordinates(ConnType icellS, std::vector<double>& 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<double> _coords;
}
}
- //================================================================================
/*!
Computes the bouding box of a given element. iP in numPol mode.
*/
- //================================================================================
-
template<class MyMeshType, class MyMatrix>
void CurveIntersector<MyMeshType,MyMatrix>::getElemBB (double* bb,
const MyMeshType& mesh,
}
}
}
+
+ /*!
+ * \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<class MyMeshType, class MyMatrix>
+ bool CurveIntersector<MyMeshType,MyMatrix>::ComputeBaryCoordsOf(double startOfSeg, double endOfSeg, double pt, double& startPos, double& endPos)
+ {
+ double deno(endOfSeg-startOfSeg);
+ startPos=(endOfSeg-pt)/deno;
+ endPos=1.-startPos;
+ return startPos>=0. && endPos>=0.;
+ }
- //================================================================================
/*! 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<class MyMeshType, class MyMatrix>
void CurveIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes (std::vector<double>& bbox,
double adjustmentEpsAbs)
}
}
- //================================================================================
/*!
* @param icellT id in target mesh in format of MyMeshType.
* @param coordsT output val that stores coordinates of the target cell
* @return true if segment is quadratic and in this case coordinates of medium node
* are placed in the middle of coordsT
*/
- //================================================================================
-
template<class MyMeshType, class MyMatrix>
- bool CurveIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates
- (ConnType icellT, std::vector<double>& coordsT)
+ bool CurveIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates(ConnType icellT, std::vector<double>& coordsT) const
{
- int nbNodesT = _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1] -
- _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)];
+ int nbNodesT(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1] - _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]);
coordsT.resize(SPACEDIM*nbNodesT);
for (ConnType iT=0; iT<nbNodesT; iT++)
{
}
return false;
}
+
+ template<class MyMeshType, class MyMatrix>
+ typename MyMeshType::MyConnType CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfTargetCellAt(ConnType icellT, ConnType nodeIdInCellT) const
+ {
+ int nbNodesT(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1] - _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]);
+ if(nodeIdInCellT>=0 && nodeIdInCellT<nbNodesT)
+ return OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::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<class MyMeshType, class MyMatrix>
- bool CurveIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinates
- (ConnType icellS, std::vector<double>& coordsS)
+ bool CurveIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinates(ConnType icellS, std::vector<double>& coordsS) const
{
- int nbNodesS = _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1] -
- _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
+ int nbNodesS = _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1] - _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
coordsS.resize(SPACEDIM*nbNodesS);
for(ConnType iS=0; iS<nbNodesS; iS++)
{
return false;
}
- //================================================================================
+ template<class MyMeshType, class MyMatrix>
+ typename MyMeshType::MyConnType CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfSourceCellAt(ConnType icellS, ConnType nodeIdInCellS) const
+ {
+ int nbNodesS(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1] - _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]);
+ if(nodeIdInCellS>=0 && nodeIdInCellS<nbNodesS)
+ return OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::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<class MyMeshType, class MyMatrix>
void CurveIntersector<MyMeshType,MyMatrix>::getDualSegments(ConnType icell,
const MyMeshType& mesh,
}
}
- //================================================================================
- /*!
- * \brief Return length of intersection of two segments
- */
- //================================================================================
-
template<class MyMeshType, class MyMatrix>
- double CurveIntersector<MyMeshType,MyMatrix>::intersectSegments(double *Coords_T,
- double *Coords_S)
+ bool CurveIntersector<MyMeshType,MyMatrix>::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
// 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
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] };
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
};
double medianSize = sqrt( medianDir[X]*medianDir[X] + medianDir[Y]*medianDir[Y] );
if ( medianSize < std::numeric_limits<double>::min() )
- return 0; // strange...
+ return false; // strange...
medianDir[X] /= medianSize, medianDir[Y] /= medianSize;
xt0 = t0[X] * medianDir[X] + t0[Y] * medianDir[Y];
xs1 = S1[X] * medianDir[X] + S1[Y] * medianDir[Y];
} // if ( SPACEDIM == 2 )
-
+ return true;
+ }
+
+ /*!
+ * \brief Return length of intersection of two segments
+ */
+ template<class MyMeshType, class MyMatrix>
+ double CurveIntersector<MyMeshType,MyMatrix>::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 );
double x1 = std::min( xt1, xs1 );
return ( x0 < x1 ) ? ( x1 - x0 ) : 0.;
}
+
+ /*!
+ * \brief Return length of intersection of two segments
+ */
+ template<class MyMeshType, class MyMatrix>
+ double CurveIntersector<MyMeshType,MyMatrix>::intersectSegments(const double *coordsT, const double *coordsS) const
+ {
+ double xs0,xs1,xt0,xt1;
+ return intersectSegmentsInternal(coordsT,coordsS,xs0,xs1,xt0,xt1);
+ }
}
--- /dev/null
+// 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 MyMeshType, class MyMatrix>
+ class CurveIntersectorP1P1PL : public CurveIntersector<MyMeshType,MyMatrix>
+ {
+ 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<ConnType>& 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
--- /dev/null
+// 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 <cassert>
+
+namespace INTERP_KERNEL
+{
+ template<class MyMeshType, class MyMatrix>
+ CurveIntersectorP1P1PL<MyMeshType,MyMatrix>::CurveIntersectorP1P1PL(const MyMeshType& meshT, const MyMeshType& meshS, double precision, double tolerance, double medianLine, int printLevel):CurveIntersector<MyMeshType,MyMatrix>(meshT, meshS, precision, tolerance, medianLine, printLevel)
+ {
+ }
+
+ template<class MyMeshType, class MyMatrix>
+ int CurveIntersectorP1P1PL<MyMeshType,MyMatrix>::getNumberOfRowsOfResMatrix() const
+ {
+ return CurveIntersector<MyMeshType,MyMatrix>::_meshT.getNumberOfNodes();
+ }
+
+ template<class MyMeshType, class MyMatrix>
+ int CurveIntersectorP1P1PL<MyMeshType,MyMatrix>::getNumberOfColsOfResMatrix() const
+ {
+ return CurveIntersector<MyMeshType,MyMatrix>::_meshS.getNumberOfNodes();
+ }
+
+ template<class MyMeshType, class MyMatrix>
+ void CurveIntersectorP1P1PL<MyMeshType,MyMatrix>::AppendValueInMatrix2(typename MyMatrix::value_type& resRow, ConnType nodeIdS0, double val0)
+ {
+ typename MyMatrix::value_type::const_iterator iterRes(resRow.find(OTT<ConnType,numPol>::indFC(nodeIdS0)));
+ if(iterRes==resRow.end())
+ {
+ resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(nodeIdS0),val0));
+ }
+ else
+ {
+ double val((*iterRes).second+val0);
+ resRow.erase(OTT<ConnType,numPol>::indFC(nodeIdS0));
+ resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(nodeIdS0),val));
+ }
+ }
+
+ template<class MyMeshType, class MyMatrix>
+ void CurveIntersectorP1P1PL<MyMeshType,MyMatrix>::AppendValueInMatrix(MyMatrix& res, ConnType nodeIdT, ConnType nodeIdS0, double val0, ConnType nodeIdS1, double val1)
+ {
+ typename MyMatrix::value_type& resRow(res[nodeIdT]);
+ AppendValueInMatrix2(resRow,nodeIdS0,val0);
+ AppendValueInMatrix2(resRow,nodeIdS1,val1);
+ }
+
+ 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))
+ throw INTERP_KERNEL::Exception("Invalid target cell detected for meshdim==1. Only SEG2 supported !");
+ assert(coordsT.size()/SPACEDIM==2);
+ 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))
+ 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.)
+ {
+ double a,b;
+ // for first
+ ConnType nodeIdT0(CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfTargetCellAt(icellT,0));
+ if(CurveIntersector<MyMeshType,MyMatrix>::ComputeBaryCoordsOf(xs0,xs1,xt0,a,b))
+ {
+ 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);
+ }
+ }
+ }
+ }
+}
+
+#endif
#include "CurveIntersectorP1P0.txx"
#include "CurveIntersectorP0P1.txx"
#include "CurveIntersectorP1P1.txx"
+#include "CurveIntersectorP1P1PL.txx"
#include "BBTree.txx"
#include <time.h>
CurveIntersector<MyMeshType,MatrixType>* intersector=0;
if(method=="P0P0")
{
- intersector = new CurveIntersectorP0P0<MyMeshType,MatrixType>
- (myMeshT, myMeshS,
- InterpolationOptions::getPrecision(),
- InterpolationOptions::getBoundingBoxAdjustmentAbs(),
- InterpolationOptions::getMedianPlane(),
- InterpolationOptions::getPrintLevel());
+ switch (InterpolationOptions::getIntersectionType())
+ {
+ case Triangulation:
+ {
+ intersector = new CurveIntersectorP0P0<MyMeshType,MatrixType>(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<MyMeshType,MatrixType>
- (myMeshT, myMeshS,
- InterpolationOptions::getPrecision(),
- InterpolationOptions::getBoundingBoxAdjustmentAbs(),
- InterpolationOptions::getMedianPlane(),
- InterpolationOptions::getPrintLevel());
+ switch (InterpolationOptions::getIntersectionType())
+ {
+ case Triangulation:
+ {
+ intersector = new CurveIntersectorP0P1<MyMeshType,MatrixType>(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<MyMeshType,MatrixType>
- (myMeshT, myMeshS,
- InterpolationOptions::getPrecision(),
- InterpolationOptions::getBoundingBoxAdjustmentAbs(),
- InterpolationOptions::getMedianPlane(),
- InterpolationOptions::getPrintLevel());
+ switch (InterpolationOptions::getIntersectionType())
+ {
+ case Triangulation:
+ {
+ intersector = 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 !");
+ }
}
else if(method=="P1P1")
{
- intersector = new CurveIntersectorP1P1<MyMeshType,MatrixType>
- (myMeshT, myMeshS,
- InterpolationOptions::getPrecision(),
- InterpolationOptions::getBoundingBoxAdjustmentAbs(),
- InterpolationOptions::getMedianPlane(),
- InterpolationOptions::getPrintLevel());
+ switch (InterpolationOptions::getIntersectionType())
+ {
+ case Triangulation:
+ intersector = 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());
+ 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\"");
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]
