Salome HOME
Fix computation height of isocel triangle with base equal zero : NaN
[tools/medcoupling.git] / src / MEDCoupling_Swig / MEDCouplingBasicsTest7.py
index 7af9a27e01def0d3cc98fd32a2a5441c620230f7..c45f74141b3e8a558b556b8bfa4ea15a985a826d 100644 (file)
@@ -1,5 +1,5 @@
 #  -*- coding: utf-8 -*-
-# Copyright (C) 2007-2021  CEA/DEN, EDF R&D
+# Copyright (C) 2007-2023  CEA, EDF
 #
 # This library is free software; you can redistribute it and/or
 # modify it under the terms of the GNU Lesser General Public
 
 
 import sys
-if sys.platform == "win32":
-    from MEDCouplingCompat import *
-else:
-    from medcoupling import *
+from medcoupling import *
 import unittest
 from math import pi,e,sqrt,cos,sin
 from datetime import datetime
@@ -452,9 +449,9 @@ class MEDCouplingBasicsTest7(unittest.TestCase):
         # non regression test in python wrapping
         rg=DataArrayInt64([0,10,29,56,75,102,121,148,167,194,213,240,259,286,305,332,351,378,397,424,443,470,489,516])
         a,b,c=DataArrayInt64([75]).splitByValueRange(rg)
-        assert(a.isEqual(DataArrayInt64([4])))
-        assert(b.isEqual(DataArrayInt64([0])))
-        assert(c.isEqual(DataArrayInt64([4])))
+        self.assertTrue(a.isEqual(DataArrayInt64([4])))
+        self.assertTrue(b.isEqual(DataArrayInt64([0])))
+        self.assertTrue(c.isEqual(DataArrayInt64([4])))
         pass
 
     def testDAIBuildExplicitArrByRanges1(self):
@@ -665,6 +662,12 @@ class MEDCouplingBasicsTest7(unittest.TestCase):
         self.assertRaises(InterpKernelException,a.indicesOfSubPart,f) # 12 in f does not exist in a
         pass
 
+    def testDAIsortToHaveConsecutivePairs(self):
+        dref=DataArrayInt64([(6, 216), (216, 218), (218, 220), (220, 222), (222, 224), (224, 226)])
+        dtest=DataArrayInt64([(6, 216), (218, 216), (224, 226), (222, 220), (218, 220), (222, 224)])
+        dtest.sortToHaveConsecutivePairs()
+        self.assertTrue(dtest.isEqual(dref))
+
     def testDAIFromLinkedListOfPairToList1(self):
         d=DataArrayInt64([(5,7),(7,3),(3,12),(12,17)])
         zeRes=DataArrayInt64([5,7,3,12,17])
@@ -992,6 +995,341 @@ class MEDCouplingBasicsTest7(unittest.TestCase):
         ]
         self.assertTrue( bm.getNodalConnectivity().isEqual(DataArrayInt(sum(conn_expected,[]))) )
 
