src/MEDSPLITTER/MEDSPLITTER_MESHCollectionDriver.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 #include <vector>
  20 #include <string>
  21 #include <map>
  22 #include <set>
  23
  24 #include <iostream>
  25 #include <fstream>
  26
  27 #include <libxml/tree.h>
  28 #include <libxml/parser.h>
  29 #include <libxml/xpath.h>
  30 #include <libxml/xpathInternals.h>
  31 #ifndef WIN32
  32 #include <sys/time.h>
  33 #else
  34 #include <time.h>
  35 #include <windows.h>
  36 #endif
  37 //Debug macros
  38 #include "MEDMEM_Utilities.hxx"
  39
  40 //MEDMEM includes
  41 #include "MEDMEM_DriversDef.hxx"
  42 #include "MEDMEM_Mesh.hxx"
  43 #include "MEDMEM_MedFileBrowser.hxx"
  44 #include "MEDMEM_Field.hxx"
  45 #include "MEDMEM_Meshing.hxx"
  46 #include "MEDMEM_CellModel.hxx"
  47 #include "MEDMEM_SkyLineArray.hxx"
  48 #include "MEDMEM_ConnectZone.hxx"
  49 #include "MEDMEM_MeshFuse.hxx"
  50 #include "MEDMEM_MedMeshDriver.hxx"
  51
  52 //MEDSPLITTER includes
  53 #include "MEDSPLITTER_Topology.hxx"
  54 #include "MEDSPLITTER_ParallelTopology.hxx"
  55 #include "MEDSPLITTER_SequentialTopology.hxx"
  56 #include "MEDSPLITTER_MESHCollectionDriver.hxx"
  57 #include "MEDSPLITTER_MESHCollection.hxx"
  58 #include "MEDSPLITTER_ParaDomainSelector.hxx"
  59
  60 using namespace MEDSPLITTER;
  61
  62 //template inclusion
  63 //#include "MEDSPLITTER_MESHCollectionDriver.H"
  64
  65
  66 MESHCollectionDriver::MESHCollectionDriver(MESHCollection* collection):_collection(collection)
  67 {
  68 }
  69
  70
  71 /*!reads a unique MED File v>=2.1
  72  * and mounts the corresponding mesh in memory
  73  *\param filename binary file
  74  *\param meshname mesh name in the MED file
  75  * */
  76 int MESHCollectionDriver::readSeq(char* filename, char* meshname)
  77 {
  78   const char* LOC = "MEDSPLITTER::MESHCollectionDriver::readSeq()";
  79   BEGIN_OF_MED(LOC);
  80
  81   _filename.resize(1);
  82   _filename[0]=string(filename);
  83   //puts the only mesh in the mesh vector
  84   //MEDMEM::MESH* new_mesh = new MEDMEM::MESH(MEDMEM::MED_DRIVER,filename, meshname);
  85   MEDMEM::MESH* new_mesh = new MEDMEM::MESH;
  86   MED_MESH_RDONLY_DRIVER meshDrv(filename,new_mesh);
  87   meshDrv.setMeshName( meshname );
  88   meshDrv.desactivateFacesComputation(); // we need not all faces
  89   meshDrv.open();
  90   meshDrv.read();
  91   meshDrv.close();
  92   (_collection->getMesh()).push_back(new_mesh);
  93
  94   _collection->setName(meshname);
  95   (_collection->getCZ()).clear();
  96   vector<int*> cellglobal,nodeglobal,faceglobal;
  97   cellglobal.resize(1);
  98   nodeglobal.resize(1);
  99   faceglobal.resize(1);
 100   cellglobal[0]=0;
 101   nodeglobal[0]=0;
 102   faceglobal[0]=0;
 103   //creation of topology from mesh
 104   //connectzone argument is 0
 105   ParallelTopology* aPT = new ParallelTopology
 106     ((_collection->getMesh()), (_collection->getCZ()), cellglobal, nodeglobal, faceglobal);
 107   _collection->setTopology(aPT);
 108   END_OF_MED(LOC);
 109   return 0;
 110 }
 111
 112 /*!
 113  * Reads the file structure to determine the list
 114  * of all the available fields
 115  *
 116  * \param field_names, vector<string> containing the field names
 117  * \param iternumber, vector<int> containing the iteration numbers
 118  * \param ordernumber, vector<int> containing the order numbers
 119  * \param types, vector<int> containing 0 for int fields and 1 for double fields
 120  *
 121  */
 122
 123 void MESHCollectionDriver::readFileStruct(vector <string>&  field_names,vector<int>& iternumber,vector <int>&  ordernumber, vector <int>& types)
 124 {
 125   const char* LOC = "MEDSPLITTER::MESHCollectionDriver::readFileStruct()";
 126   BEGIN_OF_MED(LOC);
 127
 128   const MEDMEM::MEDFILEBROWSER med_struct (_filename[0]);
 129   int nb_fields = med_struct.getNumberOfFields();
 130
 131   MESSAGE_MED("found "<<nb_fields<<" fields in file");
 132   vector<string> names = med_struct.getFieldNames();
 133   for (int ifield = 0; ifield < nb_fields; ifield++)
 134   {
 135     MEDMEM::VEC_DT_IT_ dtit=med_struct.getFieldIteration(names[ifield]);
 136     for (unsigned i = 0; i < dtit.size(); i++)
 137     {
 138       field_names.push_back(names[ifield]);
 139       iternumber.push_back(dtit[i].dt);
 140       ordernumber.push_back(dtit[i].it);
 141       types.push_back( MED_REEL64 == med_struct.getFieldType( names[ifield] ));
 142     }
 143   }
 144   END_OF_MED(LOC);
 145 }
 146
 147 //!retrieves the type of a field for a given fieldname
 148 int MESHCollectionDriver::getFieldType(const string& fieldname)
 149 {
 150   const char* LOC = "MEDSPLITTER::MESHCollectionDriver::getFieldType()";
 151   BEGIN_OF_MED(LOC);
 152   const MEDMEM::MEDFILEBROWSER med_struct (_filename[0]);
 153
 154   int ret = ( MED_REEL64 == med_struct.getFieldType( fieldname ));
 155   END_OF_MED(LOC);
 156
 157   return ret;
 158 }
 159 namespace
 160 {
 161 //================================================================================
 162 /*!
 163  * \brief Structure used in the method below
 164  */
 165 struct TJointData
 166 {
 167   char             _name[MED_NAME_SIZE+1];
 168   int              _nb_couples;
 169   med_2_3::med_int _distant_domain;
 170   med_2_3::med_geometry_type _geo_local, _geo_dist;
 171 };
 172 }
 173
 174 //================================================================================
 175 /*!
 176  * \brief Read CELL-CELL correspondences of joints with domains on other procs
 177  *  \param idomain - domain index to return correspondence for
 178  *  \param loc_domains - domians on this pocessor
 179  *  \param domain_selector - info on cell distribution among procs
 180  *  \param loc2glob_corr - out, correspondence pairs where distant ids are global
 181  */
 182 //================================================================================
 183
 184 void MESHCollectionDriver::readLoc2GlobCellConnect(int                 idomain,
 185                                                    const set<int>&     loc_domains,
 186                                                    ParaDomainSelector* domain_selector,
 187                                                    vector<int> &       loc2glob_corr)
 188 {
 189   using namespace med_2_3;
 190   med_2_3::med_err err;
 191
 192   // find joints of domains touching idomain and loaded on other processors
 193
 194   TJointData joint;
 195   list< TJointData > joints;
 196   int total_nb_couples = 0;
 197
 198   MEDMEM::MESH* loc_mesh = _collection->getMesh()[idomain];
 199   char* meshname = (char*) _meshname[idomain].c_str();
 200   char* filename = (char*) _filename[idomain].c_str();
 201   //cout << "#" << domain_selector->rank() << ": mesh - " << meshname << endl;
 202
 203   const int loc_mesh_dim = loc_mesh->getMeshDimension();
 204
 205   const MED_EN::medGeometryElement* types = loc_mesh->getTypes(MED_EN::MED_CELL);
 206   int                             nbtypes = loc_mesh->getNumberOfTypes(MED_EN::MED_CELL);
 207   const list<MED_EN::medGeometryElement>& all_types = MED_EN::meshEntities[ MED_EN::MED_CELL ];
 208
 209   med_idt fid = MEDfileOpen( filename, med_2_3::MED_ACC_RDONLY );
 210   int  njoint = MEDnSubdomainJoint(fid, meshname);
 211   for (int ijoint=0; ijoint<njoint; ijoint++)
 212   {
 213     char joint_description[MED_COMMENT_SIZE+1], name_distant[MED_NAME_SIZE+1];
 214     int nstep,nocstpncorr;
 215     err = med_2_3::MEDsubdomainJointInfo(fid,meshname, ijoint+1,
 216                                          joint._name, joint_description,
 217                                          &joint._distant_domain, name_distant,
 218                                          &nstep,&nocstpncorr);
 219     if ( err || loc_domains.count( joint._distant_domain ))
 220       continue;  // distant is on this proc
 221
 222     for (int itype=0; itype<nbtypes;itype++)
 223     {
 224       joint._geo_local = (med_geometry_type) types[itype];
 225       list<MED_EN::medGeometryElement>::const_iterator dist_type = all_types.begin();
 226       for ( ; dist_type != all_types.end(); ++dist_type )
 227       {
 228         if ( !_collection->isDimensionOK( *dist_type, loc_mesh_dim ))
 229           continue;
 230         joint._geo_dist = (med_geometry_type) *dist_type;
 231         err = MEDsubdomainCorrespondenceSize(fid, meshname, joint._name,
 232                                              MED_NO_DT, MED_NO_DT,
 233                                              med_2_3::MED_CELL, joint._geo_local,
 234                                              med_2_3::MED_CELL, joint._geo_dist,
 235                                              & joint._nb_couples );
 236         if ( !err && joint._nb_couples > 0 )
 237         {
 238           joints.push_back( joint );
 239           total_nb_couples += joint._nb_couples;
 240         }
 241       }
 242     }
 243   }
 244
 245   // read cell pairs of found joints and replace distant local ids with global ones
 246
 247   loc2glob_corr.resize( 2 * total_nb_couples );
 248   if ( total_nb_couples > 0 )
 249   {
 250     int* cell_corresp = & loc2glob_corr[0];
 251
 252     list< TJointData >::iterator jnt = joints.begin();
 253     for ( ; jnt != joints.end(); ++jnt )
 254     {
 255       // read cell couples
 256       err = MEDsubdomainCorrespondenceRd(fid, meshname, jnt->_name,
 257                                          MED_NO_DT,MED_NO_IT,
 258                                          med_2_3::MED_CELL, jnt->_geo_local,
 259                                          med_2_3::MED_CELL, jnt->_geo_dist,
 260                                          cell_corresp);
 261       if ( err ) continue;
 262
 263       // distant local ids -> global ids
 264       if ( int shift_to_global = domain_selector->getDomainShift( jnt->_distant_domain ))
 265         for ( int i = 0 ; i < jnt->_nb_couples; ++i )
 266           cell_corresp[ 2*i + 1 ] += shift_to_global;
 267
 268       cell_corresp += 2 * jnt->_nb_couples;
 269     }
 270   }
 271   MEDfileClose( fid );
 272 }
 273
 274 //================================================================================
 275 /*!
 276  * \brief Return mesh dimension from the distributed med file having been read
 277  */
 278 //================================================================================
 279
 280 int MESHCollectionDriver::readMeshDimension() const
 281 {
 282   const char* LOC = "MESHCollectionDriver::readMeshDimension(): ";
 283   if ( _filename.empty() || _meshname.empty() )
 284     throw MED_EXCEPTION( MEDMEM::STRING(LOC) << "file name or mesh name not available");
 285
 286   med_2_3::med_idt fid = med_2_3::MEDfileOpen(_filename[0].c_str(),med_2_3::MED_ACC_RDONLY);
 287   if ( fid < 0 )
 288     throw MED_EXCEPTION( MEDMEM::STRING(LOC) << "can't open file " << _filename[0]);
 289
 290   med_2_3::med_int meshDimension, spaceDimension;
 291   med_2_3::med_mesh_type meshType;
 292   char commentp3[MED_COMMENT_SIZE+1];
 293   char dtunittp3[MED_LNAME_SIZE+1];
 294   med_2_3::med_sorting_type sttp3;
 295   int nstep;
 296   med_2_3::med_axis_type axtypp3;
 297   char *t1pp3=new char[10*MED_SNAME_SIZE+1]; // preview 10-dimensional space :)
 298   char *t2pp3=new char[10*MED_SNAME_SIZE+1];
 299   int err= med_2_3::MEDmeshInfoByName(fid,_meshname[0].c_str(),
 300                                       &spaceDimension, &meshDimension, &meshType,
 301                                       commentp3,dtunittp3,&sttp3,&nstep,&axtypp3,t1pp3,t2pp3);
 302   delete [] t1pp3;
 303   delete [] t2pp3;
 304
 305   med_2_3::MEDfileClose(fid);
 306
 307   if ( err )
 308     throw MED_EXCEPTION( MEDMEM::STRING(LOC) << "mesh name is invalid");
 309
 310   return meshDimension;
 311 }
 312
 313 void MESHCollectionDriver::readSubdomain(vector<int*>& cellglobal,
 314                                          vector<int*>& faceglobal,
 315                                          vector<int*>& nodeglobal, int idomain)
 316 {
 317   const char* LOC = "MEDSPLITTER::MESHCollectionDriver::readSubdomain()";
 318   BEGIN_OF_MED(LOC);
 319   char file[256];
 320   char meshname[MED_NAME_SIZE+1];
 321
 322   strcpy(meshname,_meshname[idomain].c_str());
 323   strcpy(file,_filename[idomain].c_str());
 324   cout << "Reading "<<_meshname[idomain]<<" in "<<_filename[idomain]<<endl;
 325   MEDMEM::MESH* mesh = (_collection->getMesh())[idomain]=new MEDMEM::MESH;
 326   MED_MESH_RDONLY_DRIVER meshDrv(file,mesh);
 327   meshDrv.setMeshName( _meshname[idomain] );
 328   meshDrv.desactivateFacesComputation(); // else global face numbering becomes invalid
 329   meshDrv.open();
 330   meshDrv.read();
 331   meshDrv.close();
 332   cout <<"End of Read"<<endl;
 333   //reading MEDSPLITTER::CONNECTZONEs NODE/NODE and CELL/CELL
 334   med_2_3::med_idt fid = med_2_3::MEDfileOpen(file,med_2_3::MED_ACC_RDONLY);
 335   med_2_3::med_int njoint = med_2_3::MEDnSubdomainJoint(fid, meshname);
 336   for (int ijoint=0; ijoint<njoint; ijoint++)
 337   {
 338     int distant;
 339     char joint_description[MED_COMMENT_SIZE+1];
 340     char name[MED_NAME_SIZE+1];
 341     char name_distant[MED_NAME_SIZE+1];
 342     int nstep,nocstpncorr;
 343     int ncorr = med_2_3::MEDsubdomainJointInfo(fid,meshname, ijoint+1, name,
 344                                                joint_description,
 345                                                &distant, name_distant,&nstep,&nocstpncorr);
 346     for (int ic=0; ic<ncorr; ic++)
 347     {
 348       med_2_3::med_entity_type cor_typent_local;
 349       med_2_3::med_geometry_type cor_typgeo_local;
 350       med_2_3::med_entity_type cor_typent_dist;
 351       med_2_3::med_geometry_type cor_typgeo_dist;
 352
 353       int ncouples;
 354       med_2_3::MEDsubdomainCorrespondenceSizeInfo(fid, meshname, name, MED_NO_DT,MED_NO_IT, ic+1,
 355                                                   &cor_typent_local,  &cor_typgeo_local,
 356                                                   &cor_typent_dist, &cor_typgeo_dist,&ncouples
 357                                                   );
 358       int* node_corresp=new int[ncouples*2];
 359       if (cor_typent_local == med_2_3::MED_NODE && cor_typent_dist == med_2_3::MED_NODE)
 360       {
 361         med_2_3::MEDsubdomainCorrespondenceRd(fid, meshname, name,
 362                                               MED_NO_DT,MED_NO_IT,
 363                                               cor_typent_local, cor_typgeo_local,
 364                                               cor_typent_dist,  cor_typgeo_dist,
 365                                               node_corresp);
 366       }
 367       //constructing the connect zone and adding it to the connect zone list
 368       MEDMEM::CONNECTZONE* cz = new MEDMEM::CONNECTZONE();
 369       cz->setName(string(name));
 370       cz->setDescription(joint_description);
 371       cz->setLocalDomainNumber(idomain);
 372       cz->setDistantDomainNumber(distant);
 373       cz->setLocalMesh((_collection->getMesh())[idomain]);
 374       cz->setDistantMesh((_collection->getMesh())[distant]);
 375       cz->setNodeCorresp(node_corresp,ncouples);
 376       (_collection->getCZ()).push_back(cz);
 377
 378     }//loop on correspom_topology->nbDomain())ndances
 379   }//loop on joints
 380
 381   // Reading global numbering
 382   //
 383   int ncell=(_collection->getMesh())[idomain]->getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
 384   if (ncell>0)
 385   {
 386     int * array=new int[ncell];
 387     int offset=0;
 388     MESSAGE_MED("Reading cell global numbering for mesh "<< idomain);
 389     list<MED_EN::medGeometryElement>::const_iterator iter;
 390     char meshchar[MED_NAME_SIZE+1];
 391     strcpy(meshchar,(_collection->getMesh())[idomain]->getName().c_str());
 392     int nbtypes = (_collection->getMesh())[idomain]->getNumberOfTypes(MED_EN::MED_CELL);
 393     const MED_EN::medGeometryElement* types =(_collection->getMesh())[idomain]->getTypes(MED_EN::MED_CELL);
 394     for (int itype=0; itype<nbtypes;itype++)
 395     {
 396       MED_EN::medGeometryElement type=types[itype];
 397       if (!_collection->isDimensionOK(type,(_collection->getMesh())[idomain]->getMeshDimension()))
 398         continue;
 399       int ntype = (_collection->getMesh())[idomain]->getNumberOfElements(MED_EN::MED_CELL,type);
 400       if (ntype==0) continue;
 401       med_2_3::MEDmeshGlobalNumberRd(fid, meshname,
 402                                      MED_NO_DT, MED_NO_IT,
 403                                      med_2_3::MED_CELL, (med_2_3::med_geometry_type)type,
 404                                      array+offset );
 405       offset+=ntype;
 406     }
 407     cellglobal[idomain]=array;
 408   }
 409
 410   MESSAGE_MED("Reading node global numbering");
 411   int nnode= (_collection->getMesh())[idomain]->getNumberOfNodes();
 412   {
 413     int* array=new int[nnode];
 414     med_2_3::MEDmeshGlobalNumberRd(fid,meshname,
 415                                    MED_NO_DT, MED_NO_IT,
 416                                    med_2_3::MED_NODE, MED_POINT1,
 417                                    array);
 418     nodeglobal[idomain]=array;
 419   }
 420
 421   MESSAGE_MED("Reading face global numbering for mesh "<<idomain);
 422   MED_EN::medEntityMesh entity =
 423     (mesh->getMeshDimension()==3)?MED_EN::MED_FACE:MED_EN::MED_EDGE;
 424   int nbface=(_collection->getMesh())[idomain]->getNumberOfElements(entity,MED_EN::MED_ALL_ELEMENTS);
 425   if (nbface!=0)
 426   {
 427     int* array=new int[nbface];
 428     int offset=0;
 429     int nbtypes = mesh->getNumberOfTypes( entity );
 430     const MED_EN::medGeometryElement* types =mesh->getTypes(entity);
 431
 432     for (int itype=0; itype< nbtypes; itype++)
 433     {
 434       MED_EN::medGeometryElement type = types[itype];
 435       if (!_collection->isDimensionOK(type,mesh->getMeshDimension()-1)) continue;
 436
 437       int ntype = mesh->getNumberOfElements(entity,type);
 438       if (ntype==0) continue;
 439       med_2_3::MEDmeshGlobalNumberRd( fid, meshname,
 440                                       MED_NO_DT, MED_NO_IT,
 441                                       med_2_3::MED_CELL, (med_2_3::med_geometry_type)type,
 442                                       array+offset );
 443       offset+=ntype;
 444     }
 445     faceglobal[idomain]=array;
 446   }
 447   med_2_3::MEDfileClose(fid);
 448
 449   END_OF_MED(LOC);
 450 }
 451
 452 void MESHCollectionDriver::writeSubdomain(int idomain, int nbdomains, char* distfilename,
 453                                           ParaDomainSelector* /*domain_selector*/)
 454 {
 455   MESSAGE_MED(" Number of connect zones "<<(_collection->getCZ()).size());
 456
 457   //writing connect zones in joints
 458
 459   med_2_3::med_idt fid = med_2_3::MEDfileOpen(distfilename,med_2_3::MED_ACC_RDWR);
 460
 461   int index_joint=0;
 462
 463
 464   for (unsigned icz=0; icz<(_collection->getCZ()).size(); icz++)
 465   {
 466     if ((_collection->getCZ())[icz]->getLocalDomainNumber()==idomain)
 467     {
 468       med_2_3::med_err error;
 469       int idistant=(_collection->getCZ())[icz]->getDistantDomainNumber();
 470       char joint_name[MED_NAME_SIZE+1];
 471       sprintf(joint_name,"joint_%i",idistant+1);
 472       char desc[MED_COMMENT_SIZE+1];
 473       sprintf(desc,"connect_zone_%d",icz+1);
 474
 475       char distant_name[MED_NAME_SIZE+1];
 476       //sprintf(distant_name,"domain_%i",(_collection->getCZ())[icz]->getDistantDomainNumber());
 477
 478       //        sprintf(distant_name,(_collection->getMesh())[idistant]->getName().c_str());
 479       sprintf(distant_name,"domain_%i",idistant);
 480       char mesh_name[MED_NAME_SIZE+1];
 481
 482       strcpy (mesh_name, _collection->getMesh(idomain)->getName().c_str());
 483       SCRUTE_MED(_collection->getMesh(idomain)->getName());
 484       error = med_2_3::MEDsubdomainJointCr(fid,mesh_name, joint_name, desc,
 485                                            idistant, distant_name);
 486       if (error==-1) cout << "erreur creation de joint "<<endl;
 487
 488       /////////////////////////////////////////
 489       //writing node/node correspondency
 490       /////////////////////////////////////////
 491       int nbnodes=(_collection->getCZ())[icz]->getNodeNumber();
 492       int* node_corresp=const_cast<int*>((_collection->getCZ())[icz]->getNodeCorrespValue());
 493
 494       /* Nodes are reordered so that the ordering on the local and the distant domain
 495          correspond. The chosen order is the natural ordering on the domain
 496          with lowest proc id*/
 497       if (_collection->getSubdomainBoundaryCreates())
 498       {
 499         if (idomain<idistant)
 500           jointSort(node_corresp, nbnodes, true);
 501         else
 502           jointSort(node_corresp, nbnodes, false);
 503       }
 504       error=
 505         med_2_3::MEDsubdomainCorrespondenceWr(fid, mesh_name, joint_name,
 506                                               MED_NO_DT, MED_NO_IT,
 507                                               med_2_3::MED_NODE, MED_POINT1,
 508                                               med_2_3::MED_NODE, MED_POINT1,
 509                                               nbnodes, node_corresp);
 510       if (error==-1) cout << "erreur creation de joint "<<endl;
 511
 512       //writing cell/cell joint
 513       writeElementJoint(MED_EN::MED_CELL, icz, idomain, idistant, mesh_name,joint_name,fid);
 514       //writing face/face joint
 515       if (_collection->getSubdomainBoundaryCreates())
 516       {
 517         MED_EN::medEntityMesh constituent_entity =
 518           (_collection->getMeshDimension()==3)?MED_EN::MED_FACE:MED_EN::MED_EDGE;
 519         writeElementJoint(constituent_entity, icz, idomain, idistant, mesh_name,joint_name,fid);
 520       }
 521       index_joint++;
 522     }
 523   }
 524
 525   char meshchar[MED_NAME_SIZE+1];
 526   strcpy(meshchar,(_collection->getMesh())[idomain]->getName().c_str());
 527
 528   // Writing cell global numbering
 529   //
 530   {
 531     int ncell=_collection->getTopology()->getCellNumber(idomain);
 532     int* array=new int[ncell];
 533     _collection->getTopology()->getCellList(idomain,array);
 534     int offset=0;
 535
 536     MED_EN::MESH_ENTITIES::const_iterator currentEntity;
 537     list<MED_EN::medGeometryElement>::const_iterator iter;
 538     currentEntity  = MED_EN::meshEntities.find(MED_EN::MED_CELL);
 539
 540     int nbtypes = (_collection->getMesh())[idomain]->getNumberOfTypes(MED_EN::MED_CELL);
 541     const MED_EN::medGeometryElement* types =(_collection->getMesh())[idomain]->getTypes(MED_EN::MED_CELL);
 542     for (int itype=0; itype<nbtypes;itype++)
 543     {
 544       MED_EN::medGeometryElement type=types[itype];
 545       if (!_collection->isDimensionOK(type,(_collection->getMesh())[idomain]->getMeshDimension()))
 546         continue;
 547       int ntype = (_collection->getMesh())[idomain]->getNumberOfElements(MED_EN::MED_CELL,type);
 548       if (ntype==0) continue;
 549       med_2_3::MEDmeshGlobalNumberWr(fid,meshchar,
 550                                      MED_NO_DT, MED_NO_IT,
 551                                      med_2_3::MED_CELL, (med_2_3::med_geometry_type)type,
 552                                      ntype, array+offset);
 553       offset+=ntype;
 554     }
 555     delete[] array;
 556   }
 557
 558   MED_EN::medEntityMesh constituent_entity;
 559   if (_collection->getMeshDimension()==3)
 560     constituent_entity=MED_EN::MED_FACE;
 561   else if (_collection->getMeshDimension()==2)
 562     constituent_entity=MED_EN::MED_EDGE;
 563   else throw MEDEXCEPTION("Wrong dimension");
 564
 565
 566   //writing face global numbering
 567   {
 568     int offset=0;
 569     int nface= _collection->getTopology()->getFaceNumber(idomain);
 570     int * array=new int[nface];
 571     _collection->getTopology()->getFaceList(idomain,array);
 572
 573     int nbfacetypes = 0;
 574     const MED_EN::medGeometryElement* facetypes = 0;
 575     if ( _collection->getMesh()[idomain]->getNumberOfElements( constituent_entity, MED_ALL_ELEMENTS))
 576       {
 577         nbfacetypes = (_collection->getMesh())[idomain]->getNumberOfTypes(constituent_entity);
 578         if (nbfacetypes>0)
 579           facetypes =(_collection->getMesh())[idomain]->getTypes(constituent_entity);
 580       }
 581     for (int itype=0; itype<nbfacetypes;itype++)
 582       {
 583         MED_EN::medGeometryElement type=facetypes[itype];
 584         if (!_collection->isDimensionOK(type,_collection->getMeshDimension()-1)) continue;
 585
 586         int ntype = (_collection->getMesh())[idomain]->getNumberOfElements(constituent_entity,type);
 587         if (ntype==0) continue;
 588         med_2_3::MEDmeshGlobalNumberWr(fid,meshchar,
 589                                        MED_NO_DT, MED_NO_IT,
 590                                        med_2_3::MED_CELL, (med_2_3::med_geometry_type)type,
 591                                        ntype, array+offset );
 592         offset+=ntype;
 593       }
 594     delete[] array;
 595   }
 596
 597
 598   //writing node global numbering
 599   {
 600     int nnode=_collection->getTopology()->getNodeNumber(idomain);
 601     int* array = new int[nnode];
 602     _collection->getTopology()->getNodeList(idomain,array);
 603
 604     med_2_3::MEDmeshGlobalNumberWr(fid,meshchar,
 605                                    MED_NO_DT, MED_NO_IT,
 606                                    med_2_3::MED_NODE, MED_POINT1,
 607                                    nnode, array);
 608     delete[] array;
 609   }
 610
 611   med_2_3::MEDfileClose(fid);
 612   MESSAGE_MED("End of writing");
 613
 614 }
 615
 616 void MESHCollectionDriver::writeElementJoint(medEntityMesh entity ,
 617                                              int icz,
 618                                              int idomain,
 619                                              int idistant,
 620                                              char* mesh_name,
 621                                              char* joint_name,
 622                                              med_2_3::med_idt fid )
 623 {
 624   //////////////////////////////////////////
 625   //writing cell/cell correspondency
 626   //////////////////////////////////////////
 627   int nbcells=(_collection->getCZ())[icz]->getEntityCorrespNumber(entity,entity);
 628   const int* index = (_collection->getCZ())[icz]->getEntityCorrespIndex(entity,entity);
 629   const int* value = (_collection->getCZ())[icz]->getEntityCorrespValue(entity,entity);
 630   if ( nbcells==0 ) return; // e.g. domains have 1 common node
 631
 632   map <pair <MED_EN::medGeometryElement, MED_EN::medGeometryElement> , vector<int> > cellmap;
 633   map <MED_EN::medGeometryElement, int> local_offset;
 634   map <MED_EN::medGeometryElement, int> distant_offset;
 635
 636   //definition of the local offsets for the types present on local
 637   //and distant domains
 638   // for a mesh containing 2 triangles and 3 quads
 639   //local_offset[TRIA3]=0
 640   //local_offset[QUAD4]=2
 641
 642   int nb_types_local=(_collection->getMesh())[idomain]-> getNumberOfTypes(entity);
 643   const MED_EN::medGeometryElement* local_types = (_collection->getMesh())[idomain]->getTypes(entity);
 644   const int* local_gni = (_collection->getMesh())[idomain]-> getGlobalNumberingIndex(entity);
 645   for (int i=0; i< nb_types_local; i++)
 646   {
 647     local_offset[local_types[i]]=local_gni[i]-1;
 648   }
 649
 650   int nb_types_distant=(_collection->getMesh())[idistant]-> getNumberOfTypes(entity);
 651   const MED_EN::medGeometryElement* distant_types = (_collection->getMesh())[idistant]->getTypes(entity);
 652   const int* distant_gni = (_collection->getMesh())[idistant]-> getGlobalNumberingIndex(entity);
 653   for (int i=0; i< nb_types_distant; i++)
 654   {
 655     distant_offset[distant_types[i]]=distant_gni[i]-1;
 656   }
 657
 658   if (nb_types_local==1 && nb_types_distant==1)
 659   {
 660     MED_EN::medGeometryElement local_type =  (_collection->getMesh())[idomain]->getElementType(entity,1);
 661     MED_EN::medGeometryElement distant_type = (_collection->getMesh())[idistant]->getElementType(entity,1);
 662     vector<int> corresp;
 663     for (int i=0; i<nbcells; i++)
 664       for (int icol = index[i]-1; icol<index[i+1]-1; icol++)
 665       {
 666         corresp.push_back(i+1);
 667         corresp.push_back(value[icol]);
 668       }
 669     int size_joint = corresp.size()/2;
 670     med_2_3::MEDsubdomainCorrespondenceWr(fid, mesh_name, joint_name,
 671                                           MED_NO_DT, MED_NO_IT,
 672                                           med_2_3::MED_CELL,(med_2_3::med_geometry_type)local_type,
 673                                           med_2_3::MED_CELL,(med_2_3::med_geometry_type)distant_type,
 674                                           size_joint,&corresp[0]);
 675   }
 676   else
 677   {
 678     //classifying all the cell/cell relationships into geomtype/geomtype relationships
 679     //there exists a vector for each geomtype/geomtype pair
 680     // the vectors are stored in cellmap, a std::map with a pair<geomtype,geomtype> key
 681
 682     for (int i=0; i<nbcells; i++)
 683       for (int icol = index[i]-1; icol<index[i+1]-1; icol++)
 684       {
 685         MED_EN::medGeometryElement local_type =  (_collection->getMesh())[idomain]->getElementType(entity,i+1);
 686         MED_EN::medGeometryElement distant_type = (_collection->getMesh())[idistant]->getElementType(entity,value[icol]);
 687
 688         cellmap[make_pair(local_type, distant_type)].push_back(i+1-local_offset[local_type]);
 689         cellmap[make_pair(local_type, distant_type)].push_back(value[icol]-distant_offset[distant_type]);
 690
 691       }
 692     map <pair <MED_EN::medGeometryElement, MED_EN::medGeometryElement> , vector<int> >::const_iterator iter;
 693
 694     //going through all the (geom,geom) pairs and writing the joints
 695     for (iter= cellmap.begin(); iter != cellmap.end(); iter++)
 696     {
 697       int size= iter->second.size();
 698       int *corresp = new int[size];
 699       for (int ind=0; ind < size; ind++)
 700         corresp[ind]=(iter->second)[ind];
 701       med_2_3::med_geometry_type local_geo_elem=(med_2_3::med_geometry_type)iter->first.first;
 702       med_2_3::med_geometry_type distant_geo_elem=(med_2_3::med_geometry_type)iter->first.second;
 703       int size_joint=size/2;
 704       //med_2_3::med_err error =
 705       med_2_3::MEDsubdomainCorrespondenceWr(fid, mesh_name, joint_name,
 706                                             MED_NO_DT, MED_NO_IT,
 707                                             med_2_3::MED_CELL, local_geo_elem,
 708                                             med_2_3::MED_CELL, distant_geo_elem,
 709                                             size_joint, corresp);
 710       delete[] corresp;
 711     }
 712     // MED v 2.3.1 returns an error code when
 713     // writing a joint that is already present in the file.
 714     // Also, it returns an error code if a joint
 715     // concerns 3D elements.
 716     // Until these two items are not
 717     // changed, the following line must be commented out
 718
 719     //if (error==-1) throw MEDEXCEPTION("Error filling joint");
 720
 721
 722   }
 723 }
 724
 725 void MESHCollectionDriver::jointSort(int* elems, int nbelems, bool is_first)
 726 {
 727   //filling an ordered structure with the elem ids
 728   map <int,int> nodemap;
 729   if (is_first)
 730     for (int i=0; i<nbelems; i++)
 731       nodemap.insert(make_pair(elems[2*i],elems[2*i+1]));
 732
 733   else
 734     for (int i=0; i<nbelems; i++)
 735       nodemap.insert(make_pair(elems[2*i+1],elems[2*i]));
 736
 737   int* ptr_elems=elems;
 738
 739   //filling the vector in appropriate order
 740   for (map<int,int>::const_iterator iter=nodemap.begin(); iter!=nodemap.end(); iter++)
 741   {
 742     if (is_first)
 743     {
 744       *ptr_elems++=iter->first;
 745       *ptr_elems++=iter->second;
 746     }
 747     else
 748     {
 749       *ptr_elems++=iter->second;
 750       *ptr_elems++=iter->first;
 751     }
 752   }
 753 }