#ifndef __CURVEINTERSECTOR_HXX__
#define __CURVEINTERSECTOR_HXX__
-#include "TargetIntersector.hxx"
+#include "TargetIntersector.txx"
#include "NormalizedUnstructuredMesh.hxx"
namespace INTERP_KERNEL
double precision, double adjustmentEpsAbs, double medianLine, int printLevel);
virtual ~CurveIntersector();
void createBoundingBoxes(const MyMeshType& mesh, std::vector<double>& bbox);
- void adjustBoundingBoxes(std::vector<double>& bbox, double adjustmentEpsAbs);
static void getElemBB(double* bb, const MyMeshType& mesh, ConnType iP, ConnType nb_nodes);
static void ComputeBaryCoordsOf(double startOfSeg, double endOfSeg, double pt, double& startPos, double& endPos);
protected :
startPos = std::max(startPos,0.); startPos = std::min(startPos,1.);
endPos=1.-startPos;
}
-
- /*! Readjusts a set of bounding boxes so that they are extended
- in all dimensions for avoiding missing interesting intersections
-
- \param bbox vector containing the bounding boxes
- */
- template<class MyMeshType, class MyMatrix>
- void CurveIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes (std::vector<double>& bbox,
- double adjustmentEpsAbs)
- {
- std::size_t size = bbox.size()/(2*SPACEDIM);
- for (std::size_t i=0; i<size; i++)
- {
- for(int idim=0; idim<SPACEDIM; idim++)
- {
- bbox[i*2*SPACEDIM+2*idim ] -= adjustmentEpsAbs;
- bbox[i*2*SPACEDIM+2*idim+1] += adjustmentEpsAbs;
- }
- }
- }
/*!
* @param icellT id in target mesh in format of MyMeshType.
#ifndef __INTEGRALUNIFORMINTERSECTOR_HXX__
#define __INTEGRALUNIFORMINTERSECTOR_HXX__
-#include "TargetIntersector.hxx"
+#include "TargetIntersector.txx"
#include <cmath>
#include "NormalizedUnstructuredMesh.hxx"
#include "InterpolationOptions.hxx"
#include "MCIdType.hxx"
+#include "Intersector3D.hxx"
+#include <vector>
namespace INTERP_KERNEL
{
MyMatrixType& matrix,
const std::string& method);
INTERPKERNEL_EXPORT DuplicateFacesType retrieveDuplicateFaces() const { return _duplicate_faces; }
+
+ protected:
+ template<class MyMeshType, class MyMatrixType>
+ void performAdjustmentOfBB(Intersector3D<MyMeshType,MyMatrixType>* intersector, std::vector<double>& bbox) const
+ {
+ intersector->adjustBoundingBoxes(bbox,InterpolationOptions::getBoundingBoxAdjustment(),InterpolationOptions::getBoundingBoxAdjustmentAbs());
+ }
+
private:
SplittingPolicy _splitting_policy;
DuplicateFacesType _duplicate_faces;
+
};
}
// create BBTree structure
// - get bounding boxes
- double* bboxes = new double[6 * numSrcElems];
+ std::vector<double> bboxes(6 * numSrcElems);
ConnType* srcElemIdx = new ConnType[numSrcElems];
for(ConnType i = 0; i < numSrcElems ; ++i)
{
srcElemIdx[i] = srcElems[i]->getIndex();
}
- BBTree<3,ConnType> tree(bboxes, srcElemIdx, 0, numSrcElems, 0.);
+ // [ABN] Adjust 2D bounding box (those might be flat in the cases where the 2D surf are perfectly aligned with the axis)
+ performAdjustmentOfBB(intersector, bboxes);
+
+ BBTree<3,ConnType> tree(bboxes.data(), srcElemIdx, 0, numSrcElems, 0.);
// for each target element, get source elements with which to calculate intersection
// - calculate intersection by calling intersectCells
}
- delete [] bboxes;
delete [] srcElemIdx;
DuplicateFacesType::iterator iter;
std::vector<double> bbox;
intersector->createBoundingBoxes(myMeshS,bbox); // create the bounding boxes
- intersector->adjustBoundingBoxes(bbox, InterpolationOptions::getBoundingBoxAdjustmentAbs());
+ intersector->adjustBoundingBoxes(bbox, InterpolationOptions::getBoundingBoxAdjustment(),
+ InterpolationOptions::getBoundingBoxAdjustmentAbs());
BBTree<SPACEDIM,ConnType> my_tree(&bbox[0], 0, 0, nbMailleS);//creating the search structure
long end_filtering = clock();
#ifndef __INTERSECTOR3D_HXX__
#define __INTERSECTOR3D_HXX__
-#include "TargetIntersector.hxx"
+#include "TargetIntersector.txx"
namespace INTERP_KERNEL
{
static const NumberingPolicy numPol=MyMeshType::My_numPol;
public:
Intersector3D(const MyMeshType& targetMesh, const MyMeshType& srcMesh);
+
void getRealTargetCoordinates(ConnType icellT, std::vector<double>& coordsT) const;
void getRealSourceCoordinates(ConnType icellT, std::vector<double>& coordsT) const;
const ConnType *getStartConnOfTargetCell(ConnType icellT) const;
#ifndef __IntersectorCU_HXX__
#define __IntersectorCU_HXX__
-#include "TargetIntersector.hxx"
+#include "TargetIntersector.txx"
#include "NormalizedUnstructuredMesh.hxx"
namespace INTERP_KERNEL
#ifndef __PLANARINTERSECTOR_HXX__
#define __PLANARINTERSECTOR_HXX__
-#include "TargetIntersector.hxx"
+#include "TargetIntersector.txx"
#include "NormalizedUnstructuredMesh.hxx"
#include <map>
PlanarIntersector(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, double minDot3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel);
virtual ~PlanarIntersector();
void createBoundingBoxes(const MyMeshType& mesh, std::vector<double>& bbox);
- void adjustBoundingBoxes(std::vector<double>& bbox, double surf3DAdjustmentEps, double surf3DAdjustmentEpsAbs);
inline void getElemBB(double* bb, const MyMeshType& mesh, ConnType iP, ConnType nb_nodes);
static int Projection(double *Coords_A, double *Coords_B,
ConnType nb_NodesA, ConnType nb_NodesB, double epsilon, double md3DSurf, double minDot3DSurf, double median_plane, bool do_rotate);
}
}
- /*! Readjusts a set of bounding boxes so that they are extended
- in all dimensions for avoiding missing interesting intersections
-
- \param bbox vector containing the bounding boxes
- */
- template<class MyMeshType, class MyMatrix>
- void PlanarIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes(std::vector<double>& bbox, double surf3DAdjustmentEps, double surf3DAdjustmentEpsAbs)
- {
- /* We build the segment tree for locating possible matching intersections*/
-
- std::size_t size = bbox.size()/(2*SPACEDIM);
- for (std::size_t i=0; i<size; i++)
- {
- double max=- std::numeric_limits<double>::max();
- for(int idim=0; idim<SPACEDIM; idim++)
- {
- double Dx=bbox[i*2*SPACEDIM+1+2*idim]-bbox[i*2*SPACEDIM+2*idim];
- max=(max<Dx)?Dx:max;
- }
- for(int idim=0; idim<SPACEDIM; idim++)
- {
- bbox[i*2*SPACEDIM+2*idim ] -= surf3DAdjustmentEps*max+surf3DAdjustmentEpsAbs;
- bbox[i*2*SPACEDIM+2*idim+1] += surf3DAdjustmentEps*max+surf3DAdjustmentEpsAbs;
- }
- }
- }
-
/*!
* @param icellT id in target mesh in format of MyMeshType.
* @param coordsT output val that stores coordinates of the target cell automatically resized to the right length.
class TargetIntersector
{
public:
+ static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
typedef typename MyMeshType::MyConnType ConnType;
public:
/*!
virtual ConnType getNumberOfRowsOfResMatrix() const = 0;
virtual ConnType getNumberOfColsOfResMatrix() const = 0;
+
virtual ~TargetIntersector() { }
+ void adjustBoundingBoxes(std::vector<double>& bbox, double adjustmentEps, double adjustmentEpsAbs);
};
}
--- /dev/null
+// Copyright (C) 2007-2020 CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// Author : Anthony Geay (CEA/DEN)
+
+#ifndef __TARGETINTERSECTOR__TXX__
+#define __TARGETINTERSECTOR__TXX__
+
+#include "TargetIntersector.hxx"
+#include <limits>
+
+namespace INTERP_KERNEL
+{
+ /*! Readjusts a set of bounding boxes so that they are extended
+ in all dimensions for avoiding missing interesting intersections
+
+ @param bbox vector containing the bounding boxes
+ @param adjustmentEps relative adjustment value (a percentage of the maximal BBox dimension)
+ @param adjustmentEpsAbs absolute adjustment value (added on each side of the BBox in each dimension)
+ */
+ template<class MyMeshType, class MyMatrix>
+ void TargetIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes(std::vector<double>& bbox, double adjustmentEps, double adjustmentEpsAbs)
+ {
+ std::size_t size = bbox.size()/(2*SPACEDIM);
+ for (std::size_t i=0; i<size; i++)
+ {
+ double max=- std::numeric_limits<double>::max();
+ for(int idim=0; idim<SPACEDIM; idim++)
+ {
+ double Dx=bbox[i*2*SPACEDIM+1+2*idim]-bbox[i*2*SPACEDIM+2*idim];
+ max=(max<Dx)?Dx:max;
+ }
+ for(int idim=0; idim<SPACEDIM; idim++)
+ {
+ bbox[i*2*SPACEDIM+2*idim ] -= adjustmentEps*max+adjustmentEpsAbs;
+ bbox[i*2*SPACEDIM+2*idim+1] += adjustmentEps*max+adjustmentEpsAbs;
+ }
+ }
+ }
+
+}
+
+#endif
self.assertAlmostEqual(m[1,1],12.0,12)
self.assertAlmostEqual(m[1,2],23.0,12)
self.assertEqual(m.shape,(2,8))
- pass
+
+ def test_Interpolation2D3D_bbox_adjustment_1(self):
+ """ Interpolation 2D <-> 3D was not using bounding box adjustment.
+ In case of a 2D mesh perfectly aligned with the axis, the bounding box intersection was not working properly (flat bounding box).
+ """
+ ## Source
+ meshS = MEDCouplingUMesh('SupportOf_TEMPERATURE_OUT', 2)
+ coo = cooS = DataArrayDouble([(-0.00074999999999877595,0.00000000000000000000,0.00032540000000000005),
+ (-0.00049999999999755579,0.00025000000000140708,0.00032540000000000005),(-0.00049999999999755600,0.00000000000000000000,0.00032540000000000005),
+ (-0.00100000000000000002,0.00000000000000000000,0.00032540000000000005),(-0.00100000000000000002,0.00025000000000000543,0.00032540000000000005),
+ (-0.00075651925565617829,0.00034416831541328637,0.00032540000000000005)]) # the extra 5e-20 on Z is the true culprit :-)
+ meshS.setCoords(coo)
+ c = DataArrayInt([3, 0, 1, 2, 3, 0, 3, 4, 3, 5, 0, 4, 3, 5, 1, 0])
+ cI = DataArrayInt([0, 4, 8, 12, 16])
+ meshS.setConnectivity(c, cI)
+ meshS.checkConsistency()
+ ## Target
+ meshT = MEDCouplingUMesh('IJK_mesh', 3)
+ coo = DataArrayDouble([(-0.001,0,0.000303602),(-0.0009,0,0.000303602),(-0.0008,0,0.000303602),(-0.0007,0,0.000303602),(-0.0006,0,0.000303602),
+ (-0.0005,0,0.000303602),(-0.001,0.0005,0.000303602),(-0.0009,0.0005,0.000303602),(-0.0008,0.0005,0.000303602),(-0.0007,0.0005,0.000303602),
+ (-0.0006,0.0005,0.000303602),(-0.0005,0.0005,0.000303602),(-0.001,0,0.0003254),(-0.0009,0,0.0003254),(-0.0008,0,0.0003254),(-0.0007,0,0.0003254),
+ (-0.0006,0,0.0003254),(-0.0005,0,0.0003254),(-0.001,0.0005,0.0003254),(-0.0009,0.0005,0.0003254),(-0.0008,0.0005,0.0003254),(-0.0007,0.0005,0.0003254),
+ (-0.0006,0.0005,0.0003254),(-0.0005,0.0005,0.0003254)])
+ meshT.setCoords(coo)
+ c = DataArrayInt([18, 1, 0, 6, 7, 13, 12, 18, 19, 18, 2, 1, 7, 8, 14, 13, 19, 20, 18, 3, 2, 8, 9, 15, 14, 20, 21, 18, 4, 3, 9, 10, 16, 15, 21, 22,
+ 18, 5, 4, 10, 11, 17, 16, 22, 23])
+ cI = DataArrayInt([0, 9, 18, 27, 36, 45])
+ meshT.setConnectivity(c, cI)
+ meshT.checkConsistency()
+ ## Dummy field
+ fldSrc = MEDCouplingFieldDouble(ON_CELLS, ONE_TIME)
+ fldSrc.setMesh(meshS)
+ da = DataArrayDouble(meshS.getNumberOfCells())
+ da[:] = 50.0
+ fldSrc.setArray(da)
+ remap = MEDCouplingRemapper()
+ # remap.setBoundingBoxAdjustmentAbs(1.0e-5) # was not taken into account for 2D/3D - but we don't even need it! Default value is OK.
+ remap.prepare(meshS, meshT, "P0P0")
+ fldSrc.setNature(IntensiveMaximum)
+ fldTgt = remap.transferField(fldSrc, -1.0)
+ self.assertTrue(fldTgt.getArray().isUniform(50.0, 1e-12))
def checkMatrix(self,mat1,mat2,nbCols,eps):
self.assertEqual(len(mat1),len(mat2))
