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