1 // Copyright (C) 2007-2016 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 __INTERPOLATION3D1D_TXX__
22 #define __INTERPOLATION3D1D_TXX__
24 #include "Interpolation3D1D.hxx"
25 #include "Interpolation.txx"
26 #include "MeshElement.txx"
27 #include "PointLocator3DIntersectorP0P0.txx"
28 #include "PointLocator3DIntersectorP0P1.txx"
29 #include "PointLocator3DIntersectorP1P0.txx"
30 #include "PointLocator3DIntersectorP1P1.txx"
35 namespace INTERP_KERNEL
38 * Very similar to Interpolation3D::interpolateMeshes, except for the bounding boxes that can be
39 * adjusted in a similar fashion as in InterpolationPlanar::performAdjustmentOfBB()
41 template<class MyMeshType, class MatrixType>
42 int Interpolation3D1D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method)
44 if(InterpolationOptions::getIntersectionType() != PointLocator)
45 INTERP_KERNEL::Exception("Invalid 3D/1D intersection type specified : must be PointLocator.");
47 typedef typename MyMeshType::MyConnType ConnType;
48 // create MeshElement objects corresponding to each element of the two meshes
49 const unsigned long numSrcElems = srcMesh.getNumberOfElements();
50 const unsigned long numTargetElems = targetMesh.getNumberOfElements();
52 LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements ");
54 std::vector<MeshElement<ConnType>*> srcElems(numSrcElems);
55 std::vector<MeshElement<ConnType>*> targetElems(numTargetElems);
57 std::map<MeshElement<ConnType>*, int> indices;
59 for(unsigned long i = 0 ; i < numSrcElems ; ++i)
60 srcElems[i] = new MeshElement<ConnType>(i, srcMesh);
62 for(unsigned long i = 0 ; i < numTargetElems ; ++i)
63 targetElems[i] = new MeshElement<ConnType>(i, targetMesh);
65 Intersector3D<MyMeshType,MatrixType>* intersector=0;
66 std::string methC = InterpolationOptions::filterInterpolationMethod(method);
68 { intersector=new PointLocator3DIntersectorP0P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
70 else if(methC=="P0P1")
71 { intersector=new PointLocator3DIntersectorP0P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
73 else if(methC=="P1P0")
74 { intersector=new PointLocator3DIntersectorP1P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
76 else if(methC=="P1P1")
77 { intersector=new PointLocator3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
80 throw Exception("Invalid method chosen must be in \"P0P0\", \"P0P1\", \"P1P0\" or \"P1P1\".");
81 // create empty maps for all source elements
82 result.resize(intersector->getNumberOfRowsOfResMatrix());
84 // create BBTree structure
85 // - get bounding boxes
86 std::vector<double> bboxes(6*numSrcElems);
87 int* srcElemIdx = new int[numSrcElems];
88 for(unsigned long i = 0; i < numSrcElems ; ++i)
90 // get source bboxes in right order
91 const BoundingBox* box = srcElems[i]->getBoundingBox();
92 bboxes[6*i+0] = box->getCoordinate(BoundingBox::XMIN);
93 bboxes[6*i+1] = box->getCoordinate(BoundingBox::XMAX);
94 bboxes[6*i+2] = box->getCoordinate(BoundingBox::YMIN);
95 bboxes[6*i+3] = box->getCoordinate(BoundingBox::YMAX);
96 bboxes[6*i+4] = box->getCoordinate(BoundingBox::ZMIN);
97 bboxes[6*i+5] = box->getCoordinate(BoundingBox::ZMAX);
99 srcElemIdx[i] = srcElems[i]->getIndex();
102 adjustBoundingBoxes(bboxes);
103 const double *bboxPtr=0;
106 BBTree<3,ConnType> tree(bboxPtr, srcElemIdx, 0, numSrcElems);
108 // for each target element, get source elements with which to calculate intersection
109 // - calculate intersection by calling intersectCells
110 for(unsigned long i = 0; i < numTargetElems; ++i)
112 const BoundingBox* box = targetElems[i]->getBoundingBox();
113 const int targetIdx = targetElems[i]->getIndex();
115 // get target bbox in right order
117 targetBox[0] = box->getCoordinate(BoundingBox::XMIN);
118 targetBox[1] = box->getCoordinate(BoundingBox::XMAX);
119 targetBox[2] = box->getCoordinate(BoundingBox::YMIN);
120 targetBox[3] = box->getCoordinate(BoundingBox::YMAX);
121 targetBox[4] = box->getCoordinate(BoundingBox::ZMIN);
122 targetBox[5] = box->getCoordinate(BoundingBox::ZMAX);
124 std::vector<ConnType> intersectElems;
126 tree.getIntersectingElems(targetBox, intersectElems);
128 if ( !intersectElems.empty() )
129 intersector->intersectCells(targetIdx,intersectElems,result);
132 // free allocated memory
133 delete [] srcElemIdx;
135 int ret=intersector->getNumberOfColsOfResMatrix();
139 for(unsigned long i = 0 ; i < numSrcElems ; ++i)
143 for(unsigned long i = 0 ; i < numTargetElems ; ++i)
145 delete targetElems[i];