Salome HOME
Spread and condense on cell fields on images meshes with ghost management in 1D and 2D.
[tools/medcoupling.git] / src / MEDCoupling_Swig / MEDCouplingRemapperTest.py
index 9c463ca08e3b7710fa094977c2f84801de8e8dd5..44b0f8bd399674ba33bb3df3185dee4b7b90e099 100644 (file)
@@ -1,10 +1,10 @@
 #  -*- coding: iso-8859-1 -*-
-# Copyright (C) 2007-2013  CEA/DEN, EDF R&D
+# Copyright (C) 2007-2014  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.
+# 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
@@ -591,7 +591,7 @@ class MEDCouplingBasicsTest(unittest.TestCase):
     def testGetCrudeCSRMatrix1(self):
         """ testing CSR matrix output using numpy/scipy.
         """
-        from scipy.sparse import diags
+        from scipy.sparse import spdiags #diags
         import scipy
         from numpy import array
         arr=DataArrayDouble(3) ; arr.iota()
@@ -620,7 +620,8 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         self.assertEqual(diff.getnnz(),0)
         # IntegralGlobConstraint (division by sum of cols)
         colSum=m.sum(axis=0)
-        m_0=m*diags(array(1/colSum),[0])
+        # version 0.12.0 # m_0=m*diags(array(1/colSum),[0])
+        m_0=m*spdiags(array(1/colSum),[0],colSum.shape[1],colSum.shape[1])
         del colSum
         self.assertAlmostEqual(m_0[0,0],0.625,12)
         self.assertAlmostEqual(m_0[1,0],0.25,12)
@@ -632,7 +633,8 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         self.assertEqual(m_0.getnnz(),7)
         # ConservativeVolumic (division by sum of rows)
         rowSum=m.sum(axis=1)
-        m_1=diags(array(1/rowSum.transpose()),[0])*m
+        # version 0.12.0 # m_1=diags(array(1/rowSum.transpose()),[0])*m
+        m_1=spdiags(array(1/rowSum.transpose()),[0],rowSum.shape[0],rowSum.shape[0])*m
         del rowSum
         self.assertAlmostEqual(m_1[0,0],1.,12)
         self.assertAlmostEqual(m_1[1,0],0.4,12)
@@ -693,6 +695,96 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         delta=m0-m1t
         self.assertTrue(DataArrayDouble(delta.data).isUniform(0.,1e-17))
         pass
