src/medtool/src/MEDLoader/Swig/SauvLoaderTest.py

   1 #  -*- coding: iso-8859-1 -*-
   2 # Copyright (C) 2007-2015  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 : Edward AGAPOV (eap)
  21
  22 from MEDLoader import *
  23 import unittest, os
  24 from MEDLoaderDataForTest import MEDLoaderDataForTest
  25
  26 class SauvLoaderTest(unittest.TestCase):
  27
  28     def __getResourcesDirectory(self):
  29         med_root_dir=os.getenv("MEDTOOL_ROOT_DIR")
  30         if med_root_dir:
  31             return os.path.join( os.getenv("MEDTOOL_ROOT_DIR"), "share","resources","med")
  32         else:
  33             current_dir = os.path.dirname(os.path.realpath(__file__))
  34             return os.path.join(current_dir, "..", "..", "..", "resources")
  35         pass
  36
  37     def testSauv2Med(self):
  38         # get a file containing all types of readable piles
  39         sauvFile = os.path.join( self.__getResourcesDirectory(),"allPillesTest.sauv")
  40         self.assertTrue( os.access( sauvFile, os.F_OK))
  41
  42         # read SAUV and write MED
  43         medFile = "SauvLoaderTest.med"
  44         sr=SauvReader(sauvFile);
  45         d2=sr.loadInMEDFileDS();
  46         d2.write(medFile,0);
  47
  48         # check
  49         self.assertEqual(1,d2.getNumberOfMeshes())
  50         self.assertEqual(8+97,d2.getNumberOfFields())
  51         mm = d2.getMeshes()
  52         m = mm.getMeshAtPos(0)
  53         self.assertEqual(17,len(m.getGroupsNames()))
  54
  55         os.remove( medFile )
  56         pass
  57
  58     def testMed2Sauv(self):
  59         # read pointe.med
  60         medFile = os.path.join(self.__getResourcesDirectory(),"pointe.med")
  61         self.assertTrue( os.access( medFile, os.F_OK))
  62         pointeMed = MEDFileData.New( medFile )
  63
  64         # add 3 faces to pointeMed
  65         pointeMedMesh = pointeMed.getMeshes().getMeshAtPos(0)
  66         pointeM1D = MEDCouplingUMesh.New()
  67         pointeM1D.setCoords( pointeMedMesh.getCoords() )
  68         pointeM1D.setMeshDimension( 2 )
  69         pointeM1D.allocateCells( 3 )
  70         pointeM1D.insertNextCell( NORM_TRI3, 3, [0,1,2])
  71         pointeM1D.insertNextCell( NORM_TRI3, 3, [0,1,3])
  72         pointeM1D.insertNextCell( NORM_QUAD4, 4, [10,11,12,13])
  73         pointeM1D.finishInsertingCells()
  74         pointeMedMesh.setMeshAtLevel( -1, pointeM1D )
  75         pointeMed.getMeshes().setMeshAtPos( 0, pointeMedMesh )
  76
  77         # add a field on 2 faces to pointeMed
  78         ff1=MEDFileFieldMultiTS.New()
  79         f1=MEDCouplingFieldDouble.New(ON_GAUSS_NE,ONE_TIME)
  80         #f1.setMesh( pointeM1D )
  81         f1.setName("Field on 2 faces")
  82         d=DataArrayDouble.New()
  83         d.alloc(3+4,2)
  84         d.setInfoOnComponent(0,"sigX [MPa]")
  85         d.setInfoOnComponent(1,"sigY [GPa]")
  86         d.setValues([311,312,321,322,331,332,411,412,421,422,431,432,441,442],3+4,2)
  87         f1.setArray(d)
  88         da=DataArrayInt.New()
  89         da.setValues([0,2],2,1)
  90         da.setName("sup2")
  91         ff1.appendFieldProfile(f1,pointeMedMesh,-1,da)
  92         pointeMed.getFields().pushField( ff1 )
  93
  94         #remove fieldnodeint
  95         pointeFields = pointeMed.getFields()
  96         for i in range( pointeFields.getNumberOfFields() ):
  97             if pointeFields.getFieldAtPos(i).getName() == "fieldnodeint":
  98                 pointeFields.destroyFieldAtPos( i )
  99                 break
 100
 101         # write pointeMed to SAUV
 102         sauvFile = "SauvLoaderTest.sauv"
 103         sw=SauvWriter();
 104         sw.setMEDFileDS(pointeMed);
 105         sw.write(sauvFile);
 106
 107         # read SAUV and check
 108         sr=SauvReader.New(sauvFile);
 109         d2=sr.loadInMEDFileDS();
 110         self.assertEqual(1,d2.getNumberOfMeshes())
 111         self.assertEqual(4,d2.getNumberOfFields())
 112         m = d2.getMeshes().getMeshAtPos(0)
 113         self.assertEqual("maa1",m.getName())
 114         self.assertEqual(6,len(m.getGroupsNames()))
 115         self.assertEqual(3,m.getMeshDimension())
 116         groups = m.getGroupsNames()
 117         self.assertTrue( "groupe1" in groups )
 118         self.assertTrue( "groupe2" in groups )
 119         self.assertTrue( "groupe3" in groups )
 120         self.assertTrue( "groupe4" in groups )
 121         self.assertTrue( "groupe5" in groups )
 122         self.assertTrue( "maa1" in groups )
 123         self.assertEqual(16,m.getSizeAtLevel(0))
 124         um0 = m.getGenMeshAtLevel(0)
 125         self.assertEqual(12, um0.getNumberOfCellsWithType( NORM_TETRA4 ))
 126         self.assertEqual(2, um0.getNumberOfCellsWithType( NORM_PYRA5 ))
 127         self.assertEqual(2, um0.getNumberOfCellsWithType( NORM_HEXA8 ))
 128         um1 = m.getGenMeshAtLevel(-1)
 129         self.assertEqual(2, um1.getNumberOfCellsWithType( NORM_TRI3 ))
 130         pointeUM0 = pointeMedMesh.getGenMeshAtLevel(0)
 131         self.assertTrue(m.getCoords().isEqualWithoutConsideringStr(pointeMedMesh.getCoords(),1e-12))
 132         self.assertEqual( um0.getMeasureField(0).accumulate(0),
 133                           pointeUM0.getMeasureField(0).accumulate(0),1e-12)
 134         # check fields
 135         # fieldnodedouble
 136         fieldnodedoubleTS1 = pointeMed.getFields().getFieldWithName("fieldnodedouble")
 137         fieldnodedoubleTS2 = d2.getFields().getFieldWithName("fieldnodedouble")
 138         self.assertEqual( fieldnodedoubleTS1.getInfo(), fieldnodedoubleTS2.getInfo())
 139         self.assertEqual( fieldnodedoubleTS1.getNumberOfTS(), fieldnodedoubleTS2.getNumberOfTS())
 140         io1 = fieldnodedoubleTS1.getIterations()
 141         io2 = fieldnodedoubleTS2.getIterations()
 142         for i in range(fieldnodedoubleTS1.getNumberOfTS() ):
 143             fnd1 = fieldnodedoubleTS1.getFieldOnMeshAtLevel(ON_NODES, io1[i][0],io1[i][1],pointeUM0)
 144             fnd2 = fieldnodedoubleTS2.getFieldOnMeshAtLevel(ON_NODES, io2[i][0],io2[i][1],um0)
 145             self.assertTrue( fnd1.getArray().isEqual( fnd2.getArray(), 1e-12 ))
 146         # fieldcelldoublevector
 147         fieldnodedoubleTS1 = pointeMed.getFields().getFieldWithName("fieldcelldoublevector")
 148         fieldnodedoubleTS2 = d2.getFields().getFieldWithName("fieldcelldoublevector")
 149         self.assertEqual( fieldnodedoubleTS1.getInfo(), fieldnodedoubleTS2.getInfo())
 150         self.assertEqual( fieldnodedoubleTS1.getNumberOfTS(), fieldnodedoubleTS2.getNumberOfTS())
 151         io1 = fieldnodedoubleTS1.getIterations()
 152         io2 = fieldnodedoubleTS2.getIterations()
 153         for i in range(fieldnodedoubleTS1.getNumberOfTS() ):
 154             fnd1 = fieldnodedoubleTS1.getFieldOnMeshAtLevel(ON_CELLS, io1[i][0],io1[i][1],pointeUM0)
 155             fnd2 = fieldnodedoubleTS2.getFieldOnMeshAtLevel(ON_CELLS, io2[i][0],io2[i][1],um0)
 156             self.assertAlmostEqual( fnd1.accumulate(0), fnd2.accumulate(0) )
 157             self.assertAlmostEqual( fnd1.accumulate(1), fnd2.accumulate(1) )
 158             self.assertAlmostEqual( fnd1.accumulate(2), fnd2.accumulate(2) )
 159             # Field on 2 faces
 160             fieldOnFaces = d2.getFields().getFieldWithName(f1.getName())
 161             io1 = fieldOnFaces.getIterations()
 162             fof = fieldOnFaces.getFieldOnMeshAtLevel(f1.getTypeOfField(),io1[i][0],io1[i][1],um1)
 163             self.assertTrue( d.isEqual( fof.getArray(), 1e-12 ))
 164             pass
 165         del sr
 166         os.remove( sauvFile )
 167         pass
 168
 169     def testSauv2MedWONodeFamilyNum(self):
 170         """test for issue 0021673: [CEA 566] Bug in SauvWriter when writing meshes
 171         having no family ids on nodes."""
 172
 173         myCoords=DataArrayDouble.New([-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 ],9,2)
 174         targetConn=[0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4];
 175         targetMesh=MEDCouplingUMesh.New("BugInSauvWriter",2);
 176         targetMesh.allocateCells(5);
 177         targetMesh.insertNextCell(NORM_TRI3,3,targetConn[4:7]);
 178         targetMesh.insertNextCell(NORM_TRI3,3,targetConn[7:10]);
 179         targetMesh.insertNextCell(NORM_QUAD4,4,targetConn[0:4]);
 180         targetMesh.insertNextCell(NORM_QUAD4,4,targetConn[10:14]);
 181         targetMesh.insertNextCell(NORM_QUAD4,4,targetConn[14:18]);
 182         targetMesh.finishInsertingCells();
 183         targetMesh.setCoords(myCoords);
 184         #
 185         m=MEDFileUMesh.New()
 186         m.setMeshAtLevel(0,targetMesh)
 187         # start of bug
 188         fam=DataArrayInt.New(targetMesh.getNumberOfNodes())
 189         fam[:]=0
 190         #m.setFamilyFieldArr(1,fam)
 191         #end of bug
 192
 193         ms=MEDFileMeshes.New()
 194         ms.setMeshAtPos(0,m)
 195         meddata=MEDFileData.New()
 196         meddata.setMeshes(ms)
 197
 198         medFile = "BugInSauvWriter.sauv"
 199         sw=SauvWriter.New();
 200         sw.setMEDFileDS(meddata);
 201         sw.write(medFile);
 202
 203         os.remove( medFile )
 204         pass
 205
 206     def testSauv2MedOnPipe1D(self):
 207         """test for issue 0021745: [CEA 600] Some missing groups in mesh after reading a SAUV file with SauvReader."""
 208         sauvFile="Test_sauve_1D.sauv"
 209         # Make a sauve file with a qudratic 1D mesh
 210         m=MEDCouplingUMesh.New("pipe1D",1)
 211         m.allocateCells(2);
 212         targetConn=[0,2,1, 2,4,3]
 213         m.insertNextCell(NORM_SEG3,3,targetConn[0:3])
 214         m.insertNextCell(NORM_SEG3,3,targetConn[3:6])
 215         m.finishInsertingCells();
 216         # coords
 217         coords=[ 0.,1.,2.,4.,5. ];
 218         c=DataArrayDouble.New()
 219         c.setValues(coords,5,1)
 220         m.setCoords(c)
 221         # MEDFileUMesh
 222         mm=MEDFileUMesh.New()
 223         mm.setName(m.getName())
 224         mm.setDescription("1D mesh")
 225         mm.setCoords(c)
 226         mm.setMeshAtLevel(0,m);
 227         # MEDFileData
 228         mfd1 = MEDFileData.New()
 229         ms=MEDFileMeshes.New(); ms.setMeshAtPos(0,mm)
 230         mfd1.setMeshes(ms)
 231         # write
 232         sw=SauvWriter.New()
 233         sw.setMEDFileDS(mfd1)
 234         sw.write(sauvFile)
 235         # Check connectivity read from the sauv file
 236         sr = SauvReader.New(sauvFile)
 237         mfd2 = sr.loadInMEDFileDS()
 238         mfMesh = mfd2.getMeshes()[0]
 239         mesh = mfMesh.getMeshAtLevel(0)
 240         self.assertTrue(mesh.getNodalConnectivity().isEqual(m.getNodalConnectivity()))
 241         #
 242         del sr
 243         os.remove(sauvFile)
 244         pass
 245
 246     @unittest.skipUnless(MEDLoader.HasXDR(),"requires XDR")
 247     def testMissingGroups(self):
 248         """test for issue 0021749: [CEA 601] Some missing groups in mesh after reading a SAUV file with SauvReader."""
 249         sauvFile = os.path.join(self.__getResourcesDirectory(),"BDC-714.sauv")
 250         self.assertTrue( os.access( sauvFile, os.F_OK))
 251         name_of_group_on_cells='Slice10:ABSORBER'
 252         name_of_group_on_cells2='Slice10:00LR'
 253         sr=SauvReader.New(sauvFile)
 254         mfd2=sr.loadInMEDFileDS()
 255         mfMesh=mfd2.getMeshes()[0]
 256         #
 257         self.assertTrue(name_of_group_on_cells in mfMesh.getGroupsNames())
 258         self.assertTrue(name_of_group_on_cells2 in mfMesh.getGroupsNames())
 259         self.assertEqual(270,len(mfMesh.getGroupsNames()))
 260         #
 261         ids1=mfMesh.getGroupArr(0,name_of_group_on_cells)
 262         ids2=mfMesh.getGroupArr(0,name_of_group_on_cells2)
 263         ids1.setName("")
 264         ids2.setName("")
 265         self.assertTrue(ids1.isEqual(ids2))
 266         pass
 267
 268     def testGaussPt(self):
 269         """issue 22321: [CEA 933] Bug when reading a sauve file containing field on Gauss Pt.
 270         The problem was that a field ON_GAUSS_PT was created but no Gauss Localization
 271         was defined"""
 272
 273         # create a MEDFileData with a field ON_GAUSS_PT: 9 Gauss points, on 4 QUAD8 elements
 274         f=MEDCouplingFieldDouble(ON_GAUSS_PT)
 275         m=MEDCouplingUMesh("mesh",2) ; m.allocateCells()
 276         m.insertNextCell(NORM_QUAD8,[0,2,4,6,1,3,5,7])
 277         m.insertNextCell(NORM_QUAD8,[2,9,11,4,8,10,12,3])
 278         m.insertNextCell(NORM_QUAD8,[6,4,14,16,5,13,15,17])
 279         m.insertNextCell(NORM_QUAD8,[4,11,19,14,12,18,20,13])
 280         m.setCoords(DataArrayDouble([(0,0),(0,0.25),(0,0.5),(0.25,0.5),(0.5,0.5),(0.5,0.25),(0.5,0),(0.25,0),(0,0.75),(0,1),(0.25,1),(0.5,1),(0.5,0.75),(0.75,0.5),(1,0.5),(1,0.25),(1,0),(0.75,0),(0.75,1),(1,1),(1,0.75)],21,2))
 281         f.setMesh(m)
 282         arr=DataArrayDouble(4*9*2) ; arr.iota() ; arr.rearrange(2) ; arr.setInfoOnComponents(["YOUN []","NU []"])
 283         f.setArray(arr)
 284         refCoo=[-1,-1,1,-1,1,1,-1,1,0,-1,1,0,0,1,-1,0]
 285         gpCoo=[-0.7,-0.7,0.7,-0.7,0.7,0.7,-0.7,0.7,0,-0.7,0.7,0,0,0.7,-0.7,0,0,0]
 286         wgt=[0.3,0.3,0.3,0.3,0.4,0.4,0.4,0.4,0.7]
 287         f.setGaussLocalizationOnType(NORM_QUAD8,refCoo,gpCoo,wgt)
 288         f.setName("SIGT")
 289         f.checkCoherency()
 290         #
 291         mm=MEDFileUMesh()
 292         mm.setMeshAtLevel(0,m)
 293         mfm = MEDFileMeshes()
 294         mfm.pushMesh( mm )
 295         ff=MEDFileField1TS()
 296         ff.setFieldNoProfileSBT(f)
 297         mfmts = MEDFileFieldMultiTS()
 298         mfmts.pushBackTimeStep(ff)
 299         mff = MEDFileFields()
 300         mff.pushField( mfmts )
 301         mfd = MEDFileData.New()
 302         mfd.setFields( mff )
 303         mfd.setMeshes( mfm )
 304
 305         # convert the MED file to a SAUV file
 306         sauvFile = "SauvLoaderTest_testGaussPt.sauv"
 307         sw=SauvWriter.New();
 308         sw.setMEDFileDS(mfd);
 309         sw.write(sauvFile);
 310
 311         # convert the SAUV file back to MED
 312         sr=SauvReader.New(sauvFile);
 313         d2=sr.loadInMEDFileDS();
 314
 315         self.assertEqual( 1, d2.getNumberOfFields() )
 316         self.assertEqual( 1, d2.getNumberOfMeshes() )
 317         mfm2 = d2.getMeshes()[0]
 318         mff2 = d2.getFields()[0]
 319         m2 = mfm2.getMeshAtLevel(0)
 320         f2 = mff2.getTimeStepAtPos(0).getFieldOnMeshAtLevel(f.getTypeOfField(),0,mfm2)
 321         f2.setGaussLocalizationOnType(NORM_QUAD8,refCoo,gpCoo,wgt) # not stored in SAUV
 322         #f2.setOrder( f.getTime()[2] ) # not stored in SAUV
 323         self.assertTrue( m2.isEqual( m, 1e-12 ))
 324         self.assertTrue( f2.isEqual( f, 1e-12, 1e-12 ))
 325         del sr
 326         os.remove( sauvFile )
 327         pass
 328
 329     def testSauvWriterGroupWithOneFamily(self):
 330         """
 331         This test checks an option for sauv writing. It is requested here to copy a group from a family if a group is lying on a single family.
 332         """
 333         import re
 334         mfd=MEDLoaderDataForTest.buildAMEDFileDataWithGroupOnOneFamilyForSauv()
 335         sauvFile = "mesh.sauv"
 336         sw=SauvWriter.New()
 337         sw.setMEDFileDS(mfd)
 338         self.assertTrue(not sw.getCpyGrpIfOnASingleFamilyStatus())
 339         sw.setCpyGrpIfOnASingleFamilyStatus(True)
 340         self.assertTrue(sw.getCpyGrpIfOnASingleFamilyStatus())
 341         sw.write(sauvFile)
 342
 343         f = open(sauvFile)
 344         # String pattern for the header of the sub meshes record ("PILE" number, number of named objects, number of objects)
 345         pattern_pile= re.compile(r'\sPILE\sNUMERO\s+(?P<number>[0-9]+)NBRE\sOBJETS\sNOMMES\s+(?P<nbnamed>[0-9]+)NBRE\sOBJETS\s+(?P<nbobjects>[0-9]+)')
 346         # String pattern for a sub mesh header (cell type, number of components and three numbers)
 347         pattern_header=re.compile(r'\s+(?P<type>[0-9]+)\s+(?P<nbsubs>[0-9]+)\s+[0-9]+\s+[0-9]+\s+[0-9]+')
 348
 349         nbobjects=0
 350         line = f.readline()
 351         while(line):
 352             match_pile = pattern_pile.match(line)
 353             if match_pile:
 354                 number=int(match_pile.group("number"))
 355                 if number == 1:
 356                     nbnamed=int(match_pile.group("nbnamed"))
 357                     nbobjects=int(match_pile.group("nbobjects"))
 358                     break
 359                 pass
 360             line=f.readline()
 361             pass
 362
 363         # Skipping the objects names
 364         f.readline()
 365         # Skipping the objects ids
 366         f.readline()
 367
 368         # Looking for each sub-mesh header
 369         line = f.readline()
 370         cur_object=0
 371         while(line and cur_object < nbobjects):
 372             match_header=pattern_header.match(line)
 373             if match_header:
 374                 cell_type=int(match_header.group("type"))
 375                 nb_subs=int(match_header.group("nbsubs"))
 376                 # Looking for a compound object
 377                 if cell_type == 0:
 378                     # Testing if there is only one component
 379                     self.assertTrue(nb_subs > 1)
 380                 else:
 381                     f.readline()
 382                     f.readline()
 383                     cur_object = cur_object + 1
 384                     pass
 385                 pass
 386             line=f.readline()
 387             pass
 388         f.close()
 389         os.remove(sauvFile)
 390         pass
 391
 392     pass
 393
 394 unittest.main()