1 // Copyright (C) 2007-2021 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 // Author : Anthony Geay (CEA/DEN)
21 #ifndef __INTERPOLATION3D_TXX__
22 #define __INTERPOLATION3D_TXX__
24 #include "Interpolation3D.hxx"
25 #include "Interpolation.txx"
26 #include "MeshElement.txx"
27 #include "TransformedTriangle.hxx"
28 #include "PolyhedronIntersectorP0P0.txx"
29 #include "PointLocator3DIntersectorP0P0.txx"
30 #include "PolyhedronIntersectorP0P1.txx"
31 #include "PointLocator3DIntersectorP0P1.txx"
32 #include "PolyhedronIntersectorP1P0.txx"
33 #include "PolyhedronIntersectorP1P0Bary.txx"
34 #include "PointLocator3DIntersectorP1P0.txx"
35 #include "PolyhedronIntersectorP1P1.txx"
36 #include "PointLocator3DIntersectorP1P1.txx"
37 #include "Barycentric3DIntersectorP1P1.txx"
38 #include "MappedBarycentric3DIntersectorP1P1.txx"
40 // If defined, use recursion to traverse the binary search tree, else use the BBTree class
41 //#define USE_RECURSIVE_BBOX_FILTER
43 #ifdef USE_RECURSIVE_BBOX_FILTER
44 #include "MeshRegion.txx"
45 #include "RegionNode.hxx"
48 #else // use BBTree class
54 namespace INTERP_KERNEL
57 * Calculates the matrix of volumes of intersection between the elements of srcMesh and the elements of targetMesh.
58 * The calculation is done in two steps. First a filtering process reduces the number of pairs of elements for which the
59 * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the
60 * volume of intersection is calculated by an object of type Intersector3D for the remaining pairs, and entered into the
61 * intersection matrix.
63 * The matrix is partially sparse : it is a vector of maps of integer - double pairs.
64 * It can also be an INTERP_KERNEL::Matrix object.
65 * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
66 * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
67 * which have a non-zero intersection volume with the target element. The vector has indices running from
68 * 0 to (nb target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
69 * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
72 * @param srcMesh 3-dimensional source mesh
73 * @param targetMesh 3-dimesional target mesh, containing only tetraedra
74 * @param result matrix in which the result is stored
77 template<class MyMeshType, class MatrixType>
78 typename MyMeshType::MyConnType Interpolation3D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method)
80 typedef typename MyMeshType::MyConnType ConnType;
81 // create MeshElement objects corresponding to each element of the two meshes
82 const ConnType numSrcElems = srcMesh.getNumberOfElements();
83 const ConnType numTargetElems = targetMesh.getNumberOfElements();
85 LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements ");
87 std::vector<MeshElement<ConnType>*> srcElems(numSrcElems);
88 std::vector<MeshElement<ConnType>*> targetElems(numTargetElems);
90 std::map<MeshElement<ConnType>*, ConnType> indices;
92 for(ConnType i = 0 ; i < numSrcElems ; ++i)
93 srcElems[i] = new MeshElement<ConnType>(i, srcMesh);
95 for(ConnType i = 0 ; i < numTargetElems ; ++i)
96 targetElems[i] = new MeshElement<ConnType>(i, targetMesh);
98 Intersector3D<MyMeshType,MatrixType>* intersector=0;
99 std::string methC = InterpolationOptions::filterInterpolationMethod(method);
102 switch(InterpolationOptions::getIntersectionType())
105 intersector=new PolyhedronIntersectorP0P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy());
108 intersector=new PointLocator3DIntersectorP0P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
111 throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P0P0 interp specified : must be Triangle or PointLocator.");
114 else if(methC=="P0P1")
116 switch(InterpolationOptions::getIntersectionType())
119 intersector=new PolyhedronIntersectorP0P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy());
122 intersector=new PointLocator3DIntersectorP0P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
125 throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P0P1 interp specified : must be Triangle or PointLocator.");
128 else if(methC=="P1P0")
130 switch(InterpolationOptions::getIntersectionType())
133 intersector=new PolyhedronIntersectorP1P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy());
136 intersector=new PointLocator3DIntersectorP1P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
139 intersector=new PolyhedronIntersectorP1P0Bary<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy());
142 throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P1P0 interp specified : must be Triangle, PointLocator or Barycentric.");
145 else if(methC=="P1P1")
147 switch(InterpolationOptions::getIntersectionType())
150 intersector=new PolyhedronIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getSplittingPolicy());
153 intersector=new PointLocator3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
156 intersector=new Barycentric3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
158 case MappedBarycentric:
159 intersector=new MappedBarycentric3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
162 throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P1P1 interp specified : must be Triangle, PointLocator, Barycentric or MappedBarycentric.");
166 throw Exception("Invalid method chosen must be in \"P0P0\", \"P0P1\", \"P1P0\" or \"P1P1\".");
167 // create empty maps for all source elements
168 result.resize(intersector->getNumberOfRowsOfResMatrix());
170 #ifdef USE_RECURSIVE_BBOX_FILTER
173 * Performs a depth-first search over srcMesh, using bounding boxes to recursively eliminate the elements of targetMesh
174 * which cannot intersect smaller and smaller regions of srcMesh. At each level, each region is divided in two, forming
175 * a binary search tree with leaves consisting of only one element of the source mesh together with the elements of the
176 * target mesh that can intersect it. The recursion is implemented with a stack of RegionNodes, each one containing a
177 * source region and a target region. Each region has an associated bounding box and a vector of pointers to the elements
178 * that belong to it. Each MeshElement contains a bounding box and the global number of the corresponding element in the mesh.
181 // create initial RegionNode and fill up its source region with all the source mesh elements and
182 // its target region with all the target mesh elements whose bounding box
183 // intersects that of the source region
185 RegionNode<ConnType>* firstNode = new RegionNode<ConnType>();
187 MeshRegion<ConnType>& srcRegion = firstNode->getSrcRegion();
189 for(ConnType i = 0 ; i < numSrcElems ; ++i)
191 srcRegion.addElement(srcElems[i], srcMesh);
194 MeshRegion<ConnType>& targetRegion = firstNode->getTargetRegion();
196 for(ConnType i = 0 ; i < numTargetElems ; ++i)
198 if(!srcRegion.isDisjointWithElementBoundingBox( *(targetElems[i]) ))
200 targetRegion.addElement(targetElems[i], targetMesh);
204 // Using a stack, descend recursively, creating at each step two new RegionNodes having as source region the left and
205 // right part of the source region of the current node (created using MeshRegion::split()) and as target region all the
206 // elements of the target mesh whose bounding box intersects the corresponding part
207 // Continue until the source region contains only one element, at which point the intersection volumes are
208 // calculated with all the remaining target mesh elements and stored in the matrix if they are non-zero.
210 std::stack< RegionNode<ConnType>* > nodes;
211 nodes.push(firstNode);
213 while(!nodes.empty())
215 RegionNode<ConnType>* currNode = nodes.top();
217 LOG(4, "Popping node ");
219 if(currNode->getTargetRegion().getNumberOfElements() == 1)
222 LOG(4, " - One element");
224 MeshElement<ConnType>* targetElement = *(currNode->getTargetRegion().getBeginElements());
225 std::vector<ConnType> intersectElems;
226 for(typename std::vector< MeshElement<ConnType>* >::const_iterator iter = currNode->getSrcRegion().getBeginElements();iter != currNode->getSrcRegion().getEndElements();++iter)
227 intersectElems.push_back((*iter)->getIndex());
228 intersector->intersectCells(targetElement->getIndex(),intersectElems,result);
233 LOG(4, " - Recursion");
235 RegionNode<ConnType>* leftNode = new RegionNode<ConnType>();
236 RegionNode<ConnType>* rightNode = new RegionNode<ConnType>();
238 // split current source region
240 static BoundingBox::BoxCoord axis = BoundingBox::XMAX;
242 currNode->getTargetRegion().split(leftNode->getTargetRegion(), rightNode->getTargetRegion(), axis, targetMesh);
244 LOG(5, "After split, left target region has " << leftNode->getTargetRegion().getNumberOfElements()
245 << " elements and right target region has " << rightNode->getTargetRegion().getNumberOfElements()
248 // ugly hack to avoid problem with enum which does not start at 0
249 // I guess I ought to implement ++ for it instead ...
250 // Anyway, it basically chooses the next axis, cyclically
251 axis = (axis != BoundingBox::ZMAX) ? static_cast<BoundingBox::BoxCoord>(axis + 1) : BoundingBox::XMAX;
253 // add source elements of current node that overlap the target regions of the new nodes
254 LOG(5, " -- Adding source elements");
255 ConnType numLeftElements = 0;
256 ConnType numRightElements = 0;
257 for(typename std::vector<MeshElement<ConnType>*>::const_iterator iter = currNode->getSrcRegion().getBeginElements() ;
258 iter != currNode->getSrcRegion().getEndElements() ; ++iter)
260 LOG(6, " --- New target node");
262 if(!leftNode->getTargetRegion().isDisjointWithElementBoundingBox(**iter))
264 leftNode->getSrcRegion().addElement(*iter, srcMesh);
268 if(!rightNode->getTargetRegion().isDisjointWithElementBoundingBox(**iter))
270 rightNode->getSrcRegion().addElement(*iter, srcMesh);
276 LOG(5, "Left src region has " << numLeftElements << " elements and right src region has "
277 << numRightElements << " elements");
279 // push new nodes on stack
280 if(numLeftElements != 0)
282 nodes.push(leftNode);
289 if(numRightElements != 0)
291 nodes.push(rightNode);
299 // all nodes are deleted here
302 LOG(4, "Next iteration. Nodes left : " << nodes.size());
307 // create BBTree structure
308 // - get bounding boxes
309 double* bboxes = new double[6 * numSrcElems];
310 ConnType* srcElemIdx = new ConnType[numSrcElems];
311 for(ConnType i = 0; i < numSrcElems ; ++i)
313 // get source bboxes in right order
314 const BoundingBox* box = srcElems[i]->getBoundingBox();
315 bboxes[6*i+0] = box->getCoordinate(BoundingBox::XMIN);
316 bboxes[6*i+1] = box->getCoordinate(BoundingBox::XMAX);
317 bboxes[6*i+2] = box->getCoordinate(BoundingBox::YMIN);
318 bboxes[6*i+3] = box->getCoordinate(BoundingBox::YMAX);
319 bboxes[6*i+4] = box->getCoordinate(BoundingBox::ZMIN);
320 bboxes[6*i+5] = box->getCoordinate(BoundingBox::ZMAX);
322 // source indices have to begin with zero for BBox, I think
323 srcElemIdx[i] = srcElems[i]->getIndex();
326 BBTree<3,ConnType> tree(bboxes, srcElemIdx, 0, numSrcElems);
328 // for each target element, get source elements with which to calculate intersection
329 // - calculate intersection by calling intersectCells
330 for(ConnType i = 0; i < numTargetElems; ++i)
332 const BoundingBox* box = targetElems[i]->getBoundingBox();
333 const ConnType targetIdx = targetElems[i]->getIndex();
335 // get target bbox in right order
337 targetBox[0] = box->getCoordinate(BoundingBox::XMIN);
338 targetBox[1] = box->getCoordinate(BoundingBox::XMAX);
339 targetBox[2] = box->getCoordinate(BoundingBox::YMIN);
340 targetBox[3] = box->getCoordinate(BoundingBox::YMAX);
341 targetBox[4] = box->getCoordinate(BoundingBox::ZMIN);
342 targetBox[5] = box->getCoordinate(BoundingBox::ZMAX);
344 std::vector<ConnType> intersectElems;
346 tree.getIntersectingElems(targetBox, intersectElems);
348 if ( !intersectElems.empty() )
349 intersector->intersectCells(targetIdx,intersectElems,result);
353 delete [] srcElemIdx;
356 // free allocated memory
357 ConnType ret=intersector->getNumberOfColsOfResMatrix();
361 for(ConnType i = 0 ; i < numSrcElems ; ++i)
365 for(ConnType i = 0 ; i < numTargetElems ; ++i)
367 delete targetElems[i];