Salome HOME
22321: [CEA 933] Bug when reading a sauve file containing field on Gauss Pt
[modules/med.git] / src / MEDLoader / Swig / SauvLoaderTest.py
1 #  -*- coding: iso-8859-1 -*-
2 # Copyright (C) 2007-2013  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.
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 testSauv2Med(self):
29         # get a file containing all types of readable piles
30         self.assertTrue( os.getenv("MED_ROOT_DIR") )
31         sauvFile = os.path.join( os.getenv("MED_ROOT_DIR"), "share","salome",
32                                  "resources","med","allPillesTest.sauv")
33         self.assertTrue( os.access( sauvFile, os.F_OK))
34
35         # read SAUV and write MED
36         medFile = "SauvLoaderTest.med"
37         sr=SauvReader.New(sauvFile);
38         d2=sr.loadInMEDFileDS();
39         d2.write(medFile,0);
40
41         # check 
42         self.assertEqual(1,d2.getNumberOfMeshes())
43         self.assertEqual(8+97,d2.getNumberOfFields())
44         mm = d2.getMeshes()
45         m = mm.getMeshAtPos(0)
46         self.assertEqual(17,len(m.getGroupsNames()))
47
48         os.remove( medFile )
49         pass
50
51     def testMed2Sauv(self):
52         # read pointe.med
53         self.assertTrue( os.getenv("MED_ROOT_DIR") )
54         medFile = os.path.join( os.getenv("MED_ROOT_DIR"), "share","salome",
55                                 "resources","med","pointe.med")
56         self.assertTrue( os.access( medFile, os.F_OK))
57         pointeMed = MEDFileData.New( medFile )
58
59         # add 3 faces to pointeMed
60         pointeMedMesh = pointeMed.getMeshes().getMeshAtPos(0)
61         pointeM1D = MEDCouplingUMesh.New()
62         pointeM1D.setCoords( pointeMedMesh.getCoords() )
63         pointeM1D.setMeshDimension( 2 )
64         pointeM1D.allocateCells( 3 )
65         pointeM1D.insertNextCell( NORM_TRI3, 3, [0,1,2])
66         pointeM1D.insertNextCell( NORM_TRI3, 3, [0,1,3])
67         pointeM1D.insertNextCell( NORM_QUAD4, 4, [10,11,12,13])
68         pointeM1D.finishInsertingCells()
69         pointeMedMesh.setMeshAtLevel( -1, pointeM1D )
70         pointeMed.getMeshes().setMeshAtPos( 0, pointeMedMesh )
71
72         # add a field on 2 faces to pointeMed
73         ff1=MEDFileFieldMultiTS.New()
74         f1=MEDCouplingFieldDouble.New(ON_GAUSS_NE,ONE_TIME)
75         #f1.setMesh( pointeM1D )
76         f1.setName("Field on 2 faces")
77         d=DataArrayDouble.New()
78         d.alloc(3+4,2)
79         d.setInfoOnComponent(0,"sigX [MPa]")
80         d.setInfoOnComponent(1,"sigY [GPa]")
81         d.setValues([311,312,321,322,331,332,411,412,421,422,431,432,441,442],3+4,2)
82         f1.setArray(d)
83         da=DataArrayInt.New()
84         da.setValues([0,2],2,1)
85         da.setName("sup2")
86         ff1.appendFieldProfile(f1,pointeMedMesh,-1,da)
87         pointeMed.getFields().pushField( ff1 )
88
89         #remove fieldnodeint
90         pointeFields = pointeMed.getFields()
91         for i in range( pointeFields.getNumberOfFields() ):
92             if pointeFields.getFieldAtPos(i).getName() == "fieldnodeint":
93                 pointeFields.destroyFieldAtPos( i )
94                 break
95
96         # write pointeMed to SAUV
97         sauvFile = "SauvLoaderTest.sauv"
98         sw=SauvWriter.New();
99         sw.setMEDFileDS(pointeMed);
100         sw.write(sauvFile);
101
102         # read SAUV and check
103         sr=SauvReader.New(sauvFile);
104         d2=sr.loadInMEDFileDS();
105         self.assertEqual(1,d2.getNumberOfMeshes())
106         self.assertEqual(4,d2.getNumberOfFields())
107         m = d2.getMeshes().getMeshAtPos(0)
108         self.assertEqual("maa1",m.getName())
109         self.assertEqual(6,len(m.getGroupsNames()))
110         self.assertEqual(3,m.getMeshDimension())
111         groups = m.getGroupsNames()
112         self.assertTrue( "groupe1" in groups )
113         self.assertTrue( "groupe2" in groups )
114         self.assertTrue( "groupe3" in groups )
115         self.assertTrue( "groupe4" in groups )
116         self.assertTrue( "groupe5" in groups )
117         self.assertTrue( "maa1" in groups )
118         self.assertEqual(16,m.getSizeAtLevel(0))
119         um0 = m.getGenMeshAtLevel(0)
120         self.assertEqual(12, um0.getNumberOfCellsWithType( NORM_TETRA4 ))
121         self.assertEqual(2, um0.getNumberOfCellsWithType( NORM_PYRA5 ))
122         self.assertEqual(2, um0.getNumberOfCellsWithType( NORM_HEXA8 ))
123         um1 = m.getGenMeshAtLevel(-1)
124         self.assertEqual(2, um1.getNumberOfCellsWithType( NORM_TRI3 ))
125         pointeUM0 = pointeMedMesh.getGenMeshAtLevel(0)
126         self.assertTrue(m.getCoords().isEqualWithoutConsideringStr(pointeMedMesh.getCoords(),1e-12))
127         self.assertEqual( um0.getMeasureField(0).accumulate(0),
128                           pointeUM0.getMeasureField(0).accumulate(0),1e-12)
129         # check fields
130         # fieldnodedouble
131         fieldnodedoubleTS1 = pointeMed.getFields().getFieldWithName("fieldnodedouble")
132         fieldnodedoubleTS2 = d2.getFields().getFieldWithName("fieldnodedouble")
133         self.assertEqual( fieldnodedoubleTS1.getInfo(), fieldnodedoubleTS2.getInfo())
134         self.assertEqual( fieldnodedoubleTS1.getNumberOfTS(), fieldnodedoubleTS2.getNumberOfTS())
135         io1 = fieldnodedoubleTS1.getIterations()
136         io2 = fieldnodedoubleTS2.getIterations()
137         for i in range(fieldnodedoubleTS1.getNumberOfTS() ):
138             fnd1 = fieldnodedoubleTS1.getFieldOnMeshAtLevel(ON_NODES, io1[i][0],io1[i][1],pointeUM0)
139             fnd2 = fieldnodedoubleTS2.getFieldOnMeshAtLevel(ON_NODES, io2[i][0],io2[i][1],um0)
140             self.assertTrue( fnd1.getArray().isEqual( fnd2.getArray(), 1e-12 ))
141         # fieldcelldoublevector
142         fieldnodedoubleTS1 = pointeMed.getFields().getFieldWithName("fieldcelldoublevector")
143         fieldnodedoubleTS2 = d2.getFields().getFieldWithName("fieldcelldoublevector")
144         self.assertEqual( fieldnodedoubleTS1.getInfo(), fieldnodedoubleTS2.getInfo())
145         self.assertEqual( fieldnodedoubleTS1.getNumberOfTS(), fieldnodedoubleTS2.getNumberOfTS())
146         io1 = fieldnodedoubleTS1.getIterations()
147         io2 = fieldnodedoubleTS2.getIterations()
148         for i in range(fieldnodedoubleTS1.getNumberOfTS() ):
149             fnd1 = fieldnodedoubleTS1.getFieldOnMeshAtLevel(ON_CELLS, io1[i][0],io1[i][1],pointeUM0)
150             fnd2 = fieldnodedoubleTS2.getFieldOnMeshAtLevel(ON_CELLS, io2[i][0],io2[i][1],um0)
151             self.assertAlmostEqual( fnd1.accumulate(0), fnd2.accumulate(0) )
152             self.assertAlmostEqual( fnd1.accumulate(1), fnd2.accumulate(1) )
153             self.assertAlmostEqual( fnd1.accumulate(2), fnd2.accumulate(2) )
154             # Field on 2 faces
155             fieldOnFaces = d2.getFields().getFieldWithName(f1.getName())
156             io1 = fieldOnFaces.getIterations()
157             fof = fieldOnFaces.getFieldOnMeshAtLevel(f1.getTypeOfField(),io1[i][0],io1[i][1],um1)
158             self.assertTrue( d.isEqual( fof.getArray(), 1e-12 ))
159
160             os.remove( sauvFile )
161             pass
162         pass
163
164     def testSauv2MedWONodeFamilyNum(self):
165         """test for issue 0021673: [CEA 566] Bug in SauvWriter when writing meshes
166         having no family ids on nodes."""
167
168         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)
169         targetConn=[0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4];
170         targetMesh=MEDCouplingUMesh.New("BugInSauvWriter",2);
171         targetMesh.allocateCells(5);
172         targetMesh.insertNextCell(NORM_TRI3,3,targetConn[4:7]);
173         targetMesh.insertNextCell(NORM_TRI3,3,targetConn[7:10]);
174         targetMesh.insertNextCell(NORM_QUAD4,4,targetConn[0:4]);
175         targetMesh.insertNextCell(NORM_QUAD4,4,targetConn[10:14]);
176         targetMesh.insertNextCell(NORM_QUAD4,4,targetConn[14:18]);
177         targetMesh.finishInsertingCells();
178         targetMesh.setCoords(myCoords);
179         #
180         m=MEDFileUMesh.New()
181         m.setMeshAtLevel(0,targetMesh)
182         # start of bug
183         fam=DataArrayInt.New(targetMesh.getNumberOfNodes())
184         fam[:]=0
185         #m.setFamilyFieldArr(1,fam)
186         #end of bug
187
188         ms=MEDFileMeshes.New()
189         ms.setMeshAtPos(0,m)
190         meddata=MEDFileData.New()
191         meddata.setMeshes(ms)
192
193         medFile = "BugInSauvWriter.sauv"
194         sw=SauvWriter.New();
195         sw.setMEDFileDS(meddata);
196         sw.write(medFile);
197
198         os.remove( medFile )
199         pass
200
201     def testSauv2MedOnPipe1D(self):
202         """test for issue 0021745: [CEA 600] Some missing groups in mesh after reading a SAUV file with SauvReader."""
203         sauvFile="Test_sauve_1D.sauv"
204         # Make a sauve file with a qudratic 1D mesh
205         m=MEDCouplingUMesh.New("pipe1D",1)
206         m.allocateCells(2);
207         targetConn=[0,2,1, 2,4,3]
208         m.insertNextCell(NORM_SEG3,3,targetConn[0:3])
209         m.insertNextCell(NORM_SEG3,3,targetConn[3:6])
210         m.finishInsertingCells();
211         # coords
212         coords=[ 0.,1.,2.,4.,5. ];
213         c=DataArrayDouble.New()
214         c.setValues(coords,5,1)
215         m.setCoords(c)
216         # MEDFileUMesh
217         mm=MEDFileUMesh.New()
218         mm.setName(m.getName())
219         mm.setDescription("1D mesh")
220         mm.setCoords(c)
221         mm.setMeshAtLevel(0,m);
222         # MEDFileData
223         mfd1 = MEDFileData.New()
224         ms=MEDFileMeshes.New(); ms.setMeshAtPos(0,mm)
225         mfd1.setMeshes(ms)
226         # write
227         sw=SauvWriter.New()
228         sw.setMEDFileDS(mfd1)
229         sw.write(sauvFile)
230         # Check connectivity read from the sauv file
231         sr = SauvReader.New(sauvFile)
232         mfd2 = sr.loadInMEDFileDS()
233         mfMesh = mfd2.getMeshes()[0]
234         mesh = mfMesh.getMeshAtLevel(0)
235         self.assertTrue(mesh.getNodalConnectivity().isEqual(m.getNodalConnectivity()))
236         #
237         os.remove(sauvFile)
238         pass
239
240     def testMissingGroups(self):
241         """test for issue 0021749: [CEA 601] Some missing groups in mesh after reading a SAUV file with SauvReader."""
242         self.assertTrue( os.getenv("MED_ROOT_DIR") )
243         sauvFile = os.path.join( os.getenv("MED_ROOT_DIR"), "share","salome",
244                                  "resources","med","BDC-714.sauv")
245         self.assertTrue( os.access( sauvFile, os.F_OK))
246         name_of_group_on_cells='Slice10:ABSORBER'
247         name_of_group_on_cells2='Slice10:00LR'
248         sr=SauvReader.New(sauvFile)
249         mfd2=sr.loadInMEDFileDS()
250         mfMesh=mfd2.getMeshes()[0]
251         #
252         self.assertTrue(name_of_group_on_cells in mfMesh.getGroupsNames())
253         self.assertTrue(name_of_group_on_cells2 in mfMesh.getGroupsNames())
254         self.assertEqual(270,len(mfMesh.getGroupsNames()))
255         #
256         ids1=mfMesh.getGroupArr(0,name_of_group_on_cells)
257         ids2=mfMesh.getGroupArr(0,name_of_group_on_cells2)
258         ids1.setName("")
259         ids2.setName("")
260         self.assertTrue(ids1.isEqual(ids2))
261         pass
262
263     def testGaussPt(self):
264         """issue 22321: [CEA 933] Bug when reading a sauve file containing field on Gauss Pt.
265         The problem was that a field ON_GAUSS_PT was created but no Gauss Localization
266         was defined"""
267
268         # create a MEDFileData with a field ON_GAUSS_PT: 9 Gauss points, on 4 QUAD8 elements
269         f=MEDCouplingFieldDouble(ON_GAUSS_PT)
270         m=MEDCouplingUMesh("mesh",2) ; m.allocateCells()
271         m.insertNextCell(NORM_QUAD8,[0,2,4,6,1,3,5,7])
272         m.insertNextCell(NORM_QUAD8,[2,9,11,4,8,10,12,3])
273         m.insertNextCell(NORM_QUAD8,[6,4,14,16,5,13,15,17])
274         m.insertNextCell(NORM_QUAD8,[4,11,19,14,12,18,20,13])
275         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))
276         f.setMesh(m)
277         arr=DataArrayDouble(4*9*2) ; arr.iota() ; arr.rearrange(2) ; arr.setInfoOnComponents(["YOUN []","NU []"])
278         f.setArray(arr)
279         refCoo=[-1,-1,1,-1,1,1,-1,1,0,-1,1,0,0,1,-1,0]
280         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]
281         wgt=[0.3,0.3,0.3,0.3,0.4,0.4,0.4,0.4,0.7]
282         f.setGaussLocalizationOnType(NORM_QUAD8,refCoo,gpCoo,wgt)
283         f.setName("SIGT")
284         f.checkCoherency()
285         #
286         mm=MEDFileUMesh()
287         mm.setMeshAtLevel(0,m)
288         mfm = MEDFileMeshes()
289         mfm.pushMesh( mm )
290         ff=MEDFileField1TS()
291         ff.setFieldNoProfileSBT(f)
292         mfmts = MEDFileFieldMultiTS()
293         mfmts.pushBackTimeStep(ff)
294         mff = MEDFileFields()
295         mff.pushField( mfmts )
296         mfd = MEDFileData.New()
297         mfd.setFields( mff )
298         mfd.setMeshes( mfm )
299
300         # convert the MED file to a SAUV file
301         sauvFile = "SauvLoaderTest::testGaussPt.sauv"
302         sw=SauvWriter.New();
303         sw.setMEDFileDS(mfd);
304         sw.write(sauvFile);
305
306         # convert the SAUV file back to MED
307         sr=SauvReader.New(sauvFile);
308         d2=sr.loadInMEDFileDS();
309
310         self.assertEqual( 1, d2.getNumberOfFields() )
311         self.assertEqual( 1, d2.getNumberOfMeshes() )
312         mfm2 = d2.getMeshes()[0]
313         mff2 = d2.getFields()[0]
314         m2 = mfm2.getMeshAtLevel(0)
315         f2 = mff2.getTimeStepAtPos(0).getFieldOnMeshAtLevel(f.getTypeOfField(),0,mfm2)
316         f2.setGaussLocalizationOnType(NORM_QUAD8,refCoo,gpCoo,wgt) # not stored in SAUV
317         #f2.setOrder( f.getTime()[2] ) # not stored in SAUV
318         self.assertTrue( m2.isEqual( m, 1e-12 ))
319         self.assertTrue( f2.isEqual( f, 1e-12, 1e-12 ))
320
321         os.remove( sauvFile )
322
323
324         pass
325     pass
326
327 unittest.main()