1 // Copyright (C) 2007-2019 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 #ifndef __INTERPOLATION2D3D_TXX__
20 #define __INTERPOLATION2D3D_TXX__
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"
39 namespace INTERP_KERNEL
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.
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].
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
62 template<class MyMeshType, class MyMatrixType>
63 typename MyMeshType::MyConnType Interpolation2D3D::interpolateMeshes(const MyMeshType& srcMesh,
64 const MyMeshType& targetMesh,
66 const std::string& method)
68 typedef typename MyMeshType::MyConnType ConnType;
69 // create MeshElement objects corresponding to each element of the two meshes
70 const ConnType numSrcElems = srcMesh.getNumberOfElements();
71 const ConnType numTargetElems = targetMesh.getNumberOfElements();
73 LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements ");
75 std::vector<MeshElement<ConnType>*> srcElems(numSrcElems);
76 std::vector<MeshElement<ConnType>*> targetElems(numTargetElems);
78 std::map<MeshElement<ConnType>*, ConnType> indices;
79 DuplicateFacesType intersectFaces;
81 for(ConnType i = 0 ; i < numSrcElems ; ++i)
82 srcElems[i] = new MeshElement<ConnType>(i, srcMesh);
84 for(ConnType i = 0 ; i < numTargetElems ; ++i)
85 targetElems[i] = new MeshElement<ConnType>(i, targetMesh);
87 Intersector3D<MyMeshType,MyMatrixType>* intersector=0;
88 std::string methC = InterpolationOptions::filterInterpolationMethod(method);
89 const double dimCaracteristic = CalculateCharacteristicSizeOfMeshes(srcMesh, targetMesh, InterpolationOptions::getPrintLevel());
92 switch(InterpolationOptions::getIntersectionType())
95 intersector=new Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrixType>(targetMesh,
100 getSplittingPolicy());
103 throw INTERP_KERNEL::Exception("Invalid 2D to 3D intersection type for P0P0 interp specified : must be Triangulation.");
107 throw Exception("Invalid method chosen must be in \"P0P0\".");
108 // create empty maps for all source elements
109 matrix.resize(intersector->getNumberOfRowsOfResMatrix());
111 // create BBTree structure
112 // - get bounding boxes
113 double* bboxes = new double[6 * numSrcElems];
114 ConnType* srcElemIdx = new ConnType[numSrcElems];
115 for(ConnType i = 0; i < numSrcElems ; ++i)
117 // get source bboxes in right order
118 const BoundingBox* box = srcElems[i]->getBoundingBox();
119 bboxes[6*i+0] = box->getCoordinate(BoundingBox::XMIN);
120 bboxes[6*i+1] = box->getCoordinate(BoundingBox::XMAX);
121 bboxes[6*i+2] = box->getCoordinate(BoundingBox::YMIN);
122 bboxes[6*i+3] = box->getCoordinate(BoundingBox::YMAX);
123 bboxes[6*i+4] = box->getCoordinate(BoundingBox::ZMIN);
124 bboxes[6*i+5] = box->getCoordinate(BoundingBox::ZMAX);
126 // source indices have to begin with zero for BBox, I think
127 srcElemIdx[i] = srcElems[i]->getIndex();
130 BBTree<3,ConnType> tree(bboxes, srcElemIdx, 0, numSrcElems, 0.);
132 // for each target element, get source elements with which to calculate intersection
133 // - calculate intersection by calling intersectCells
134 for(ConnType i = 0; i < numTargetElems; ++i)
136 const BoundingBox* box = targetElems[i]->getBoundingBox();
137 const ConnType targetIdx = targetElems[i]->getIndex();
139 // get target bbox in right order
141 targetBox[0] = box->getCoordinate(BoundingBox::XMIN);
142 targetBox[1] = box->getCoordinate(BoundingBox::XMAX);
143 targetBox[2] = box->getCoordinate(BoundingBox::YMIN);
144 targetBox[3] = box->getCoordinate(BoundingBox::YMAX);
145 targetBox[4] = box->getCoordinate(BoundingBox::ZMIN);
146 targetBox[5] = box->getCoordinate(BoundingBox::ZMAX);
148 std::vector<ConnType> intersectElems;
150 tree.getIntersectingElems(targetBox, intersectElems);
152 if ( !intersectElems.empty() )
153 intersector->intersectCells(targetIdx, intersectElems, matrix);
158 delete [] srcElemIdx;
160 DuplicateFacesType::iterator iter;
161 for (iter = intersectFaces.begin(); iter != intersectFaces.end(); ++iter)
163 if (iter->second.size() > 1)
165 _duplicate_faces.insert(std::make_pair(iter->first, iter->second));
169 // free allocated memory
170 ConnType ret=intersector->getNumberOfColsOfResMatrix();
174 for(ConnType i = 0 ; i < numSrcElems ; ++i)
178 for(ConnType i = 0 ; i < numTargetElems ; ++i)
180 delete targetElems[i];