From: ageay Date: Fri, 26 Aug 2011 16:21:16 +0000 (+0000) Subject: Cross dimension intersectors. X-Git-Tag: V6_main_FINAL~969 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=793c5142dc1bcf90d00c8f83ea0ee4a83df14397;p=tools%2Fmedcoupling.git Cross dimension intersectors. --- 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..62e921708 --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation2D1D.txx @@ -0,0 +1,155 @@ +// 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 +{ + + /** \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/Interpolation3D2D.cxx b/src/INTERP_KERNEL/Interpolation3D2D.cxx new file mode 100644 index 000000000..d29392342 --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation3D2D.cxx @@ -0,0 +1,40 @@ +// 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 +// + +#include "Interpolation3D2D.hxx" + +namespace INTERP_KERNEL +{ + /** + * \defgroup interpolation3D Interpolation3D + * \class Interpolation3D + * \brief Class used to calculate the volumes of intersection between the elements of two 3D meshes. + * + */ + /** + * Default constructor + * + */ + Interpolation3D2D::Interpolation3D2D() + { + } + Interpolation3D2D::Interpolation3D2D(const InterpolationOptions& io):Interpolation(io) + { + } +} diff --git a/src/INTERP_KERNEL/Interpolation3D2D.hxx b/src/INTERP_KERNEL/Interpolation3D2D.hxx new file mode 100644 index 000000000..f68a3b672 --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation3D2D.hxx @@ -0,0 +1,54 @@ +// 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 __INTERPOLATION3D2D_HXX__ +#define __INTERPOLATION3D2D_HXX__ + +#include +#include + +#include "Interpolation.hxx" +#include "NormalizedUnstructuredMesh.hxx" +#include "InterpolationOptions.hxx" + +namespace INTERP_KERNEL +{ + class Interpolation3D2D : public Interpolation + { + public: + typedef std::map > DuplicateFacesType; + + Interpolation3D2D(); + Interpolation3D2D(const InterpolationOptions& io); + template + int interpolateMeshes(const MyMeshType& srcMesh, + const MyMeshType& targetMesh, + MyMatrixType& matrix, + const char *method); + DuplicateFacesType retrieveDuplicateFaces() const + { + return _duplicate_faces; + } + private: + SplittingPolicy _splitting_policy; + DuplicateFacesType _duplicate_faces; + }; +} + +#endif diff --git a/src/INTERP_KERNEL/Interpolation3D2D.txx b/src/INTERP_KERNEL/Interpolation3D2D.txx new file mode 100644 index 000000000..27ddb7e8f --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation3D2D.txx @@ -0,0 +1,187 @@ +// 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 __INTERPOLATION3D2D_TXX__ +#define __INTERPOLATION3D2D_TXX__ + +#include "Interpolation3D2D.hxx" +#include "Interpolation.txx" +#include "MeshElement.txx" +#include "TransformedTriangle.hxx" +#include "Polyhedron3D2DIntersectorP0P0.txx" +#include "PointLocator3DIntersectorP0P0.txx" +#include "PolyhedronIntersectorP0P1.txx" +#include "PointLocator3DIntersectorP0P1.txx" +#include "PolyhedronIntersectorP1P0.txx" +#include "PolyhedronIntersectorP1P0Bary.txx" +#include "PointLocator3DIntersectorP1P0.txx" +#include "PolyhedronIntersectorP1P1.txx" +#include "PointLocator3DIntersectorP1P1.txx" +#include "Log.hxx" + +#include "BBTree.txx" + +namespace INTERP_KERNEL +{ + /** + * Calculates the matrix of volumes of intersection between the elements of srcMesh and the elements of targetMesh. + * The calculation is done 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 Intersector3D 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. + * It can also be an INTERP_KERNEL::Matrix object. + * 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 (nb 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 srcMesh 3-dimensional source mesh + * @param targetMesh 3-dimesional target mesh, containing only tetraedra + * @param matrix matrix in which the result is stored + * + */ + template + int Interpolation3D2D::interpolateMeshes(const MyMeshType& srcMesh, + const MyMeshType& targetMesh, + MyMatrixType& matrix, + const char *method) + { + typedef typename MyMeshType::MyConnType ConnType; + // create MeshElement objects corresponding to each element of the two meshes + const unsigned long numSrcElems = srcMesh.getNumberOfElements(); + const unsigned long numTargetElems = targetMesh.getNumberOfElements(); + + LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements "); + + std::vector*> srcElems(numSrcElems); + std::vector*> targetElems(numTargetElems); + + std::map*, int> indices; + DuplicateFacesType intersectFaces; + + for(unsigned long i = 0 ; i < numSrcElems ; ++i) + srcElems[i] = new MeshElement(i, srcMesh); + + for(unsigned long i = 0 ; i < numTargetElems ; ++i) + targetElems[i] = new MeshElement(i, targetMesh); + + Intersector3D* intersector=0; + std::string methC = InterpolationOptions::filterInterpolationMethod(method); + const double dimCaracteristic = CalculateCharacteristicSizeOfMeshes(srcMesh, targetMesh, InterpolationOptions::getPrintLevel()); + if(methC=="P0P0") + { + switch(InterpolationOptions::getIntersectionType()) + { + case Triangulation: + intersector=new Polyhedron3D2DIntersectorP0P0(targetMesh, + srcMesh, + dimCaracteristic, + getPrecision(), + intersectFaces, + getSplittingPolicy()); + break; + default: + throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P0P0 interp specified : must be Triangulation."); + } + } + else + throw Exception("Invalid method choosed must be in \"P0P0\"."); + // create empty maps for all source elements + matrix.resize(intersector->getNumberOfRowsOfResMatrix()); + + // create BBTree structure + // - get bounding boxes + double* bboxes = new double[6 * numSrcElems]; + int* srcElemIdx = new int[numSrcElems]; + for(unsigned long i = 0; i < numSrcElems ; ++i) + { + // get source bboxes in right order + const BoundingBox* box = srcElems[i]->getBoundingBox(); + bboxes[6*i+0] = box->getCoordinate(BoundingBox::XMIN); + bboxes[6*i+1] = box->getCoordinate(BoundingBox::XMAX); + bboxes[6*i+2] = box->getCoordinate(BoundingBox::YMIN); + bboxes[6*i+3] = box->getCoordinate(BoundingBox::YMAX); + bboxes[6*i+4] = box->getCoordinate(BoundingBox::ZMIN); + bboxes[6*i+5] = box->getCoordinate(BoundingBox::ZMAX); + + // source indices have to begin with zero for BBox, I think + srcElemIdx[i] = srcElems[i]->getIndex(); + } + + BBTree<3,ConnType> tree(bboxes, srcElemIdx, 0, numSrcElems, 0.); + + // for each target element, get source elements with which to calculate intersection + // - calculate intersection by calling intersectCells + for(unsigned long i = 0; i < numTargetElems; ++i) + { + const BoundingBox* box = targetElems[i]->getBoundingBox(); + const int targetIdx = targetElems[i]->getIndex(); + + // get target bbox in right order + double targetBox[6]; + targetBox[0] = box->getCoordinate(BoundingBox::XMIN); + targetBox[1] = box->getCoordinate(BoundingBox::XMAX); + targetBox[2] = box->getCoordinate(BoundingBox::YMIN); + targetBox[3] = box->getCoordinate(BoundingBox::YMAX); + targetBox[4] = box->getCoordinate(BoundingBox::ZMIN); + targetBox[5] = box->getCoordinate(BoundingBox::ZMAX); + + std::vector intersectElems; + + tree.getIntersectingElems(targetBox, intersectElems); + + if ( !intersectElems.empty() ) + intersector->intersectCells(targetIdx, intersectElems, matrix); + + } + + delete [] bboxes; + delete [] srcElemIdx; + + DuplicateFacesType::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)); + } + } + + // free allocated memory + int ret=intersector->getNumberOfColsOfResMatrix(); + + delete intersector; + + for(unsigned long i = 0 ; i < numSrcElems ; ++i) + { + delete srcElems[i]; + } + for(unsigned long i = 0 ; i < numTargetElems ; ++i) + { + delete targetElems[i]; + } + 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 diff --git a/src/INTERP_KERNEL/PlanarIntersector.hxx b/src/INTERP_KERNEL/PlanarIntersector.hxx index 12a3455ad..395ff8b39 100644 --- a/src/INTERP_KERNEL/PlanarIntersector.hxx +++ b/src/INTERP_KERNEL/PlanarIntersector.hxx @@ -23,6 +23,9 @@ #include "TargetIntersector.hxx" #include "NormalizedUnstructuredMesh.hxx" +#include +#include + namespace INTERP_KERNEL { class TranslationRotationMatrix; @@ -35,6 +38,7 @@ namespace INTERP_KERNEL static const int MESHDIM=MyMeshType::MY_MESHDIM; typedef typename MyMeshType::MyConnType ConnType; static const NumberingPolicy numPol=MyMeshType::My_numPol; + typedef typename std::map > DuplicateFacesType; public: //! \addtogroup InterpKerGrpIntPlan @{ PlanarIntersector(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel); @@ -45,6 +49,10 @@ namespace INTERP_KERNEL inline void getElemBB(double* bb, const MyMeshType& mesh, ConnType iP, ConnType nb_nodes); static int projection(double *Coords_A, double *Coords_B, int nb_NodesA, int nb_NodesB, double epsilon, double md3DSurf, double median_plane, bool do_rotate); + virtual const DuplicateFacesType* getIntersectFaces() const + { + return NULL; + } protected : int projectionThis(double *Coords_A, double *Coords_B, int nb_NodesA, int nb_NodesB); void getRealTargetCoordinates(ConnType icellT, std::vector& coordsT); diff --git a/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.hxx b/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.hxx new file mode 100644 index 000000000..c94c51cb0 --- /dev/null +++ b/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.hxx @@ -0,0 +1,79 @@ +// 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 __POLYHEDRON3D2DINTERSECTORP0P0_HXX__ +#define __POLYHEDRON3D2DINTERSECTORP0P0_HXX__ + +#include "Intersector3DP0P0.hxx" +#include "SplitterTetra.hxx" +#include "NormalizedUnstructuredMesh.hxx" + +namespace INTERP_KERNEL +{ + + + /** + * \brief Class responsible for calculating intersection between a hexahedron target element and + * the source elements. + * + */ + template + class Polyhedron3D2DIntersectorP0P0 : public Intersector3DP0P0 + { + typedef typename std::map > DuplicateFacesType; + + 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: + + Polyhedron3D2DIntersectorP0P0(const MyMeshType& targetMesh, + const MyMeshType& srcMesh, + const double dimCaracteristic, + const double precision, + DuplicateFacesType& intersectFaces, + SplittingPolicy policy = PLANAR_FACE_5); + + ~Polyhedron3D2DIntersectorP0P0(); + + void intersectCells(ConnType targetCell, + const std::vector& srcCells, + MyMatrixType& matrix); + + private: + void releaseArrays(); + private: + /// pointers to the SplitterTetra objects representing the tetrahedra + /// that result from the splitting of the hexahedron target cell + std::vector< SplitterTetra* > _tetra; + + SplitterTetra2 _split; + + double _dim_caracteristic; + double _precision; + + DuplicateFacesType& _intersect_faces; + + }; +} + +#endif diff --git a/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx b/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx new file mode 100644 index 000000000..12a17bd45 --- /dev/null +++ b/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx @@ -0,0 +1,174 @@ +// 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 __POLYHEDRON3D2DINTERSECTORP0P0_TXX__ +#define __POLYHEDRON3D2DINTERSECTORP0P0_TXX__ + +#include "Polyhedron3D2DIntersectorP0P0.hxx" +#include "Intersector3DP0P0.txx" +#include "MeshUtils.hxx" + +#include "SplitterTetra.txx" + +namespace INTERP_KERNEL +{ + + /** + * Constructor creating object from target cell global number + * The constructor first calculates the necessary nodes, + * (depending on the splitting policy) and then splits the hexahedron into + * tetrahedra, placing these in the internal vector _tetra. + * + * @param targetMesh mesh containing the target elements + * @param srcMesh mesh containing the source elements + * @param policy splitting policy to be used + */ + template + Polyhedron3D2DIntersectorP0P0::Polyhedron3D2DIntersectorP0P0(const MyMeshType& targetMesh, + const MyMeshType& srcMesh, + const double dimCaracteristic, + const double precision, + DuplicateFacesType& intersectFaces, + SplittingPolicy policy) + : Intersector3DP0P0(targetMesh,srcMesh), + _split(targetMesh,srcMesh,policy), + _dim_caracteristic(dimCaracteristic), + _precision(precision), + _intersect_faces(intersectFaces) + { + } + + /** + * Destructor. + * Liberates the SplitterTetra objects and potential sub-node points that have been allocated. + * + */ + template + Polyhedron3D2DIntersectorP0P0::~Polyhedron3D2DIntersectorP0P0() + { + releaseArrays(); + } + + template + void Polyhedron3D2DIntersectorP0P0::releaseArrays() + { + for(typename std::vector< SplitterTetra* >::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter) + delete *iter; + _split.releaseArrays(); + _tetra.clear(); + } + + /** + * Calculates the volume of intersection of an element in the source mesh and the target element + * represented by the object. + * The calculation is performed by calling the corresponding method for + * each SplitterTetra object created by the splitting. + * + * @param targetCell in C mode. + * @param srcCells in C mode. + * + */ + template + void Polyhedron3D2DIntersectorP0P0::intersectCells(ConnType targetCell, + const std::vector& srcCells, + MyMatrixType& matrix) + { + int nbOfNodesT=Intersector3D::_target_mesh.getNumberOfNodesOfElement(OTT::indFC(targetCell)); + releaseArrays(); + _split.splitTargetCell(targetCell,nbOfNodesT,_tetra); + + for(typename std::vector::const_iterator iterCellS=srcCells.begin();iterCellS!=srcCells.end();iterCellS++) + { + double surface = 0.; + std::multiset listOfTetraFacesTreated; + std::set listOfTetraFacesColinear; + + // calculate the coordinates of the nodes + const NumberingPolicy numPol=MyMeshType::My_numPol; + typename MyMeshType::MyConnType cellSrc = *iterCellS; + int cellSrcIdx = OTT::indFC(cellSrc); + NormalizedCellType normCellType=Intersector3D::_src_mesh.getTypeOfElement(cellSrcIdx); + const CellModel& cellModelCell=CellModel::GetCellModel(normCellType); + const MyMeshType& _src_mesh = Intersector3D::_src_mesh; + unsigned nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(cellSrcIdx) : cellModelCell.getNumberOfNodes(); + int *polyNodes=new int[nbOfNodes4Type]; + double **polyCoords = new double*[nbOfNodes4Type]; + for(int i = 0;i<(int)nbOfNodes4Type;++i) + { + // we could store mapping local -> global numbers too, but not sure it is worth it + const int globalNodeNum = getGlobalNumberOfNode(i, OTT::indFC(*iterCellS), _src_mesh); + polyNodes[i]=globalNodeNum; + polyCoords[i] = (double*)_src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum; + } + + for(typename std::vector*>::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter) + surface += (*iter)->intersectSourceFace(normCellType, + nbOfNodes4Type, + polyNodes, + polyCoords, + _dim_caracteristic, + _precision, + listOfTetraFacesTreated, + listOfTetraFacesColinear); + + if(surface!=0.) { + + matrix[targetCell].insert(std::make_pair(cellSrcIdx, surface)); + + bool isSrcFaceColinearWithFaceOfTetraTargetCell = false; + std::set::iterator iter; + for (iter = listOfTetraFacesColinear.begin(); iter != listOfTetraFacesColinear.end(); ++iter) + { + if (listOfTetraFacesTreated.count(*iter) != 1) + { + isSrcFaceColinearWithFaceOfTetraTargetCell = false; + break; + } + else + { + isSrcFaceColinearWithFaceOfTetraTargetCell = true; + } + } + + if (isSrcFaceColinearWithFaceOfTetraTargetCell) + { + DuplicateFacesType::iterator intersectFacesIter = _intersect_faces.find(cellSrcIdx); + if (intersectFacesIter != _intersect_faces.end()) + { + intersectFacesIter->second.insert(targetCell); + } + else + { + std::set targetCellSet; + targetCellSet.insert(targetCell); + _intersect_faces.insert(std::make_pair(cellSrcIdx, targetCellSet)); + } + + } + + } + + delete[] polyNodes; + delete[] polyCoords; + + } + _split.releaseArrays(); + } +} + +#endif