src/INTERP_KERNEL/Interpolation1D0D.txx

   1 // Copyright (C) 2018  CEA/DEN, EDF R&D
   2 //
   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.
   7 //
   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.
  12 //
  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
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (EDF R&D)
  20
  21 #ifndef __INTERPOLATION1D0D_TXX__
  22 #define __INTERPOLATION1D0D_TXX__
  23
  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"
  31 #include "Log.hxx"
  32
  33 #include "BBTree.txx"
  34
  35 #include "InterpKernelAssert.hxx"
  36
  37 namespace INTERP_KERNEL
  38 {
  39   /**
  40    *  Very similar to Interpolation3D::interpolateMeshes, except for the bounding boxes that can be
  41    *  adjusted in a similar fashion as in InterpolationPlanar::performAdjustmentOfBB()
  42    **/
  43   template<class MyMeshType, class MatrixType>
  44   int Interpolation1D0D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method)
  45   {
  46 #if __cplusplus >= 201103L
  47     constexpr int SPACEDIM=MyMeshType::MY_SPACEDIM;
  48     using ConnType=typename MyMeshType::MyConnType;
  49     IKAssert(SPACEDIM==3);
  50
  51     if(InterpolationOptions::getIntersectionType() != PointLocator)
  52       INTERP_KERNEL::Exception("Invalid 1D/0D intersection type specified : must be PointLocator.");
  53
  54     std::string methC ( InterpolationOptions::filterInterpolationMethod(method) );
  55     if(methC!="P1P1")
  56       throw Exception("Invalid method chosen must be in \"P1P1\".");
  57
  58     const double epsilon(getPrecision());
  59     // create MeshElement objects corresponding to each element of the two meshes
  60     const unsigned long numSrcElems(srcMesh.getNumberOfElements()), numTargetElems(targetMesh.getNumberOfElements());
  61
  62     LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements ");
  63
  64     std::vector<MeshElement<ConnType>*> srcElems(numSrcElems);
  65
  66     std::map<MeshElement<ConnType>*, int> indices;
  67
  68     for(unsigned long i = 0 ; i < numSrcElems ; ++i)
  69       srcElems[i] = new MeshElement<ConnType>(i, srcMesh);
  70
  71     // create empty maps for all source elements
  72     result.resize(targetMesh.getNumberOfNodes());
  73
  74     // create BBTree structure
  75     // - get bounding boxes
  76     std::vector<double> bboxes(2*SPACEDIM*numSrcElems);
  77     int* srcElemIdx = new int[numSrcElems];
  78     for(unsigned long i = 0; i < numSrcElems ; ++i)
  79       {
  80         // get source bboxes in right order
  81         srcElems[i]->getBoundingBox()->toCompactData(bboxes.data()+6*i);
  82         srcElemIdx[i] = srcElems[i]->getIndex();
  83       }
  84
  85     adjustBoundingBoxes(bboxes);
  86     const double *bboxPtr(nullptr);
  87     if(numSrcElems>0)
  88       bboxPtr=&bboxes[0];
  89     BBTree<SPACEDIM,ConnType> tree(bboxPtr, srcElemIdx, 0, numSrcElems);
  90     const ConnType *trgConnPtr(targetMesh.getConnectivityPtr()),*trgConnIPtr(targetMesh.getConnectivityIndexPtr());
  91     const ConnType *srcConnPtr(srcMesh.getConnectivityPtr()),*srcConnIPtr(srcMesh.getConnectivityIndexPtr());
  92     const double *trgCooPtr(targetMesh.getCoordinatesPtr()),*srcCooPtr(srcMesh.getCoordinatesPtr());
  93     for(unsigned long i = 0; i < numTargetElems; ++i)
  94       {
  95         IKAssert(trgConnIPtr[i+1]==i+1 && trgConnIPtr[i]==i);
  96         std::vector<ConnType> srcSegCondidates;
  97         const double *trgCellPosition(trgCooPtr+SPACEDIM*trgConnPtr[i]);
  98         typename MatrixType::value_type& resRow(result[trgConnPtr[i]]);
  99         tree.getElementsAroundPoint(trgCellPosition, srcSegCondidates);
 100         for(auto srcSeg: srcSegCondidates)
 101           {
 102             IKAssertMsg(srcConnIPtr[srcSeg+1]==2*(srcSeg+1) && srcConnIPtr[srcSeg]==2*srcSeg,"Only implemented for linear 1D source");
 103             double bc0(0.),bc1(0.);
 104             ConnType srcNode0(srcConnPtr[2*srcSeg]),srcNode1(srcConnPtr[2*srcSeg+1]);
 105             if(IsPointOn3DSeg(srcCooPtr+SPACEDIM*srcNode0,srcCooPtr+SPACEDIM*srcNode1,trgCellPosition,epsilon,bc0,bc1))
 106               {
 107                 resRow.insert(std::make_pair(srcNode0,bc0));
 108                 resRow.insert(std::make_pair(srcNode1,bc1));
 109                 continue;
 110               }
 111           }
 112       }
 113     delete [] srcElemIdx;
 114     for(unsigned long i = 0 ; i < numSrcElems ; ++i)
 115       delete srcElems[i];
 116     return srcMesh.getNumberOfNodes();
 117   }
 118 #else
 119   throw INTERP_KERNEL::Exception("Go buying a C++11 compiler :)");
 120 #endif
 121 }
 122
 123 #endif