From: eap Date: Thu, 8 Apr 2010 13:57:04 +0000 (+0000) Subject: EHPOC X-Git-Tag: V5_1_main_FINAL~143 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=59082d07300172a63f6938c7f4058dddb5e859bf;p=tools%2Fmedcoupling.git EHPOC * Interpolation of 1D meshes has been implemented --- diff --git a/src/INTERP_KERNEL/CurveIntersector.hxx b/src/INTERP_KERNEL/CurveIntersector.hxx new file mode 100644 index 000000000..d139d366a --- /dev/null +++ b/src/INTERP_KERNEL/CurveIntersector.hxx @@ -0,0 +1,73 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __CURVEINTERSECTOR_HXX__ +#define __CURVEINTERSECTOR_HXX__ + +#include "TargetIntersector.hxx" +#include "NormalizedUnstructuredMesh.hxx" + +namespace INTERP_KERNEL +{ + template + class CurveIntersector : public TargetIntersector + { + 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; + public: + CurveIntersector(const MyMeshType& meshT, const MyMeshType& meshS, + double precision, double adjustmentEpsAbs, double medianLine, int printLevel); + virtual ~CurveIntersector(); + 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); + + protected : + bool getRealTargetCoordinates(ConnType icellT, std::vector& coordsT); + bool getRealSourceCoordinates(ConnType icellS, std::vector& coordsS); + double intersectSegments(double *Coords_T, double *Coords_S); + + struct TDualSegment + { + std::vector _coords; + int _nodeId; // in mesh mode + }; + static void getDualSegments(ConnType icell, + const MyMeshType& mesh, + std::vector& segments); + + protected: + const ConnType *_connectT; + const ConnType *_connectS; + const double *_coordsT; + const double *_coordsS; + const ConnType *_connIndexT; + const ConnType *_connIndexS; + const MyMeshType& _meshT; + const MyMeshType& _meshS; + double _tolerance; + double _precision; + double _median_line; + int _print_level; + }; +} + +#endif diff --git a/src/INTERP_KERNEL/CurveIntersector.txx b/src/INTERP_KERNEL/CurveIntersector.txx new file mode 100644 index 000000000..469c4ab8d --- /dev/null +++ b/src/INTERP_KERNEL/CurveIntersector.txx @@ -0,0 +1,396 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __CURVEINTERSECTOR_TXX__ +#define __CURVEINTERSECTOR_TXX__ + +#include "CurveIntersector.hxx" + +#include + +namespace INTERP_KERNEL +{ + template + CurveIntersector + ::CurveIntersector(const MyMeshType& meshT, const MyMeshType& meshS, + double precision, double tolerance, double medianLine, int printLevel): + _meshT(meshT), + _meshS(meshS), + _tolerance(tolerance), + _precision(precision), + _median_line(medianLine), + _print_level(printLevel) + { + if ( SPACEDIM != 1 && SPACEDIM != 2 ) + throw Exception("CurveIntersector(): space dimension of mesh must be 1 or 2"); + if ( MESHDIM != 1 ) + throw Exception("CurveIntersector(): mesh dimension must be 1"); + + _connectT = meshT.getConnectivityPtr(); + _connectS = meshS.getConnectivityPtr(); + _connIndexT = meshT.getConnectivityIndexPtr(); + _connIndexS = meshS.getConnectivityIndexPtr(); + _coordsT = meshT.getCoordinatesPtr(); + _coordsS = meshS.getCoordinatesPtr(); + } + + template + CurveIntersector::~CurveIntersector() + { + } + + //================================================================================ + /*! + \brief creates the bounding boxes for all the cells of mesh \a mesh + + \param mesh structure pointing to the mesh + \param bbox vector containing the bounding boxes + */ + //================================================================================ + + template + void CurveIntersector::createBoundingBoxes (const MyMeshType& mesh, + std::vector& bbox) + { + long nbelems = mesh.getNumberOfElements(); + bbox.resize(2*SPACEDIM* nbelems); + const double* coords = mesh.getCoordinatesPtr(); + const ConnType* conn = mesh.getConnectivityPtr(); + const ConnType* conn_index = mesh.getConnectivityIndexPtr(); + int ibox=0; + for(long icell=0; icell::max(); + bbox[2*SPACEDIM*ibox+2*idim+1] = -std::numeric_limits::max(); + } + //updating the bounding box with each node of the element + for (int j=0; j + ::coo2C(conn[OTT::conn2C(conn_index[icell]+j)]); + for(int idim=0; idim x ) ? bbox[ibox*2*SPACEDIM + 2*idim+1] : x; + } + } + ibox++; + } + } + + //================================================================================ + /*! + Computes the bouding box of a given element. iP in numPol mode. + */ + //================================================================================ + + template + void CurveIntersector::getElemBB (double* bb, + const MyMeshType& mesh, + ConnType iP, + ConnType nb_nodes) + { + const double* coords = mesh.getCoordinatesPtr(); + const ConnType* conn_index = mesh.getConnectivityIndexPtr(); + const ConnType* conn = mesh.getConnectivityPtr(); + //initializing bounding box limits + for(int idim=0; idim::max(); + bb[2*idim+1] = -std::numeric_limits::max(); + } + + for (ConnType i=0; i::coo2C(conn[OTT::conn2C(conn_index[OTT::ind2C(iP)]+i)])); + for(int idim=0; idimbb[2*idim+1]) ? x : bb[2*idim+1]; + } + } + } + + //================================================================================ + /*! 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) + { + long size = bbox.size()/(2*SPACEDIM); + for (int i=0; i + bool CurveIntersector::getRealTargetCoordinates + (ConnType icellT, std::vector& coordsT) + { + int nbNodesT = _connIndexT[OTT::ind2C(icellT)+1] - + _connIndexT[OTT::ind2C(icellT)]; + coordsT.resize(SPACEDIM*nbNodesT); + for (ConnType iT=0; iT::coo2C(_connectT[OTT::conn2C(_connIndexT[OTT::ind2C(icellT)]+iT)])+idim]; + } + } + if ( nbNodesT > 2 ) + { + for(int idim=0; idim + bool CurveIntersector::getRealSourceCoordinates + (ConnType icellS, std::vector& coordsS) + { + int nbNodesS = _connIndexS[OTT::ind2C(icellS)+1] - + _connIndexS[OTT::ind2C(icellS)]; + coordsS.resize(SPACEDIM*nbNodesS); + for(ConnType iS=0; iS::coo2C(_connectS[OTT::conn2C(_connIndexS[OTT::ind2C(icellS)]+iS)])+idim]; + } + } + if ( nbNodesS > 2 ) + { + for(int idim=0; idim + void CurveIntersector::getDualSegments(ConnType icell, + const MyMeshType& mesh, + std::vector& segments) + { + // get coordinates of cell nodes + int nbNodes; + std::vector ncoords; + std::vector nodeIds; + { + const ConnType *connect = mesh.getConnectivityPtr(); + const ConnType *connIndex = mesh.getConnectivityIndexPtr(); + const double *coords = mesh.getCoordinatesPtr(); + + nbNodes = connIndex[icell+1] - connIndex[icell]; + + ncoords.resize(SPACEDIM*nbNodes); + nodeIds.resize(nbNodes); + + for(int i=0; i::conn2C(connIndex[OTT::ind2C(icell)]+i)]; + ncoords[SPACEDIM*i+idim] = coords[SPACEDIM*OTT::coo2C(nodeIds[i])+idim]; + } + if ( nbNodes > 2 ) // quadratic segment, put medium node in the middle + { + for(int idim=0; idim + double CurveIntersector::intersectSegments(double *Coords_T, + double *Coords_S) + { + double xt0 = Coords_T[0], xt1 = Coords_T[1]; + double xs0 = Coords_S[0], xs1 = Coords_S[1]; + if ( SPACEDIM == 2 ) + { + // Pass 2D->1D + + enum { X=0, Y }; + + // check if two segments overlap in 2D within tolerance + + double* t0 = Coords_T; + double* t1 = Coords_T + 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 + t01[X] /= tSize, t01[Y] /= tSize; // normalize t01 + + double* s0 = Coords_S; + double* s1 = Coords_S + 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 + double nt01_x_t0s1 = t0s1[X] * t01[Y] - t0s1[Y] * t01[X]; // t0s1 dot norm of t01 + double dist_ts0 = fabs( nt01_x_t0s0 ); // dist from tgt seg to s0 + double dist_ts1 = fabs( nt01_x_t0s1 ); // dist from tgt seg to s1 + 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 + + double S0[2] = { s0[X], s0[Y] }; + double S1[2] = { s1[X], s1[Y] }; + if ( s0_out_of_tol ) // put s0 within tolerance + { + double t = _tolerance * nt01_x_t0s0 / dist_ts0; // signed tolerance + double r = ( nt01_x_t0s0 - t ) / ( nt01_x_t0s0 - nt01_x_t0s1 ); + S0[X] = s0[X] * ( 1.-r ) + s1[X] * r; + S0[Y] = s0[Y] * ( 1.-r ) + s1[Y] * r; + } + if ( s1_out_of_tol ) // put s1 within tolerance + { + double t = _tolerance * nt01_x_t0s1 / dist_ts1; // signed tolerance + double r = ( nt01_x_t0s1 - t ) / ( nt01_x_t0s1 - nt01_x_t0s0 ); + S1[X] = s1[X] * ( 1.-r ) + s0[X] * r; + S1[Y] = s1[Y] * ( 1.-r ) + s0[Y] * r; + } + + // project tgt and src segments to median line + + 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 + s01[X] /= sSize, s01[Y] /= sSize; // normalize s01 + + // make t01 and s01 codirected + double t01_x_s01 = t01[X] * s01[X] + t01[Y] * s01[Y]; // t01 dot s01 + if ( t01_x_s01 < 0 ) + s01[X] = -s01[X], s01[Y] = -s01[Y]; + + double medianDir[2] = { + t01[X] * ( 1.-_median_line) + s01[X] * _median_line, + t01[Y] * ( 1.-_median_line) + s01[Y] * _median_line + }; + double medianSize = sqrt( medianDir[X]*medianDir[X] + medianDir[Y]*medianDir[Y] ); + if ( medianSize < std::numeric_limits::min() ) + return 0; // strange... + medianDir[X] /= medianSize, medianDir[Y] /= medianSize; + + xt0 = t0[X] * medianDir[X] + t0[Y] * medianDir[Y]; + xt1 = t1[X] * medianDir[X] + t1[Y] * medianDir[Y]; + xs0 = S0[X] * medianDir[X] + S0[Y] * medianDir[Y]; + xs1 = S1[X] * medianDir[X] + S1[Y] * medianDir[Y]; + + } // if ( SPACEDIM == 2 ) + + if ( xt0 > xt1 ) std::swap( xt0, xt1 ); + if ( xs0 > xs1 ) std::swap( xs0, xs1 ); + + double x0 = std::max( xt0, xs0 ); + double x1 = std::min( xt1, xs1 ); + return ( x0 < x1 ) ? ( x1 - x0 ) : 0.; + } + +} + +#endif diff --git a/src/INTERP_KERNEL/CurveIntersectorP0P0.hxx b/src/INTERP_KERNEL/CurveIntersectorP0P0.hxx new file mode 100644 index 000000000..3638f49b2 --- /dev/null +++ b/src/INTERP_KERNEL/CurveIntersectorP0P0.hxx @@ -0,0 +1,46 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __CURVEINTERSECTORP0P0_HXX__ +#define __CURVEINTERSECTORP0P0_HXX__ + +#include "CurveIntersector.hxx" + +namespace INTERP_KERNEL +{ + template + class CurveIntersectorP0P0 : 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; + + CurveIntersectorP0P0(const MyMeshType& meshT, const MyMeshType& meshS, + double precision, double tolerance, + double medianLine, int printLevel); + public: + int getNumberOfRowsOfResMatrix() const; + int getNumberOfColsOfResMatrix() const; + void intersectCells(ConnType icellT, + const std::vector& icellsS, MyMatrix& res); + }; +} + +#endif diff --git a/src/INTERP_KERNEL/CurveIntersectorP0P0.txx b/src/INTERP_KERNEL/CurveIntersectorP0P0.txx new file mode 100644 index 000000000..d3feea029 --- /dev/null +++ b/src/INTERP_KERNEL/CurveIntersectorP0P0.txx @@ -0,0 +1,77 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +#ifndef __CURVEINTERSECTORP0P0_TXX__ +#define __CURVEINTERSECTORP0P0_TXX__ + +#include "CurveIntersectorP0P0.hxx" +#include "CurveIntersector.txx" + +#define BASE_INTERSECTOR CurveIntersector + +namespace INTERP_KERNEL +{ + template + CurveIntersectorP0P0 + ::CurveIntersectorP0P0(const MyMeshType& meshT, const MyMeshType& meshS, + double precision, double tolerance, + double medianLine, int printLevel): + BASE_INTERSECTOR(meshT, meshS, precision, tolerance, medianLine, printLevel) + { + } + + template + int CurveIntersectorP0P0 + ::getNumberOfRowsOfResMatrix() const + { + return BASE_INTERSECTOR::_meshT.getNumberOfElements(); + } + + template + int CurveIntersectorP0P0 + ::getNumberOfColsOfResMatrix() const + { + return BASE_INTERSECTOR::_meshS.getNumberOfElements(); + } + + template + void CurveIntersectorP0P0 + ::intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res) + { + typename MyMatrix::value_type& resRow = res[icellT]; + std::vector coordsT; + int t, nbSegT = 1 + BASE_INTERSECTOR::getRealTargetCoordinates(icellT,coordsT); + for ( t = 0; t < nbSegT; ++t ) + for(typename std::vector::const_iterator + iter=icellsS.begin(); iter!=icellsS.end(); iter++) + { + int iS = *iter; + std::vector coordsS; + int s, nbSegS = 1 + BASE_INTERSECTOR::getRealSourceCoordinates(iS,coordsS); + for ( s = 0; s < nbSegS; ++s ) + { + double surf = BASE_INTERSECTOR::intersectSegments(coordsT.data() + t*SPACEDIM, + coordsS.data() + s*SPACEDIM); + if(surf!=0.) + resRow.insert(std::make_pair(OTT::indFC(iS),surf)); + } + } + } +} +#undef BASE_INTERSECTOR + +#endif diff --git a/src/INTERP_KERNEL/CurveIntersectorP0P1.hxx b/src/INTERP_KERNEL/CurveIntersectorP0P1.hxx new file mode 100644 index 000000000..021002d21 --- /dev/null +++ b/src/INTERP_KERNEL/CurveIntersectorP0P1.hxx @@ -0,0 +1,46 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __CURVEINTERSECTORP0P1_HXX__ +#define __CURVEINTERSECTORP0P1_HXX__ + +#include "CurveIntersector.hxx" + +namespace INTERP_KERNEL +{ + template + class CurveIntersectorP0P1 : 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; + + CurveIntersectorP0P1(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; + }; +} + +#endif diff --git a/src/INTERP_KERNEL/CurveIntersectorP0P1.txx b/src/INTERP_KERNEL/CurveIntersectorP0P1.txx new file mode 100644 index 000000000..6ff2ac993 --- /dev/null +++ b/src/INTERP_KERNEL/CurveIntersectorP0P1.txx @@ -0,0 +1,95 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +#ifndef __CurveIntersectorP0P1_TXX__ +#define __CurveIntersectorP0P1_TXX__ + +#include "CurveIntersectorP0P1.hxx" +#include "CurveIntersector.txx" + +#define BASE_INTERSECTOR CurveIntersector + +namespace INTERP_KERNEL +{ + template + CurveIntersectorP0P1 + ::CurveIntersectorP0P1(const MyMeshType& meshT, const MyMeshType& meshS, + double precision, double tolerance, + double medianLine, int printLevel): + BASE_INTERSECTOR(meshT, meshS, precision, tolerance, medianLine, printLevel) + { + } + + template + int CurveIntersectorP0P1 + ::getNumberOfRowsOfResMatrix() const + { + return BASE_INTERSECTOR::_meshT.getNumberOfNodes(); + } + + template + int CurveIntersectorP0P1 + ::getNumberOfColsOfResMatrix() const + { + return BASE_INTERSECTOR::_meshS.getNumberOfElements(); + } + + //================================================================================ + /*! + * \brief Project from segments to nodes + */ + //================================================================================ + + template + void CurveIntersectorP0P1 + ::intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res) + { + std::vector segmentsT; + BASE_INTERSECTOR::getDualSegments( icellT, BASE_INTERSECTOR::_meshT, segmentsT); + for ( int t = 0; t < segmentsT.size(); ++t ) + { + typename MyMatrix::value_type& resRow = res[ OTT::ind2C( segmentsT[t]._nodeId )]; + for(typename std::vector::const_iterator + iter=icellsS.begin(); iter!=icellsS.end(); iter++) + { + int iS = *iter; + std::vector coordsS; + int s, nbSegS = 1 + BASE_INTERSECTOR::getRealSourceCoordinates(iS,coordsS); + for ( s = 0; s < nbSegS; ++s ) + { + double surf = BASE_INTERSECTOR::intersectSegments(segmentsT[t]._coords.data(), + coordsS.data() + s*SPACEDIM); + if(surf!=0.) + { + typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT::indFC(iS)); + if(iterRes==resRow.end()) + resRow.insert(std::make_pair(OTT::indFC(iS),surf)); + else + { + surf+=(*iterRes).second; + resRow.erase(OTT::indFC(iS)); + resRow.insert(std::make_pair(OTT::indFC(iS),surf)); + } + } + } + } + } + } +} +#undef BASE_INTERSECTOR + +#endif diff --git a/src/INTERP_KERNEL/CurveIntersectorP1P0.hxx b/src/INTERP_KERNEL/CurveIntersectorP1P0.hxx new file mode 100644 index 000000000..11fdc4766 --- /dev/null +++ b/src/INTERP_KERNEL/CurveIntersectorP1P0.hxx @@ -0,0 +1,44 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __CURVEINTERSECTORP1P0_HXX__ +#define __CURVEINTERSECTORP1P0_HXX__ + +#include "CurveIntersector.hxx" + +namespace INTERP_KERNEL +{ + template + class CurveIntersectorP1P0 : 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; + + CurveIntersectorP1P0(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; + }; +} + +#endif diff --git a/src/INTERP_KERNEL/CurveIntersectorP1P0.txx b/src/INTERP_KERNEL/CurveIntersectorP1P0.txx new file mode 100644 index 000000000..49008d088 --- /dev/null +++ b/src/INTERP_KERNEL/CurveIntersectorP1P0.txx @@ -0,0 +1,95 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +#ifndef __CurveIntersectorP1P0_TXX__ +#define __CurveIntersectorP1P0_TXX__ + +#include "CurveIntersectorP1P0.hxx" +#include "CurveIntersector.txx" + +#define BASE_INTERSECTOR CurveIntersector + +namespace INTERP_KERNEL +{ + template + CurveIntersectorP1P0 + ::CurveIntersectorP1P0(const MyMeshType& meshT, const MyMeshType& meshS, + double precision, double tolerance, + double medianLine, int printLevel): + BASE_INTERSECTOR (meshT, meshS, precision, tolerance, medianLine, printLevel) + { + } + + template + int CurveIntersectorP1P0 + ::getNumberOfRowsOfResMatrix() const + { + return BASE_INTERSECTOR::_meshT.getNumberOfElements(); + } + + template + int CurveIntersectorP1P0 + ::getNumberOfColsOfResMatrix() const + { + return BASE_INTERSECTOR::_meshS.getNumberOfNodes(); + } + + //================================================================================ + /*! + * \brief Project from source nodes to target segments + */ + //================================================================================ + + template + void CurveIntersectorP1P0 + ::intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res) + { + typename MyMatrix::value_type& resRow = res[ icellT ]; + std::vector segmentsS; + std::vector coordsT; + int t, nbSegT = 1 + BASE_INTERSECTOR::getRealTargetCoordinates(icellT,coordsT); + for ( t = 0; t < nbSegT; ++t ) + for(typename std::vector::const_iterator + iter=icellsS.begin(); iter!=icellsS.end(); iter++) + { + int iS = *iter; + BASE_INTERSECTOR::getDualSegments( OTT::ind2C(iS), + BASE_INTERSECTOR::_meshS, segmentsS); + for ( int s = 0; s < segmentsS.size(); ++s ) + { + double surf = BASE_INTERSECTOR::intersectSegments(segmentsS[s]._coords.data(), + coordsT.data() + t*SPACEDIM); + if(surf!=0.) + { + int nS = segmentsS[s]._nodeId; + typename MyMatrix::value_type::const_iterator iterRes=resRow.find(nS); + if(iterRes==resRow.end()) + resRow.insert(std::make_pair(nS,surf)); + else + { + surf+=(*iterRes).second; + resRow.erase(nS); + resRow.insert(std::make_pair(nS,surf)); + } + } + } + } + } +} +#undef BASE_INTERSECTOR + +#endif diff --git a/src/INTERP_KERNEL/CurveIntersectorP1P1.hxx b/src/INTERP_KERNEL/CurveIntersectorP1P1.hxx new file mode 100644 index 000000000..067e51290 --- /dev/null +++ b/src/INTERP_KERNEL/CurveIntersectorP1P1.hxx @@ -0,0 +1,46 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __CURVEINTERSECTORP1P1_HXX__ +#define __CURVEINTERSECTORP1P1_HXX__ + +#include "CurveIntersector.hxx" + +namespace INTERP_KERNEL +{ + template + class CurveIntersectorP1P1 : 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; + + CurveIntersectorP1P1(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; + }; +} + +#endif diff --git a/src/INTERP_KERNEL/CurveIntersectorP1P1.txx b/src/INTERP_KERNEL/CurveIntersectorP1P1.txx new file mode 100644 index 000000000..fa39b2e5f --- /dev/null +++ b/src/INTERP_KERNEL/CurveIntersectorP1P1.txx @@ -0,0 +1,90 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +#ifndef __CurveIntersectorP1P1_TXX__ +#define __CurveIntersectorP1P1_TXX__ + +#include "CurveIntersectorP1P1.hxx" +#include "CurveIntersector.txx" + +#define BASE_INTERSECTOR CurveIntersector + +namespace INTERP_KERNEL +{ + template + CurveIntersectorP1P1 + ::CurveIntersectorP1P1(const MyMeshType& meshT, const MyMeshType& meshS, + double precision, double tolerance, + double medianLine, int printLevel): + BASE_INTERSECTOR (meshT, meshS, precision, tolerance, medianLine, printLevel) + { + } + + template + int CurveIntersectorP1P1 + ::getNumberOfRowsOfResMatrix() const + { + return BASE_INTERSECTOR::_meshT.getNumberOfNodes(); + } + + template + int CurveIntersectorP1P1 + ::getNumberOfColsOfResMatrix() const + { + return BASE_INTERSECTOR::_meshS.getNumberOfNodes(); + } + + template + void CurveIntersectorP1P1 + ::intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res) + { + std::vector segmentsT, segmentsS; + BASE_INTERSECTOR::getDualSegments( icellT, BASE_INTERSECTOR::_meshT, segmentsT); + for ( int t = 0; t < segmentsT.size(); ++t ) + { + typename MyMatrix::value_type& resRow = res[ OTT::ind2C( segmentsT[t]._nodeId )]; + for(typename std::vector::const_iterator + iter=icellsS.begin(); iter!=icellsS.end(); iter++) + { + int iS = *iter; + BASE_INTERSECTOR::getDualSegments( OTT::ind2C(iS), + BASE_INTERSECTOR::_meshS, segmentsS); + for ( int s = 0; s < segmentsS.size(); ++s ) + { + double surf = BASE_INTERSECTOR::intersectSegments(segmentsT[t]._coords.data(), + segmentsS[s]._coords.data()); + if(surf!=0.) + { + int nS = segmentsS[s]._nodeId; + typename MyMatrix::value_type::const_iterator iterRes=resRow.find(nS); + if(iterRes==resRow.end()) + resRow.insert(std::make_pair(nS,surf)); + else + { + surf+=(*iterRes).second; + resRow.erase(nS); + resRow.insert(std::make_pair(nS,surf)); + } + } + } + } + } + } +} +#undef BASE_INTERSECTOR + +#endif diff --git a/src/INTERP_KERNEL/Interpolation1D.hxx b/src/INTERP_KERNEL/Interpolation1D.hxx new file mode 100755 index 000000000..f07a22cc9 --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation1D.hxx @@ -0,0 +1,34 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __INTERPOLATION1D_HXX__ +#define __INTERPOLATION1D_HXX__ + +#include "InterpolationCurve.hxx" + +namespace INTERP_KERNEL +{ + class Interpolation1D : public InterpolationCurve + { + public: + Interpolation1D() { } + Interpolation1D(const InterpolationOptions& io):InterpolationCurve(io) {} + }; +} + +#endif diff --git a/src/INTERP_KERNEL/Interpolation1D.txx b/src/INTERP_KERNEL/Interpolation1D.txx new file mode 100644 index 000000000..b25b77c4c --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation1D.txx @@ -0,0 +1,26 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __INTERPOLATION1D_TXX__ +#define __INTERPOLATION1D_TXX__ + +#include "Interpolation1D.hxx" + +#include "InterpolationCurve.txx" + +#endif diff --git a/src/INTERP_KERNEL/Interpolation2DCurve.hxx b/src/INTERP_KERNEL/Interpolation2DCurve.hxx new file mode 100644 index 000000000..bff59b26b --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation2DCurve.hxx @@ -0,0 +1,37 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __INTERPOLATION2DCURVE_HXX__ +#define __INTERPOLATION2DCURVE_HXX__ + +#include "InterpolationCurve.hxx" +#include "InterpolationOptions.hxx" + +namespace INTERP_KERNEL +{ + class Interpolation2DCurve : public InterpolationCurve + { + public: + Interpolation2DCurve(); + Interpolation2DCurve(const InterpolationOptions& io); + // geometric precision, intersection tolerance, coice of the median line, + void setOptions(double precision, double tolerance, double medianLine); + }; +} + +#endif diff --git a/src/INTERP_KERNEL/Interpolation2DCurve.txx b/src/INTERP_KERNEL/Interpolation2DCurve.txx new file mode 100644 index 000000000..04bd48728 --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation2DCurve.txx @@ -0,0 +1,63 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __INTERPOLATION2DCURVE_TXX__ +#define __INTERPOLATION2DCURVE_TXX__ + +#include "Interpolation2DCurve.hxx" +#include "InterpolationCurve.txx" + +namespace INTERP_KERNEL +{ + Interpolation2DCurve::Interpolation2DCurve() + { + // to have non-zero default thickness of target element + InterpolationOptions::setBoundingBoxAdjustmentAbs( InterpolationOptions::getPrecision() ); + } + + Interpolation2DCurve::Interpolation2DCurve + (const InterpolationOptions& io):InterpolationCurve(io) + { + // to have non-zero default thickness of target element + InterpolationOptions::setBoundingBoxAdjustmentAbs( InterpolationOptions::getPrecision() ); + } + + /** + * \brief Function used to set the options for the intersection calculation + * \details The following options can be modified: + * -# Precision: Level of precision of the computations. + * - Values: positive real number. + * - Default: 1.0E-12. + * -# Tolerance: Thickness of target element. + * - Values: positive real number. + * - Default: 1.0E-12. + * -# Median line: Position of the median line where both segments will be projected. + * - Values: real number between 0.0 and 1.0. + * - Default: 0.5 + */ + void Interpolation2DCurve::setOptions (double precision, + double tolerance, + double medianLine) + { + InterpolationOptions::setPrecision(precision); + InterpolationOptions::setBoundingBoxAdjustmentAbs(tolerance); + InterpolationOptions::setMedianPlane(medianLine); + } +} + +#endif diff --git a/src/INTERP_KERNEL/InterpolationCurve.hxx b/src/INTERP_KERNEL/InterpolationCurve.hxx new file mode 100644 index 000000000..fe3706e93 --- /dev/null +++ b/src/INTERP_KERNEL/InterpolationCurve.hxx @@ -0,0 +1,42 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __INTERPOLATIONCURVE_HXX__ +#define __INTERPOLATIONCURVE_HXX__ + +#include "Interpolation.hxx" +#include "InterpolationOptions.hxx" + +namespace INTERP_KERNEL +{ + template + class InterpolationCurve : public Interpolation< InterpolationCurve > + { + public: + InterpolationCurve(); + InterpolationCurve(const InterpolationOptions & io); + + // Main function to interpolate + template + int interpolateMeshes(const MyMeshType& mesh1, const MyMeshType& mesh2, + MatrixType& result, const char *method); + + }; +} + +#endif diff --git a/src/INTERP_KERNEL/InterpolationCurve.txx b/src/INTERP_KERNEL/InterpolationCurve.txx new file mode 100644 index 000000000..d7181ce61 --- /dev/null +++ b/src/INTERP_KERNEL/InterpolationCurve.txx @@ -0,0 +1,171 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __INTERPOLATIONCURVE_TXX__ +#define __INTERPOLATIONCURVE_TXX__ + +#include "InterpolationCurve.hxx" +#include "InterpolationOptions.hxx" +#include "CurveIntersectorP0P0.txx" +#include "CurveIntersectorP1P0.txx" +#include "CurveIntersectorP0P1.txx" +#include "CurveIntersectorP1P1.txx" +#include "BBTree.txx" + +#include + +namespace INTERP_KERNEL +{ + /** + * \defgroup interpolationCurve InterpolationCurve + * + * \class InterpolationCurve + * \brief Class used to compute the coefficients of the interpolation matrix between + * two local meshes in two dimensions. + */ + template + InterpolationCurve::InterpolationCurve() + { + } + + template + InterpolationCurve::InterpolationCurve (const InterpolationOptions& io) + :Interpolation< InterpolationCurve >(io) + { + } + + /** \brief Main function to interpolate 1D meshes. + \details The algorithm proceeds in two steps: first a filtering process reduces the number of pairs of elements for which the + * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the + * volume of intersection is calculated by an object of type IntersectorPlanar for the remaining pairs, and entered into the + * intersection matrix. + * + * The matrix is partially sparse : it is a vector of maps of integer - double pairs. + * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless + * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements + * which have a non-zero intersection volume with the target element. The vector has indices running from + * 0 to (#target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however, + * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j]. + * + + * @param myMeshS Planar source mesh + * @Param myMeshT Planar target mesh + * @return vector containing for each element i of the source mesh, a map giving for each element j + * of the target mesh which i intersects, the area of the intersection + * + */ + template + template + int InterpolationCurve::interpolateMeshes (const MyMeshType& myMeshS, + const MyMeshType& myMeshT, + MatrixType& result, + const char * method) + { + static const int SPACEDIM=MyMeshType::MY_SPACEDIM; + typedef typename MyMeshType::MyConnType ConnType; + static const NumberingPolicy numPol = MyMeshType::My_numPol; + + long global_start = clock(); + int counter=0; + + long nbMailleS = myMeshS.getNumberOfElements(); + long nbMailleT = myMeshT.getNumberOfElements(); + + CurveIntersector* intersector=0; + std::string meth(method); + if(meth=="P0P0") { + intersector = new CurveIntersectorP0P0 + (myMeshT, myMeshS, + InterpolationOptions::getPrecision(), + InterpolationOptions::getBoundingBoxAdjustmentAbs(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getPrintLevel()); + } + else if(meth=="P0P1") { + intersector = new CurveIntersectorP0P1 + (myMeshT, myMeshS, + InterpolationOptions::getPrecision(), + InterpolationOptions::getBoundingBoxAdjustmentAbs(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getPrintLevel()); + } + else if(meth=="P1P0") { + intersector = new CurveIntersectorP1P0 + (myMeshT, myMeshS, + InterpolationOptions::getPrecision(), + InterpolationOptions::getBoundingBoxAdjustmentAbs(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getPrintLevel()); + } + else if(meth=="P1P1") { + intersector = new CurveIntersectorP1P1 + (myMeshT, myMeshS, + InterpolationOptions::getPrecision(), + InterpolationOptions::getBoundingBoxAdjustmentAbs(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getPrintLevel()); + } + else + throw INTERP_KERNEL::Exception("Invalid method specified ! Must be in : \"P0P0\" \"P0P1\" \"P1P0\" or \"P1P1\""); + /****************************************************************/ + /* Create a search tree based on the bounding boxes */ + /* Instanciate the intersector and initialise the result vector */ + /****************************************************************/ + + long start_filtering=clock(); + + std::vector bbox; + intersector->createBoundingBoxes(myMeshS,bbox); // create the bounding boxes + intersector->adjustBoundingBoxes(bbox, InterpolationOptions::getBoundingBoxAdjustmentAbs()); + BBTree my_tree(&bbox[0], 0, 0, nbMailleS);//creating the search structure + + long end_filtering = clock(); + + result.resize(intersector->getNumberOfRowsOfResMatrix());//on initialise. + + /****************************************************/ + /* Loop on the target cells - core of the algorithm */ + /****************************************************/ + long start_intersection = clock(); + const ConnType *connIndxT = myMeshT.getConnectivityIndexPtr(); + for(int iT=0; iT intersecting_elems; + double bb[2*SPACEDIM]; + intersector->getElemBB(bb,myMeshT,OTT::indFC(iT),nb_nodesT); + my_tree.getIntersectingElems(bb, intersecting_elems); + intersector->intersectCells(iT,intersecting_elems,result); + counter += intersecting_elems.size(); + } + int ret = intersector->getNumberOfColsOfResMatrix(); + delete intersector; + + if (InterpolationOptions::getPrintLevel() >= 1) { + long end_intersection=clock(); + std::cout << "Filtering time= " << end_filtering-start_filtering << std::endl; + std::cout << "Intersection time= " << end_intersection-start_intersection << std::endl; + long global_end =clock(); + std::cout << "Number of computed intersections = " << counter << std::endl; + std::cout << "Global time= " << global_end - global_start << std::endl; + } + return ret; + } +} + +#endif diff --git a/src/INTERP_KERNEL/Makefile.am b/src/INTERP_KERNEL/Makefile.am index fd748819a..b2f2bbb0b 100644 --- a/src/INTERP_KERNEL/Makefile.am +++ b/src/INTERP_KERNEL/Makefile.am @@ -141,7 +141,20 @@ VTKNormalizedUnstructuredMesh.txx \ VectorUtils.hxx \ VolSurfFormulae.hxx \ VolSurfUser.hxx \ -VolSurfUser.txx +VolSurfUser.txx \ +CurveIntersector.hxx \ +CurveIntersector.txx \ +CurveIntersectorP0P0.hxx \ +CurveIntersectorP0P0.txx \ +CurveIntersectorP0P1.hxx \ +CurveIntersectorP1P0.hxx \ +CurveIntersectorP1P1.hxx \ +Interpolation1D.hxx \ +Interpolation1D.txx \ +Interpolation2DCurve.hxx \ +Interpolation2DCurve.txx \ +InterpolationCurve.hxx \ +InterpolationCurve.txx EXTRA_DIST += \