+
+    @unittest.skipUnless(MEDCouplingHasNumPyBindings() and MEDCouplingHasSciPyBindings(),"requires numpy AND scipy")
+    def testNonConformWithRemapper_1(self):
+        coo=DataArrayDouble([-0.396700000780411,-0.134843245350081,-0.0361311386958691,-0.407550009429364,-0.13484324535008,-0.0361311386958923,-0.396700000780411,-0.132191446077668,-0.0448729493559049,-0.407550009429364,-0.132191446077666,-0.0448729493559254,-0.396700000780411,-0.128973582738749,-0.0534226071577727,-0.407550009429364,-0.128973582738747,-0.0534226071577904,-0.396700000780411,-0.128348829636458,-0.0346583696473619,-0.407550009429364,-0.128348829636457,-0.0346583696473822,-0.396700000780411,-0.125874740261886,-0.0430683597970123,-0.407550009429364,-0.125874740261885,-0.0430683597970302,-0.396700000780411,-0.122905344829122,-0.051310216195766,-0.407550009429364,-0.12290534482912,-0.0513102161957814],12,3)
+        conn=DataArrayInt([2,9,3,11,2,3,5,11,2,8,9,11,2,10,8,11,2,5,4,11,2,4,10,11,3,0,1,6,3,1,7,6,3,2,0,6,3,8,2,6,3,7,9,6,3,9,8,6])
+        m=MEDCoupling1SGTUMesh("mesh",NORM_TETRA4)
+        m.setNodalConnectivity(conn)
+        m.setCoords(coo)
+        # m is ready
+        m1,d,di,rd,rdi=m.buildUnstructured().buildDescendingConnectivity()
+        rdi2=rdi.deltaShiftIndex()
+        cellIds=rdi2.getIdsEqual(1)
+        skinAndNonConformCells=m1[cellIds]
+        skinAndNonConformCells.zipCoords() # at this point skinAndNonConformCells contains non conform cells and skin cells. Now trying to split them in two parts.
+        #
+        rem=MEDCouplingRemapper()
+        rem.setMaxDistance3DSurfIntersect(1e-12)
+        rem.setMinDotBtwPlane3DSurfIntersect(0.99)# this line is important it is to tell to remapper to select only cells with very close orientation 
+        rem.prepare(skinAndNonConformCells,skinAndNonConformCells,"P0P0")
+        mat=rem.getCrudeCSRMatrix()
+        indptr=DataArrayInt(mat.indptr)
+        indptr2=indptr.deltaShiftIndex()
+        cellIdsOfNonConformCells=indptr2.getIdsNotEqual(1)
+        cellIdsOfSkin=indptr2.getIdsEqual(1)
+        self.assertTrue(cellIdsOfSkin.isEqual(DataArrayInt([1,2,3,5,6,7,8,9,10,11,12,13,14,15,16,17,19,20,21,23])))
+        self.assertTrue(cellIdsOfNonConformCells.isEqual(DataArrayInt([0,4,18,22])))
+        pass
+
+    def test3D1DOnP1P0_1(self):
+        """ This test focused on P1P0 interpolation with a source with meshDim=1 spaceDim=3 and a target with meshDim=3.
+        This test has revealed a bug in remapper. A reverse matrix is computed so a reverse method should be given in input.
+        """
+        target=MEDCouplingCMesh()
+        arrX=DataArrayDouble([0,1]) ; arrY=DataArrayDouble([0,1]) ; arrZ=DataArrayDouble(11) ; arrZ.iota()
+        target.setCoords(arrX,arrY,arrZ)
+        target=target.buildUnstructured() ; target.setName("TargetSecondaire")
+        #
+        sourceCoo=DataArrayDouble([(0.5,0.5,0.1),(0.5,0.5,1.2),(0.5,0.5,1.6),(0.5,0.5,1.8),(0.5,0.5,2.43),(0.5,0.5,2.55),(0.5,0.5,4.1),(0.5,0.5,4.4),(0.5,0.5,4.9),(0.5,0.5,5.1),(0.5,0.5,7.6),(0.5,0.5,7.7),(0.5,0.5,8.2),(0.5,0.5,8.4),(0.5,0.5,8.6),(0.5,0.5,8.8),(0.5,0.5,9.2),(0.5,0.5,9.6),(0.5,0.5,11.5)])
+        source=MEDCoupling1SGTUMesh("SourcePrimaire",NORM_SEG2)
+        source.setCoords(sourceCoo)
+        source.allocateCells()
+        for i in xrange(len(sourceCoo)-1):
+            source.insertNextCell([i,i+1])
+            pass
+        source=source.buildUnstructured()
+        fsource=MEDCouplingFieldDouble(ON_NODES) ; fsource.setName("field")
+        fsource.setMesh(source)
+        arr=DataArrayDouble(len(sourceCoo)) ; arr.iota(0.7) ; arr*=arr
+        fsource.setArray(arr)
+        fsource.setNature(ConservativeVolumic)
+        #
+        rem=MEDCouplingRemapper()
+        rem.setIntersectionType(PointLocator)
+        rem.prepare(source,target,"P1P0")
+        f2Test=rem.transferField(fsource,-27)
+        self.assertEqual(f2Test.getName(),fsource.getName())
+        self.assertEqual(f2Test.getMesh().getHiddenCppPointer(),target.getHiddenCppPointer())
+        expArr=DataArrayDouble([0.49,7.956666666666667,27.29,-27,59.95666666666667,94.09,-27,125.69,202.89,296.09])
+        self.assertTrue(f2Test.getArray().isEqual(expArr,1e-12))
+        f2Test=rem.reverseTransferField(f2Test,-36)
+        self.assertEqual(f2Test.getName(),fsource.getName())
+        self.assertEqual(f2Test.getMesh().getHiddenCppPointer(),source.getHiddenCppPointer())
+        expArr2=DataArrayDouble([0.49,7.956666666666667,7.956666666666667,7.956666666666667,27.29,27.29,59.95666666666667,59.95666666666667,59.95666666666667,94.09,125.69,125.69,202.89,202.89,202.89,202.89,296.09,296.09,-36.])
+        self.assertTrue(f2Test.getArray().isEqual(expArr2,1e-12))
+        pass
+
+    def testRemapperAMR1(self):
+        """ This test is the origin of the ref values for MEDCouplingBasicsTest.testAMR2"""
+        coarse=DataArrayDouble(35) ; coarse.iota(0) #X=5,Y=7
+        fine=DataArrayDouble(3*2*4*4) ; fine.iota(0) #X=3,Y=2 refined by 4
+        MEDCouplingIMesh.CondenseFineToCoarse([5,7],fine,[(1,4),(2,4)],[4,4],coarse)
+        #
+        m=MEDCouplingCartesianAMRMesh("mesh",2,[6,8],[0.,0.],[1.,1.])
+        trgMesh=m.buildUnstructured()
+        m.addPatch([(1,4),(2,4)],[4,4])
+        srcMesh=m[0].getMesh().buildUnstructured()
+        srcField=MEDCouplingFieldDouble(ON_CELLS)
+        fine2=DataArrayDouble(3*2*4*4) ; fine2.iota(0) ; srcField.setArray(fine2)
+        srcField.setMesh(srcMesh) ; srcField.setNature(Integral)
+        #
+        trgField=MEDCouplingFieldDouble(ON_CELLS)
+        coarse2=DataArrayDouble(35) ; coarse2.iota(0) ; trgField.setArray(coarse2)
+        trgField.setMesh(trgMesh) ; trgField.setNature(Integral)
+        #
+        rem=MEDCouplingRemapper()
+        rem.prepare(srcMesh,trgMesh,"P0P0")
+        rem.partialTransfer(srcField,trgField)
+        #
+        self.assertTrue(coarse.isEqual(trgField.getArray(),1e-12))
+        pass
     
     def build2DSourceMesh_1(self):
         sourceCoords=[-0.3,-0.3, 0.7,-0.3, -0.3,0.7, 0.7,0.7]