Salome HOME
b6195f6df1abfefeea187b429a92e2c50e308779
[modules/med.git] / src / medtool / src / MEDLoader / Test / SauvLoaderTest.cxx
1 // Copyright (C) 2007-2015  CEA/DEN, EDF R&D
2 //
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.
7 //
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.
12 //
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
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19
20 #include "SauvLoaderTest.hxx"
21
22 #include "SauvReader.hxx"
23 #include "SauvWriter.hxx"
24 #include "MEDFileData.hxx"
25 #include "MEDCouplingFieldDouble.hxx"
26 #include "MEDCouplingMemArray.hxx"
27
28 #ifdef WIN32
29 # include <windows.h>
30 #else
31 # include <unistd.h>
32 #endif
33
34 #include <vector>
35 #include <string>
36
37 using namespace ParaMEDMEM;
38
39 void SauvLoaderTest::testSauv2Med()
40 {
41   // read a file containing all types of readable piles
42   std::string file = getResourceFile("allPillesTest.sauv");
43   MEDCouplingAutoRefCountObjectPtr<SauvReader> sr=SauvReader::New(file.c_str());
44   MEDCouplingAutoRefCountObjectPtr<MEDFileData> d2=sr->loadInMEDFileDS();
45   // write MED
46   d2->write("allPillesTest.med",0);
47   // check 
48   CPPUNIT_ASSERT_EQUAL(1,d2->getNumberOfMeshes());
49   CPPUNIT_ASSERT_EQUAL(8+97,d2->getNumberOfFields());
50   MEDFileMesh * m = d2->getMeshes()->getMeshAtPos(0);
51   CPPUNIT_ASSERT_EQUAL(17,int(m->getGroupsNames().size()));
52 }
53
54 void SauvLoaderTest::testMed2SauvOnAMeshWithVoidFamily()
55 {
56   // Create a mesh with 2 quads.
57   const int spaceDim = 2;
58   const int nbOfNodes = 6;
59   double coords[nbOfNodes*spaceDim] = {0,0, 1,0, 1,1, 0,1, 2,0, 2,1};
60   int conn[8]={0,1,2,3, 1,4,5,2};
61   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> mesh2d=MEDCouplingUMesh::New("Mesh",spaceDim);
62   mesh2d->allocateCells(2);
63   mesh2d->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);
64   mesh2d->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+4);
65   mesh2d->finishInsertingCells();
66   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> myCoords=DataArrayDouble::New();
67   myCoords->alloc(nbOfNodes,spaceDim);
68   std::copy(coords,coords+nbOfNodes*spaceDim,myCoords->getPointer());
69   mesh2d->setCoords(myCoords);
70
71   // create a MedFileUMesh
72   MEDCouplingAutoRefCountObjectPtr<MEDFileUMesh> m= MEDFileUMesh::New();
73   m->setMeshAtLevel(0,mesh2d);
74
75   // Create families and groups
76
77   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> fam = DataArrayInt::New();
78   fam->alloc(2,1);
79   int elemsFams[2] = {-2,-3};
80   std::copy(elemsFams,elemsFams+2,fam->getPointer());
81
82   m->setFamilyFieldArr(0,fam);
83
84   std::map<std::string,int> theFamilies;
85   theFamilies["FAM_-1"]=-1;
86   theFamilies["FAM_-2"]=-2;
87   theFamilies["FAM_-3"]=-3;
88
89   std::map<std::string, std::vector<std::string> > theGroups;
90   theGroups["Group1"].push_back("FAM_-2");
91   theGroups["Group2"].push_back("FAM_-3");
92   theGroups["Grouptot"].push_back("FAM_-1");
93   theGroups["Grouptot"].push_back("FAM_-2");
94   theGroups["Grouptot"].push_back("FAM_-3");
95   m->setFamilyInfo(theFamilies);
96   m->setGroupInfo(theGroups);
97
98   // write to MED for visual check
99   //const char* medFile = "mesh_with_void_family.med";
100   //m->write(medFile, 2);
101
102   // write to SAUV
103   const char* sauvFile = "mesh_with_void_family.sauv";
104   MEDCouplingAutoRefCountObjectPtr<MEDFileData> medData = MEDFileData::New();
105   MEDCouplingAutoRefCountObjectPtr<MEDFileMeshes> medMeshes = MEDFileMeshes::New();
106   MEDCouplingAutoRefCountObjectPtr<SauvWriter> sw=SauvWriter::New();
107   medMeshes->setMeshAtPos(0, m);
108   medData->setMeshes(medMeshes);
109   sw->setMEDFileDS(medData);
110   sw->write(sauvFile);
111
112   // read SAUV and check groups
113   MEDCouplingAutoRefCountObjectPtr<SauvReader> sr=SauvReader::New(sauvFile);
114   MEDCouplingAutoRefCountObjectPtr<MEDFileData> d2=sr->loadInMEDFileDS();
115   MEDFileUMesh* m2 = static_cast<MEDFileUMesh*>( d2->getMeshes()->getMeshAtPos(0) );
116   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> group1 = m2->getGroup(0, "Group1");
117   CPPUNIT_ASSERT_EQUAL(1,(int)group1->getNumberOfCells());
118   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> group2 = m2->getGroup(0, "Group2");
119   CPPUNIT_ASSERT_EQUAL(1,(int)group2->getNumberOfCells());
120   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> grptot = m2->getGroup(0, "Grouptot");
121   CPPUNIT_ASSERT_EQUAL(2,(int)grptot->getNumberOfCells());
122 }
123
124 void SauvLoaderTest::testSauv2MedOnA3SubsField()
125 {
126   // read SAUV
127   std::string sauvFile = getResourceFile("portico_3subs.sauv");
128   MEDCouplingAutoRefCountObjectPtr<SauvReader> sr=SauvReader::New(sauvFile.c_str());
129   MEDCouplingAutoRefCountObjectPtr<MEDFileData> d2=sr->loadInMEDFileDS();
130   // check mesh
131   MEDFileUMesh* m2 = static_cast<MEDFileUMesh*>(d2->getMeshes()->getMeshAtPos(0));
132   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> mesh1d = m2->getMeshAtLevel(0);
133   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> length1dField = mesh1d->getMeasureField(0);
134   std::cout << "Length of 1d elements: " << length1dField->accumulate(0) << std::endl;
135   CPPUNIT_ASSERT_DOUBLES_EQUAL(3, length1dField->accumulate(0), 1e-12);
136   // check field
137   MEDCouplingAutoRefCountObjectPtr<MEDFileFieldMultiTS> field =
138     dynamic_cast<MEDFileFieldMultiTS *>(d2->getFields()->getFieldWithName("CHAM1D"));
139   std::cout << "Number of components in field: " << field->getInfo().size() << std::endl;
140   CPPUNIT_ASSERT_EQUAL(6,(int)field->getInfo().size());
141   std::vector< std::pair<int,int> > timesteps = field->getIterations();
142
143   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> field1d =
144     field->getFieldOnMeshAtLevel(ON_GAUSS_NE, timesteps[0].first, timesteps[0].second, 0, m2);
145
146   // Check first component of the field
147   // 2 gauss points per element => 12 values
148   double values[12] = {
149       -7.687500000000e-03,
150       -7.687500000000e-03,
151       -4.562500000000e-03,
152       -4.562500000000e-03,
153       -8.208333333333e-03,
154       -8.208333333333e-03,
155       -6.125000000000e-03,
156       -6.125000000000e-03,
157       -4.041666666666e-03,
158       -4.041666666666e-03,
159       -6.111413346910e-07,
160       -6.111413346910e-07};
161
162   for (int i=0; i < field1d->getNumberOfTuples(); i++)
163   {
164     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[i], field1d->getIJ(i, 0), 1e-12 );
165   }
166 }
167
168 void SauvLoaderTest::testMed2Sauv()
169 {
170   // read pointe.med
171   std::string file = getResourceFile("pointe.med");
172   MEDCouplingAutoRefCountObjectPtr<MEDFileData> pointeMed=MEDFileData::New(file.c_str());
173
174   // add 3 faces to pointeMed
175   MEDFileUMesh* pointeMedMesh = static_cast<MEDFileUMesh*>(pointeMed->getMeshes()->getMeshAtPos(0));
176   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> pointeM1D = MEDCouplingUMesh::New();
177   DataArrayDouble     *coords = pointeMedMesh->getCoords();
178   pointeM1D->setCoords( coords );
179   pointeM1D->setMeshDimension( 2 );
180   pointeM1D->allocateCells( 3 );
181   int conn[]=
182     {
183       0,1,2, 0,1,3, 10,11,12,13
184     };
185   pointeM1D->insertNextCell( INTERP_KERNEL::NORM_TRI3, 3, conn);
186   pointeM1D->insertNextCell( INTERP_KERNEL::NORM_TRI3, 3, conn+3);
187   pointeM1D->insertNextCell( INTERP_KERNEL::NORM_QUAD4, 4, conn+6);
188   pointeM1D->finishInsertingCells();
189   pointeMedMesh->setMeshAtLevel( -1, pointeM1D );
190   pointeMed->getMeshes()->setMeshAtPos( 0, pointeMedMesh );
191
192   // add a field on 2 faces to pointeMed
193   MEDCouplingAutoRefCountObjectPtr<MEDFileFieldMultiTS> ff1=MEDFileFieldMultiTS::New();
194   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> f1=MEDCouplingFieldDouble::New(ON_GAUSS_NE,ONE_TIME);
195   f1->setMesh( pointeM1D );
196   f1->setName("Field on 2 faces");
197   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> d=DataArrayDouble::New();
198   d->alloc(3+4,2);
199   d->setInfoOnComponent(0,"sigX [MPa]");
200   d->setInfoOnComponent(1,"sigY [GPa]");
201   double vals[2*(3+4)] =
202     {
203       311,312,321,322,331,332,411,412,421,422,431,432,441,442
204     };
205   std::copy(vals,vals+d->getNbOfElems(),d->getPointer());
206   f1->setArray(d);
207   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> da=DataArrayInt::New();
208   int ids[] =
209     {
210       0,2
211     };
212   da->alloc(2,1);
213   std::copy(ids,ids+da->getNbOfElems(),da->getPointer());
214   da->setName("sup2");
215   ff1->appendFieldProfile(f1,pointeMedMesh,-1,da);
216   pointeMed->getFields()->pushField( ff1 );
217
218   // remove "fieldnodeint"
219   MEDFileFields* pointeFields = pointeMed->getFields();
220   for ( int i = 0; i < pointeFields->getNumberOfFields(); ++i )
221     {
222       MEDCouplingAutoRefCountObjectPtr<MEDFileAnyTypeFieldMultiTS> ts = pointeFields->getFieldAtPos(i);
223       if ( std::string("fieldnodeint") == ts->getName())
224         {
225           pointeFields->destroyFieldAtPos( i );
226           break;
227         }
228     }
229   // write pointeMed to SAUV
230   const char* sauvFile = "pointe.sauv";
231   MEDCouplingAutoRefCountObjectPtr<SauvWriter> sw=SauvWriter::New();
232   sw->setMEDFileDS(pointeMed);
233   sw->write(sauvFile);
234
235   // read SAUV and check
236   MEDCouplingAutoRefCountObjectPtr<SauvReader> sr=SauvReader::New(sauvFile);
237   MEDCouplingAutoRefCountObjectPtr<MEDFileData> d2=sr->loadInMEDFileDS();
238   CPPUNIT_ASSERT_EQUAL(1,d2->getNumberOfMeshes());
239   CPPUNIT_ASSERT_EQUAL(4,d2->getNumberOfFields());
240   MEDFileUMesh * m = static_cast<MEDFileUMesh*>( d2->getMeshes()->getMeshAtPos(0) );
241   CPPUNIT_ASSERT_EQUAL(std::string("maa1"),std::string(m->getName() ));
242   CPPUNIT_ASSERT_EQUAL(3,m->getMeshDimension());
243   std::vector<std::string > groups = m->getGroupsNames();
244   CPPUNIT_ASSERT_EQUAL(6,(int)groups.size());
245   CPPUNIT_ASSERT( std::find(groups.begin(),groups.end(),"groupe1") != groups.end() );
246   CPPUNIT_ASSERT( std::find(groups.begin(),groups.end(),"groupe2") != groups.end() );
247   CPPUNIT_ASSERT( std::find(groups.begin(),groups.end(),"groupe3") != groups.end() );
248   CPPUNIT_ASSERT( std::find(groups.begin(),groups.end(),"groupe4") != groups.end() );
249   CPPUNIT_ASSERT( std::find(groups.begin(),groups.end(),"groupe5") != groups.end() );
250   CPPUNIT_ASSERT( std::find(groups.begin(),groups.end(),"maa1") != groups.end() );
251   CPPUNIT_ASSERT_EQUAL(16,m->getSizeAtLevel(0));
252   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> um0 = m->getGenMeshAtLevel(0);
253   CPPUNIT_ASSERT_EQUAL(12, um0->getNumberOfCellsWithType( INTERP_KERNEL::NORM_TETRA4 ));
254   CPPUNIT_ASSERT_EQUAL(2,  um0->getNumberOfCellsWithType( INTERP_KERNEL::NORM_PYRA5 ));
255   CPPUNIT_ASSERT_EQUAL(2,  um0->getNumberOfCellsWithType( INTERP_KERNEL::NORM_HEXA8 ));
256   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> um1 = m->getGenMeshAtLevel(-1);
257   CPPUNIT_ASSERT_EQUAL(2, um1->getNumberOfCellsWithType( INTERP_KERNEL::NORM_TRI3 ));
258   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> pointeUM0 =
259     static_cast<MEDCouplingUMesh*>( pointeMedMesh->getGenMeshAtLevel(0));
260   DataArrayDouble *coo = m->getCoords();
261   DataArrayDouble *pointeCoo = pointeMedMesh->getCoords();
262   CPPUNIT_ASSERT(coo->isEqualWithoutConsideringStr(*pointeCoo,1e-12));
263   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol = um0->getMeasureField(0);
264   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> pointeVol = pointeUM0->getMeasureField(0);
265   CPPUNIT_ASSERT_DOUBLES_EQUAL( vol->accumulate(0), pointeVol->accumulate(0),1e-12);
266   // check fields
267   // fieldnodedouble
268   MEDCouplingAutoRefCountObjectPtr<MEDFileFieldMultiTS> fieldnodedoubleTS1 =
269     dynamic_cast<MEDFileFieldMultiTS *>(pointeMed->getFields()->getFieldWithName("fieldnodedouble"));
270   MEDCouplingAutoRefCountObjectPtr<MEDFileFieldMultiTS> fieldnodedoubleTS2 =
271     dynamic_cast<MEDFileFieldMultiTS *>(d2->getFields()->getFieldWithName("fieldnodedouble"));
272   CPPUNIT_ASSERT_EQUAL( fieldnodedoubleTS1->getInfo().size(), fieldnodedoubleTS2->getInfo().size());
273   for ( size_t i = 0; i < fieldnodedoubleTS1->getInfo().size(); ++i )
274     CPPUNIT_ASSERT_EQUAL( fieldnodedoubleTS1->getInfo()[i], fieldnodedoubleTS2->getInfo()[i]);
275   CPPUNIT_ASSERT_EQUAL( fieldnodedoubleTS1->getNumberOfTS(), fieldnodedoubleTS2->getNumberOfTS());
276   std::vector< std::pair<int,int> > io1 = fieldnodedoubleTS1->getIterations();
277   std::vector< std::pair<int,int> > io2 = fieldnodedoubleTS2->getIterations();
278   for ( int i =0; i < fieldnodedoubleTS1->getNumberOfTS(); ++i )
279     {
280       MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> fnd1 =
281         fieldnodedoubleTS1->getFieldOnMeshAtLevel(ON_NODES, io1[i].first,io1[i].second,pointeUM0);
282       MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> fnd2 =
283         fieldnodedoubleTS2->getFieldOnMeshAtLevel(ON_NODES, io2[i].first,io2[i].second,um0);
284       CPPUNIT_ASSERT( fnd1->getArray()->isEqual( *fnd2->getArray(), 1e-12 ));
285     }
286   // fieldcelldoublevector
287   MEDCouplingAutoRefCountObjectPtr<MEDFileFieldMultiTS> fieldcelldoublevectorTS1 =
288     dynamic_cast<MEDFileFieldMultiTS *>(pointeMed->getFields()->getFieldWithName("fieldcelldoublevector"));
289   MEDCouplingAutoRefCountObjectPtr<MEDFileFieldMultiTS> fieldcelldoublevectorTS2 =
290     dynamic_cast<MEDFileFieldMultiTS *>(d2->getFields()->getFieldWithName("fieldcelldoublevector"));
291   CPPUNIT_ASSERT_EQUAL( fieldcelldoublevectorTS1->getInfo().size(), fieldcelldoublevectorTS2->getInfo().size());
292   for ( size_t i = 0; i < fieldcelldoublevectorTS1->getInfo().size(); ++i )
293     CPPUNIT_ASSERT_EQUAL( fieldcelldoublevectorTS1->getInfo()[i], fieldcelldoublevectorTS2->getInfo()[i]);
294   CPPUNIT_ASSERT_EQUAL( fieldcelldoublevectorTS1->getNumberOfTS(), fieldcelldoublevectorTS2->getNumberOfTS());
295   io1 = fieldcelldoublevectorTS1->getIterations();
296   io2 = fieldcelldoublevectorTS2->getIterations();
297   for ( int i =0; i < fieldcelldoublevectorTS1->getNumberOfTS(); ++i )
298     {
299       MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> fnd1 =
300         fieldcelldoublevectorTS1->getFieldOnMeshAtLevel(ON_CELLS, io1[i].first,io1[i].second,pointeUM0);
301       MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> fnd2 =
302         fieldcelldoublevectorTS2->getFieldOnMeshAtLevel(ON_CELLS, io2[i].first,io2[i].second,um0);
303       CPPUNIT_ASSERT_DOUBLES_EQUAL( fnd1->accumulate(0), fnd2->accumulate(0), 1e-12 );
304       CPPUNIT_ASSERT_DOUBLES_EQUAL( fnd1->accumulate(1), fnd2->accumulate(1), 1e-12 );
305       CPPUNIT_ASSERT_DOUBLES_EQUAL( fnd1->accumulate(2), fnd2->accumulate(2), 1e-12 );
306     }
307   // "Field on 2 faces"
308   MEDCouplingAutoRefCountObjectPtr<MEDFileFieldMultiTS> fieldOnFaces =
309     dynamic_cast<MEDFileFieldMultiTS *>(d2->getFields()->getFieldWithName(f1->getName().c_str()));
310   io1 = fieldOnFaces->getIterations();
311   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> fof =
312     fieldOnFaces->getFieldOnMeshAtLevel(f1->getTypeOfField(),io1[0].first,io1[0].second,um1);
313   CPPUNIT_ASSERT( d->isEqual( *fof->getArray(), 1e-12 ));
314 }
315
316 void SauvLoaderTest::tearDown()
317 {
318   const int nbFilesToRemove = 3;
319   const char* fileToRemove[nbFilesToRemove] = { "allPillesTest.med", "pointe.sauv", "mesh_with_void_family.sauv" };
320   for ( int i = 0; i < nbFilesToRemove; ++i )
321     {
322 #ifdef WIN32
323       if (GetFileAttributes(fileToRemove[i]) != INVALID_FILE_ATTRIBUTES)
324 #else
325         if (access(fileToRemove[i], F_OK) == 0)
326 #endif
327       remove(fileToRemove[i]);
328   }
329 }
330
331 std::string SauvLoaderTest::getResourceFile( const std::string& filename )
332 {
333   std::string resourceFile = "";
334
335   if ( getenv("top_srcdir") ) {
336     // we are in 'make test' step
337     resourceFile = getenv("top_srcdir");
338     resourceFile += "/resources/";
339   }
340   else if ( getenv("MED_ROOT_DIR") ) {
341     // use MED_ROOT_DIR env.var
342     resourceFile = getenv("MED_ROOT_DIR");
343     resourceFile += "/share/salome/resources/med/";
344   }
345   resourceFile += filename;
346 #ifdef WIN32
347   std::string fixedpath = resourceFile;
348   for ( int i=0; i < fixedpath.length(); ++i )
349     if (fixedpath[i] == '/')
350       fixedpath[i] = '\\';
351   return fixedpath;
352 #endif
353   return resourceFile;
354 }