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 (EDF R&D)
23 #include "Interpolation3D1D.hxx"
24 #include "Interpolation.txx"
25 #include "MeshElement.txx"
26 #include "PointLocator3DIntersectorP0P0.txx"
27 #include "PointLocator3DIntersectorP0P1.txx"
28 #include "PointLocator3DIntersectorP1P0.txx"
29 #include "PointLocator3DIntersectorP1P1.txx"
36 namespace INTERP_KERNEL
39 * Very similar to Interpolation3D::interpolateMeshes, except for the bounding boxes that can be
40 * adjusted in a similar fashion as in InterpolationPlanar::performAdjustmentOfBB()
42 template<class MyMeshType, class MatrixType>
43 typename MyMeshType::MyConnType Interpolation3D1D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method)
45 if(InterpolationOptions::getIntersectionType() != PointLocator)
46 INTERP_KERNEL::Exception("Invalid 3D/1D-0D intersection type specified : must be PointLocator.");
48 typedef typename MyMeshType::MyConnType ConnType;
49 // create MeshElement objects corresponding to each element of the two meshes
50 const ConnType numSrcElems = srcMesh.getNumberOfElements();
51 const ConnType numTargetElems = targetMesh.getNumberOfElements();
53 LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements ");
55 std::vector< std::unique_ptr< MeshElement<ConnType> > > srcElems(numSrcElems);
56 std::vector< std::unique_ptr< MeshElement<ConnType> > > targetElems(numTargetElems);
58 std::map<MeshElement<ConnType>*, ConnType> indices;
60 for(ConnType i = 0 ; i < numSrcElems ; ++i)
61 srcElems[i].reset( new MeshElement<ConnType>(i, srcMesh) );
63 for(ConnType i = 0 ; i < numTargetElems ; ++i)
64 targetElems[i].reset( new MeshElement<ConnType>(i, targetMesh) );
66 std::unique_ptr< Intersector3D<MyMeshType,MatrixType> > intersector;
67 std::string methC = InterpolationOptions::filterInterpolationMethod(method);
69 { intersector.reset( new PointLocator3DIntersectorP0P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
71 else if(methC=="P0P1")
72 { intersector.reset( new PointLocator3DIntersectorP0P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
74 else if(methC=="P1P0")
75 { intersector.reset( new PointLocator3DIntersectorP1P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
77 else if(methC=="P1P1")
78 { intersector.reset( new PointLocator3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
81 throw Exception("Invalid method chosen must be in \"P0P0\", \"P0P1\", \"P1P0\" or \"P1P1\".");
82 // create empty maps for all source elements
83 result.resize(intersector->getNumberOfRowsOfResMatrix());
85 // create BBTree structure
86 // - get bounding boxes
87 std::vector<double> bboxes(6*numSrcElems);
88 std::unique_ptr<ConnType[]> srcElemIdx{ new ConnType[numSrcElems] };
89 for(ConnType i = 0; i < numSrcElems ; ++i)
91 // get source bboxes in right order
92 const BoundingBox* box = srcElems[i]->getBoundingBox();
93 bboxes[6*i+0] = box->getCoordinate(BoundingBox::XMIN);
94 bboxes[6*i+1] = box->getCoordinate(BoundingBox::XMAX);
95 bboxes[6*i+2] = box->getCoordinate(BoundingBox::YMIN);
96 bboxes[6*i+3] = box->getCoordinate(BoundingBox::YMAX);
97 bboxes[6*i+4] = box->getCoordinate(BoundingBox::ZMIN);
98 bboxes[6*i+5] = box->getCoordinate(BoundingBox::ZMAX);
100 srcElemIdx[i] = srcElems[i]->getIndex();
103 adjustBoundingBoxes(bboxes);
104 const double *bboxPtr = nullptr;
106 bboxPtr=bboxes.data();
107 BBTree<3,ConnType> tree(bboxPtr, srcElemIdx.get(), 0, numSrcElems);
109 // for each target element, get source elements with which to calculate intersection
110 // - calculate intersection by calling intersectCells
111 for(ConnType i = 0; i < numTargetElems; ++i)
113 const BoundingBox* box = targetElems[i]->getBoundingBox();
114 const ConnType targetIdx = targetElems[i]->getIndex();
116 // get target bbox in right order
118 targetBox[0] = box->getCoordinate(BoundingBox::XMIN);
119 targetBox[1] = box->getCoordinate(BoundingBox::XMAX);
120 targetBox[2] = box->getCoordinate(BoundingBox::YMIN);
121 targetBox[3] = box->getCoordinate(BoundingBox::YMAX);
122 targetBox[4] = box->getCoordinate(BoundingBox::ZMIN);
123 targetBox[5] = box->getCoordinate(BoundingBox::ZMAX);
125 std::vector<ConnType> intersectElems;
127 tree.getIntersectingElems(targetBox, intersectElems);
129 if ( !intersectElems.empty() )
130 intersector->intersectCells(targetIdx,intersectElems,result);
132 return intersector->getNumberOfColsOfResMatrix();