X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FINTERP_KERNEL%2FInterpolation3D1D.txx;h=35e01c71c6605fd1c939791023d515cdbb73b04f;hb=1b746b38f3cdeae6654a9501f37fde5e56e59288;hp=15d100d1035e7726be928520ac02df8454126131;hpb=4e3de127a60f9ffb63abae1d17ed01f3f6711fd2;p=tools%2Fmedcoupling.git

diff --git a/src/INTERP_KERNEL/Interpolation3D1D.txx b/src/INTERP_KERNEL/Interpolation3D1D.txx
old mode 100644
new mode 100755
index 15d100d10..35e01c71c
--- a/src/INTERP_KERNEL/Interpolation3D1D.txx
+++ b/src/INTERP_KERNEL/Interpolation3D1D.txx
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2019  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2024  CEA, EDF
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -16,10 +16,9 @@
 //
 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
-// Author : Anthony Geay (CEA/DEN)
+// Author : Anthony Geay (EDF R&D)
 
-#ifndef __INTERPOLATION3D1D_TXX__
-#define __INTERPOLATION3D1D_TXX__
+#pragma once
 
 #include "Interpolation3D1D.hxx"
 #include "Interpolation.txx"
@@ -30,7 +29,7 @@
 #include "PointLocator3DIntersectorP1P1.txx"
 #include "Log.hxx"
 
