1 // Copyright (C) 2007-2016 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 #include "SauvLoaderTest.hxx"
22 #include "SauvReader.hxx"
23 #include "SauvWriter.hxx"
24 #include "MEDFileData.hxx"
25 #include "MEDCouplingFieldDouble.hxx"
26 #include "MEDCouplingMemArray.hxx"
27 #include "TestInterpKernelUtils.hxx" // getResourceFile()
38 using namespace MEDCoupling;
40 void SauvLoaderTest::testSauv2Med()
42 // read a file containing all types of readable piles
43 std::string file = INTERP_TEST::getResourceFile("allPillesTest.sauv", 3);
44 MCAuto<SauvReader> sr=SauvReader::New(file.c_str());
45 MCAuto<MEDFileData> d2=sr->loadInMEDFileDS();
47 d2->write("allPillesTest.med",0);
49 CPPUNIT_ASSERT_EQUAL(1,d2->getNumberOfMeshes());
50 CPPUNIT_ASSERT_EQUAL(8+97,d2->getNumberOfFields());
51 MEDFileMesh * m = d2->getMeshes()->getMeshAtPos(0);
52 CPPUNIT_ASSERT_EQUAL(17,int(m->getGroupsNames().size()));
55 void SauvLoaderTest::testMed2SauvOnAMeshWithVoidFamily()
57 // Create a mesh with 2 quads.
58 const int spaceDim = 2;
59 const int nbOfNodes = 6;
60 double coords[nbOfNodes*spaceDim] = {0,0, 1,0, 1,1, 0,1, 2,0, 2,1};
61 int conn[8]={0,1,2,3, 1,4,5,2};
62 MCAuto<MEDCouplingUMesh> mesh2d=MEDCouplingUMesh::New("Mesh",spaceDim);
63 mesh2d->allocateCells(2);
64 mesh2d->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);
65 mesh2d->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+4);
66 mesh2d->finishInsertingCells();
67 MCAuto<DataArrayDouble> myCoords=DataArrayDouble::New();
68 myCoords->alloc(nbOfNodes,spaceDim);
69 std::copy(coords,coords+nbOfNodes*spaceDim,myCoords->getPointer());
70 mesh2d->setCoords(myCoords);
72 // create a MedFileUMesh
73 MCAuto<MEDFileUMesh> m= MEDFileUMesh::New();
74 m->setMeshAtLevel(0,mesh2d);
76 // Create families and groups
78 MCAuto<DataArrayInt> fam = DataArrayInt::New();
80 int elemsFams[2] = {-2,-3};
81 std::copy(elemsFams,elemsFams+2,fam->getPointer());
83 m->setFamilyFieldArr(0,fam);
85 std::map<std::string,int> theFamilies;
86 theFamilies["FAM_-1"]=-1;
87 theFamilies["FAM_-2"]=-2;
88 theFamilies["FAM_-3"]=-3;
90 std::map<std::string, std::vector<std::string> > theGroups;
91 theGroups["Group1"].push_back("FAM_-2");
92 theGroups["Group2"].push_back("FAM_-3");
93 theGroups["Grouptot"].push_back("FAM_-1");
94 theGroups["Grouptot"].push_back("FAM_-2");
95 theGroups["Grouptot"].push_back("FAM_-3");
96 m->setFamilyInfo(theFamilies);
97 m->setGroupInfo(theGroups);
99 // write to MED for visual check
100 //const char* medFile = "mesh_with_void_family.med";
101 //m->write(medFile, 2);
104 const char* sauvFile = "mesh_with_void_family.sauv";
105 MCAuto<MEDFileData> medData = MEDFileData::New();
106 MCAuto<MEDFileMeshes> medMeshes = MEDFileMeshes::New();
107 MCAuto<SauvWriter> sw=SauvWriter::New();
108 medMeshes->setMeshAtPos(0, m);
109 medData->setMeshes(medMeshes);
110 sw->setMEDFileDS(medData);
113 // read SAUV and check groups
114 MCAuto<SauvReader> sr=SauvReader::New(sauvFile);
115 MCAuto<MEDFileData> d2=sr->loadInMEDFileDS();
116 MEDFileUMesh* m2 = static_cast<MEDFileUMesh*>( d2->getMeshes()->getMeshAtPos(0) );
117 MCAuto<MEDCouplingUMesh> group1 = m2->getGroup(0, "Group1");
118 CPPUNIT_ASSERT_EQUAL(1,(int)group1->getNumberOfCells());
119 MCAuto<MEDCouplingUMesh> group2 = m2->getGroup(0, "Group2");
120 CPPUNIT_ASSERT_EQUAL(1,(int)group2->getNumberOfCells());
121 MCAuto<MEDCouplingUMesh> grptot = m2->getGroup(0, "Grouptot");
122 CPPUNIT_ASSERT_EQUAL(2,(int)grptot->getNumberOfCells());
125 void SauvLoaderTest::testSauv2MedOnA3SubsField()
128 std::string sauvFile = INTERP_TEST::getResourceFile("portico_3subs.sauv", 3);
129 MCAuto<SauvReader> sr=SauvReader::New(sauvFile.c_str());
130 MCAuto<MEDFileData> d2=sr->loadInMEDFileDS();
132 MEDFileUMesh* m2 = static_cast<MEDFileUMesh*>(d2->getMeshes()->getMeshAtPos(0));
133 MCAuto<MEDCouplingUMesh> mesh1d = m2->getMeshAtLevel(0);
134 MCAuto<MEDCouplingFieldDouble> length1dField = mesh1d->getMeasureField(0);
135 std::cout << "Length of 1d elements: " << length1dField->accumulate(0) << std::endl;
136 CPPUNIT_ASSERT_DOUBLES_EQUAL(3, length1dField->accumulate(0), 1e-12);
138 MCAuto<MEDFileFieldMultiTS> field =
139 dynamic_cast<MEDFileFieldMultiTS *>(d2->getFields()->getFieldWithName("CHAM1D"));
140 std::cout << "Number of components in field: " << field->getInfo().size() << std::endl;
141 CPPUNIT_ASSERT_EQUAL(6,(int)field->getInfo().size());
142 std::vector< std::pair<int,int> > timesteps = field->getIterations();
144 MCAuto<MEDCouplingFieldDouble> field1d =
145 field->getFieldOnMeshAtLevel(ON_GAUSS_NE, timesteps[0].first, timesteps[0].second, 0, m2);
147 // Check first component of the field
148 // 2 gauss points per element => 12 values
149 double values[12] = {
161 -6.111413346910e-07};
163 for (int i=0; i < field1d->getNumberOfTuples(); i++)
165 CPPUNIT_ASSERT_DOUBLES_EQUAL( values[i], field1d->getIJ(i, 0), 1e-12 );
169 void SauvLoaderTest::testMed2Sauv()
172 std::string file = INTERP_TEST::getResourceFile("pointe.med", 3);
173 MCAuto<MEDFileData> pointeMed=MEDFileData::New(file.c_str());
175 // add 3 faces to pointeMed
176 MEDFileUMesh* pointeMedMesh = static_cast<MEDFileUMesh*>(pointeMed->getMeshes()->getMeshAtPos(0));
177 MCAuto<MEDCouplingUMesh> pointeM1D = MEDCouplingUMesh::New();
178 DataArrayDouble *coords = pointeMedMesh->getCoords();
179 pointeM1D->setCoords( coords );
180 pointeM1D->setMeshDimension( 2 );
181 pointeM1D->allocateCells( 3 );
184 0,1,2, 0,1,3, 10,11,12,13
186 pointeM1D->insertNextCell( INTERP_KERNEL::NORM_TRI3, 3, conn);
187 pointeM1D->insertNextCell( INTERP_KERNEL::NORM_TRI3, 3, conn+3);
188 pointeM1D->insertNextCell( INTERP_KERNEL::NORM_QUAD4, 4, conn+6);
189 pointeM1D->finishInsertingCells();
190 pointeMedMesh->setMeshAtLevel( -1, pointeM1D );
191 pointeMed->getMeshes()->setMeshAtPos( 0, pointeMedMesh );
193 // add a field on 2 faces to pointeMed
194 MCAuto<MEDFileFieldMultiTS> ff1=MEDFileFieldMultiTS::New();
195 MCAuto<MEDCouplingFieldDouble> f1=MEDCouplingFieldDouble::New(ON_GAUSS_NE,ONE_TIME);
196 f1->setMesh( pointeM1D );
197 f1->setName("Field on 2 faces");
198 MCAuto<DataArrayDouble> d=DataArrayDouble::New();
200 d->setInfoOnComponent(0,"sigX [MPa]");
201 d->setInfoOnComponent(1,"sigY [GPa]");
202 double vals[2*(3+4)] =
204 311,312,321,322,331,332,411,412,421,422,431,432,441,442
206 std::copy(vals,vals+d->getNbOfElems(),d->getPointer());
208 MCAuto<DataArrayInt> da=DataArrayInt::New();
214 std::copy(ids,ids+da->getNbOfElems(),da->getPointer());
216 ff1->appendFieldProfile(f1,pointeMedMesh,-1,da);
217 pointeMed->getFields()->pushField( ff1 );
219 // remove "fieldnodeint"
220 MEDFileFields* pointeFields = pointeMed->getFields();
221 for ( int i = 0; i < pointeFields->getNumberOfFields(); ++i )
223 MCAuto<MEDFileAnyTypeFieldMultiTS> ts = pointeFields->getFieldAtPos(i);
224 if ( std::string("fieldnodeint") == ts->getName())
226 pointeFields->destroyFieldAtPos( i );
230 // write pointeMed to SAUV
231 const char* sauvFile = "pointe.sauv";
232 MCAuto<SauvWriter> sw=SauvWriter::New();
233 sw->setMEDFileDS(pointeMed);
236 // read SAUV and check
237 MCAuto<SauvReader> sr=SauvReader::New(sauvFile);
238 MCAuto<MEDFileData> d2=sr->loadInMEDFileDS();
239 CPPUNIT_ASSERT_EQUAL(1,d2->getNumberOfMeshes());
240 CPPUNIT_ASSERT_EQUAL(4,d2->getNumberOfFields());
241 MEDFileUMesh * m = static_cast<MEDFileUMesh*>( d2->getMeshes()->getMeshAtPos(0) );
242 CPPUNIT_ASSERT_EQUAL(std::string("maa1"),std::string(m->getName() ));
243 CPPUNIT_ASSERT_EQUAL(3,m->getMeshDimension());
244 std::vector<std::string > groups = m->getGroupsNames();
245 CPPUNIT_ASSERT_EQUAL(6,(int)groups.size());
246 CPPUNIT_ASSERT( std::find(groups.begin(),groups.end(),"groupe1") != groups.end() );
247 CPPUNIT_ASSERT( std::find(groups.begin(),groups.end(),"groupe2") != groups.end() );
248 CPPUNIT_ASSERT( std::find(groups.begin(),groups.end(),"groupe3") != groups.end() );
249 CPPUNIT_ASSERT( std::find(groups.begin(),groups.end(),"groupe4") != groups.end() );
250 CPPUNIT_ASSERT( std::find(groups.begin(),groups.end(),"groupe5") != groups.end() );
251 CPPUNIT_ASSERT( std::find(groups.begin(),groups.end(),"maa1") != groups.end() );
252 CPPUNIT_ASSERT_EQUAL(16,m->getSizeAtLevel(0));
253 MCAuto<MEDCouplingMesh> um0 = m->getMeshAtLevel(0);
254 CPPUNIT_ASSERT_EQUAL(12, um0->getNumberOfCellsWithType( INTERP_KERNEL::NORM_TETRA4 ));
255 CPPUNIT_ASSERT_EQUAL(2, um0->getNumberOfCellsWithType( INTERP_KERNEL::NORM_PYRA5 ));
256 CPPUNIT_ASSERT_EQUAL(2, um0->getNumberOfCellsWithType( INTERP_KERNEL::NORM_HEXA8 ));
257 MCAuto<MEDCouplingMesh> um1 = m->getMeshAtLevel(-1);
258 CPPUNIT_ASSERT_EQUAL(1, um1->getNumberOfCellsWithType( INTERP_KERNEL::NORM_TRI3 ));
259 //CPPUNIT_ASSERT_EQUAL(2, um1->getNumberOfCellsWithType( INTERP_KERNEL::NORM_TRI3 ));
260 MCAuto<MEDCouplingUMesh> pointeUM0 =
261 static_cast<MEDCouplingUMesh*>( pointeMedMesh->getMeshAtLevel(0));
262 DataArrayDouble *coo = m->getCoords();
263 DataArrayDouble *pointeCoo = pointeMedMesh->getCoords();
264 CPPUNIT_ASSERT(coo->isEqualWithoutConsideringStr(*pointeCoo,1e-12));
265 MCAuto<MEDCouplingFieldDouble> vol = um0->getMeasureField(0);
266 MCAuto<MEDCouplingFieldDouble> pointeVol = pointeUM0->getMeasureField(0);
267 CPPUNIT_ASSERT_DOUBLES_EQUAL( vol->accumulate(0), pointeVol->accumulate(0),1e-12);
270 MCAuto<MEDFileFieldMultiTS> fieldnodedoubleTS1 =
271 dynamic_cast<MEDFileFieldMultiTS *>(pointeMed->getFields()->getFieldWithName("fieldnodedouble"));
272 MCAuto<MEDFileFieldMultiTS> fieldnodedoubleTS2 =
273 dynamic_cast<MEDFileFieldMultiTS *>(d2->getFields()->getFieldWithName("fieldnodedouble"));
274 CPPUNIT_ASSERT_EQUAL( fieldnodedoubleTS1->getInfo().size(), fieldnodedoubleTS2->getInfo().size());
275 for ( size_t i = 0; i < fieldnodedoubleTS1->getInfo().size(); ++i )
276 CPPUNIT_ASSERT_EQUAL( fieldnodedoubleTS1->getInfo()[i], fieldnodedoubleTS2->getInfo()[i]);
277 CPPUNIT_ASSERT_EQUAL( fieldnodedoubleTS1->getNumberOfTS(), fieldnodedoubleTS2->getNumberOfTS());
278 std::vector< std::pair<int,int> > io1 = fieldnodedoubleTS1->getIterations();
279 std::vector< std::pair<int,int> > io2 = fieldnodedoubleTS2->getIterations();
280 for ( int i =0; i < fieldnodedoubleTS1->getNumberOfTS(); ++i )
282 MCAuto<MEDCouplingFieldDouble> fnd1 =
283 fieldnodedoubleTS1->getFieldOnMeshAtLevel(ON_NODES, io1[i].first,io1[i].second,pointeUM0);
284 MCAuto<MEDCouplingFieldDouble> fnd2 =
285 fieldnodedoubleTS2->getFieldOnMeshAtLevel(ON_NODES, io2[i].first,io2[i].second,um0);
286 CPPUNIT_ASSERT( fnd1->getArray()->isEqual( *fnd2->getArray(), 1e-12 ));
288 // fieldcelldoublevector
289 MCAuto<MEDFileFieldMultiTS> fieldcelldoublevectorTS1 =
290 dynamic_cast<MEDFileFieldMultiTS *>(pointeMed->getFields()->getFieldWithName("fieldcelldoublevector"));
291 MCAuto<MEDFileFieldMultiTS> fieldcelldoublevectorTS2 =
292 dynamic_cast<MEDFileFieldMultiTS *>(d2->getFields()->getFieldWithName("fieldcelldoublevector"));
293 CPPUNIT_ASSERT_EQUAL( fieldcelldoublevectorTS1->getInfo().size(), fieldcelldoublevectorTS2->getInfo().size());
294 for ( size_t i = 0; i < fieldcelldoublevectorTS1->getInfo().size(); ++i )
295 CPPUNIT_ASSERT_EQUAL( fieldcelldoublevectorTS1->getInfo()[i], fieldcelldoublevectorTS2->getInfo()[i]);
296 CPPUNIT_ASSERT_EQUAL( fieldcelldoublevectorTS1->getNumberOfTS(), fieldcelldoublevectorTS2->getNumberOfTS());
297 io1 = fieldcelldoublevectorTS1->getIterations();
298 io2 = fieldcelldoublevectorTS2->getIterations();
299 for ( int i =0; i < fieldcelldoublevectorTS1->getNumberOfTS(); ++i )
301 MCAuto<MEDCouplingFieldDouble> fnd1 =
302 fieldcelldoublevectorTS1->getFieldOnMeshAtLevel(ON_CELLS, io1[i].first,io1[i].second,pointeUM0);
303 MCAuto<MEDCouplingFieldDouble> fnd2 =
304 fieldcelldoublevectorTS2->getFieldOnMeshAtLevel(ON_CELLS, io2[i].first,io2[i].second,um0);
305 CPPUNIT_ASSERT_DOUBLES_EQUAL( fnd1->accumulate(0), fnd2->accumulate(0), 1e-12 );
306 CPPUNIT_ASSERT_DOUBLES_EQUAL( fnd1->accumulate(1), fnd2->accumulate(1), 1e-12 );
307 CPPUNIT_ASSERT_DOUBLES_EQUAL( fnd1->accumulate(2), fnd2->accumulate(2), 1e-12 );
309 // "Field on 2 faces"
310 MCAuto<MEDFileFieldMultiTS> fieldOnFaces =
311 dynamic_cast<MEDFileFieldMultiTS *>(d2->getFields()->getFieldWithName(f1->getName().c_str()));
312 io1 = fieldOnFaces->getIterations();
313 MCAuto<MEDCouplingFieldDouble> fof =
314 fieldOnFaces->getFieldOnMeshAtLevel(f1->getTypeOfField(),io1[0].first,io1[0].second,um1);
315 CPPUNIT_ASSERT( d->isEqual( *fof->getArray(), 1e-12 ));
318 void SauvLoaderTest::testCellsWithLingNames()
320 // test IMP 3285: [CEA 1778] SauvReader: only keep the meshes named in the table MED_MAIL
321 std::string file = INTERP_TEST::getResourceFile("test_MED_MAIL.sauv", 3);
322 MCAuto<SauvReader> sr=SauvReader::New(file.c_str());
323 MCAuto<MEDFileData> d2=sr->loadInMEDFileDS();
324 // check that the mesh contains
325 // - Nombre de noeuds : 74
326 // - Nombre de mailles de type MED_TRIA3 : 6
327 // - Nombre de mailles de type MED_QUAD4 : 43
328 // - Nombre de mailles de type MED_HEXA8 : 24
329 // - Nombre de mailles de type MED_PENTA6 : 3
330 MEDFileUMesh* m = static_cast<MEDFileUMesh*>( d2->getMeshes()->getMeshAtPos(0));
331 CPPUNIT_ASSERT_EQUAL(6, m->getNumberOfCellsWithType( INTERP_KERNEL::NORM_TRI3 ));
332 CPPUNIT_ASSERT_EQUAL(43, m->getNumberOfCellsWithType( INTERP_KERNEL::NORM_QUAD4 ));
333 CPPUNIT_ASSERT_EQUAL(24, m->getNumberOfCellsWithType( INTERP_KERNEL::NORM_HEXA8 ));
334 CPPUNIT_ASSERT_EQUAL(3, m->getNumberOfCellsWithType( INTERP_KERNEL::NORM_PENTA6 ));
337 void SauvLoaderTest::tearDown()
339 const int nbFilesToRemove = 3;
340 const char* fileToRemove[nbFilesToRemove] = { "allPillesTest.med", "pointe.sauv", "mesh_with_void_family.sauv" };
341 for ( int i = 0; i < nbFilesToRemove; ++i )
344 if (GetFileAttributes(fileToRemove[i]) != INVALID_FILE_ATTRIBUTES)
346 if (access(fileToRemove[i], F_OK) == 0)
348 remove(fileToRemove[i]);