Salome HOME
Small fix.
[tools/medcoupling.git] / src / MEDPartitioner / Test / MEDPARTITIONERTest.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 "MEDPARTITIONERTest.hxx"
21
22 #include "MEDPARTITIONER_MeshCollection.hxx"
23 #include "MEDPARTITIONER_ParallelTopology.hxx"
24 #include "MEDPARTITIONER_ParaDomainSelector.hxx"
25 #include "MEDPARTITIONER_Utils.hxx"
26
27 #include "CellModel.hxx"
28 #include "MEDFileMesh.hxx"
29 #include "MEDLoader.hxx"
30 #include "MEDLoaderBase.hxx"
31 #include "MEDCouplingUMesh.hxx"
32 #include "MEDCouplingExtrudedMesh.hxx"
33 #include "MEDCouplingFieldDouble.hxx"
34 #include "MEDCouplingMemArray.hxx"
35 #include "MEDCouplingMultiFields.hxx"
36
37 #include <cppunit/TestAssert.h>
38
39 #include <sstream>
40 #include <fstream>
41 #include <cmath>
42 #include <list>
43 #include <stdexcept>
44 #include <cstdlib>
45 #include <vector>
46 #ifdef WIN32
47 #include<direct.h>
48 #define getcwd _getcwd
49 #else
50 #include <unistd.h>
51 #endif
52
53 #ifdef HAVE_MPI
54 #include <mpi.h>
55 #endif
56
57 using namespace std;
58 using namespace ParaMEDMEM;
59 using namespace MEDPARTITIONER;
60
61 void MEDPARTITIONERTest::setSize(int ni, int nj, int nk)
62 {
63   this->_ni=ni;  //nb of hexa9
64   this->_nj=nj;
65   this->_nk=nk;
66   this->_ntot=_ni*_nj*_nk;
67   string ijk=IntToStr(ni)+"x"+IntToStr(nj)+"x"+IntToStr(nk);
68   this->_file_name="tmp_testMesh_"+ijk+".med";
69   this->_file_name_with_faces="tmp_testMeshWithFaces_"+ijk+".med";
70   string ij=IntToStr(ni)+"x"+IntToStr(nj);
71   this->_file_name2="tmp_testMesh_"+ij+".med";
72   this->_mesh_name="testMesh";
73 }
74
75 void MEDPARTITIONERTest::setSmallSize()
76 {
77   setSize(2,3,5); //nb of hexa9
78 }
79
80 void MEDPARTITIONERTest::setMedianSize()
81 {
82   setSize(20,30,50); //nb of hexa9
83 }
84
85 void MEDPARTITIONERTest::setbigSize()
86 {
87   setSize(200,300,500); //nb of hexa9
88 }
89
90 std::string MEDPARTITIONERTest::getPartitionerExe() const
91 {
92   std::string execName;
93   if ( getenv("MEDCOUPLING_ROOT_DIR") )
94     {
95       execName=getenv("MEDCOUPLING_ROOT_DIR");  //.../INSTALL/MED
96       execName+="/bin/medpartitioner";
97     }
98   else
99     {
100       execName = getcwd(NULL, 0);
101       execName += "/../../MEDPartitioner/medpartitioner";
102       if (! std::ifstream(execName.c_str()))
103         CPPUNIT_FAIL("Can't find medpartitioner, please set MEDCOUPLING_ROOT_DIR");
104     }
105   return execName;
106 }
107
108 // ============================================================================
109 /*!
110  *  Set up the environment called at every CPPUNIT_TEST ()
111  */
112 // ============================================================================
113 void MEDPARTITIONERTest::setUp()
114 {
115   this->_verbose=0;
116 #if defined(HAVE_MPI)
117   if (MyGlobals::_Rank==-1)  //do once only
118     {
119       MPI_Init(0,0);
120       MPI_Comm_size(MPI_COMM_WORLD, &MyGlobals::_World_Size);
121       MPI_Comm_rank(MPI_COMM_WORLD, &MyGlobals::_Rank);
122     }
123 #else
124   //sequential : no MPI
125   MyGlobals::_World_Size=1;
126   MyGlobals::_Rank=0;
127 #endif
128
129   if (_verbose>10)
130     {
131 #if defined(HAVE_MPI)
132       cout<<"\ndefined(HAVE_MPI)"<<endl;
133 #else
134       cout<<"\nNOT defined(HAVE_MPI)"<<endl;
135 #endif
136 #if defined(MED_ENABLE_PARMETIS)
137       cout<<"defined(MED_ENABLE_PARMETIS)"<<endl;
138 #else
139       cout<<"NOT defined(MED_ENABLE_PARMETIS)"<<endl;
140 #endif
141 #if defined(MED_ENABLE_METIS)
142       cout<<"defined(MED_ENABLE_METIS)"<<endl;
143 #else
144       cout<<"NOT defined(MED_ENABLE_METIS)"<<endl;
145 #endif
146 #if defined(MED_ENABLE_SCOTCH)
147       cout<<"defined(MED_ENABLE_SCOTCH)"<<endl;
148 #else
149       cout<<"NOT defined(MED_ENABLE_SCOTCH)"<<endl;
150 #endif
151     }
152 }
153
154 // ============================================================================
155 /*!
156  *  - delete
157  */
158 // ============================================================================
159 void MEDPARTITIONERTest::tearDown()
160 {
161 }
162
163 ParaMEDMEM::MEDCouplingUMesh * MEDPARTITIONERTest::buildCUBE3DMesh()
164 //only hexa8
165 {
166   vector<int> conn;
167   vector<double> coor;
168   for (int k=0; k<=_nk; k++)
169     for (int j=0; j<=_nj; j++)
170       for (int i=0; i<=_ni; i++)
171         {
172           coor.push_back(i+.1);
173           coor.push_back(j+.2);
174           coor.push_back(k+.3);
175         }
176   int ii;
177   for (int k=0; k<_nk; k++)
178     for (int j=0; j<_nj; j++)
179       for (int i=0; i<_ni; i++)
180         {
181           ii=i + j*(_ni+1) + k*(_ni+1)*(_nj+1);
182           conn.push_back(ii);
183           conn.push_back(ii+1);
184           ii=ii + _ni + 2 ;
185           conn.push_back(ii);
186           conn.push_back(ii-1);
187
188           ii=i + j*(_ni+1) + (k+1)*(_ni+1)*(_nj+1);
189           conn.push_back(ii);
190           conn.push_back(ii+1);
191           ii=ii + _ni + 2 ;
192           conn.push_back(ii);
193           conn.push_back(ii-1);
194         }
195
196   /*
197   if (_verbose)  //only for debug
198     {
199       cout<< "\nnb coor " << (_ni+1)*(_nj+1)*(_nk+1)*3 << " " << coor.size() << endl;
200       for (int i=0; i<(int)coor.size(); i++)
201         cout << coor[i] << " ";
202       cout << endl;
203       cout << "\nnb conn " << (_ni)*(_nj)*(_nk)*8 << " " << conn.size() << endl;
204       for (int i=0; i<(int)conn.size(); i=i+8)
205         {
206           for (int j=0; j<8; j++)
207             cout << conn[i+j] << " ";
208           cout << endl;
209         }
210       cout << endl;
211     }
212   */
213
214   MEDCouplingUMesh *mesh=MEDCouplingUMesh::New();
215   mesh->setMeshDimension(3);
216   int nbc=conn.size()/8; //nb of cells
217   int nbv=coor.size()/3; //nb of vertices
218   mesh->allocateCells(nbc);
219   for(int i=0; i<nbc; i++)
220     {
221       int onehexa[8];
222       std::copy(conn.begin()+i*8,conn.begin()+(i+1)*8,onehexa);
223       if (false) //(_verbose)
224         {
225           for (int j=0; j<8; j++) cout<<onehexa[j]<<" ";
226           cout<<endl;
227         }
228       mesh->insertNextCell(INTERP_KERNEL::NORM_HEXA8,8,onehexa);
229     }
230   mesh->finishInsertingCells();
231   DataArrayDouble *myCoords=DataArrayDouble::New();
232   myCoords->alloc(nbv,3);
233   std::copy(coor.begin(),coor.end(),myCoords->getPointer());
234   mesh->setCoords(myCoords);
235   mesh->setName(_mesh_name.c_str());
236   myCoords->decrRef();
237   mesh->checkCoherency();
238   return mesh;
239 }
240
241 ParaMEDMEM::MEDCouplingUMesh * MEDPARTITIONERTest::buildCARRE3DMesh()
242 //only quad4 in oblique (k=j)
243 {
244   vector<int> conn;
245   vector<double> coor;
246   for (int j=0; j<=_nj; j++)
247     for (int i=0; i<=_ni; i++)
248       {
249         int k=j;
250         coor.push_back(i+.1);
251         coor.push_back(j+.2);
252         coor.push_back(k+.3);
253       }
254   int ii;
255   int k=0;
256   for (int j=0; j<_nj; j++)
257     for (int i=0; i<_ni; i++)
258       {
259         ii=i + j*(_ni+1) + k*(_ni+1)*(_nj+1);
260         conn.push_back(ii);
261         conn.push_back(ii+1);
262         ii=ii + _ni + 2 ;
263         conn.push_back(ii);
264         conn.push_back(ii-1);
265       }
266
267   if (false) //(_verbose)
268     {
269       cout<<"\nnb coor "<<(_ni+1)*(_nj+1)*3<<" "<<coor.size()<<endl;
270       for (int i=0; i<(int)coor.size(); i++)
271         cout << coor[i] << " ";
272       cout<<endl;
273       cout<<"\nnb conn "<<(_ni)*(_nj)*4<<" "<<conn.size()<<endl;
274       for (int i=0; i<(int)conn.size(); i=i+4)
275         {
276           for (int j=0; j<4; j++) cout<<conn[i+j]<<" ";
277           cout<<endl;
278         }
279       cout<<endl;
280     }
281
282   MEDCouplingUMesh *mesh=MEDCouplingUMesh::New();
283   mesh->setMeshDimension(2);
284   int nbc=conn.size()/4; //nb of cells
285   int nbv=coor.size()/3; //nb of vertices
286   mesh->allocateCells(nbc);
287   for(int i=0; i<nbc; i++)
288     {
289       int onequa[4];
290       std::copy(conn.begin()+i*4,conn.begin()+(i+1)*4,onequa);
291       if (false) //(_verbose)
292         {
293           for (int j=0; j<4; j++) cout<<onequa[j]<<" ";
294           cout<<endl;
295         }
296       mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,onequa);
297     }
298   mesh->finishInsertingCells();
299   DataArrayDouble *myCoords=DataArrayDouble::New();
300   myCoords->alloc(nbv,3);
301   std::copy(coor.begin(),coor.end(),myCoords->getPointer());
302   mesh->setCoords(myCoords);
303   mesh->setName(_mesh_name.c_str());
304   myCoords->decrRef();
305   mesh->checkCoherency();
306   return mesh;
307 }
308
309 ParaMEDMEM::MEDCouplingUMesh * MEDPARTITIONERTest::buildFACE3DMesh()
310 //only quad4 on a global face of the CUBE3D (k=0)
311 {
312   vector<int> conn;
313   vector<double> coor;
314   for (int j=0; j<=_nj; j++)
315     for (int i=0; i<=_ni; i++)
316       {
317         int k=0;
318         coor.push_back(i+.1);
319         coor.push_back(j+.2);
320         coor.push_back(k+.3);
321       }
322   int ii;
323   int k=0;
324   for (int j=0; j<_nj; j++)
325     for (int i=0; i<_ni; i++)
326       {
327         ii=i + j*(_ni+1) + k*(_ni+1)*(_nj+1);
328         conn.push_back(ii);
329         conn.push_back(ii+1);
330         ii=ii + _ni + 2 ;
331         conn.push_back(ii);
332         conn.push_back(ii-1);
333       }
334
335   if (false) //(_verbose)
336     {
337       cout<<"\nnb coor "<<(_ni+1)*(_nj+1)*3<<" "<<coor.size()<<endl;
338       for (int i=0; i<(int)coor.size(); i++)
339         cout << coor[i] << " ";
340       cout<<endl;
341       cout<<"\nnb conn "<<(_ni)*(_nj)*4<<" "<<conn.size()<<endl;
342       for (int i=0; i<(int)conn.size(); i=i+4)
343         {
344           for (int j=0; j<4; j++)
345             cout << conn[i+j] << " ";
346           cout << endl;
347         }
348       cout << endl;
349     }
350
351   MEDCouplingUMesh *mesh=MEDCouplingUMesh::New();
352   mesh->setMeshDimension(2);
353   int nbc=conn.size()/4; //nb of cells
354   int nbv=coor.size()/3; //nb of vertices
355   mesh->allocateCells(nbc);
356   for(int i=0; i<nbc; i++)
357     {
358       int onequa[4];
359       std::copy(conn.begin()+i*4,conn.begin()+(i+1)*4,onequa);
360       if (false) //(_verbose)
361         {
362           for (int j=0; j<4; j++) cout<<onequa[j]<<" ";
363           cout<<endl;
364         }
365       mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,onequa);
366     }
367   mesh->finishInsertingCells();
368   DataArrayDouble *myCoords=DataArrayDouble::New();
369   myCoords->alloc(nbv,3);
370   std::copy(coor.begin(),coor.end(),myCoords->getPointer());
371   mesh->setCoords(myCoords);
372   mesh->setName(_mesh_name.c_str());
373   myCoords->decrRef();
374   mesh->checkCoherency();
375   return mesh;
376 }
377
378 MEDCouplingFieldDouble * MEDPARTITIONERTest::buildVecFieldOnCells(string myfileName)
379 {
380   //int ni=2,nj=3,nk=5; //nb of hexa9
381   vector<double> field;
382   for (int k=0; k<_nk; k++)
383     for (int j=0; j<_nj; j++)
384       for (int i=0; i<_ni; i++)
385         {
386           field.push_back(i+.1);
387           field.push_back(j+.2);
388           field.push_back(k+.3);
389         }
390
391   MEDCouplingUMesh *mesh=MEDLoader::ReadUMeshFromFile(myfileName.c_str(),_mesh_name.c_str(),0);
392   int nbOfCells=mesh->getNumberOfCells();
393   MEDCouplingFieldDouble *f1=MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME);
394   f1->setName("VectorFieldOnCells");
395   f1->setDescription("DescriptionOfFieldOnCells"); //not saved in file?
396   f1->setMesh(mesh);
397   DataArrayDouble *myField=DataArrayDouble::New();
398   myField->alloc(nbOfCells,3);
399   std::copy(field.begin(),field.end(),myField->getPointer());
400   f1->setArray(myField);
401   myField->setInfoOnComponent(0,"vx");
402   myField->setInfoOnComponent(1,"vy");
403   myField->setInfoOnComponent(2,"vz");
404   myField->decrRef();
405   f1->setTime(2.,0,1);
406   f1->checkCoherency();
407   mesh->decrRef();
408   return f1;
409 }
410
411 MEDCouplingFieldDouble * MEDPARTITIONERTest::buildVecFieldOnNodes()
412 {
413   //int ni=2,nj=3,nk=5; //nb of hexa9
414   vector<double> field;
415   for (int k=0; k<=_nk; k++)
416     for (int j=0; j<=_nj; j++)
417       for (int i=0; i<=_ni; i++)
418         {
419           field.push_back(i+.1);
420           field.push_back(j+.2);
421           field.push_back(k+.3);
422         }
423
424   MEDCouplingUMesh *mesh=MEDLoader::ReadUMeshFromFile(_file_name.c_str(),_mesh_name.c_str(),0);
425   int nbOfNodes=mesh->getNumberOfNodes();
426   MEDCouplingFieldDouble *f1=MEDCouplingFieldDouble::New(ON_NODES,ONE_TIME);
427   f1->setName("VectorFieldOnNodes");
428   f1->setDescription("DescriptionOfFieldOnNodes"); //not saved in file?
429   f1->setMesh(mesh);
430   DataArrayDouble *myField=DataArrayDouble::New();
431   myField->alloc(nbOfNodes,3);
432   std::copy(field.begin(),field.end(),myField->getPointer());
433   f1->setArray(myField);
434   myField->setInfoOnComponent(0,"vx");
435   myField->setInfoOnComponent(1,"vy");
436   myField->setInfoOnComponent(2,"vz");
437   myField->decrRef();
438   f1->setTime(2.,0,1);
439   f1->checkCoherency();
440   mesh->decrRef();
441   return f1;
442 }
443
444
445 void MEDPARTITIONERTest::createTestMeshWithoutField()
446 {
447   {
448     MEDCouplingUMesh * mesh = buildCUBE3DMesh();
449     MEDLoader::WriteUMesh(_file_name.c_str(),mesh,true);
450     if (_verbose) cout<<endl<<_file_name<<" created"<<endl;
451     if (_ntot<1000000) //too long
452       {
453         MEDCouplingUMesh *mesh_rw=MEDLoader::ReadUMeshFromFile(_file_name.c_str(),mesh->getName().c_str(),0);
454         if (_verbose) cout<<_file_name<<" reread"<<endl;
455         CPPUNIT_ASSERT(mesh->isEqual(mesh_rw,1e-12));
456         mesh_rw->decrRef();
457       }
458     mesh->decrRef();
459   }
460
461   {
462     vector<const ParaMEDMEM::MEDCouplingUMesh*> meshes;
463     MEDCouplingUMesh * mesh1 = buildCUBE3DMesh();
464     MEDCouplingUMesh * mesh2 = buildFACE3DMesh();
465     mesh1->setName("testMesh");
466     mesh2->setName("theFaces");
467     mesh2->tryToShareSameCoordsPermute(*mesh1, 1e-9);
468     mesh2->checkCoherency();
469     mesh1->checkCoherency();
470     meshes.push_back(mesh1);
471     meshes.push_back(mesh2);
472     MEDLoader::WriteUMeshes(_file_name_with_faces.c_str(), meshes, true);
473
474     ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(_file_name_with_faces.c_str(), mesh1->getName().c_str());
475     DataArrayInt* FacesFam=DataArrayInt::New();
476     FacesFam->alloc(mfm->getSizeAtLevel(-1),1);
477     FacesFam->fillWithValue(-1);
478     DataArrayInt* CellsFam=DataArrayInt::New();
479     CellsFam->alloc(mfm->getSizeAtLevel(0),1);
480     CellsFam->fillWithValue(1);
481     mfm->setFamilyFieldArr(-1,FacesFam);
482     mfm->setFamilyFieldArr(0,CellsFam);
483     map<string,int> theFamilies;
484     theFamilies["FAMILLE_ZERO"]=0;
485     theFamilies["FamilyFaces"]=-1;
486     theFamilies["FamilyCells"]=1;
487     map<string, vector<string> > theGroups;
488     theGroups["GroupFaces"].push_back("FamilyFaces");
489     theGroups["GroupCells"].push_back("FamilyCells");
490     mfm->setFamilyInfo(theFamilies);
491     mfm->setGroupInfo(theGroups);
492     mfm->write(_file_name_with_faces.c_str(),0);
493     FacesFam->decrRef();
494     CellsFam->decrRef();
495
496     /*ce truc marche pas!
497       ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(_file_name_with_faces.c_str(), mesh1->getName());
498       vector<const ParaMEDMEM::MEDCouplingUMesh*> ms;
499       ms.push_back(mesh2);
500       mfm->setGroupsFromScratch(-1, ms);
501       mfm->write(_file_name_with_faces.c_str(),0);
502     */
503
504     if (_verbose) cout<<endl<<_file_name_with_faces<<" created"<<endl;
505     if (_ntot<1000000) //too long
506       {
507         MEDCouplingUMesh *mesh_rw=MEDLoader::ReadUMeshFromFile(_file_name_with_faces.c_str(),mesh1->getName().c_str(),0);
508         if (_verbose) cout<<_file_name_with_faces<<" reread"<<endl;
509         CPPUNIT_ASSERT(mesh1->isEqual(mesh_rw,1e-12));
510         mesh_rw->decrRef();
511       }
512     mesh1->decrRef();
513     mesh2->decrRef();
514     mfm->decrRef();
515   }
516
517   {
518     MEDCouplingUMesh * mesh = buildCARRE3DMesh();
519     MEDLoader::WriteUMesh(_file_name2.c_str(),mesh,true);
520     if (_verbose) cout<<endl<<_file_name2<<" created"<<endl;
521     MEDCouplingUMesh *mesh_rw=MEDLoader::ReadUMeshFromFile(_file_name2.c_str(),mesh->getName().c_str(),0);
522     if (_verbose) cout<<_file_name2<<" reread"<<endl;
523     CPPUNIT_ASSERT(mesh->isEqual(mesh_rw,1e-12));
524     mesh_rw->decrRef();
525     mesh->decrRef();
526   }
527 }
528
529 /*
530 create a set of nbx*nby*nbz files mesh of ni*ny*nz cells
531 */
532 void MEDPARTITIONERTest::createHugeTestMesh(int ni, int nj, int nk, int nbx, int nby, int nbz, int nbTarget)
533 {
534   setSize(ni,nj,nk);
535   _nb_target_huge=nbTarget;
536   MEDCouplingUMesh * mesh = buildCUBE3DMesh();
537   //int nbx=1, nby=1, nbz=2;
538   std::vector< double > cooDep,cooFin;
539   mesh->getCoordinatesOfNode(0, cooDep);
540   mesh->getCoordinatesOfNode(mesh->getNumberOfNodes()-1, cooFin);
541   //cout<<endl<<cooDep[0]<<" "<<cooDep[1]<<" "<<cooDep[2]<<endl;
542   //cout<<cooFin[0]<<" "<<cooFin[1]<<" "<<cooFin[2]<<endl;
543
544   string tagXml="\
545 <?xml version=\"1.0\" encoding=\"UTF-8\"?>\n \
546 <root>\n \
547   <version maj=\"2\" min=\"3\" ver=\"1\"/>\n \
548   <description what=\"\" when=\"YYMMDDHHmm\"/>\n \
549   <content>\n \
550     <mesh name=\"testMesh\"/>\n \
551   </content>\n \
552   <splitting>\n \
553     <subdomain number=\"$subdomainNumber\"/>\n \
554     <global_numbering present=\"no\"/>\n \
555   </splitting>\n \
556   <files>\n$tagSubfile \
557   </files>\n \
558   <mapping>\n$tagMesh \
559   </mapping>\n \
560 </root>\n";
561
562   string tagSubfiles, tagSubfile="\
563     <subfile id=\"$xyz\">\n \
564       <name>$fileName</name>\n \
565       <machine>localhost</machine>\n \
566     </subfile>\n";
567   string tagMeshes, tagMesh="\
568     <mesh name=\"testMesh\">\n \
569       <chunk subdomain=\"$xyz\">\n \
570         <name>testMesh</name>\n \
571       </chunk>\n \
572     </mesh>\n";
573
574   int xyz=1;
575   string sxyz;
576   DataArrayDouble* coordsInit=mesh->getCoords()->deepCpy();
577   double* ptrInit=coordsInit->getPointer();
578   double deltax=cooFin[0]-cooDep[0];
579   double deltay=cooFin[1]-cooDep[1];
580   double deltaz=cooFin[2]-cooDep[2];
581
582   double dz=0.;
583   for (int z=0; z<nbz; z++)
584     {
585       double dy=0.;
586       for (int y=0; y<nby; y++)
587         {
588           double dx=0.;
589           for (int x=0; x<nbx; x++)
590             {
591               string fileName;
592               sxyz=IntToStr(xyz);
593               fileName="tmp_testMeshHuge_"+IntToStr(_ni)+"x"+IntToStr(_nj)+"x"+IntToStr(_nk)+"_"+sxyz+".med";
594
595               DataArrayDouble* coords=mesh->getCoords();
596               //int nbOfComp=coords->getNumberOfComponents();  //be 3D
597               int nbOfTuple=coords->getNumberOfTuples();
598               double* ptr=coords->getPointer();
599               double* ptrini=ptrInit;
600               for (int i=0; i<nbOfTuple; i++)
601                 {
602                   *ptr=(*ptrini)+dx; ptr++; ptrini++; //be 3D
603                   *ptr=(*ptrini)+dy; ptr++; ptrini++;
604                   *ptr=(*ptrini)+dz; ptr++; ptrini++;
605                 }
606
607               MEDLoader::WriteUMesh(fileName.c_str(),mesh,true);
608
609               tagSubfiles+=tagSubfile;
610               tagSubfiles.replace(tagSubfiles.find("$xyz"),4,sxyz);
611               tagSubfiles.replace(tagSubfiles.find("$fileName"),9,fileName);
612
613               tagMeshes+=tagMesh;
614               tagMeshes.replace(tagMeshes.find("$xyz"),4,sxyz);
615               xyz++;
616               dx+=deltax;
617             }
618           dy+=deltay;
619         }
620       dz+=deltaz;
621     }
622   coordsInit->decrRef();
623
624   tagXml.replace(tagXml.find("$subdomainNumber"),16,sxyz);
625   tagXml.replace(tagXml.find("$tagSubfile"),11,tagSubfiles);
626   tagXml.replace(tagXml.find("$tagMesh"),8,tagMeshes);
627
628   string nameFileXml;
629   _file_name_huge_xml="tmp_testMeshHuge_"+IntToStr(_ni)+"x"+IntToStr(_nj)+"x"+IntToStr(_nk)+"_"+sxyz+".xml";
630   std::ofstream f(_file_name_huge_xml.c_str());
631   f<<tagXml;
632   f.close();
633   //cout<<"\n"<<tagXml<<endl;
634   if (_verbose)
635     cout<<endl<<nameFileXml<<" created"<<endl;
636   mesh->decrRef();
637 }
638
639 void MEDPARTITIONERTest::createTestMeshWithVecFieldOnCells()
640 {
641   {
642     string name=_file_name;
643     MEDCouplingFieldDouble *f1=buildVecFieldOnCells(name);
644     name.replace(name.find(".med"),4,"_WithVecFieldOnCells.med");
645     MEDLoader::WriteField(name.c_str(),f1,true);
646     f1->setTime(3.,1,1);  //time,it,order
647     f1->applyFunc("x/2.");
648     MEDLoader::WriteField(name.c_str(),f1,false);
649     if (_verbose) cout<<endl<<name<<" created"<<endl;
650     if (_ntot<1000000) //too long
651       {
652         MEDCouplingFieldDouble *f2=MEDLoader::ReadFieldCell(name.c_str(),f1->getMesh()->getName().c_str(),0,f1->getName().c_str(),0,1);
653         //DataArrayDouble *res=f2->getArray();
654         if (_verbose) cout<<name<<" reread"<<endl;
655         //CPPUNIT_ASSERT(f1->isEqual(f2,1e-12,1e-12));
656         f2->decrRef();
657       }
658     f1->decrRef();
659   }
660   {
661     string name=_file_name;
662     MEDCouplingFieldDouble *f1=buildVecFieldOnCells(name);
663     name.replace(name.find(".med"),4,"_WithVecFieldOnGaussNe.med");
664     MEDCouplingFieldDouble *f3=MEDCouplingFieldDouble::New(ON_GAUSS_NE,ONE_TIME);
665     f3->setMesh(f1->getMesh());
666     //cout<<"\nNumberOfMeshPlacesExpected "<<f3->getNumberOfMeshPlacesExpected()<<" "
667     //                                     /*<<getNumberOfTuples(f1->getMesh())<<" "*/
668     //                                     <<f3->getMesh()->getNumberOfNodes()<<" "
669     //                                     <<f3->getMesh()->getNumberOfCells()<<endl;
670     f3->setName("MyFieldOnGaussNE");
671     f3->setDescription("MyDescriptionNE");
672     DataArrayDouble *array=DataArrayDouble::New();
673     //int nb=f1->getMesh()->getNumberOfNodes();
674
675     /*8 pt de gauss by cell
676       int nb=f3->getMesh()->getNumberOfCells()*8;
677       array->alloc(nb,2);
678       double *ptr=array->getPointer();
679       for (int i=0; i<nb*2; i=i+2) {ptr[i]=(double)(i/8) ; ptr[i]=2.*(double)(i/8);}
680     */
681
682     //more nbptgauss=8 by default needs set MEDCouplingFieldDiscretizationPerCell
683     //theory: (may be) http://www.code-aster.org/V2/doc/v9/fr/man_r/r3/r3.06.03.pdf
684     int nbptgauss=8; //nb pt de gauss by cell
685     int nbcell=f3->getMesh()->getNumberOfCells();
686     int nb=nbcell*nbptgauss;
687     int nbcomp=2;
688     array->alloc(nb,nbcomp);
689     double *ptr=array->getPointer();
690     int ii=0;
691     for (int i=0; i<nbcell; i++)
692       for (int j=0; j<nbptgauss; j++)
693         for (int k=0; k<nbcomp; k++)
694           {
695             //123.4 for 12th cell,3rd component, 4th gausspoint
696             ptr[ii]=(double)((i+1)*10+(k+1))+((double)(j+1))/10.;
697             ii++;
698           }
699     array->setInfoOnComponent(0,"vGx");
700     array->setInfoOnComponent(1,"vGy");
701     f3->setTime(4.,5,6);
702     f3->setArray(array);
703     array->decrRef();
704     MEDLoader::WriteField(name.c_str(),f3,true);
705     if (_verbose) cout<<endl<<name<<" created"<<endl;
706     f3->checkCoherency();
707     f1->decrRef();
708     if (_ntot<1000000) //too long
709       {
710         MEDCouplingFieldDouble* f4=MEDLoader::ReadField(ON_GAUSS_NE,
711                                                         name.c_str(), f3->getMesh()->getName().c_str(), 0, "MyFieldOnGaussNE", 5, 6);
712         if (_verbose) cout<<"MyFieldOnGaussNE reread"<<endl;
713         f4->decrRef();
714       }
715     f3->decrRef();
716   }
717   {
718     string name=_file_name_with_faces;
719     MEDCouplingFieldDouble *f1=buildVecFieldOnCells(name);
720     name.replace(name.find(".med"),4,"_WithVecFieldOnCells.med");
721     MEDLoader::WriteField(name.c_str(),f1,true);
722     if (_verbose) cout<<endl<<name<<" created"<<endl;
723     if (_ntot<1000000) //too long
724       {
725         MEDCouplingFieldDouble *f2=MEDLoader::ReadFieldCell(name.c_str(),f1->getMesh()->getName().c_str(),0,f1->getName().c_str(),0,1);
726         if (_verbose) cout<<name<<" reread"<<endl;
727         //CPPUNIT_ASSERT(f1->isEqual(f2,1e-12,1e-12)); assertion failed!!
728         f2->decrRef();
729       }
730     f1->decrRef();
731   }
732 }
733
734 void MEDPARTITIONERTest::createTestMeshWithVecFieldOnNodes()
735 {
736   MEDCouplingFieldDouble *f1=buildVecFieldOnNodes();
737   string name=_file_name;
738   name.replace(name.find(".med"),4,"_WithVecFieldOnNodes.med");
739   MEDLoader::WriteField(name.c_str(),f1,true);
740   if (_verbose) cout<<endl<<name<<" created"<<endl;
741   if (_ntot<1000000) //too long
742     {
743       MEDCouplingFieldDouble *f2=MEDLoader::ReadFieldNode(name.c_str(),f1->getMesh()->getName().c_str(),0,f1->getName().c_str(),0,1);
744       if (_verbose) cout<<name<<" reread"<<endl;
745       //CPPUNIT_ASSERT(f1->isEqual(f2,1e-12,1e-12)); assertion failed!!
746       f2->decrRef();
747     }
748   f1->decrRef();
749 }
750
751 void MEDPARTITIONERTest::verifyTestMeshWithVecFieldOnNodes()
752 {
753   string name=_file_name;
754   name.replace(name.find(".med"),4,"_WithVecFieldOnNodes.med");
755   MEDCouplingUMesh * m=MEDLoader::ReadUMeshFromFile(name.c_str(),_mesh_name.c_str(),0);
756   std::set<INTERP_KERNEL::NormalizedCellType> types(m->getAllGeoTypes());
757   if (_verbose)
758     {
759       cout<<"\n types in "<<name<<" : ";
760       //for (std::set<INTERP_KERNEL::NormalizedCellType>::iterator t=types.begin(); t!=types.end(); ++t) cout<<" "<<*t;
761       for (std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator t=types.begin(); t!=types.end(); ++t)
762         {
763           //INTERP_KERNEL::CellModel essai=INTERP_KERNEL::CellModel::GetCellModel(*t);
764           cout<<" "<<(INTERP_KERNEL::CellModel::GetCellModel(*t)).getRepr();
765         }
766       cout<<endl;
767     }
768   m->decrRef();
769
770   MEDFileUMesh * mf = MEDFileUMesh::New(_file_name.c_str(),_mesh_name.c_str(),-1,-1);
771   vector<int> lev;
772   lev=mf->getNonEmptyLevels();
773   if (_verbose)
774     {
775       cout<<" levels in "<<name<<" : ";
776       for (vector<int>::iterator l=lev.begin(); l!=lev.end(); ++l) cout<<" "<<*l;
777       cout<<endl;
778     }
779   mf->decrRef();
780 }
781
782 void MEDPARTITIONERTest::createTestMeshes()
783 {
784   createTestMeshWithoutField();
785   createTestMeshWithVecFieldOnCells();
786   createTestMeshWithVecFieldOnNodes();
787 }
788
789 void MEDPARTITIONERTest::deleteTestMeshes()
790 {
791   string cmd="rm *tmp_testMesh*";
792   if (_verbose) cout<<endl<<cmd<<endl;
793   system(cmd.c_str());  //may be not if debug
794 }
795
796 void MEDPARTITIONERTest::testMeshCollectionSingle()
797 {
798   setSmallSize();
799   createTestMeshes();
800   MyGlobals::_World_Size=1;
801   MyGlobals::_Rank=0;
802   string fileName=_file_name_with_faces;
803   MEDPARTITIONER::ParaDomainSelector parallelizer(false);
804   MEDPARTITIONER::MeshCollection collection(fileName,parallelizer);
805   CPPUNIT_ASSERT(collection.isParallelMode());
806   CPPUNIT_ASSERT_EQUAL(3, collection.getMeshDimension());
807   CPPUNIT_ASSERT(collection.getName()=="testMesh");
808   CPPUNIT_ASSERT_EQUAL(1,collection.getNbOfLocalMeshes());
809   CPPUNIT_ASSERT_EQUAL(1,collection.getNbOfGlobalMeshes());
810   CPPUNIT_ASSERT_EQUAL(_ni*_nj*_nk,collection.getNbOfLocalCells());
811   CPPUNIT_ASSERT_EQUAL(_ni*_nj,collection.getNbOfLocalFaces());
812 }
813
814 void MEDPARTITIONERTest::testMeshCollectionXml()
815 {
816   setSmallSize();
817   createHugeTestMesh(_ni, _nj, _nk, 2, 2, 2, 32); //xml but not so huge
818   string fileName=_file_name_huge_xml;
819   MEDPARTITIONER::ParaDomainSelector parallelizer(false);
820   MEDPARTITIONER::MeshCollection collection(fileName,parallelizer);
821   CPPUNIT_ASSERT(collection.isParallelMode());
822   CPPUNIT_ASSERT_EQUAL(3, collection.getMeshDimension());
823   CPPUNIT_ASSERT(collection.getName()=="testMesh");
824   CPPUNIT_ASSERT_EQUAL(8,collection.getNbOfLocalMeshes());
825   CPPUNIT_ASSERT_EQUAL(8,collection.getNbOfGlobalMeshes());
826   CPPUNIT_ASSERT_EQUAL(_ni*_nj*_nk*8,collection.getNbOfLocalCells());
827   CPPUNIT_ASSERT_EQUAL(0,collection.getNbOfLocalFaces());
828 }
829
830
831 //#################for metis
832
833
834
835 #if defined(MED_ENABLE_METIS)
836 void MEDPARTITIONERTest::testMeshCollectionSinglePartitionMetis()
837 {
838   setSmallSize();
839   createTestMeshes();
840   //MyGlobals::_Verbose=500;
841   string fileName=_file_name_with_faces;
842   int ndomains=2;
843   bool split_family=false;
844   bool empty_groups=false;
845   MEDPARTITIONER::ParaDomainSelector parallelizer(false);
846   MEDPARTITIONER::MeshCollection collection(fileName,parallelizer);
847
848   MEDPARTITIONER::ParallelTopology* aPT = (MEDPARTITIONER::ParallelTopology*) collection.getTopology();
849   aPT->setGlobalNumerotationDefault(collection.getParaDomainSelector());
850   //Creating the graph and partitioning it
851   auto_ptr< MEDPARTITIONER::Topology > new_topo;
852   new_topo.reset( collection.createPartition(ndomains,MEDPARTITIONER::Graph::METIS) );
853   //Creating a new mesh collection from the partitioning
854   MEDPARTITIONER::MeshCollection new_collection(collection,new_topo.get(),split_family,empty_groups);
855
856   //example to create files
857   //MyGlobals::_General_Informations.clear();
858   //MyGlobals::_General_Informations.push_back(SerializeFromString("finalMeshName=Merge"));
859   //if (MyGlobals::_Verbose>100) cout << "generalInformations : \n"<<ReprVectorOfString(MyGlobals::_General_Informations);
860   //new_collection.write("ttmp")
861
862   CPPUNIT_ASSERT(new_collection.isParallelMode());
863   CPPUNIT_ASSERT_EQUAL(3, new_collection.getMeshDimension());
864   CPPUNIT_ASSERT(new_collection.getName()==collection.getName());
865   CPPUNIT_ASSERT_EQUAL(ndomains,new_collection.getNbOfLocalMeshes());
866   CPPUNIT_ASSERT_EQUAL(ndomains,new_collection.getNbOfGlobalMeshes());
867   CPPUNIT_ASSERT_EQUAL(collection.getNbOfLocalCells(),new_collection.getNbOfLocalCells());
868   CPPUNIT_ASSERT_EQUAL(collection.getNbOfLocalFaces(),new_collection.getNbOfLocalFaces());
869 }
870
871 void MEDPARTITIONERTest::testMeshCollectionComplexPartitionMetis()
872 {
873   setSmallSize();
874   createHugeTestMesh(_ni, _nj, _nk, 2, 2, 2, 32); //xml on 2*2*2 meshes but not so huge
875   string fileName=_file_name_huge_xml;
876   bool split_family=false;
877   bool empty_groups=false;
878   MEDPARTITIONER::ParaDomainSelector parallelizer(false);
879   MEDPARTITIONER::MeshCollection collection(fileName,parallelizer);
880
881   MEDPARTITIONER::ParallelTopology* aPT = (MEDPARTITIONER::ParallelTopology*) collection.getTopology();
882   aPT->setGlobalNumerotationDefault(collection.getParaDomainSelector());
883
884   for (int ndomains=2 ; ndomains<=16 ; ndomains++)
885     {
886       //Creating the graph and partitioning it
887       auto_ptr< MEDPARTITIONER::Topology > new_topo;
888       new_topo.reset( collection.createPartition(ndomains,MEDPARTITIONER::Graph::METIS) );
889       //Creating a new mesh collection from the partitioning
890       MEDPARTITIONER::MeshCollection new_collection(collection,new_topo.get(),split_family,empty_groups);
891
892       CPPUNIT_ASSERT_EQUAL(ndomains,new_collection.getNbOfLocalMeshes());
893       CPPUNIT_ASSERT_EQUAL(ndomains,new_collection.getNbOfGlobalMeshes());
894       CPPUNIT_ASSERT_EQUAL(collection.getNbOfLocalCells(),new_collection.getNbOfLocalCells());
895       CPPUNIT_ASSERT_EQUAL(collection.getNbOfLocalFaces(),new_collection.getNbOfLocalFaces());
896     }
897 }
898
899 void MEDPARTITIONERTest::testMetisSmallSize()
900 {
901   //#if !defined(HAVE_MPI)
902   setSmallSize();
903   createTestMeshes();
904   std::string MetisOrScotch("metis");
905   launchMetisOrScotchMedpartitionerOnTestMeshes(MetisOrScotch);
906   verifyMetisOrScotchMedpartitionerOnSmallSizeForMesh(MetisOrScotch);
907   verifyMetisOrScotchMedpartitionerOnSmallSizeForFieldOnCells(MetisOrScotch);
908   verifyMetisOrScotchMedpartitionerOnSmallSizeForFieldOnGaussNe(MetisOrScotch);
909   //#endif
910 }
911 #endif
912
913
914 //#################for scotch
915
916
917 #if defined(MED_ENABLE_SCOTCH)
918 void MEDPARTITIONERTest::testMeshCollectionSinglePartitionScotch()
919 {
920   setSmallSize();
921   createTestMeshes();
922   //MyGlobals::_Verbose=500;
923   string fileName=_file_name_with_faces;
924   int ndomains=2;
925   bool split_family=false;
926   bool empty_groups=false;
927   MEDPARTITIONER::ParaDomainSelector parallelizer(false);
928   MEDPARTITIONER::MeshCollection collection(fileName,parallelizer);
929
930   MEDPARTITIONER::ParallelTopology* aPT = (MEDPARTITIONER::ParallelTopology*) collection.getTopology();
931   aPT->setGlobalNumerotationDefault(collection.getParaDomainSelector());
932   //Creating the graph and partitioning it
933   auto_ptr< MEDPARTITIONER::Topology > new_topo;
934   new_topo.reset( collection.createPartition(ndomains,MEDPARTITIONER::Graph::SCOTCH) );
935   //Creating a new mesh collection from the partitioning
936   MEDPARTITIONER::MeshCollection new_collection(collection,new_topo.get(),split_family,empty_groups);
937
938   //example to create files
939   //MyGlobals::_General_Informations.clear();
940   //MyGlobals::_General_Informations.push_back(SerializeFromString("finalMeshName=Merge"));
941   //if (MyGlobals::_Verbose>100) cout << "generalInformations : \n"<<ReprVectorOfString(MyGlobals::_General_Informations);
942   //new_collection.write("ttmp")
943
944   CPPUNIT_ASSERT(new_collection.isParallelMode());
945   CPPUNIT_ASSERT_EQUAL(3, new_collection.getMeshDimension());
946   CPPUNIT_ASSERT(new_collection.getName()==collection.getName());
947   CPPUNIT_ASSERT_EQUAL(ndomains,new_collection.getNbOfLocalMeshes());
948   CPPUNIT_ASSERT_EQUAL(ndomains,new_collection.getNbOfGlobalMeshes());
949   CPPUNIT_ASSERT_EQUAL(collection.getNbOfLocalCells(),new_collection.getNbOfLocalCells());
950   CPPUNIT_ASSERT_EQUAL(collection.getNbOfLocalFaces(),new_collection.getNbOfLocalFaces());
951 }
952
953 void MEDPARTITIONERTest::testMeshCollectionComplexPartitionScotch()
954 {
955   setSmallSize();
956   createHugeTestMesh(_ni, _nj, _nk, 2, 2, 2, 32); //xml on 2*2*2 meshes but not so huge
957   string fileName=_file_name_huge_xml;
958   bool split_family=false;
959   bool empty_groups=false;
960   MEDPARTITIONER::ParaDomainSelector parallelizer(false);
961   MEDPARTITIONER::MeshCollection collection(fileName,parallelizer);
962
963   MEDPARTITIONER::ParallelTopology* aPT = (MEDPARTITIONER::ParallelTopology*) collection.getTopology();
964   aPT->setGlobalNumerotationDefault(collection.getParaDomainSelector());
965
966   for (int ndomains=2 ; ndomains<=16 ; ndomains++)
967     {
968       //Creating the graph and partitioning it
969       auto_ptr< MEDPARTITIONER::Topology > new_topo;
970       new_topo.reset( collection.createPartition(ndomains,MEDPARTITIONER::Graph::SCOTCH) );
971       //Creating a new mesh collection from the partitioning
972       MEDPARTITIONER::MeshCollection new_collection(collection,new_topo.get(),split_family,empty_groups);
973
974       CPPUNIT_ASSERT_EQUAL(ndomains,new_collection.getNbOfLocalMeshes());
975       CPPUNIT_ASSERT_EQUAL(ndomains,new_collection.getNbOfGlobalMeshes());
976       CPPUNIT_ASSERT_EQUAL(collection.getNbOfLocalCells(),new_collection.getNbOfLocalCells());
977       CPPUNIT_ASSERT_EQUAL(collection.getNbOfLocalFaces(),new_collection.getNbOfLocalFaces());
978     }
979 }
980
981 void MEDPARTITIONERTest::testScotchSmallSize()
982 {
983   //#if !defined(HAVE_MPI)
984   setSmallSize();
985   createTestMeshes();
986   std::string MetisOrScotch("scotch");
987   launchMetisOrScotchMedpartitionerOnTestMeshes(MetisOrScotch);
988   verifyMetisOrScotchMedpartitionerOnSmallSizeForMesh(MetisOrScotch);
989   verifyMetisOrScotchMedpartitionerOnSmallSizeForFieldOnCells(MetisOrScotch);
990   verifyMetisOrScotchMedpartitionerOnSmallSizeForFieldOnGaussNe(MetisOrScotch);
991   //#endif
992 }
993 #endif
994
995 void MEDPARTITIONERTest::launchMetisOrScotchMedpartitionerOnTestMeshes(std::string MetisOrScotch)
996 {
997   int res;
998   string cmd,execName,sourceName,targetName;
999
1000   execName=getPartitionerExe();
1001
1002   cmd="which "+execName+" 2>/dev/null 1>/dev/null";  //no trace
1003   res=system(cmd.c_str());
1004   CPPUNIT_ASSERT_EQUAL_MESSAGE(execName + " - INVALID PATH TO medpartitioner", 0, res);
1005
1006   cmd=execName+" --ndomains=2 --split-method="+MetisOrScotch;  //on same proc
1007   sourceName=_file_name;
1008   targetName=_file_name;
1009   targetName.replace(targetName.find(".med"),4,"_partitionedTo2_");
1010   cmd+=" --input-file="+sourceName+" --output-file="+targetName+" --verbose="+IntToStr(_verbose);
1011   if (_verbose) cout<<endl<<cmd<<endl;
1012   res=system(cmd.c_str());
1013   CPPUNIT_ASSERT_EQUAL(0, res);
1014
1015   cmd=execName+" --ndomains=5 --split-method="+MetisOrScotch; //on less proc
1016   sourceName=_file_name;
1017   targetName=_file_name;
1018   targetName.replace(targetName.find(".med"),4,"_partitionedTo5_");
1019   cmd+=" --input-file="+sourceName+" --output-file="+targetName+" --verbose="+IntToStr(_verbose);
1020   if (_verbose) cout<<endl<<cmd<<endl;
1021   res=system(cmd.c_str());
1022   CPPUNIT_ASSERT_EQUAL(0, res);
1023
1024   cmd=execName+" --ndomains=1 --split-method="+MetisOrScotch;  //on 1 proc
1025   sourceName=targetName+".xml";
1026   targetName=_file_name;
1027   targetName.replace(targetName.find(".med"),4,"_remergedFrom5_");
1028   cmd+=" --input-file="+sourceName+" --output-file="+targetName+" --verbose="+IntToStr(_verbose);
1029   if (_verbose) cout<<endl<<cmd<<endl;
1030   res=system(cmd.c_str());
1031   CPPUNIT_ASSERT_EQUAL(0, res);
1032
1033   cmd=execName+" --ndomains=1 --split-method="+MetisOrScotch;  //on more proc
1034   //sourceName=targetName+".xml";
1035   targetName=_file_name;
1036   targetName.replace(targetName.find(".med"),4,"_remergedFrom5_");
1037   cmd+=" --input-file="+sourceName+" --output-file="+targetName+" --verbose="+IntToStr(_verbose);
1038   if (_verbose) cout<<endl<<cmd<<endl;
1039   res=system(cmd.c_str());
1040   CPPUNIT_ASSERT_EQUAL(0, res);
1041 }
1042
1043 void MEDPARTITIONERTest::verifyMetisOrScotchMedpartitionerOnSmallSizeForMesh(std::string MetisOrScotch)
1044 {
1045   int res;
1046   string fileName,cmd,execName,sourceName,targetName,input;
1047   execName=getPartitionerExe();
1048   fileName=_file_name_with_faces;
1049
1050   ParaMEDMEM::MEDFileUMesh* initialMesh=ParaMEDMEM::MEDFileUMesh::New(fileName.c_str(),_mesh_name.c_str());
1051   ParaMEDMEM::MEDCouplingUMesh* cellMesh=initialMesh->getLevel0Mesh(false);
1052   ParaMEDMEM::MEDCouplingUMesh* faceMesh=initialMesh->getLevelM1Mesh(false);
1053
1054   cmd=execName+" --ndomains=5 --split-method="+MetisOrScotch;  //on same proc
1055   sourceName=fileName;
1056   targetName=fileName;
1057   targetName.replace(targetName.find(".med"),4,"_partitionedTo5_");
1058   cmd+=" --input-file="+sourceName+" --output-file="+targetName+" --verbose="+IntToStr(_verbose);
1059   if (_verbose) cout<<endl<<cmd<<endl;
1060   res=system(cmd.c_str());
1061   CPPUNIT_ASSERT_EQUAL(0, res);
1062   input=targetName+".xml";
1063
1064   MEDPARTITIONER::ParaDomainSelector parallelizer(false);
1065   MEDPARTITIONER::MeshCollection collection(input,parallelizer);
1066   CPPUNIT_ASSERT_EQUAL(3, collection.getMeshDimension());
1067   std::vector<ParaMEDMEM::MEDCouplingUMesh*>cellMeshes=collection.getMesh();
1068   CPPUNIT_ASSERT_EQUAL(5, (int) cellMeshes.size());
1069   int nbcells=0;
1070   for (std::size_t i = 0; i < cellMeshes.size(); i++)
1071     nbcells+=cellMeshes[i]->getNumberOfCells();
1072   CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), nbcells);
1073
1074   std::vector<ParaMEDMEM::MEDCouplingUMesh*>faceMeshes=collection.getFaceMesh();
1075   CPPUNIT_ASSERT_EQUAL(5, (int) faceMeshes.size());
1076   int nbfaces=0;
1077   for (std::size_t i=0; i < faceMeshes.size(); i++)
1078     nbfaces+=faceMeshes[i]->getNumberOfCells();
1079   CPPUNIT_ASSERT_EQUAL(faceMesh->getNumberOfCells(), nbfaces);
1080
1081   //merge split meshes and test equality
1082   cmd=execName+" --ndomains=1 --split-method="+MetisOrScotch;  //on same proc
1083   sourceName=targetName+".xml";
1084   targetName=fileName;
1085   targetName.replace(targetName.find(".med"),4,"_remergedFrom5_");
1086   cmd+=" --input-file="+sourceName+" --output-file="+targetName+" --verbose="+IntToStr(_verbose);
1087   if (_verbose) cout<<endl<<cmd<<endl;
1088   res=system(cmd.c_str());
1089   CPPUNIT_ASSERT_EQUAL(0, res);
1090
1091   string refusedName=targetName+"1.med";
1092   ParaMEDMEM::MEDFileUMesh* refusedMesh=ParaMEDMEM::MEDFileUMesh::New(refusedName.c_str(),_mesh_name.c_str());
1093   ParaMEDMEM::MEDCouplingUMesh* refusedCellMesh=refusedMesh->getLevel0Mesh(false);
1094   ParaMEDMEM::MEDCouplingUMesh* refusedFaceMesh=refusedMesh->getLevelM1Mesh(false);
1095
1096   CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), refusedCellMesh->getNumberOfCells());
1097   CPPUNIT_ASSERT_EQUAL(faceMesh->getNumberOfCells(), refusedFaceMesh->getNumberOfCells());
1098
1099   /*not the good job
1100     ParaMEDMEM::MEDCouplingMesh* mergeCell=cellMesh->mergeMyselfWith(refusedCellMesh);
1101     CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), mergeCell->getNumberOfCells());
1102
1103     ParaMEDMEM::MEDCouplingMesh* mergeFace=faceMesh->mergeMyselfWith(refusedFaceMesh);
1104     CPPUNIT_ASSERT_EQUAL(faceMesh->getNumberOfCells(), mergeFace->getNumberOfCells());
1105
1106     CPPUNIT_ASSERT(faceMesh->isEqual(refusedFaceMesh,1e-12));
1107   */
1108
1109   std::vector<const MEDCouplingUMesh *> meshes;
1110   std::vector<DataArrayInt *> corr;
1111   meshes.push_back(cellMesh);
1112   refusedCellMesh->tryToShareSameCoordsPermute(*cellMesh, 1e-9);
1113   meshes.push_back(refusedCellMesh);
1114   MEDCouplingUMesh* fusedCell=MEDCouplingUMesh::FuseUMeshesOnSameCoords(meshes,0,corr);
1115   CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), fusedCell->getNumberOfCells());
1116
1117   meshes.resize(0);
1118   for (std::size_t i = 0; i < corr.size(); i++)
1119     corr[i]->decrRef();
1120   corr.resize(0);
1121   meshes.push_back(faceMesh);
1122   refusedFaceMesh->tryToShareSameCoordsPermute(*faceMesh, 1e-9);
1123   meshes.push_back(refusedFaceMesh);
1124   MEDCouplingUMesh* fusedFace=MEDCouplingUMesh::FuseUMeshesOnSameCoords(meshes,0,corr);
1125   CPPUNIT_ASSERT_EQUAL(faceMesh->getNumberOfCells(), fusedFace->getNumberOfCells());
1126
1127   for (std::size_t i = 0; i < corr.size(); i++)
1128     corr[i]->decrRef();
1129   fusedFace->decrRef();
1130   refusedFaceMesh->decrRef();
1131   faceMesh->decrRef();
1132   fusedCell->decrRef();
1133   refusedCellMesh->decrRef();
1134   refusedMesh->decrRef();
1135   cellMesh->decrRef();
1136   initialMesh->decrRef();
1137   //done in ~collection
1138   //for (int i = 0; i < faceMeshes.size(); i++) faceMeshes[i]->decrRef();
1139   //for (int i = 0; i < cellMeshes.size(); i++) cellMeshes[i]->decrRef();
1140 }
1141
1142 void MEDPARTITIONERTest::verifyMetisOrScotchMedpartitionerOnSmallSizeForFieldOnCells(std::string MetisOrScotch)
1143 {
1144   int res;
1145   string fileName,cmd,execName,sourceName,targetName,input;
1146   execName=getPartitionerExe();
1147   fileName=_file_name;
1148   fileName.replace(fileName.find(".med"),4,"_WithVecFieldOnCells.med");
1149
1150   ParaMEDMEM::MEDFileUMesh* initialMesh=ParaMEDMEM::MEDFileUMesh::New(fileName.c_str(),_mesh_name.c_str());
1151   ParaMEDMEM::MEDCouplingUMesh* cellMesh=initialMesh->getLevel0Mesh(false);
1152
1153   cmd=execName+" --ndomains=5 --split-method="+MetisOrScotch;  //on same proc
1154   sourceName=fileName;
1155   targetName=fileName;
1156   targetName.replace(targetName.find(".med"),4,"_partitionedTo5_");
1157   cmd+=" --input-file="+sourceName+" --output-file="+targetName+" --verbose="+IntToStr(_verbose);
1158   if (_verbose) cout<<endl<<cmd<<endl;
1159   res=system(cmd.c_str());
1160   CPPUNIT_ASSERT_EQUAL(0, res);
1161   input=targetName+".xml";
1162
1163   //merge split meshes and test equality
1164   cmd=execName+" --ndomains=1 --split-method="+MetisOrScotch;  //on same proc
1165   sourceName=targetName+".xml";
1166   targetName=fileName;
1167   targetName.replace(targetName.find(".med"),4,"_remergedFrom5_");
1168   cmd+=" --input-file="+sourceName+" --output-file="+targetName+" --verbose="+IntToStr(_verbose);
1169   if (_verbose) cout<<endl<<cmd<<endl;
1170   res=system(cmd.c_str());
1171   CPPUNIT_ASSERT_EQUAL(0, res);
1172
1173   string refusedName=targetName+"1.med";
1174   ParaMEDMEM::MEDFileUMesh* refusedMesh=ParaMEDMEM::MEDFileUMesh::New(refusedName.c_str(),_mesh_name.c_str());
1175   ParaMEDMEM::MEDCouplingUMesh* refusedCellMesh=refusedMesh->getLevel0Mesh(false);
1176
1177   CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), refusedCellMesh->getNumberOfCells());
1178
1179   std::vector<const MEDCouplingUMesh *> meshes;
1180   std::vector<DataArrayInt *> corr;
1181   meshes.push_back(cellMesh);
1182   refusedCellMesh->tryToShareSameCoordsPermute(*cellMesh, 1e-9);
1183   meshes.push_back(refusedCellMesh);
1184   MEDCouplingUMesh* fusedCell=MEDCouplingUMesh::FuseUMeshesOnSameCoords(meshes,0,corr);
1185   CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), fusedCell->getNumberOfCells());
1186
1187   MEDCouplingFieldDouble* field1=MEDLoader::ReadFieldCell(fileName.c_str(),initialMesh->getName().c_str(),0,"VectorFieldOnCells",0,1);
1188   MEDCouplingFieldDouble* field2=MEDLoader::ReadFieldCell(refusedName.c_str(),refusedCellMesh->getName().c_str(),0,"VectorFieldOnCells",0,1);
1189
1190   int nbcells=corr[1]->getNumberOfTuples();
1191   CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), nbcells);
1192   //use corr to test equality of field
1193   DataArrayDouble* f1=field1->getArray();
1194   DataArrayDouble* f2=field2->getArray();
1195   if (_verbose>300)
1196     {
1197       cout<<"\nf1 : "<<f1->reprZip();
1198       cout<<"\nf2 : "<<f2->reprZip(); //field2->advancedRepradvancedRepr();
1199       for (std::size_t i = 0; i < corr.size(); i++)
1200         cout << "\ncorr " << i << " : " << corr[i]->reprZip();
1201
1202     }
1203   int nbequal=0;
1204   int nbcomp=field1->getNumberOfComponents();
1205   double* p1=f1->getPointer();
1206   double* p2=f2->getPointer();
1207   int* pc=corr[1]->getPointer();
1208   for (int i = 0; i < nbcells; i++)
1209     {
1210       int i1=pc[i]*nbcomp;
1211       int i2=i*nbcomp;
1212       for (int j = 0; j < nbcomp; j++)
1213         {
1214           if (p1[i1+j]==p2[i2+j]) nbequal++;
1215           //cout<<" "<<p1[i1+j]<<"="<<p2[i2+j];
1216         }
1217     }
1218   CPPUNIT_ASSERT_EQUAL(nbcells*nbcomp, nbequal);
1219
1220   for (std::size_t i = 0; i < corr.size(); i++)
1221     corr[i]->decrRef();
1222   field1->decrRef();
1223   field2->decrRef();
1224   fusedCell->decrRef();
1225   refusedMesh->decrRef();
1226   refusedCellMesh->decrRef();
1227   cellMesh->decrRef();
1228   initialMesh->decrRef();
1229 }
1230
1231 void MEDPARTITIONERTest::verifyMetisOrScotchMedpartitionerOnSmallSizeForFieldOnGaussNe(std::string MetisOrScotch)
1232 {
1233   int res;
1234   string fileName,cmd,execName,sourceName,targetName,input;
1235   execName=getPartitionerExe();
1236   fileName=_file_name;
1237   fileName.replace(fileName.find(".med"),4,"_WithVecFieldOnGaussNe.med");
1238
1239   ParaMEDMEM::MEDFileUMesh* initialMesh=ParaMEDMEM::MEDFileUMesh::New(fileName.c_str(),_mesh_name.c_str());
1240   ParaMEDMEM::MEDCouplingUMesh* cellMesh=initialMesh->getLevel0Mesh(false);
1241
1242   cmd=execName+" --ndomains=5 --split-method="+MetisOrScotch;  //on same proc
1243   sourceName=fileName;
1244   targetName=fileName;
1245   targetName.replace(targetName.find(".med"),4,"_partitionedTo5_");
1246   cmd+=" --input-file="+sourceName+" --output-file="+targetName+" --verbose="+IntToStr(_verbose);
1247   if (_verbose) cout<<endl<<cmd<<endl;
1248   res=system(cmd.c_str());
1249   CPPUNIT_ASSERT_EQUAL(0, res);
1250   input=targetName+".xml";
1251
1252   //merge split meshes and test equality
1253   cmd=execName+" --ndomains=1 --split-method="+MetisOrScotch;  //on same proc
1254   sourceName=targetName+".xml";
1255   targetName=fileName;
1256   targetName.replace(targetName.find(".med"),4,"_remergedFrom5_");
1257   cmd+=" --input-file="+sourceName+" --output-file="+targetName+" --verbose="+IntToStr(_verbose);
1258   if (_verbose) cout<<endl<<cmd<<endl;
1259   res=system(cmd.c_str());
1260   CPPUNIT_ASSERT_EQUAL(0, res);
1261
1262   string refusedName=targetName+"1.med";
1263   ParaMEDMEM::MEDFileUMesh* refusedMesh=ParaMEDMEM::MEDFileUMesh::New(refusedName.c_str(),_mesh_name.c_str());
1264   ParaMEDMEM::MEDCouplingUMesh* refusedCellMesh=refusedMesh->getLevel0Mesh(false);
1265
1266   CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), refusedCellMesh->getNumberOfCells());
1267
1268   std::vector<const MEDCouplingUMesh *> meshes;
1269   std::vector<DataArrayInt *> corr;
1270   meshes.push_back(cellMesh);
1271   refusedCellMesh->tryToShareSameCoordsPermute(*cellMesh, 1e-9);
1272   meshes.push_back(refusedCellMesh);
1273   MEDCouplingUMesh* fusedCell=MEDCouplingUMesh::FuseUMeshesOnSameCoords(meshes,0,corr);
1274   CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), fusedCell->getNumberOfCells());
1275
1276   MEDCouplingFieldDouble* field1=MEDLoader::ReadField(ON_GAUSS_NE,fileName.c_str(),initialMesh->getName().c_str(),0,"MyFieldOnGaussNE",5,6);
1277   MEDCouplingFieldDouble* field2=MEDLoader::ReadField(ON_GAUSS_NE,refusedName.c_str(),refusedCellMesh->getName().c_str(),0,"MyFieldOnGaussNE",5,6);
1278
1279   int nbcells=corr[1]->getNumberOfTuples();
1280   CPPUNIT_ASSERT_EQUAL(cellMesh->getNumberOfCells(), nbcells);
1281   //use corr to test equality of field
1282   DataArrayDouble* f1=field1->getArray();
1283   DataArrayDouble* f2=field2->getArray();
1284   if (_verbose>300)
1285     {
1286       cout << "\nf1 : " << f1->reprZip(); //123.4 for 12th cell,3rd component, 4th gausspoint
1287       cout << "\nf2 : " << f2->reprZip(); //field2->advancedRepradvancedRepr();
1288       for (std::size_t i = 0; i < corr.size(); i++)
1289         cout << "\ncorr " << i << " : " << corr[i]->reprZip();
1290
1291     }
1292   int nbequal=0;
1293   int nbptgauss=8;
1294   int nbcomp=field1->getNumberOfComponents();
1295   double* p1=f1->getPointer();
1296   double* p2=f2->getPointer();
1297   int* pc=corr[1]->getPointer();
1298   for (int i = 0; i < nbcells; i++)
1299     {
1300       int i1=pc[i]*nbcomp*nbptgauss;
1301       int i2=i*nbcomp*nbptgauss;
1302       for (int j = 0; j < nbcomp*nbptgauss; j++)
1303         {
1304           if (p1[i1+j]==p2[i2+j]) nbequal++;
1305           //cout<<" "<<p1[i1+j]<<"="<<p2[i2+j];
1306         }
1307     }
1308   CPPUNIT_ASSERT_EQUAL(nbcells*nbcomp*nbptgauss, nbequal);
1309
1310   for (std::size_t i = 0; i < corr.size(); i++)
1311     corr[i]->decrRef();
1312   field1->decrRef();
1313   field2->decrRef();
1314   fusedCell->decrRef();
1315   refusedMesh->decrRef();
1316   refusedCellMesh->decrRef();
1317   cellMesh->decrRef();
1318   initialMesh->decrRef();
1319 }
1320
1321 //================================================================================
1322 /*!
1323  * \brief Test for 0021756: [CEA 602] MEDPartitioner improvements
1324  */
1325 //================================================================================
1326
1327 void MEDPARTITIONERTest::testCreateBoundaryFaces2D()
1328 {
1329   // Fixed complains are:
1330   // - 2D is not available
1331   // - groups and family handling is bugged (probably due to bug in the handling
1332   //   of arrayTo in castIntField())
1333   // - creates boundary faces option is not handled
1334
1335   // Create a 2D mesh in a file
1336
1337   const char fileName[] = "tmp_testCreateBoundaryFaces2D.med";
1338
1339   const int idFam1 = 3, idFam2 = 2;
1340   int nbFam1, nbFam2, nbc;
1341   {
1342     const int nbX = 20, nbY = 15;
1343     vector<int> conn;
1344     vector<double> coor;
1345     for (int j=0; j<=nbY; j++)
1346       for (int i=0; i<=nbX; i++)
1347         {
1348           coor.push_back(i+.1);
1349           coor.push_back(j+.2);
1350         }
1351     int ii;
1352     for (int j=0; j<nbY; j++)
1353       for (int i=0; i<nbX; i++)
1354         {
1355           ii=i + j*(nbX+1);
1356           conn.push_back(ii);
1357           conn.push_back(ii+1);
1358           ii=ii + nbX + 2 ;
1359           conn.push_back(ii);
1360           conn.push_back(ii-1);
1361         }
1362     MEDCouplingUMesh *mesh=MEDCouplingUMesh::New();
1363     mesh->setMeshDimension(2);
1364
1365     nbc=conn.size()/4; //nb of cells
1366     mesh->allocateCells(nbc);
1367     int* pConn = &conn[0];
1368     for(int i=0; i<nbc; i++, pConn+=4)
1369       mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,pConn);
1370     mesh->finishInsertingCells();
1371
1372     int nbv=coor.size()/2; //nb of vertices
1373     DataArrayDouble *myCoords=DataArrayDouble::New();
1374     myCoords->useArray( &coor[0], /*ownership=*/false, CPP_DEALLOC, nbv, 2 );
1375     mesh->setCoords(myCoords);
1376     mesh->setName("FacesIn2D");
1377     myCoords->decrRef();
1378     mesh->checkCoherency();
1379
1380     // groups of cells
1381     DataArrayInt* cellsFam=DataArrayInt::New();
1382     cellsFam->alloc(nbc,1);
1383     nbFam1 = nbc/3, nbFam2 = nbc/2;
1384     int iE = 0;
1385     for ( int i = 0; i < nbFam1; ++i ) cellsFam->getPointer()[ iE++ ] = idFam1;
1386     for ( int i = 0; i < nbFam2; ++i ) cellsFam->getPointer()[ iE++ ] = idFam2;
1387     for (            ; iE < nbc;     ) cellsFam->getPointer()[ iE++ ] = 0;
1388     map<string,int> theFamilies;
1389     theFamilies["FAMILLE_ZERO"]=0;
1390     theFamilies["Family1"     ]=idFam1;
1391     theFamilies["Family2"     ]=idFam2;
1392     map<string, vector<string> > theGroups;
1393     theGroups["Group1"].push_back("Family1");
1394     theGroups["Group2"].push_back("Family2");
1395
1396     // write mesh
1397     MEDFileUMesh * fileMesh = MEDFileUMesh::New();
1398     fileMesh->setMeshAtLevel(0, mesh);
1399     fileMesh->setFamilyInfo(theFamilies);
1400     fileMesh->setGroupInfo(theGroups);
1401     fileMesh->setFamilyFieldArr(0, cellsFam);
1402     fileMesh->write(fileName,2);
1403
1404     cellsFam->decrRef();
1405     mesh    ->decrRef();
1406     fileMesh->decrRef();
1407
1408   } // mesh creation
1409
1410   // Partition the mesh into 4 parts
1411
1412   const int ndomains = 4;
1413   ParaDomainSelector parallelizer(false);
1414   MeshCollection collection(fileName,parallelizer);
1415   ParallelTopology* aPT = (ParallelTopology*) collection.getTopology();
1416   aPT->setGlobalNumerotationDefault(collection.getParaDomainSelector());
1417
1418   std::auto_ptr< Topology > new_topo;
1419 #if defined(MED_ENABLE_METIS) || defined(MED_ENABLE_PARMETIS)
1420   new_topo.reset( collection.createPartition(ndomains,Graph::METIS) );
1421 #endif
1422 #if defined(MED_ENABLE_SCOTCH)
1423   if ( !new_topo.get() )
1424     new_topo.reset( collection.createPartition(ndomains,Graph::SCOTCH) );
1425 #endif
1426   if ( !new_topo.get() )
1427     return;
1428
1429   // Check that "2D is available"
1430
1431   const char xmlName[] = "tmp_testCreateBoundaryFaces2D";
1432   {
1433     MyGlobals::_Creates_Boundary_Faces = true;
1434     MeshCollection new_collection(collection,new_topo.get());
1435
1436     CPPUNIT_ASSERT_EQUAL(ndomains,new_collection.getNbOfLocalMeshes());
1437     CPPUNIT_ASSERT_EQUAL(ndomains,new_collection.getNbOfGlobalMeshes());
1438     CPPUNIT_ASSERT_EQUAL(collection.getNbOfLocalCells(),new_collection.getNbOfLocalCells());
1439     CPPUNIT_ASSERT_EQUAL(0,collection.getNbOfLocalFaces());
1440     CPPUNIT_ASSERT      (new_collection.getNbOfLocalFaces() > 0 );
1441
1442     MyGlobals::_General_Informations.clear();
1443     MyGlobals::_General_Informations.push_back(SerializeFromString("finalMeshName=2D"));
1444     new_collection.write( xmlName );
1445   }
1446
1447   // Check that "groups and family handling is NOT bugged"
1448
1449   MeshCollection new_collection(std::string(xmlName)+".xml");
1450   std::map< int, int > famId2nb; // count total nb of cells in divided families
1451   std::map< int, int >::iterator id2nn;
1452   {
1453     const std::vector<ParaMEDMEM::DataArrayInt*>& famIdsVec = new_collection.getCellFamilyIds();
1454     for ( size_t i = 0; i < famIdsVec.size(); ++i )
1455       {
1456         ParaMEDMEM::DataArrayInt* famIdsArr = famIdsVec[i];
1457         for ( int j = famIdsArr->getNbOfElems()-1; j >= 0; --j )
1458           {
1459             id2nn = famId2nb.insert( make_pair( famIdsArr->getPointer()[j], 0 )).first;
1460             id2nn->second++;
1461           }
1462       }
1463   }
1464   CPPUNIT_ASSERT_EQUAL( 3, (int) famId2nb.size() ); // 3 fams/groups in all
1465   CPPUNIT_ASSERT_EQUAL( 1, (int) famId2nb.count( 0      ));
1466   CPPUNIT_ASSERT_EQUAL( 1, (int) famId2nb.count( idFam1 ));
1467   CPPUNIT_ASSERT_EQUAL( 1, (int) famId2nb.count( idFam2 ));
1468   CPPUNIT_ASSERT_EQUAL( nbFam1, famId2nb[ idFam1 ]);
1469   CPPUNIT_ASSERT_EQUAL( nbFam2, famId2nb[ idFam2 ]);
1470   CPPUNIT_ASSERT_EQUAL( nbc - nbFam1 - nbFam2, famId2nb[ 0 ]);
1471
1472   // Check that "creates boundary faces option is handled"
1473
1474   famId2nb.clear();
1475   const std::vector<ParaMEDMEM::DataArrayInt*>& famIdsVec = new_collection.getFaceFamilyIds();
1476   for ( size_t i = 0; i < famIdsVec.size(); ++i )
1477     {
1478       ParaMEDMEM::DataArrayInt* famIdsArr = famIdsVec[i];
1479       for ( int j = famIdsArr->getNbOfElems()-1; j >= 0; --j )
1480         {
1481           id2nn = famId2nb.insert( make_pair( famIdsArr->getPointer()[j], 0 )).first;
1482           id2nn->second++;
1483         }
1484     }
1485
1486   CPPUNIT_ASSERT( !famId2nb.empty() );
1487
1488   // for each "JOINT_n_p_..." group there must be "JOINT_p_n_..." group
1489   // of the same size
1490   std::map<std::string,int>& famName2id = new_collection.getFamilyInfo();
1491   std::map<std::string,int>::iterator na2id = famName2id.begin(), na2id2;
1492   std::set< int > okFamIds;
1493   okFamIds.insert(0);
1494   for ( ; na2id != famName2id.end(); ++na2id )
1495     {
1496       if ( okFamIds.count( na2id->second ) || na2id->first[0] != 'J')
1497         continue;
1498       na2id2 = na2id;
1499       bool groupOK = false;
1500       while ( !groupOK && ++na2id2 != famName2id.end() )
1501         groupOK = ( na2id2->first.find_first_not_of( na2id->first ) == std::string::npos );
1502
1503       CPPUNIT_ASSERT( groupOK );
1504       CPPUNIT_ASSERT( na2id->second != na2id2->second);
1505       CPPUNIT_ASSERT_EQUAL( 1, (int) famId2nb.count( na2id2->second ));
1506       CPPUNIT_ASSERT_EQUAL( 1, (int) famId2nb.count( na2id->second ));
1507       CPPUNIT_ASSERT_EQUAL( (int) famId2nb[ na2id2->second ],
1508                             (int) famId2nb[ na2id->second ]);
1509       okFamIds.insert( na2id2->second );
1510     }
1511 }