X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FINTERP_KERNEL%2FInterpolation3D.txx;h=67327f861a5b9765dacc630a018ba83611910156;hb=be4c3bb042d5426fbbe54378b9d7b35173ab27ef;hp=3597aa1377a8632a8d737e4b25ba682db5b23900;hpb=48782c06022ca2caa36f849cb5a29ea4fe2aaa83;p=tools%2Fmedcoupling.git diff --git a/src/INTERP_KERNEL/Interpolation3D.txx b/src/INTERP_KERNEL/Interpolation3D.txx index 3597aa137..67327f861 100644 --- a/src/INTERP_KERNEL/Interpolation3D.txx +++ b/src/INTERP_KERNEL/Interpolation3D.txx @@ -1,33 +1,43 @@ -// Copyright (C) 2007-2008 CEA/DEN, EDF R&D +// Copyright (C) 2007-2013 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 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. +// 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 +// 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 +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // +// Author : Anthony Geay (CEA/DEN) + #ifndef __INTERPOLATION3D_TXX__ #define __INTERPOLATION3D_TXX__ #include "Interpolation3D.hxx" +#include "Interpolation.txx" #include "MeshElement.txx" #include "TransformedTriangle.hxx" -#include "PolyhedronIntersector.txx" +#include "PolyhedronIntersectorP0P0.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 "Barycentric3DIntersectorP1P1.txx" #include "Log.hxx" /// If defined, use recursion to traverse the binary search tree, else use the BBTree class -#define USE_RECURSIVE_BBOX_FILTER +//#define USE_RECURSIVE_BBOX_FILTER #ifdef USE_RECURSIVE_BBOX_FILTER #include "MeshRegion.txx" @@ -42,23 +52,6 @@ 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 - * - */ - Interpolation3D::Interpolation3D() - { - } - Interpolation3D::Interpolation3D(const InterpolationOptions& io):Interpolation(io) - { - } - /** * 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 @@ -74,14 +67,14 @@ namespace INTERP_KERNEL * 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 result matrix in which the result is stored * */ template - int Interpolation3D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const char *method) + int Interpolation3D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method) { typedef typename MyMeshType::MyConnType ConnType; // create MeshElement objects corresponding to each element of the two meshes @@ -92,30 +85,86 @@ namespace INTERP_KERNEL std::vector*> srcElems(numSrcElems); std::vector*> targetElems(numTargetElems); - + std::map*, int> indices; - + 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(method); - if(method=="P0P0") - intersector=new PolyhedronIntersector(targetMesh, srcMesh, getSplittingPolicy()); - else if(method=="P0P1") - intersector=new PolyhedronIntersectorP0P1(targetMesh, srcMesh, getSplittingPolicy()); - else if(method=="P1P0") - intersector=new PolyhedronIntersectorP1P0(targetMesh, srcMesh, getSplittingPolicy()); + std::string methC = InterpolationOptions::filterInterpolationMethod(method); + if(methC=="P0P0") + { + switch(InterpolationOptions::getIntersectionType()) + { + case Triangulation: + intersector=new PolyhedronIntersectorP0P0(targetMesh, srcMesh, getSplittingPolicy()); + break; + case PointLocator: + intersector=new PointLocator3DIntersectorP0P0(targetMesh, srcMesh, getPrecision()); + break; + default: + throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P0P0 interp specified : must be Triangle or PointLocator."); + } + } + else if(methC=="P0P1") + { + switch(InterpolationOptions::getIntersectionType()) + { + case Triangulation: + intersector=new PolyhedronIntersectorP0P1(targetMesh, srcMesh, getSplittingPolicy()); + break; + case PointLocator: + intersector=new PointLocator3DIntersectorP0P1(targetMesh, srcMesh, getPrecision()); + break; + default: + throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P0P1 interp specified : must be Triangle or PointLocator."); + } + } + else if(methC=="P1P0") + { + switch(InterpolationOptions::getIntersectionType()) + { + case Triangulation: + intersector=new PolyhedronIntersectorP1P0(targetMesh, srcMesh, getSplittingPolicy()); + break; + case PointLocator: + intersector=new PointLocator3DIntersectorP1P0(targetMesh, srcMesh, getPrecision()); + break; + case Barycentric: + intersector=new PolyhedronIntersectorP1P0Bary(targetMesh, srcMesh, getSplittingPolicy()); + break; + default: + throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P1P0 interp specified : must be Triangle, PointLocator or Barycentric."); + } + } + else if(methC=="P1P1") + { + switch(InterpolationOptions::getIntersectionType()) + { + case Triangulation: + intersector=new PolyhedronIntersectorP1P1(targetMesh, srcMesh, getSplittingPolicy()); + break; + case PointLocator: + intersector=new PointLocator3DIntersectorP1P1(targetMesh, srcMesh, getPrecision()); + break; + case Barycentric: + intersector=new Barycentric3DIntersectorP1P1(targetMesh, srcMesh, getPrecision()); + break; + default: + throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P1P1 interp specified : must be Triangle or PointLocator."); + } + } else - throw Exception("Invalid method choosed must be in \"P0P0\", \"P0P1\"."); + throw Exception("Invalid method choosed must be in \"P0P0\", \"P0P1\", \"P1P0\" or \"P1P1\"."); // create empty maps for all source elements result.resize(intersector->getNumberOfRowsOfResMatrix()); #ifdef USE_RECURSIVE_BBOX_FILTER - + /* * Performs a depth-first search over srcMesh, using bounding boxes to recursively eliminate the elements of targetMesh * which cannot intersect smaller and smaller regions of srcMesh. At each level, each region is divided in two, forming @@ -130,7 +179,7 @@ namespace INTERP_KERNEL // intersects that of the source region RegionNode* firstNode = new RegionNode(); - + MeshRegion& srcRegion = firstNode->getSrcRegion(); for(unsigned long i = 0 ; i < numSrcElems ; ++i) @@ -181,11 +230,11 @@ namespace INTERP_KERNEL RegionNode* leftNode = new RegionNode(); RegionNode* rightNode = new RegionNode(); - + // split current source region //} decide on axis static BoundingBox::BoxCoord axis = BoundingBox::XMAX; - + currNode->getTargetRegion().split(leftNode->getTargetRegion(), rightNode->getTargetRegion(), axis, targetMesh); LOG(5, "After split, left target region has " << leftNode->getTargetRegion().getNumberOfElements() @@ -205,19 +254,19 @@ namespace INTERP_KERNEL iter != currNode->getSrcRegion().getEndElements() ; ++iter) { LOG(6, " --- New target node"); - + if(!leftNode->getTargetRegion().isDisjointWithElementBoundingBox(**iter)) { leftNode->getSrcRegion().addElement(*iter, srcMesh); ++numLeftElements; } - + if(!rightNode->getTargetRegion().isDisjointWithElementBoundingBox(**iter)) { rightNode->getSrcRegion().addElement(*iter, srcMesh); ++numRightElements; } - + } LOG(5, "Left src region has " << numLeftElements << " elements and right src region has " @@ -242,7 +291,7 @@ namespace INTERP_KERNEL delete rightNode; } } - + // all nodes are deleted here delete currNode; @@ -250,11 +299,11 @@ namespace INTERP_KERNEL } #else // Use BBTree - + // create BBTree structure // - get bounding boxes - double bboxes[6 * numSrcElems]; - int srcElemIdx[numSrcElems]; + 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 @@ -269,9 +318,9 @@ namespace INTERP_KERNEL // source indices have to begin with zero for BBox, I think srcElemIdx[i] = srcElems[i]->getIndex(); } - + BBTree<3,ConnType> tree(bboxes, srcElemIdx, 0, numSrcElems); - + // 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) @@ -292,9 +341,13 @@ namespace INTERP_KERNEL tree.getIntersectingElems(targetBox, intersectElems); - intersector->intersectCells(targetIdx,intersectElems,result); + if ( !intersectElems.empty() ) + intersector->intersectCells(targetIdx,intersectElems,result); } - + + delete [] bboxes; + delete [] srcElemIdx; + #endif // free allocated memory int ret=intersector->getNumberOfColsOfResMatrix(); @@ -312,7 +365,6 @@ namespace INTERP_KERNEL return ret; } - } #endif