src/Plugins/MEDReader/plugin/Test/testMEDReader22.py

   1 #  -*- coding: iso-8859-1 -*-
   2 # Copyright (C) 2020-2022  CEA/DEN, EDF R&D
   3 #
   4 # This library is free software; you can redistribute it and/or
   5 # modify it under the terms of the GNU Lesser General Public
   6 # License as published by the Free Software Foundation; either
   7 # version 2.1 of the License, or (at your option) any later version.
   8 #
   9 # This library is distributed in the hope that it will be useful,
  10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 # Lesser General Public License for more details.
  13 #
  14 # You should have received a copy of the GNU Lesser General Public
  15 # License along with this library; if not, write to the Free Software
  16 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 #
  18 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19 #
  20 # Author : Anthony Geay (EDF R&D)
  21
  22
  23 __doc__ = """
  24 Test of GroupsAsMultiBlocks filter in the MEDReader plugin.
  25 """
  26
  27 from paraview.simple import *
  28 paraview.simple._DisableFirstRenderCameraReset()
  29 from vtk.util import numpy_support
  30 import numpy as np
  31 import medcoupling as mc
  32 from MEDReaderHelper import WriteInTmpDir,RetriveBaseLine
  33
  34 def MyAssert(clue):
  35     if not clue:
  36         raise RuntimeError("Assertion failed !")
  37
  38 def generateCase(fname):
  39     arr = mc.DataArrayDouble([0,1,2,3,4,5])
  40     m = mc.MEDCouplingCMesh()
  41     m.setCoords(arr,arr)
  42     m = m.buildUnstructured()
  43     m.changeSpaceDimension(3,0.)
  44     m.setName("mesh")
  45     mm = mc.MEDFileUMesh()
  46     mm[0] = m
  47     grp_BottomLeft = mc.DataArrayInt([0,1,5,6]) ; grp_BottomLeft.setName("BottomLeft")
  48     grp_BottomRight = mc.DataArrayInt([3,4,8,9]) ; grp_BottomRight.setName("BottomRight")
  49     grp_TopLeft = mc.DataArrayInt([15,16,20,21]) ; grp_TopLeft.setName("TopLeft")
  50     grp_TopRight = mc.DataArrayInt([18,19,23,24]) ; grp_TopRight.setName("TopRight")
  51     mm.setGroupsAtLevel(0,[grp_BottomLeft,grp_BottomRight,grp_TopLeft,grp_TopRight])
  52     mm.write(fname,2)
  53     f = mc.MEDCouplingFieldDouble(mc.ON_CELLS)
  54     f.setMesh(m)
  55     f.setArray(m.computeCellCenterOfMass().magnitude())
  56     f.setName("field")
  57     mc.WriteFieldUsingAlreadyWrittenMesh(fname,f)
  58     f2 = mc.MEDCouplingFieldDouble(mc.ON_NODES)
  59     f2.setMesh(m)
  60     f2.setArray(m.getCoords().magnitude())
  61     f2.setName("field2")
  62     mc.WriteFieldUsingAlreadyWrittenMesh(fname,f2)
  63     return f,f2,grp_BottomLeft,grp_BottomRight,grp_TopLeft,grp_TopRight
  64
  65 @WriteInTmpDir
  66 def test():
  67     fname = "testMEDReader22.med"
  68     f,f2,grp_BottomLeft,grp_BottomRight,grp_TopLeft,grp_TopRight = generateCase(fname)
  69     reader = MEDReader(FileName=fname)
  70     reader.AllArrays = ['TS0/mesh/ComSup0/field@@][@@P0','TS0/mesh/ComSup0/field2@@][@@P1','TS0/mesh/ComSup0/mesh@@][@@P0']
  71     reader.AllTimeSteps = ['0000']
  72
  73     groupsNames = GroupsNames(Input=reader)
  74     groupsNames.UpdatePipeline()
  75     MyAssert( servermanager.Fetch(groupsNames).GetNumberOfBlocks() == 1 )
  76     gn = servermanager.Fetch(groupsNames).GetBlock(0)
  77     blockIDS = numpy_support.vtk_to_numpy(gn.GetColumnByName("Block ID"))
  78     MyAssert(np.all(blockIDS==np.array([0,1,2,3],dtype=np.int32)))
  79     grpNames = gn.GetColumnByName("Group Name")
  80     MyAssert(grpNames.GetNumberOfTuples()==4)
  81     MyAssert([grpNames.GetValue(i) for i in range(4)] == ['BottomLeft','BottomRight','TopLeft','TopRight'])
  82
  83     groupsAsMultiBlocks = GroupsAsMultiBlocks(Input=reader)
  84     groupsAsMultiBlocks.UpdatePipeline()
  85     blocks = servermanager.Fetch(groupsAsMultiBlocks)
  86     MyAssert(blocks.GetNumberOfBlocks()==4)
  87     # test of bottom left
  88     ds_bl = blocks.GetBlock(0)
  89     MyAssert(ds_bl.GetNumberOfCells()==4)
  90     bl_ref_conn = np.array([ 1,  0,  6,  7,  2,  1,  7,  8,  7,  6, 12, 13,  8,  7, 13, 14],dtype=np.int32)
  91     MyAssert(ds_bl.GetCellData().GetNumberOfArrays() == 3 )# 3 for field, mesh and FamilyIdCell
  92     MyAssert(np.all( bl_ref_conn == numpy_support.vtk_to_numpy( ds_bl.GetCells().GetConnectivityArray() )) )
  93     MyAssert( mc.DataArrayDouble(numpy_support.vtk_to_numpy( ds_bl.GetCellData().GetArray("field") )).isEqual(f.getArray()[grp_BottomLeft],1e-12) )
  94     # test of bottom right
  95     ds_br = blocks.GetBlock(1)
  96     MyAssert(ds_br.GetNumberOfCells()==4)
  97     br_ref_conn = np.array([ 4,  3,  9, 10,  5,  4, 10, 11, 10,  9, 15, 16, 11, 10, 16, 17],dtype=np.int32)
  98     MyAssert(np.all( br_ref_conn == numpy_support.vtk_to_numpy( ds_br.GetCells().GetConnectivityArray() )) )
  99     MyAssert( mc.DataArrayDouble(numpy_support.vtk_to_numpy( ds_br.GetCellData().GetArray("field") )).isEqual(f.getArray()[grp_BottomRight],1e-12) )
 100     # test of top left
 101     ds_tl = blocks.GetBlock(2)
 102     MyAssert(ds_tl.GetNumberOfCells()==4)
 103     tl_ref_conn = np.array([ 19, 18, 24, 25, 20, 19, 25, 26, 25, 24, 30, 31, 26, 25, 31, 32],dtype=np.int32)
 104     MyAssert(np.all( tl_ref_conn == numpy_support.vtk_to_numpy( ds_tl.GetCells().GetConnectivityArray() )) )
 105     MyAssert( mc.DataArrayDouble(numpy_support.vtk_to_numpy( ds_tl.GetCellData().GetArray("field") )).isEqual(f.getArray()[grp_TopLeft],1e-12) )
 106     # test of top right
 107     ds_tr = blocks.GetBlock(3)
 108     MyAssert(ds_tr.GetNumberOfCells()==4)
 109     tr_ref_conn = np.array([ 22, 21, 27, 28, 23, 22, 28, 29, 28, 27, 33, 34, 29, 28, 34, 35],dtype=np.int32)
 110     MyAssert(np.all( tr_ref_conn == numpy_support.vtk_to_numpy( ds_tr.GetCells().GetConnectivityArray() )) )
 111     MyAssert( mc.DataArrayDouble(numpy_support.vtk_to_numpy( ds_tr.GetCellData().GetArray("field") )).isEqual(f.getArray()[grp_TopRight],1e-12) )
 112     #
 113     for ds in [ds_bl,ds_br,ds_tl,ds_tr]:
 114         MyAssert(ds.GetPointData().GetNumberOfArrays() == 1 )# 1 for field2
 115         MyAssert(ds.GetNumberOfPoints()==36) # for performance reasons all blocks share the same input coordinates
 116         MyAssert(mc.DataArrayDouble( numpy_support.vtk_to_numpy( ds.GetPointData().GetArray("field2") ) ).isEqual(f2.getArray(),1e-12) )
 117
 118 if __name__ == "__main__":
 119     test()
 120