src/INTERP_KERNEL/Interpolation2D3D.txx

   1 // Copyright (C) 2007-2022  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 #ifndef __INTERPOLATION2D3D_TXX__
  20 #define __INTERPOLATION2D3D_TXX__
  21
  22 #include "Interpolation2D3D.hxx"
  23 #include "Interpolation.txx"
  24 #include "MeshElement.txx"
  25 #include "TransformedTriangle.hxx"
  26 #include "Polyhedron3D2DIntersectorP0P0.txx"
  27 #include "PointLocator3DIntersectorP0P0.txx"
  28 #include "PolyhedronIntersectorP0P1.txx"
  29 #include "PointLocator3DIntersectorP0P1.txx"
  30 #include "PolyhedronIntersectorP1P0.txx"
  31 #include "PolyhedronIntersectorP1P0Bary.txx"
  32 #include "PointLocator3DIntersectorP1P0.txx"
  33 #include "PolyhedronIntersectorP1P1.txx"
  34 #include "PointLocator3DIntersectorP1P1.txx"
  35 #include "InterpolationHelper.txx"
  36
  37 #include "BBTree.txx"
  38
  39 namespace INTERP_KERNEL
  40 {
  41   /**
  42    * Calculates the matrix of volumes of intersection between the elements of srcMesh and the elements of targetMesh.
  43    * The calculation is done in two steps. First a filtering process reduces the number of pairs of elements for which the
  44    * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the
  45    * volume of intersection is calculated by an object of type Intersector3D for the remaining pairs, and entered into the
  46    * intersection matrix.
  47    *
  48    * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
  49    * It can also be an INTERP_KERNEL::Matrix object.
  50    * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
  51    * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
  52    * which have a non-zero intersection volume with the target element. The vector has indices running from
  53    * 0 to (nb target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
  54    * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
  55    *
  56
  57    * @param srcMesh     3DSurf source mesh (meshDim=2,spaceDim=3)
  58    * @param targetMesh  3D target mesh, containing only tetraedra
  59    * @param matrix      matrix in which the result is stored
  60    *
  61    */
  62   template<class MyMeshType, class MyMatrixType>
  63   typename MyMeshType::MyConnType Interpolation2D3D::interpolateMeshes(const MyMeshType& srcMesh,
  64                                            const MyMeshType& targetMesh,
  65                                            MyMatrixType& matrix,
  66                                            const std::string& method)
  67   {
  68     typedef typename MyMeshType::MyConnType ConnType;
  69     // create MeshElement objects corresponding to each element of the two meshes
  70     const ConnType numTargetElems = targetMesh.getNumberOfElements();
  71
  72     LOG(2, "Target mesh has " << numTargetElems << " elements ");
  73
  74     DuplicateFacesType intersectFaces;
  75
  76     std::unique_ptr< Intersector3D<MyMeshType,MyMatrixType> > intersector;
  77     std::string methC = InterpolationOptions::filterInterpolationMethod(method);
  78     const double dimCaracteristic = CalculateCharacteristicSizeOfMeshes(srcMesh, targetMesh, InterpolationOptions::getPrintLevel());
  79     if(methC=="P0P0")
  80       {
  81         switch(InterpolationOptions::getIntersectionType())
  82           {
  83           case Triangulation:
  84             intersector.reset( new Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrixType>(targetMesh,
  85                                                                                    srcMesh,
  86                                                                                    dimCaracteristic,
  87                                                                                    getPrecision(),
  88                                                                                    intersectFaces,
  89                                                                                    getSplittingPolicy()) );
  90             break;
  91           default:
  92             throw INTERP_KERNEL::Exception("Invalid 2D to 3D intersection type for P0P0 interp specified : must be Triangulation.");
  93           }
  94       }
  95     else
  96       throw Exception("Invalid method chosen must be in \"P0P0\".");
  97     // create empty maps for all source elements
  98     matrix.resize(intersector->getNumberOfRowsOfResMatrix());
  99
 100     // create BBTree structure
 101     // [ABN] Adjust 2D bounding box (those might be flat in the cases where the 2D surf are perfectly aligned with the axis)
 102     BBTreeStandAlone<3,ConnType> tree( BuildBBTreeWithAdjustment(srcMesh,[this,&intersector](double *bbox, typename MyMeshType::MyConnType sz){ this->performAdjustmentOfBB(intersector.get(),bbox,sz); }) );
 103
 104     // for each target element, get source elements with which to calculate intersection
 105     // - calculate intersection by calling intersectCells
 106     for(ConnType i = 0; i < numTargetElems; ++i)
 107       {
 108         MeshElement<ConnType> trgMeshElem(i, targetMesh);
 109
 110         const BoundingBox *box = trgMeshElem.getBoundingBox();
 111
 112         // get target bbox in right order
 113         double targetBox[6];
 114         box->fillInXMinXmaxYminYmaxZminZmaxFormat(targetBox);
 115
 116         std::vector<ConnType> intersectElems;
 117
 118         tree.getIntersectingElems(targetBox, intersectElems);
 119
 120         if ( !intersectElems.empty() )
 121           intersector->intersectCells(i, intersectElems, matrix);
 122       }
 123
 124     DuplicateFacesType::iterator iter;
 125     for (iter = intersectFaces.begin(); iter != intersectFaces.end(); ++iter)
 126       {
 127         if (iter->second.size() > 1)
 128           {
 129             _duplicate_faces.insert(std::make_pair(iter->first, iter->second));
 130           }
 131       }
 132
 133     // free allocated memory
 134     ConnType ret=intersector->getNumberOfColsOfResMatrix();
 135
 136     return ret;
 137
 138   }
 139 }
 140
 141 #endif