X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FMEDCoupling_Swig%2FMEDCouplingBasicsTest7.py;h=4ce8081f213b86d2aea742c5c866ac6393ca3ba1;hb=1b746b38f3cdeae6654a9501f37fde5e56e59288;hp=6b3955a860bc823c87e0f2e7756ae3f645b8a1ce;hpb=b307fa3ee9c6d9e08082e2ccc832b28a17fd6d2c;p=tools%2Fmedcoupling.git diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py index 6b3955a86..4ce8081f2 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py @@ -1,5 +1,5 @@ # -*- coding: utf-8 -*- -# Copyright (C) 2007-2021 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 @@ -20,10 +20,7 @@ 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]) @@ -764,7 +767,7 @@ class MEDCouplingBasicsTest7(unittest.TestCase): arr1 = DataArrayInt([0,1,1,2,2,3,4,4,5,5,5,11]) self.assertTrue(DataArrayInt.FindPermutationFromFirstToSecondDuplicate(arr0,arr1).isEqual(DataArrayInt([8,5,3,1,6,9,4,2,0,11,10,7]))) self.assertTrue(DataArrayInt.FindPermutationFromFirstToSecondDuplicate(arr1,arr0).isEqual(DataArrayInt([8,3,7,2,6,1,4,11,0,5,10,9]))) - + def testDAIIndexOfSameConsecutiveValueGroups(self): arr = DataArrayInt([0,1,1,2,2,3,4,4,5,5,5,11]) self.assertTrue(arr.indexOfSameConsecutiveValueGroups().isEqual(DataArrayInt([0,1,3,5,6,8,11,12]))) @@ -778,14 +781,14 @@ class MEDCouplingBasicsTest7(unittest.TestCase): self.assertTrue(sk2.getValuesArray().isEqual(arr)) self.assertTrue(sk2.getIndexArray().isEqual(DataArrayInt([0,13,16,37,84]))) - def testSkyLineUniqueNotSortedByPack(self): + def testSkyLineUniqueNotSortedByPack(self): arrI = DataArrayInt([0,3,9,15,18,24,36,48,54]) arr = DataArrayInt([1,4,5,0,4,5,2,5,6,3,6,7,1,5,6,2,6,7,0,1,5,5,8,9,0,1,4,6,9,10,1,2,4,6,8,9,2,3,5,7,9,10,1,2,5,7,10,11,2,3,6,6,10,11]) sk = MEDCouplingSkyLineArray(arrI,arr) sk2 = sk.uniqueNotSortedByPack() self.assertTrue(sk2.getIndexArray().isEqual(DataArrayInt([0,3,8,13,16,21,29,37,42]))) self.assertTrue(sk2.getValuesArray().isEqual(DataArrayInt([1,4,5,0,2,4,5,6,1,3,5,6,7,2,6,7,0,1,5,8,9,0,1,2,4,6,8,9,10,1,2,3,5,7,9,10,11,2,3,6,10,11]))) - + def testSkyLineAggregatePacks1(self): arr = DataArrayDouble(3) ; arr.iota() m = MEDCouplingCMesh() ; m.setCoords(arr,arr) ; m = m.buildUnstructured() @@ -884,7 +887,7 @@ class MEDCouplingBasicsTest7(unittest.TestCase): self.assertTrue(lsk.getIndexArray().isEqual(DataArrayInt([0, 3, 7, 9]))) self.assertTrue(rsk.getValuesArray().isEqual(DataArrayInt([11, 12, 13, 14, 15, 16, 17, 18, 19]))) self.assertTrue(rsk.getIndexArray().isEqual(DataArrayInt([0, 3, 7, 9]))) - + def testPenta18GaussNE(self): conn = [1,0,2,4,3,5,6,7,8,9,13,14,11,10,15,12,17,16] coo = DataArrayDouble([(27.237499999999997, 9.8, 0.0), (26.974999999999994, 9.8, 0.0), (27.111517409545634, 9.532083869948877, 0.0), (27.237499999999997, 9.8, 0.5000000000000001), (26.974999999999994, 9.8, 0.5000000000000002), (27.111517409545634, 9.532083869948877, 0.5), (27.106249999999996, 9.8, 0.0), (27.17450870477282, 9.666041934974439, 0.0), (27.04325870477281, 9.666041934974439, 0.0), (27.106249999999996, 9.8, 0.5000000000000001), (27.237499999999997, 9.8, 0.25), (26.974999999999994, 9.8, 0.2500000000000001), (27.106249999999996, 9.8, 0.2500000000000001), (27.174508704772816, 9.666041934974439, 0.5), (27.043258704772814, 9.666041934974439, 0.5000000000000001), (27.111517409545634, 9.532083869948877, 0.25), (27.043258704772818, 9.666041934974436, 0.25000000000000006), (27.174508704772816, 9.666041934974436, 0.25)]) @@ -956,7 +959,518 @@ class MEDCouplingBasicsTest7(unittest.TestCase): self.assertEqual(arr3.getInfoOnComponents(),comps) self.assertTrue(arr3.isEqual(arr)) - pass + def testComputeMeshCenterOfMass0(self): + #2D + arr = DataArrayDouble(5) ; arr.iota() + m = MEDCouplingCMesh() ; m.setCoords(arr,arr) ; m=m.buildUnstructured() + self.assertTrue( m.computeMeshCenterOfMass().isEqual(DataArrayDouble([2,2],1,2),1e-12) ) + #3D + m = MEDCouplingCMesh() ; m.setCoords(arr,arr,arr) ; m=m.buildUnstructured() + self.assertTrue( m.computeMeshCenterOfMass().isEqual(DataArrayDouble([2,2,2],1,3),1e-12) ) + + def testBugPenta15_0(self): + """ + Non regression test from Roberto Da Via pointing error in connectivity of 5th sub face of penta15 cell. + """ + coo=DataArrayDouble([ + (0,1,1),(0,0,1),(1,0,1), + (0,1,0),(0,0,0),(1,0,0), + (0,0.5,1),(0.5,0,1),(0.5,0.5,1), + (0,0.5,0),(0.5,0,0),(0.5,0.5,0), + (0,1,0.5),(0,0,0.5),(1,0,0.5) + ]) + + m = MEDCouplingUMesh("penta15",3) + m.setCoords(coo) + m.allocateCells() + m.insertNextCell(NORM_PENTA15,list(range(15))) + bm = m.buildBoundaryMesh(True) + bm.writeVTK("boundary.vtu") + conn_expected = [ + [6, 0, 1, 2, 6, 7, 8], + [6, 3, 5, 4, 11, 10, 9], + [8, 0, 3, 4, 1, 12, 9, 13, 6], + [8, 1, 4, 5, 2, 13, 10, 14, 7], + [8, 2, 5, 3, 0, 14, 11, 12, 8] # old = [8, 2, 4, 5, 0, 14, 11, 12, 8] + ] + 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 testComputeTriangleHeight1(self): + m = MEDCouplingUMesh("mesh",2) + m.setCoords(DataArrayDouble([(0,0,0),(0,0,0),(10,0,0)])) + m.allocateCells() + m.insertNextCell(NORM_TRI3, [0,1,2]) + m = MEDCoupling1SGTUMesh(m) + res = m.computeTriangleHeight() + expected = DataArrayDouble([(10,0,0)]) + self.assertTrue( res.isEqual(expected,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]))) + + def testGetCellContainingPointOnPolyhedronWithPlanarFace(self): + """ + See CEA spns #40783 + In case of polyhedron with a face defined by several colinear points, + we must use other non colinear points to be able to define a face from these three points + to define the relative position of the test point to this face + """ + eps = 1.0e-12 + coo = DataArrayDouble( [ (0.176, 0.1125, 1.05), + (0.176, 0.120375, 1.05), + (0.176, 0.120375, 1.0), + (0.176, 0.1125, 1.0), + (0.176000000000000018, 0.12825, 1.05), + (0.176000000000000018, 0.12825, 1.0), + (0.207, 0.1125, 1.05), + (0.207, 0.1125, 1.0), + (0.207, 0.12825, 1.05), + (0.207, 0.12825, 1.0)] ) + + m = MEDCouplingUMesh("Mesh",3) + m.setCoords(coo) + m.allocateCells() + # put -1 to separate faces connectivity + # substract -1 from mdump table ids + m.insertNextCell(NORM_POLYHED,[0, 1, 2, 3, -1, + 1, 4, 5, 2, -1, + 6, 7, 9, 8, -1, + 3, 7, 6, 0, -1, + 9, 5, 4, 8, -1, + 3, 2, 5, 9, 7, -1, # PB in this order + #7, 3, 2, 5, 9, -1, # OK in this order + 1, 0, 6, 8, 4]) + + # test point inside the box + pt_above = (0.2, 0.12, 1.07) + pt_below = (0.2, 0.12, 0.9) + pt_inside = (0.2, 0.12, 1.025) + pts = DataArrayDouble([pt_above, pt_below, pt_inside]) + a,b = m.getCellsContainingPoints(pts, eps) + self.assertTrue(a.isEqual(DataArrayInt([0]))) + # only the third point is inside + self.assertTrue(b.isEqual(DataArrayInt([0,0,0,1]))) + + # rotate the mesh to see if getCellsContainingPoints works + # even if point is not inside bounding box + center=coo[0] + vector=[1.,0.,0.] + m.rotate(center,vector,-pi/4.); + + # test 3 points: above, below and inside + pt_above = (0.19, 0.09, 1.04) + pt_below = (0.19, 0.11, 1.02) + pt_inside = (0.19, 0.10, 1.02) + pts_rotated = DataArrayDouble([pt_above, pt_below, pt_inside]) + + a,b = m.getCellsContainingPoints(pts_rotated, eps) + self.assertTrue(a.isEqual(DataArrayInt([0]))) + # only the third point is inside + self.assertTrue(b.isEqual(DataArrayInt([0,0,0,1]))) + + def testGetCellContainingPointOnPolyhedronWithPlanarFaceWithManyNodes(self): + """ + Similar test with many colinear nodes on the planar face + """ + eps = 1.0e-12 + coo = DataArrayDouble( [(0.176000000000000018, 0.120375, 1.0 ), + (0.176000000000000018, 0.128250, 1.0 ), + (0.176000000000000018, 0.136125, 1.0 ), + (0.176000000000000018, 0.144, 1.0 ), + (0.176000000000000018, 0.151875, 1.0 ), + (0.176000000000000018, 0.159750, 1.0 ), + (0.176000000000000018, 0.167625, 1.0 ), + (0.176000000000000018, 0.1755, 1.0 ), + (0.176000000000000018, 0.183375, 1.0 ), + (0.176000000000000018, 0.191250, 1.0 ), + (0.176000000000000018, 0.199125, 1.0 ), + (0.176, 0.207, 1.0 ), + (0.207, 0.207, 1.0 ), + (0.176, 0.1125, 1.0 ), + (0.207, 0.1125, 1.0 ), + (0.176, 0.120375, 1.05), + (0.176000000000000018, 0.128250, 1.05), + (0.176000000000000018, 0.136125, 1.05), + (0.176000000000000018, 0.144, 1.05), + (0.176000000000000018, 0.151875, 1.05), + (0.176000000000000018, 0.159750, 1.05), + (0.176000000000000018, 0.167625, 1.05), + (0.176000000000000018, 0.1755, 1.05), + (0.176000000000000018, 0.183375, 1.05), + (0.176000000000000018, 0.191250, 1.05), + (0.176000000000000018, 0.199125, 1.05), + (0.176, 0.207, 1.05), + (0.207, 0.207, 1.05), + (0.176, 0.1125, 1.05), + (0.207, 0.1125, 1.05)]) + + m = MEDCouplingUMesh("Mesh",3) + m.setCoords(coo) + m.allocateCells() + # put -1 to separate faces connectivity + # substract -1 from mdump table ids + m.insertNextCell(NORM_POLYHED, + [13, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, -1, #1 + 29, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 28, -1, #2 + 14, 29, 28, 13, -1, #3 + 11, 26, 27, 12, -1, #4 + 12, 27, 29, 14, -1, #5 + 13, 28, 15, 0, -1, #6 + 0, 15, 16, 1, -1, #7 + 1, 16, 17, 2, -1, #8 + 2, 17, 18, 3, -1, #9 + 3, 18, 19, 4, -1, #10 + 4, 19, 20, 5, -1, #11 + 5, 20, 21, 6, -1, #12 + 6, 21, 22, 7, -1, #13 + 7, 22, 23, 8, -1, #14 + 8, 23, 24, 9, -1, #15 + 9, 24, 25, 10, -1, #16 + 10, 25, 26, 11] ) + + ##### See CEA 40783: error with polyhedrons (box split by on edge on its face) + pt_above = (0.1915, 0.15975, 1.07) + pt_below = (0.1915, 0.15975, 0.9) + pt_inside = (0.1915, 0.15975, 1.025) + pts = DataArrayDouble([pt_above, pt_below, pt_inside]) + a,b = m.getCellsContainingPoints(pts,eps) + self.assertTrue(a.isEqual(DataArrayInt([0]))) + # only the third point is inside + self.assertTrue(b.isEqual(DataArrayInt([0,0,0,1]))) + if __name__ == '__main__': unittest.main()