-#include "BBTree.txx"
+#include "InterpolationHelper.txx"
 
 namespace INTERP_KERNEL
 {
@@ -39,114 +38,56 @@ namespace INTERP_KERNEL
    *  adjusted in a similar fashion as in InterpolationPlanar::performAdjustmentOfBB()
    **/
   template<class MyMeshType, class MatrixType>
-  int Interpolation3D1D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method)
+  typename MyMeshType::MyConnType Interpolation3D1D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const std::string& method)
   {
     if(InterpolationOptions::getIntersectionType() != PointLocator)
-      INTERP_KERNEL::Exception("Invalid 3D/1D intersection type specified : must be PointLocator.");
+      INTERP_KERNEL::Exception("Invalid 3D/1D-0D intersection type specified : must be PointLocator.");
 
     typedef typename MyMeshType::MyConnType ConnType;
     // create MeshElement objects corresponding to each element of the two meshes
-    const unsigned long numSrcElems = srcMesh.getNumberOfElements();
-    const unsigned long numTargetElems = targetMesh.getNumberOfElements();
+    const ConnType numTargetElems = targetMesh.getNumberOfElements();
 
-    LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements ");
+    LOG(2, "Target mesh has " << numTargetElems << " elements ");
 
-    std::vector<MeshElement<ConnType>*> srcElems(numSrcElems);
-    std::vector<MeshElement<ConnType>*> targetElems(numTargetElems);
-
-    std::map<MeshElement<ConnType>*, int> indices;
-
-    for(unsigned long i = 0 ; i < numSrcElems ; ++i)
-      srcElems[i] = new MeshElement<ConnType>(i, srcMesh);       
-
-    for(unsigned long i = 0 ; i < numTargetElems ; ++i)
-      targetElems[i] = new MeshElement<ConnType>(i, targetMesh);
-
-    Intersector3D<MyMeshType,MatrixType>* intersector=0;
+    std::unique_ptr< Intersector3D<MyMeshType,MatrixType> > intersector;
     std::string methC = InterpolationOptions::filterInterpolationMethod(method);
     if(methC=="P0P0")
-      { intersector=new PointLocator3DIntersectorP0P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
+      { intersector.reset( new PointLocator3DIntersectorP0P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
       }
     else if(methC=="P0P1")
-      {  intersector=new PointLocator3DIntersectorP0P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
+      {  intersector.reset( new PointLocator3DIntersectorP0P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
       }
     else if(methC=="P1P0")
-      {  intersector=new PointLocator3DIntersectorP1P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
+      {  intersector.reset( new PointLocator3DIntersectorP1P0<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
       }
     else if(methC=="P1P1")
-      {  intersector=new PointLocator3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
+      {  intersector.reset( new PointLocator3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision()) );
       }
     else
       throw Exception("Invalid method chosen must be in \"P0P0\", \"P0P1\", \"P1P0\" or \"P1P1\".");
     // create empty maps for all source elements
     result.resize(intersector->getNumberOfRowsOfResMatrix());
 
-    // create BBTree structure
-    // - get bounding boxes
-    std::vector<double> bboxes(6*numSrcElems);
-    int* srcElemIdx = new int[numSrcElems];
-    for(unsigned long i = 0; i < numSrcElems ; ++i)
-      {
-        // get source bboxes in right order
-        const BoundingBox* box = srcElems[i]->getBoundingBox();
-        bboxes[6*i+0] = box->getCoordinate(BoundingBox::XMIN);
-        bboxes[6*i+1] = box->getCoordinate(BoundingBox::XMAX);
-        bboxes[6*i+2] = box->getCoordinate(BoundingBox::YMIN);
-        bboxes[6*i+3] = box->getCoordinate(BoundingBox::YMAX);
-        bboxes[6*i+4] = box->getCoordinate(BoundingBox::ZMIN);
-        bboxes[6*i+5] = box->getCoordinate(BoundingBox::ZMAX);
-
-        srcElemIdx[i] = srcElems[i]->getIndex();
-      }
-
-    adjustBoundingBoxes(bboxes);
-    const double *bboxPtr=0;
-    if(numSrcElems>0)
-      bboxPtr=&bboxes[0];
-    BBTree<3,ConnType> tree(bboxPtr, srcElemIdx, 0, numSrcElems);
+    BBTreeStandAlone<3,ConnType> tree( BuildBBTreeWithAdjustment(srcMesh,[this](double *bbox, typename MyMeshType::MyConnType sz){ this->adjustBoundingBoxes(bbox,sz); }) );
 
     // for each target element, get source elements with which to calculate intersection
     // - calculate intersection by calling intersectCells
-    for(unsigned long i = 0; i < numTargetElems; ++i)
+    for(ConnType i = 0; i < numTargetElems; ++i)
       {
-        const BoundingBox* box = targetElems[i]->getBoundingBox();
-        const int targetIdx = targetElems[i]->getIndex();
-
+        MeshElement<ConnType> trgMeshElem(i, targetMesh);
+        
+        const BoundingBox* box = trgMeshElem.getBoundingBox();
         // get target bbox in right order
         double targetBox[6];
-        targetBox[0] = box->getCoordinate(BoundingBox::XMIN);
-        targetBox[1] = box->getCoordinate(BoundingBox::XMAX);
-        targetBox[2] = box->getCoordinate(BoundingBox::YMIN);
-        targetBox[3] = box->getCoordinate(BoundingBox::YMAX);
-        targetBox[4] = box->getCoordinate(BoundingBox::ZMIN);
-        targetBox[5] = box->getCoordinate(BoundingBox::ZMAX);
+        box->fillInXMinXmaxYminYmaxZminZmaxFormat(targetBox);
 
         std::vector<ConnType> intersectElems;
 
         tree.getIntersectingElems(targetBox, intersectElems);
 
         if ( !intersectElems.empty() )
-          intersector->intersectCells(targetIdx,intersectElems,result);
+          intersector->intersectCells(i,intersectElems,result);
       }
-
-    // free allocated memory
-    delete [] srcElemIdx;
-
-    int ret=intersector->getNumberOfColsOfResMatrix();
-
-    delete intersector;
-
-    for(unsigned long i = 0 ; i < numSrcElems ; ++i)
-      {
-        delete srcElems[i];
-      }
-    for(unsigned long i = 0 ; i < numTargetElems ; ++i)
-      {
-        delete targetElems[i];
-      }
-    return ret;
-
+    return intersector->getNumberOfColsOfResMatrix();
   }
 }
-
-#endif