From: bph Date: Wed, 5 Oct 2011 08:55:31 +0000 (+0000) Subject: Prestation INTERP_KERNEL D. POIZAT X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=c36c2f327fdeadc4031f9a88c453a1146ce9e018;p=tools%2Fmedcoupling.git Prestation INTERP_KERNEL D. POIZAT --- diff --git a/src/INTERP_KERNEL/Interpolation2D1D.hxx b/src/INTERP_KERNEL/Interpolation2D1D.hxx new file mode 100644 index 000000000..522322cb1 --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation2D1D.hxx @@ -0,0 +1,53 @@ +// Copyright (C) 2007-2011 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 __INTERPOLATION2D1D_HXX__ +#define __INTERPOLATION2D1D_HXX__ + +#include "Interpolation.hxx" +#include "Planar2D1DIntersectorP0P0.hxx" +#include "NormalizedUnstructuredMesh.hxx" +#include "InterpolationOptions.hxx" + +namespace INTERP_KERNEL +{ + class Interpolation2D1D : public Interpolation + { + public: + typedef std::map > DuplicateFacesType; + + Interpolation2D1D() { setOrientation(2); } + Interpolation2D1D(const InterpolationOptions& io):Interpolation(io) { } + public: + + // Main function to interpolate triangular and quadratic meshes + template + int interpolateMeshes(const MyMeshType& meshS, const MyMeshType& meshT, MatrixType& result, const char *method); + DuplicateFacesType retrieveDuplicateFaces() const + { + return _duplicate_faces; + } + private: + DuplicateFacesType _duplicate_faces; + private: + double _dim_caracteristic; + }; +} + +#endif diff --git a/src/INTERP_KERNEL/Interpolation2D1D.txx b/src/INTERP_KERNEL/Interpolation2D1D.txx new file mode 100644 index 000000000..e8b682eef --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation2D1D.txx @@ -0,0 +1,162 @@ +// Copyright (C) 2007-2011 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 __INTERPOLATION2D1D_TXX__ +#define __INTERPOLATION2D1D_TXX__ + +#include "Interpolation2D1D.hxx" + +namespace INTERP_KERNEL +{ + /** + * \defgroup interpolation2D1D Interpolation2D1D + * \class Interpolation2D1D + * \brief Class used to calculate the segments of intersection between the elements of a 2D linear source mesh + * and a 2D surface target mesh. + * + */ + + /** \brief Main function to interpolate triangular or quadrangular 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 + int Interpolation2D1D::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; + /***********************************************************/ + /* Check both meshes are made of triangles and quadrangles */ + /***********************************************************/ + + long nbMailleS=myMeshS.getNumberOfElements(); + + /**************************************************/ + /* Search the characteristic size of the meshes */ + /**************************************************/ + + int printLevel = InterpolationOptions::getPrintLevel(); + _dim_caracteristic = CalculateCharacteristicSizeOfMeshes(myMeshS, myMeshT, printLevel); + if (printLevel>=1) + { + std::cout << "Interpolation2D1D::computation of the intersections" << std::endl; + } + + PlanarIntersector* intersector=0; + std::string meth = InterpolationOptions::filterInterpolationMethod(method); + if(meth=="P0P0") + { + switch (InterpolationOptions::getIntersectionType()) + { + case Geometric2D: + intersector=new Geometric2DIntersector(myMeshT, myMeshS, _dim_caracteristic, + InterpolationOptions::getMaxDistance3DSurfIntersect(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getPrecision(), + InterpolationOptions::getOrientation()); + break; + default: + throw INTERP_KERNEL::Exception("Invalid intersection type ! Must be : Geometric2D"); + } + } + else + throw INTERP_KERNEL::Exception("Invalid method specified or intersection type ! Must be : \"P0P0\""); + + /****************************************************************/ + /* 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 + const double *bboxPtr=0; + if(nbMailleS>0) + bboxPtr=&bbox[0]; + BBTree my_tree(bboxPtr, 0, 0,nbMailleS, -getPrecision());//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(); + long nbelem_type=myMeshT.getNumberOfElements(); + 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(); + intersecting_elems.clear(); + } + int ret=intersector->getNumberOfColsOfResMatrix(); + + const DuplicateFacesType& intersectFaces = *intersector->getIntersectFaces(); + DuplicateFacesType::const_iterator iter; + for (iter = intersectFaces.begin(); iter != intersectFaces.end(); ++iter) + { + if (iter->second.size() > 1) + { + _duplicate_faces.insert(std::make_pair(iter->first, iter->second)); + } + } + + 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/Planar2D1DIntersectorP0P0.hxx b/src/INTERP_KERNEL/Planar2D1DIntersectorP0P0.hxx new file mode 100644 index 000000000..3a79ad04c --- /dev/null +++ b/src/INTERP_KERNEL/Planar2D1DIntersectorP0P0.hxx @@ -0,0 +1,59 @@ +// Copyright (C) 2007-2011 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 __PLANAR2D1DINTERSECTORP0P0_HXX__ +#define __PLANAR2D1DINTERSECTORP0P0_HXX__ + +#include "PlanarIntersector.hxx" + +namespace INTERP_KERNEL +{ + template + class Planar2D1DIntersectorP0P0 : public PlanarIntersector + { + 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; + protected: + Planar2D1DIntersectorP0P0(const MyMeshType& meshT, const MyMeshType& meshS, + double dimCaracteristic, double precision, double md3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel); + public: + int getNumberOfRowsOfResMatrix() const; + int getNumberOfColsOfResMatrix() const; + const typename PlanarIntersector::DuplicateFacesType* getIntersectFaces() const + { + return &_intersect_faces; + } + void intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res); + /*! + * Contrary to intersectCells method here icellS and icellT are \b not in \b C mode but in mode of MyMeshType. + */ + double intersectGeometry1D(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS, + bool& isColinear) + { return asLeaf().intersectGeometry1D(icellT,icellS,nbNodesT,nbNodesS, isColinear); } + protected: + ConcreteP0P0Intersector& asLeaf() { return static_cast(*this); } + private: + typename PlanarIntersector::DuplicateFacesType _intersect_faces; + }; +} + +#endif diff --git a/src/INTERP_KERNEL/Planar2D1DIntersectorP0P0.txx b/src/INTERP_KERNEL/Planar2D1DIntersectorP0P0.txx new file mode 100644 index 000000000..d0cd94647 --- /dev/null +++ b/src/INTERP_KERNEL/Planar2D1DIntersectorP0P0.txx @@ -0,0 +1,80 @@ +// Copyright (C) 2007-2011 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 __PLANAR2D1DINTERSECTORP0P0_TXX__ +#define __PLANAR2D1DINTERSECTORP0P0_TXX__ + +#include "Planar2D1DIntersectorP0P0.hxx" + +namespace INTERP_KERNEL +{ + template + Planar2D1DIntersectorP0P0::Planar2D1DIntersectorP0P0(const MyMeshType& meshT, const MyMeshType& meshS, + double dimCaracteristic, double precision, double md3DSurf, double medianPlane, + bool doRotate, int orientation, int printLevel): + PlanarIntersector(meshT,meshS,dimCaracteristic,precision,md3DSurf,medianPlane,doRotate,orientation,printLevel) + { + } + + template + int Planar2D1DIntersectorP0P0::getNumberOfRowsOfResMatrix() const + { + return PlanarIntersector::_meshT.getNumberOfElements(); + } + + template + int Planar2D1DIntersectorP0P0::getNumberOfColsOfResMatrix() const + { + return PlanarIntersector::_meshS.getNumberOfElements(); + } + + template + void Planar2D1DIntersectorP0P0::intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res) + { + int nbNodesT=PlanarIntersector::_connIndexT[icellT+1]-PlanarIntersector::_connIndexT[icellT]; + typename MyMatrix::value_type& resRow=res[icellT]; + for(typename std::vector::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++) + { + int iS=*iter; + int nbNodesS=PlanarIntersector::_connIndexS[iS+1]-PlanarIntersector::_connIndexS[iS]; + bool isColinear = false; + double surf=intersectGeometry1D(OTT::indFC(icellT),OTT::indFC(iS), + nbNodesT,nbNodesS, isColinear); + surf=PlanarIntersector::getValueRegardingOption(surf); + if(surf!=0.) + resRow.insert(std::make_pair(OTT::indFC(iS),surf)); + + if (isColinear) + { + typename PlanarIntersector::DuplicateFacesType::iterator intersectFacesIter = _intersect_faces.find(iS); + if (intersectFacesIter != _intersect_faces.end()) + { + intersectFacesIter->second.insert(icellT); + } + else + { + std::set targetCellSet; + targetCellSet.insert(icellT); + _intersect_faces.insert(std::make_pair(iS, targetCellSet)); + } + } + } + } +} + +#endif