1 // Copyright (C) 2018-2024 CEA, EDF
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)
21 #ifndef __INTERPOLATION1D0D_TXX__
22 #define __INTERPOLATION1D0D_TXX__
24 #include "Interpolation1D0D.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 #include "InterpKernelAssert.hxx"
37 namespace INTERP_KERNEL
40 * Very similar to Interpolation3D::interpolateMeshes, except for the bounding boxes that can be
41 * adjusted in a similar fashion as in InterpolationPlanar::performAdjustmentOfBB()
43 template<class MyMeshType, class MatrixType>
44 typename MyMeshType::MyConnType Interpolation1D0D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method)
46 constexpr int SPACEDIM=MyMeshType::MY_SPACEDIM;
47 using ConnType=typename MyMeshType::MyConnType;
48 IKAssert(SPACEDIM==3);
50 if(InterpolationOptions::getIntersectionType() != PointLocator)
51 INTERP_KERNEL::Exception("Invalid 1D/0D intersection type specified : must be PointLocator.");
53 std::string methC ( InterpolationOptions::filterInterpolationMethod(method) );
55 throw Exception("Invalid method chosen must be in \"P1P1\".");
57 const double epsilon(getPrecision());
58 // create MeshElement objects corresponding to each element of the two meshes
59 const ConnType numSrcElems(srcMesh.getNumberOfElements()), numTargetElems(targetMesh.getNumberOfElements());
61 LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements ");
63 std::vector<MeshElement<ConnType>*> srcElems(numSrcElems);
65 std::map<MeshElement<ConnType>*, int> indices;
67 for(ConnType i = 0 ; i < numSrcElems ; ++i)
68 srcElems[i] = new MeshElement<ConnType>(i, srcMesh);
70 // create empty maps for all source elements
71 result.resize(targetMesh.getNumberOfNodes());
73 // create BBTree structure
74 // - get bounding boxes
75 std::vector<double> bboxes(2*SPACEDIM*numSrcElems);
76 for(ConnType i = 0; i < numSrcElems ; ++i)
78 // get source bboxes in right order
79 srcElems[i]->getBoundingBox()->toCompactData(bboxes.data()+6*i);
82 adjustBoundingBoxes(bboxes);
83 const double *bboxPtr(nullptr);
86 BBTree<SPACEDIM,ConnType> tree(bboxPtr, nullptr, 0, numSrcElems);
87 const ConnType *trgConnPtr(targetMesh.getConnectivityPtr()),*trgConnIPtr(targetMesh.getConnectivityIndexPtr());
88 const ConnType *srcConnPtr(srcMesh.getConnectivityPtr()),*srcConnIPtr(srcMesh.getConnectivityIndexPtr());
89 const double *trgCooPtr(targetMesh.getCoordinatesPtr()),*srcCooPtr(srcMesh.getCoordinatesPtr());
90 for(ConnType i = 0; i < numTargetElems; ++i)
92 IKAssert(trgConnIPtr[i+1]==i+1 && trgConnIPtr[i]==i);
93 std::vector<ConnType> srcSegCondidates;
94 const double *trgCellPosition(trgCooPtr+SPACEDIM*trgConnPtr[i]);
95 typename MatrixType::value_type& resRow(result[trgConnPtr[i]]);
96 tree.getElementsAroundPoint(trgCellPosition, srcSegCondidates);
97 for(auto srcSeg: srcSegCondidates)
99 IKAssertMsg(srcConnIPtr[srcSeg+1]==2*(srcSeg+1) && srcConnIPtr[srcSeg]==2*srcSeg,"Only implemented for linear 1D source");
100 double bc0(0.),bc1(0.);
101 ConnType srcNode0(srcConnPtr[2*srcSeg]),srcNode1(srcConnPtr[2*srcSeg+1]);
102 if(IsPointOn3DSeg(srcCooPtr+SPACEDIM*srcNode0,srcCooPtr+SPACEDIM*srcNode1,trgCellPosition,epsilon,bc0,bc1))
104 resRow.insert(std::make_pair(srcNode0,bc0));
105 resRow.insert(std::make_pair(srcNode1,bc1));
110 for(ConnType i = 0 ; i < numSrcElems ; ++i)
112 return srcMesh.getNumberOfNodes();