+    def testBugWithPolyhedInterpWithMoreThan255Nodes(self):
+        """
+        [EDF25207] : Check interpolation containing polyhedron with more than 255 nodes is OK at bbox computation stage
+        """
+        n = 8
+        arr = DataArrayDouble(n) ; arr.iota()
+        m = MEDCouplingCMesh()
+        m.setCoords(arr,arr,arr)
+        m = m.buildUnstructured()
+        skin = m.computeSkin()
+        skin.zipCoords()
+        # check that skin contains more than 2**8-1 node to reveal bug
+        self.assertTrue(skin.getNumberOfNodes()>255)
+        # Build a single polyhedron cell from skin
+        skin1 = MEDCoupling1SGTUMesh(skin)
+        conn = skin1.getNodalConnectivity()
+        conn.rearrange( MEDCouplingUMesh.GetNumberOfNodesOfGeometricType(skin1.getCellModelEnum()) )
+        connPolyhed = conn.changeNbOfComponents(MEDCouplingUMesh.GetNumberOfNodesOfGeometricType(skin1.getCellModelEnum())+1,-1)
+        connPolyhed.rearrange(1)
+        connPolyhed.popBackSilent()
+        meshSinglePolyhed = MEDCouplingUMesh("",3)
+        meshSinglePolyhed.allocateCells()
+        meshSinglePolyhed.insertNextCell(NORM_POLYHED,connPolyhed.getValues())
+        meshSinglePolyhed.setCoords( skin1.getCoords() )
+
+        rem = MEDCouplingRemapper()
+        rem.prepare(meshSinglePolyhed,m,"P0P0")
+        res = rem.getCrudeMatrix()
+        self.assertTrue( all([len(elt)==1 for elt in res]) )
+        self.assertTrue( all([elt[0]>0.99 and elt[0]<1.01 for elt in res]) )
+
+    @unittest.skipUnless(MEDCouplingHasNumPyBindings(),"requires numpy")
+    def testShapeFuncAndDerivative0(self):
+        """
+        Test values returned by MEDCoupling on HEXA27 element of shape function and its derivatives.
+        See https://www.code-aster.org/V2/doc/v10/fr/man_r/r3/r3.01.01.pdf
+        """
+        import numpy as np
+        
+        ref_coords_hexa27_med = [[-1.0, -1.0, -1.0], [-1.0, 1.0, -1.0], [1.0, 1.0, -1.0], [1.0, -1.0, -1.0], [-1.0, -1.0, 1.0], [-1.0, 1.0, 1.0], [1.0, 1.0, 1.0], [1.0, -1.0, 1.0], [-1.0, 0.0, -1.0], [0.0, 1.0, -1.0], [1.0, 0.0, -1.0], [0.0, -1.0, -1.0], [-1.0, 0.0, 1.0], [0.0, 1.0, 1.0], [1.0, 0.0, 1.0], [0.0, -1.0, 1.0], [-1.0, -1.0, 0.0], [-1.0, 1.0, 0.0], [1.0, 1.0, 0.0], [1.0, -1.0, 0.0], [0.0, 0.0, -1.0], [-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]]
+
+        def coor2index(coor):
+            zeMap = {-1.0 : 0, 0.0 : 2 , 1.0 : 1}
+            return zeMap[coor]
+
+        vcoor2index = np.vectorize( coor2index )
+        node2ijk_hexa27_med = vcoor2index( np.array(ref_coords_hexa27_med) )
+
+        def N_1d_quad(x):
+            return np.array([-0.5*x*(1-x), 0.5*x*(x+1), 1.-x*x])
+
+        def N_3d_hexa27(x, i, j, k):
+            return N_1d_quad(x[0])[i]*N_1d_quad(x[1])[j]*N_1d_quad(x[2])[k]
+
+        def N_hexa27(x):
+            return np.array([N_3d_hexa27(x, *node2ijk_hexa27_med[node,:]) for node in range(27)])
+
+        # Implementing shape function derivatives
+        def diff_N_1d_quad(x):
+            return np.array([x-0.5, x+0.5, -2.*x])
+
+        def diff_N_3d_hexa27(x, i, j, k):
+            return np.array([diff_N_1d_quad(x[0])[i]*N_1d_quad(x[1])[j]     *N_1d_quad(x[2])[k],
+                            N_1d_quad(x[0])[i]     *diff_N_1d_quad(x[1])[j]*N_1d_quad(x[2])[k],
+                            N_1d_quad(x[0])[i]     *N_1d_quad(x[1])[j]     *diff_N_1d_quad(x[2])[k]])
+
+        def diff_N_hexa27(x):
+            return np.array([diff_N_3d_hexa27(x, *node2ijk_hexa27_med[node,:]) for node in range(27)])
+        # computation of ref values
+        posInRefCoord = [-0.85685375,-0.90643355,-0.90796825]
+        ref = N_hexa27( np.array(posInRefCoord) )
+        ref2 = diff_N_hexa27( np.array(posInRefCoord) )
+        # computation using MEDCoupling
+        gl = MEDCouplingGaussLocalization(NORM_HEXA27,sum(ref_coords_hexa27_med,[]),posInRefCoord,[1])
+        mcShapeFunc = gl.getShapeFunctionValues()
+        mcShapeFunc.rearrange(1)
+        self.assertTrue( mcShapeFunc.isEqual(DataArrayDouble(ref),1e-12) )
+
+        mvDevOfShapeFunc = gl.getDerivativeOfShapeFunctionValues()
+        mvDevOfShapeFunc.rearrange(1)
+        ref2_mc = DataArrayDouble(ref2) ; ref2_mc.rearrange(1)
+        self.assertTrue( mvDevOfShapeFunc.isEqual( ref2_mc, 1e-12) )
+
+    def testShapeFuncAndDerivative1(self):
+        """
+        This test focus
+        """
+        def GetShapeFunc(ref_coord,vec):
+            gl3 = MEDCouplingGaussLocalization(gt,sum(ref_coord,[]), vec, [1])
+            funVal = gl3.getShapeFunctionValues()
+            funVal.rearrange(1)
+            return funVal
+
+        def GetDerivative(ref_coord,vec):
+            gl3 = MEDCouplingGaussLocalization(gt,sum(ref_coord,[]), vec, [1])
+            funVal = gl3.getDerivativeOfShapeFunctionValues()
+            return funVal
+        vec = [-0.85685375,-0.90643355,-0.90796825]
+        eps = 1e-6
+        # 3D cells
+        for gt in [NORM_TETRA4,NORM_TETRA10,NORM_HEXA8,NORM_PENTA6,NORM_PYRA5,NORM_PYRA13,NORM_PENTA15,NORM_PENTA6,NORM_PENTA18,NORM_HEXA20,NORM_HEXA27]: # type of cell for which derivatives are implemented
+            ref_coord = [list(elt) for elt in MEDCouplingGaussLocalization.GetDefaultReferenceCoordinatesOf(gt).getValuesAsTuple()]
+
+            der_computed = GetDerivative(ref_coord,vec)
+            der_computed.rearrange(3)
+
+            der_deduced = ( GetShapeFunc(ref_coord,[vec[0]+eps,vec[1],vec[2]])-GetShapeFunc(ref_coord,vec) ) / eps
+            delta_X = der_computed[:,0]-der_deduced
+            delta_X.abs()
+            self.assertTrue(delta_X.findIdsNotInRange(-1e-4,+1e-4).empty())
+
+            der_deduced = ( GetShapeFunc(ref_coord,[vec[0],vec[1]+eps,vec[2]])-GetShapeFunc(ref_coord,vec) ) / eps
+            delta_Y = der_computed[:,1]-der_deduced
+            delta_Y.abs()
+            self.assertTrue(delta_Y.findIdsNotInRange(-1e-5,+1e-5).empty())
+
+            der_deduced = ( GetShapeFunc(ref_coord,[vec[0],vec[1],vec[2]+eps])-GetShapeFunc(ref_coord,vec) ) / eps
+            delta_Z = der_computed[:,2]-der_deduced
+            delta_Z.abs()
+            self.assertTrue(delta_Z.findIdsNotInRange(-1e-5,+1e-5).empty())
+
+        for gt,ref_coord  in [(NORM_TETRA4,[[0.0, 1.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]]),(NORM_TETRA10,[[0.0, 1.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 0.5], [0.0, 0.5, 0.5], [0.5, 0.5, 0.0], [0.5, 0.0, 0.0], [0.5, 0.0, 0.5]]),(NORM_HEXA8,[[-1.0, -1.0, -1.0], [-1.0, 1.0, -1.0], [1.0, 1.0, -1.0], [1.0, -1.0, -1.0], [-1.0, -1.0, 1.0], [-1.0, 1.0, 1.0], [1.0, 1.0, 1.0], [1.0, -1.0, 1.0]]),(NORM_HEXA8,[[-1.0, 1.0, 0.0], [-1.0, -1.0, 0.0], [1.0, -1.0, 0.0], [1.0, 1.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]),(NORM_HEXA8,[[-1.0, -1.0, 0.0], [-1.0, 1.0, 0.0], [1.0, 1.0, 0.0], [1.0, -1.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]),(NORM_PENTA6,[[-1.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [-1.0, -0.0, 1.0], [1.0, 1.0, 0.0], [1.0, 0.0, 0.0], [1.0, 0.0, 1.0]]),(NORM_PENTA6,[[-1.0, 1.0, 0.0], [-1.0, -1.0, 0.0], [1.0, -1.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]),(NORM_PENTA6,[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]),(NORM_PYRA5,[[1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]),(NORM_PYRA13, [[1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [0.5, -0.5, 0.0], [-0.5, -0.5, 0.0], [-0.5, 0.5, 0.0], [0.5, 0.5, 0.0], [0.5, 0.0, 0.5], [0.0, -0.5, 0.5], [-0.5, 0.0, 0.5], [0.0, 0.5, 0.5]]),(NORM_PENTA15,[[-1.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [-1.0, -0.0, 1.0], [1.0, 1.0, 0.0], [1.0, 0.0, 0.0], [1.0, 0.0, 1.0], [-1.0, 0.5, 0.0], [-1.0, 0.0, 0.5], [-1.0, 0.5, 0.5], [1.0, 0.5, 0.0], [1.0, 0.0, 0.5], [1.0, 0.5, 0.5], [0.0, 1.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]),(NORM_PENTA18,[[-1.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [-1.0, -0.0, 1.0], [1.0, 1.0, 0.0], [1.0, 0.0, 0.0], [1.0, 0.0, 1.0], [-1.0, 0.5, 0.0], [-1.0, 0.0, 0.5], [-1.0, 0.5, 0.5], [1.0, 0.5, 0.0], [1.0, 0.0, 0.5], [1.0, 0.5, 0.5], [0.0, 1.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.5, 0.0], [0.0, 0.0, 0.5], [0.0, 0.5, 0.5]]),(NORM_HEXA20,[[-1.0, -1.0, -1.0], [-1.0, 1.0, -1.0], [1.0, 1.0, -1.0], [1.0, -1.0, -1.0], [-1.0, -1.0, 1.0], [-1.0, 1.0, 1.0], [1.0, 1.0, 1.0], [1.0, -1.0, 1.0], [-1.0, 0.0, -1.0], [0.0, 1.0, -1.0], [1.0, 0.0, -1.0], [0.0, -1.0, -1.0], [-1.0, 0.0, 1.0], [0.0, 1.0, 1.0], [1.0, 0.0, 1.0], [0.0, -1.0, 1.0], [-1.0, -1.0, 0.0], [-1.0, 1.0, 0.0], [1.0, 1.0, 0.0], [1.0, -1.0, 0.0]])]: # type of cell for which derivatives are implemented
+            
+            der_computed = GetDerivative(ref_coord,vec)
+            der_computed.rearrange(3)
+
+            der_deduced = ( GetShapeFunc(ref_coord,[vec[0]+eps,vec[1],vec[2]])-GetShapeFunc(ref_coord,vec) ) / eps
+            delta_X = der_computed[:,0]-der_deduced
+            delta_X.abs()
+            self.assertTrue(delta_X.findIdsNotInRange(-1e-4,+1e-4).empty())
+
+            der_deduced = ( GetShapeFunc(ref_coord,[vec[0],vec[1]+eps,vec[2]])-GetShapeFunc(ref_coord,vec) ) / eps
+            delta_Y = der_computed[:,1]-der_deduced
+            delta_Y.abs()
+            self.assertTrue(delta_Y.findIdsNotInRange(-1e-5,+1e-5).empty())
+
+            der_deduced = ( GetShapeFunc(ref_coord,[vec[0],vec[1],vec[2]+eps])-GetShapeFunc(ref_coord,vec) ) / eps
+            delta_Z = der_computed[:,2]-der_deduced
+            delta_Z.abs()
+            self.assertTrue(delta_Z.findIdsNotInRange(-1e-5,+1e-5).empty())
+
+        # 2D cells
+        vec = [0.64,0.2]
+
+        for gt in [NORM_QUAD4,NORM_QUAD8,NORM_QUAD9,NORM_TRI3,NORM_TRI6,NORM_TRI7]:
+            ref_coord = [list(elt) for elt in MEDCouplingGaussLocalization.GetDefaultReferenceCoordinatesOf(gt).getValuesAsTuple()]
+
+            der_computed = GetDerivative(ref_coord,vec)
+            der_computed.rearrange(2)
+
+            der_deduced = ( GetShapeFunc(ref_coord,[vec[0]+eps,vec[1]])-GetShapeFunc(ref_coord,vec) ) / eps
+            delta_X = der_computed[:,0]-der_deduced
+            delta_X.abs()
+            self.assertTrue(delta_X.findIdsNotInRange(-1e-5,+1e-5).empty())
+
+            der_deduced = ( GetShapeFunc(ref_coord,[vec[0],vec[1]+eps])-GetShapeFunc(ref_coord,vec) ) / eps
+            delta_Y = der_computed[:,1]-der_deduced
+            delta_Y.abs()
+            self.assertTrue(delta_Y.findIdsNotInRange(-1e-4,+1e-4).empty())
+
+        # B version of TRI6, QUAD4 and QUAD8
+        for gt,ref_coord in [(NORM_TRI3,[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]),(NORM_TRI6,[[0., 0.], [1., 0.], [0., 1.], [0.5, 0.], [0.5, 0.5], [0., 0.5]]),
+            (NORM_QUAD4,[[-1., -1.], [1., -1.], [1., 1.], [-1., 1.]]),(NORM_QUAD4,[[-1., -1.], [-1., 1.], [1., 1.], [1., -1.]]),(NORM_QUAD4,[[-1., 0.], [1., 0.], [0., 0.], [0., 0.]]),(NORM_QUAD8,[[-1., -1.], [1., -1.], [1., 1.], [-1., 1.], [0., -1.], [1., 0.], [0., 1.], [-1., 0.]])]:
+            der_computed = GetDerivative(ref_coord,vec)
+            der_computed.rearrange(2)
+
+            der_deduced = ( GetShapeFunc(ref_coord,[vec[0]+eps,vec[1]])-GetShapeFunc(ref_coord,vec) ) / eps
+            delta_X = der_computed[:,0]-der_deduced
+            delta_X.abs()
+            self.assertTrue(delta_X.findIdsNotInRange(-1e-5,+1e-5).empty())
+
+            der_deduced = ( GetShapeFunc(ref_coord,[vec[0],vec[1]+eps])-GetShapeFunc(ref_coord,vec) ) / eps
+            delta_Y = der_computed[:,1]-der_deduced
+            delta_Y.abs()
+            self.assertTrue(delta_Y.findIdsNotInRange(-1e-5,+1e-5).empty())
+            
+        # 1D cells
+        vec = [0.64]
+
+        for gt in [NORM_SEG2,NORM_SEG3,NORM_SEG4]:
+            ref_coord = [list(elt) for elt in MEDCouplingGaussLocalization.GetDefaultReferenceCoordinatesOf(gt).getValuesAsTuple()]
+
+            der_computed = GetDerivative(ref_coord,vec)
+            der_computed.rearrange(1)
+
+            der_deduced = ( GetShapeFunc(ref_coord,[vec[0]+eps])-GetShapeFunc(ref_coord,vec) ) / eps
+            delta_X = der_computed[:,0]-der_deduced
+            delta_X.abs()
+            self.assertTrue(delta_X.findIdsNotInRange(-1e-5,+1e-5).empty())
+        
+        #B version of SEG2
+        for gt,ref_coord in [(NORM_SEG2,[[0.], [1.]])]:
+            der_computed = GetDerivative(ref_coord,vec)
+            der_computed.rearrange(1)
+
+            der_deduced = ( GetShapeFunc(ref_coord,[vec[0]+eps])-GetShapeFunc(ref_coord,vec) ) / eps
+            delta_X = der_computed[:,0]-der_deduced
+            delta_X.abs()
+            self.assertTrue(delta_X.findIdsNotInRange(-1e-5,+1e-5).empty())
+        
+
+    def testComputeTriangleHeight0(self):
+        arr = DataArrayDouble([0,1])
+        m = MEDCouplingCMesh() ; m.setCoords(arr,arr)
+        m = m.buildUnstructured() ; m.simplexize(0) ; m = MEDCoupling1SGTUMesh(m)
+        res = m.computeTriangleHeight()
+        expected = DataArrayDouble([(1.0, 1.0, sqrt(2)/2.0), (sqrt(2)/2.0, 1.0, 1.0)])
+        self.assertTrue( res.isEqual(expected,1e-12) )
+        m.changeSpaceDimension(3,100)
+        res2 = m.computeTriangleHeight()
+        self.assertTrue( res2.isEqual(expected,1e-12) )
+        expected2 = DataArrayDouble([sqrt(2)/2.0, sqrt(2)/2.0])
+        self.assertTrue( res2.minPerTuple().isEqual(expected2,1e-12) )
+
+    def testDAILocateComponentId0(self):
+        arr = DataArrayInt( [(0, 1, 2), (3, 4, 5), (6, 2, 3), (7, 8, 9), (9, 0, 10), (11, 12, 13), (14, 5, 11), (15, 16, 17)] )
+        valToSearchIntoTuples = DataArrayInt( [1, 4, 6, 8, 10, 12, 14, 16, 17] )
+        tupleIdHint = DataArrayInt( [0, 1, 2, 3, 4, 5, 6, 7, 7] )
+        ret = arr.locateComponentId( valToSearchIntoTuples, tupleIdHint )
+        self.assertTrue( ret.isEqual(DataArrayInt([1, 1, 0, 1, 2, 1, 0, 1, 2]) ) )
+        pass
+
+    def testMeasureOnGaussPtMeshDimNotEqualSpaceDim0(self):
+        """
+        [EDF26877] : This test focuses on computation of measure field on field on Gauss Point in the special case where SpaceDim
+        are not eqaul to the meshDim.
+        """
+        seg2 = MEDCouplingUMesh("mesh",1)
+        seg2.setCoords(DataArrayDouble([(3,3),(4,4)]))
+        seg2.allocateCells()
+        seg2.insertNextCell(NORM_SEG2,[0,1])
+        fff=MEDCouplingFieldDouble.New(ON_GAUSS_PT) ; fff.setName("CH1RB") ; fff.setNature(IntensiveMaximum)
+        fff.setMesh(seg2)
+        fff.setGaussLocalizationOnCells([0], [0.,1.], [0.333333333333333], [1.0])
+        disc = fff.getDiscretization()
+        # spaceDim = 2 meshDim = 1
+        self.assertTrue( disc.getMeasureField(seg2,True).getArray().isEqual(DataArrayDouble([sqrt(2.0)]),1e-10) )
+        # spaceDim = 3 meshDim = 1
+        seg2.setCoords(DataArrayDouble([(3,3,3),(4,4,4)]))
+        self.assertTrue( disc.getMeasureField(seg2,True).getArray().isEqual(DataArrayDouble([sqrt(3.0)]),1e-10) )
+        # spaceDim = 3 meshDim = 2
+        tri = MEDCouplingUMesh("mesh",2)
+        tri.setCoords( DataArrayDouble([(0,0,0),(1,1,0),(2,2,2)]) )
+        tri.allocateCells()
+        tri.insertNextCell(NORM_TRI3,[0,1,2])
+        fff=MEDCouplingFieldDouble.New(ON_GAUSS_PT) ; fff.setName("CH1RB") ; fff.setNature(IntensiveMaximum)
+        fff.setMesh(tri)
+        fff.setGaussLocalizationOnCells(list(range(0, 1)), [0., 0., 1., 0., 0., 1.], [0.3333333333333333, 0.3333333333333333], [0.5])
+        disc = fff.getDiscretization()
+        self.assertTrue( disc.getMeasureField(tri,True).getArray().isEqual( tri.getMeasureField(True).getArray(), 1e-10) )
+        pass
+
+    def testUMeshExplodeMeshTo(self):
+        """
+        [EDF27988] : implementation of reduceToCells implies implementation of MEDCouplingUMesh.explodeMeshTo
+        """
+        arr = DataArrayDouble(5) ; arr.iota()
+        m = MEDCouplingCMesh() ; m.setCoords(arr,arr,arr)
+        m = m.buildUnstructured()
+        m1 = m[::2] ; m2 = m[1::2]
+        m1.simplexize(PLANAR_FACE_5)
+        m = MEDCouplingUMesh.MergeUMeshesOnSameCoords([m1,m2])
+        mE1 = m.explodeMeshTo(-1)
+        ref1 = m.buildDescendingConnectivity()
+        mE2 = m.explodeMeshTo(-2)
+        ref2 = m.explode3DMeshTo1D()
+        mE3 = m.explodeMeshTo(-3)
+        self.assertTrue( len(mE1) ==5 )
+        self.assertTrue( mE1[0].getNodalConnectivity().isEqual(ref1[0].getNodalConnectivity()) )
+        self.assertTrue( mE1[0].getNodalConnectivityIndex().isEqual(ref1[0].getNodalConnectivityIndex()) )
+        self.assertTrue( mE1[0].getCoords().getHiddenCppPointer() == m.getCoords().getHiddenCppPointer() )
+        for i in range(1,5):
+            self.assertTrue( mE1[i].isEqual(ref1[i]) )
+        #
+        self.assertTrue( len(mE2) ==5 )
+        self.assertTrue( mE2[0].getNodalConnectivity().isEqual(ref2[0].getNodalConnectivity()) )
+        self.assertTrue( mE2[0].getNodalConnectivityIndex().isEqual(ref2[0].getNodalConnectivityIndex()) )
+        self.assertTrue( mE2[0].getCoords().getHiddenCppPointer() == m.getCoords().getHiddenCppPointer() )
+        for i in range(1,5):
+            self.assertTrue( mE2[i].isEqual(ref2[i]) )
+        #
+        self.assertTrue( mE3[0].getMeshDimension() == 0 )
+        self.assertTrue( mE3[0].getNumberOfCells() == mE3[0].getNumberOfNodes() )
+        a,b = m.getReverseNodalConnectivity()
+        self.assertTrue( mE3[3].isEqual(a) and mE3[4].isEqual(b) )
+        ref3_2 = (m.getNodalConnectivityIndex().deltaShiftIndex()-1) ; ref3_2.computeOffsetsFull()
+        self.assertTrue( ref3_2.isEqual(mE3[2]) )
+        tmp = m.getNodalConnectivityIndex().deepCopy() ; tmp.popBackSilent() ; tmp = tmp.buildComplement( len(m.getNodalConnectivity()) ) ; ref3_1 = m.getNodalConnectivity()[tmp]
+        self.assertTrue( ref3_1.isEqual(mE3[1]) )
+        #
+        cellsInPolyh = [37,160]
+        polyh = m[cellsInPolyh]
+        polyh.convertAllToPoly()
+        m[cellsInPolyh] = polyh
+        pE3 = m.explodeMeshTo(-3)
+        self.assertTrue( pE3[0].getMeshDimension() == 0 )
+        self.assertTrue( pE3[0].getNumberOfCells() == pE3[0].getNumberOfNodes() )
+        a,b = m.getReverseNodalConnectivity()
+        self.assertTrue( pE3[3].isEqual(a) and pE3[4].isEqual(b) )
+        self.assertTrue( pE3[2].isEqual(mE3[2]) ) # indexed arrays are the same
+
+        ref_a,ref_b = DataArrayInt.ExtractFromIndexedArrays( DataArrayInt(cellsInPolyh).buildComplement(m.getNumberOfCells()), mE3[1], mE3[2] )
+        a,b = DataArrayInt.ExtractFromIndexedArrays( DataArrayInt(cellsInPolyh).buildComplement(m.getNumberOfCells()), pE3[1], pE3[2] )
+        self.assertTrue( ref_a.isEqual(a) )
+        self.assertTrue( ref_b.isEqual(b) )
+        for cell in cellsInPolyh:
+            ref_c,ref_d = DataArrayInt.ExtractFromIndexedArrays( cell, mE3[1], mE3[2] ) ; ref_c.sort()
+            c,d = DataArrayInt.ExtractFromIndexedArrays( cell, pE3[1], pE3[2] )
+            self.assertTrue( ref_c.isEqual(c) )
+            self.assertTrue( ref_d.isEqual(d) )
+    
+    def testGetCellContainingPointRelativeEps(self):
+        """
+        See EDF27860 : This test checks that detection of point inside a cell works by normalizing cell around origin with factor equal to the max delta of bbox along axis X, Y or Z.
+        """
+        # in this test cell is vuluntary far from origin {15260.775604514516, 11197.646906189217, 14187.820484060947}
+        # and caracteritic size is ~ 1500
+        coo = DataArrayDouble( [(14724.199858870656, 11928.888084722483, 14442.32726944039), (14788.407409534622, 11992.60694822231, 14453.86181555231), (15572.175148726046, 10798.586790270576, 14471.54225356788), (15643.898717334796, 10853.094666047728, 14477.233802854305), (15005.31495255754, 11573.261110174888, 13933.313698681504), (15070.29423166349, 11636.377758513776, 13946.650959030132), (15797.351350158377, 10466.40572765595, 13965.524190108257), (15869.808770928525, 10519.99285973948, 13972.419352086607), (15273.866774426142, 11216.458197414971, 13433.169979717744), (15340.421031616577, 11277.882145177837, 13446.53598386297), (16013.382514001762, 10132.795887638129, 13465.184281842226), (16086.979064572806, 10184.802292369684, 13472.147425473782)] )
+        m = MEDCouplingUMesh("",3)
+        m.setCoords(coo)
+        m.allocateCells()
+        m.insertNextCell(NORM_TETRA4,[0,5,4,6])
+        m.insertNextCell(NORM_TETRA4,[4,5,9,7])
+
+        ##### See EDF2760 pt is outside cell 0 (6e-4) and 1 (8e-4)
+        pt = DataArrayDouble([(15263.41200205526, 11314.957094727113, 13950.0)])
+        a,b = m.getCellsContainingPoints(pt,1e-3)
+        self.assertTrue(a.isEqual(DataArrayInt([0,1])))
+        self.assertTrue(b.isEqual(DataArrayInt([0,2])))
+
+        # by shifting pt by 10 along Z pt in only inside cell # 0
+        pt += [0,0,10]
+        a1,b1 = m.getCellsContainingPoints(pt,1e-3)
+        self.assertTrue(a1.isEqual(DataArrayInt([0])))
+        self.assertTrue(b1.isEqual(DataArrayInt([0,1])))
 
     pass