--- /dev/null
+// 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<Interpolation2D1D>
+ {
+ public:
+ typedef std::map<int,std::set<int> > DuplicateFacesType;
+
+ Interpolation2D1D() { setOrientation(2); }
+ Interpolation2D1D(const InterpolationOptions& io):Interpolation<Interpolation2D1D>(io) { }
+ public:
+
+ // Main function to interpolate triangular and quadratic meshes
+ template<class MyMeshType, class MatrixType>
+ 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
--- /dev/null
+// 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<class MyMeshType, class MatrixType>
+ 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<MyMeshType,MatrixType>* intersector=0;
+ std::string meth = InterpolationOptions::filterInterpolationMethod(method);
+ if(meth=="P0P0")
+ {
+ switch (InterpolationOptions::getIntersectionType())
+ {
+ case Geometric2D:
+ intersector=new Geometric2DIntersector<MyMeshType,MatrixType,Planar2D1DIntersectorP0P0>(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<double> bbox;
+ intersector->createBoundingBoxes(myMeshS,bbox); // create the bounding boxes
+ const double *bboxPtr=0;
+ if(nbMailleS>0)
+ bboxPtr=&bbox[0];
+ BBTree<SPACEDIM,ConnType> 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<nbelem_type; iT++)
+ {
+ int nb_nodesT=connIndxT[iT+1]-connIndxT[iT];
+ std::vector<int> intersecting_elems;
+ double bb[2*SPACEDIM];
+ intersector->getElemBB(bb,myMeshT,OTT<ConnType,numPol>::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
--- /dev/null
+// 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<Interpolation3D2D>(io)
+ {
+ }
+}
--- /dev/null
+// 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 <set>
+#include <map>
+
+#include "Interpolation.hxx"
+#include "NormalizedUnstructuredMesh.hxx"
+#include "InterpolationOptions.hxx"
+
+namespace INTERP_KERNEL
+{
+ class Interpolation3D2D : public Interpolation<Interpolation3D2D>
+ {
+ public:
+ typedef std::map<int,std::set<int> > DuplicateFacesType;
+
+ Interpolation3D2D();
+ Interpolation3D2D(const InterpolationOptions& io);
+ template<class MyMeshType, class MyMatrixType>
+ 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
--- /dev/null
+// 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<class MyMeshType, class MyMatrixType>
+ 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<MeshElement<ConnType>*> srcElems(numSrcElems);
+ std::vector<MeshElement<ConnType>*> targetElems(numTargetElems);
+
+ std::map<MeshElement<ConnType>*, int> indices;
+ DuplicateFacesType intersectFaces;
+
+ for(unsigned long i = 0 ; i < numSrcElems ; ++i)
+ srcElems[i] = new MeshElement<ConnType>(i, srcMesh);
+
+ for(unsigned long i = 0 ; i < numTargetElems ; ++i)
+ targetElems[i] = new MeshElement<ConnType>(i, targetMesh);
+
+ Intersector3D<MyMeshType,MyMatrixType>* 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<MyMeshType,MyMatrixType>(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<ConnType> 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
--- /dev/null
+// 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 MyMeshType, class MyMatrix, class ConcreteP0P0Intersector>
+ class Planar2D1DIntersectorP0P0 : public PlanarIntersector<MyMeshType,MyMatrix>
+ {
+ public:
+ static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
+ static const int MESHDIM=MyMeshType::MY_MESHDIM;
+ typedef typename MyMeshType::MyConnType ConnType;
+ static const NumberingPolicy numPol=MyMeshType::My_numPol;
+ 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<MyMeshType,MyMatrix>::DuplicateFacesType* getIntersectFaces() const
+ {
+ return &_intersect_faces;
+ }
+ void intersectCells(ConnType icellT, const std::vector<ConnType>& 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<ConcreteP0P0Intersector&>(*this); }
+ private:
+ typename PlanarIntersector<MyMeshType,MyMatrix>::DuplicateFacesType _intersect_faces;
+ };
+}
+
+#endif
--- /dev/null
+// 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<class MyMeshType, class MyMatrix, class ConcreteP0P0Intersector>
+ Planar2D1DIntersectorP0P0<MyMeshType,MyMatrix,ConcreteP0P0Intersector>::Planar2D1DIntersectorP0P0(const MyMeshType& meshT, const MyMeshType& meshS,
+ double dimCaracteristic, double precision, double md3DSurf, double medianPlane,
+ bool doRotate, int orientation, int printLevel):
+ PlanarIntersector<MyMeshType,MyMatrix>(meshT,meshS,dimCaracteristic,precision,md3DSurf,medianPlane,doRotate,orientation,printLevel)
+ {
+ }
+
+ template<class MyMeshType, class MyMatrix, class ConcreteP0P0Intersector>
+ int Planar2D1DIntersectorP0P0<MyMeshType,MyMatrix,ConcreteP0P0Intersector>::getNumberOfRowsOfResMatrix() const
+ {
+ return PlanarIntersector<MyMeshType,MyMatrix>::_meshT.getNumberOfElements();
+ }
+
+ template<class MyMeshType, class MyMatrix, class ConcreteP0P0Intersector>
+ int Planar2D1DIntersectorP0P0<MyMeshType,MyMatrix,ConcreteP0P0Intersector>::getNumberOfColsOfResMatrix() const
+ {
+ return PlanarIntersector<MyMeshType,MyMatrix>::_meshS.getNumberOfElements();
+ }
+
+ template<class MyMeshType, class MyMatrix, class ConcreteP0P0Intersector>
+ void Planar2D1DIntersectorP0P0<MyMeshType,MyMatrix,ConcreteP0P0Intersector>::intersectCells(ConnType icellT, const std::vector<ConnType>& icellsS, MyMatrix& res)
+ {
+ int nbNodesT=PlanarIntersector<MyMeshType,MyMatrix>::_connIndexT[icellT+1]-PlanarIntersector<MyMeshType,MyMatrix>::_connIndexT[icellT];
+ typename MyMatrix::value_type& resRow=res[icellT];
+ for(typename std::vector<ConnType>::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++)
+ {
+ int iS=*iter;
+ int nbNodesS=PlanarIntersector<MyMeshType,MyMatrix>::_connIndexS[iS+1]-PlanarIntersector<MyMeshType,MyMatrix>::_connIndexS[iS];
+ bool isColinear = false;
+ double surf=intersectGeometry1D(OTT<ConnType,numPol>::indFC(icellT),OTT<ConnType,numPol>::indFC(iS),
+ nbNodesT,nbNodesS, isColinear);
+ surf=PlanarIntersector<MyMeshType,MyMatrix>::getValueRegardingOption(surf);
+ if(surf!=0.)
+ resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(iS),surf));
+
+ if (isColinear)
+ {
+ typename PlanarIntersector<MyMeshType,MyMatrix>::DuplicateFacesType::iterator intersectFacesIter = _intersect_faces.find(iS);
+ if (intersectFacesIter != _intersect_faces.end())
+ {
+ intersectFacesIter->second.insert(icellT);
+ }
+ else
+ {
+ std::set<int> targetCellSet;
+ targetCellSet.insert(icellT);
+ _intersect_faces.insert(std::make_pair(iS, targetCellSet));
+ }
+ }
+ }
+ }
+}
+
+#endif
#include "TargetIntersector.hxx"
#include "NormalizedUnstructuredMesh.hxx"
+#include <map>
+#include <set>
+
namespace INTERP_KERNEL
{
class TranslationRotationMatrix;
static const int MESHDIM=MyMeshType::MY_MESHDIM;
typedef typename MyMeshType::MyConnType ConnType;
static const NumberingPolicy numPol=MyMeshType::My_numPol;
+ typedef typename std::map<int,std::set<int> > 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);
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<double>& coordsT);
--- /dev/null
+// 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 MyMeshType, class MyMatrixType>
+ class Polyhedron3D2DIntersectorP0P0 : public Intersector3DP0P0<MyMeshType,MyMatrixType>
+ {
+ typedef typename std::map<int,std::set<int> > 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<ConnType>& 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<MyMeshType>* > _tetra;
+
+ SplitterTetra2<MyMeshType> _split;
+
+ double _dim_caracteristic;
+ double _precision;
+
+ DuplicateFacesType& _intersect_faces;
+
+ };
+}
+
+#endif
--- /dev/null
+// 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<class MyMeshType, class MyMatrixType>
+ Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrixType>::Polyhedron3D2DIntersectorP0P0(const MyMeshType& targetMesh,
+ const MyMeshType& srcMesh,
+ const double dimCaracteristic,
+ const double precision,
+ DuplicateFacesType& intersectFaces,
+ SplittingPolicy policy)
+ : Intersector3DP0P0<MyMeshType,MyMatrixType>(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<class MyMeshType, class MyMatrixType>
+ Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrixType>::~Polyhedron3D2DIntersectorP0P0()
+ {
+ releaseArrays();
+ }
+
+ template<class MyMeshType, class MyMatrixType>
+ void Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrixType>::releaseArrays()
+ {
+ for(typename std::vector< SplitterTetra<MyMeshType>* >::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<class MyMeshType, class MyMatrixType>
+ void Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrixType>::intersectCells(ConnType targetCell,
+ const std::vector<ConnType>& srcCells,
+ MyMatrixType& matrix)
+ {
+ int nbOfNodesT=Intersector3D<MyMeshType,MyMatrixType>::_target_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(targetCell));
+ releaseArrays();
+ _split.splitTargetCell(targetCell,nbOfNodesT,_tetra);
+
+ for(typename std::vector<ConnType>::const_iterator iterCellS=srcCells.begin();iterCellS!=srcCells.end();iterCellS++)
+ {
+ double surface = 0.;
+ std::multiset<TriangleFaceKey> listOfTetraFacesTreated;
+ std::set<TriangleFaceKey> listOfTetraFacesColinear;
+
+ // calculate the coordinates of the nodes
+ const NumberingPolicy numPol=MyMeshType::My_numPol;
+ typename MyMeshType::MyConnType cellSrc = *iterCellS;
+ int cellSrcIdx = OTT<ConnType,numPol>::indFC(cellSrc);
+ NormalizedCellType normCellType=Intersector3D<MyMeshType,MyMatrixType>::_src_mesh.getTypeOfElement(cellSrcIdx);
+ const CellModel& cellModelCell=CellModel::GetCellModel(normCellType);
+ const MyMeshType& _src_mesh = Intersector3D<MyMeshType,MyMatrixType>::_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<ConnType,numPol>::indFC(*iterCellS), _src_mesh);
+ polyNodes[i]=globalNodeNum;
+ polyCoords[i] = (double*)_src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum;
+ }
+
+ for(typename std::vector<SplitterTetra<MyMeshType>*>::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<TriangleFaceKey>::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<int> targetCellSet;
+ targetCellSet.insert(targetCell);
+ _intersect_faces.insert(std::make_pair(cellSrcIdx, targetCellSet));
+ }
+
+ }
+
+ }
+
+ delete[] polyNodes;
+ delete[] polyCoords;
+
+ }
+ _split.releaseArrays();
+ }
+}
+
+#endif
