Back to index

salome-med  6.5.0
MEDSPLITTER_MESHCollectionDriver.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 #include <vector>
00020 #include <string>
00021 #include <map>
00022 #include <set>
00023 
00024 #include <iostream>
00025 #include <fstream>
00026 
00027 #include <libxml/tree.h>
00028 #include <libxml/parser.h>
00029 #include <libxml/xpath.h>
00030 #include <libxml/xpathInternals.h>
00031 #ifndef WIN32
00032 #include <sys/time.h>
00033 #else
00034 #include <time.h>
00035 #include <windows.h>
00036 #endif
00037 //Debug macros
00038 #include "MEDMEM_Utilities.hxx"
00039 
00040 //MEDMEM includes
00041 #include "MEDMEM_DriversDef.hxx"
00042 #include "MEDMEM_Mesh.hxx"
00043 #include "MEDMEM_MedFileBrowser.hxx"
00044 #include "MEDMEM_Field.hxx"
00045 #include "MEDMEM_Meshing.hxx"
00046 #include "MEDMEM_CellModel.hxx"
00047 #include "MEDMEM_SkyLineArray.hxx"
00048 #include "MEDMEM_ConnectZone.hxx"
00049 #include "MEDMEM_MeshFuse.hxx"
00050 #include "MEDMEM_MedMeshDriver.hxx"
00051 
00052 //MEDSPLITTER includes
00053 #include "MEDSPLITTER_Topology.hxx"
00054 #include "MEDSPLITTER_ParallelTopology.hxx"
00055 #include "MEDSPLITTER_SequentialTopology.hxx"
00056 #include "MEDSPLITTER_MESHCollectionDriver.hxx"
00057 #include "MEDSPLITTER_MESHCollection.hxx"
00058 #include "MEDSPLITTER_ParaDomainSelector.hxx"
00059 
00060 using namespace MEDSPLITTER;
00061 
00062 //template inclusion
00063 //#include "MEDSPLITTER_MESHCollectionDriver.H"
00064 
00065 
00066 MESHCollectionDriver::MESHCollectionDriver(MESHCollection* collection):_collection(collection)
00067 {
00068 }
00069 
00070 
00076 int MESHCollectionDriver::readSeq(char* filename, char* meshname)
00077 {
00078   const char* LOC = "MEDSPLITTER::MESHCollectionDriver::readSeq()";
00079   BEGIN_OF_MED(LOC);
00080 
00081   _filename.resize(1);
00082   _filename[0]=string(filename);
00083   //puts the only mesh in the mesh vector
00084   //MEDMEM::MESH* new_mesh = new MEDMEM::MESH(MEDMEM::MED_DRIVER,filename, meshname);
00085   MEDMEM::MESH* new_mesh = new MEDMEM::MESH;
00086   MED_MESH_RDONLY_DRIVER meshDrv(filename,new_mesh);
00087   meshDrv.setMeshName( meshname );
00088   meshDrv.desactivateFacesComputation(); // we need not all faces
00089   meshDrv.open();
00090   meshDrv.read();
00091   meshDrv.close();
00092   (_collection->getMesh()).push_back(new_mesh);
00093 
00094   _collection->setName(meshname);
00095   (_collection->getCZ()).clear();
00096   vector<int*> cellglobal,nodeglobal,faceglobal;
00097   cellglobal.resize(1);
00098   nodeglobal.resize(1);
00099   faceglobal.resize(1);
00100   cellglobal[0]=0;
00101   nodeglobal[0]=0;
00102   faceglobal[0]=0;
00103   //creation of topology from mesh 
00104   //connectzone argument is 0
00105   ParallelTopology* aPT = new ParallelTopology
00106     ((_collection->getMesh()), (_collection->getCZ()), cellglobal, nodeglobal, faceglobal);
00107   _collection->setTopology(aPT);
00108   END_OF_MED(LOC);
00109   return 0;
00110 }
00111 
00123 void MESHCollectionDriver::readFileStruct(vector <string>&  field_names,vector<int>& iternumber,vector <int>&  ordernumber, vector <int>& types)
00124 {
00125   const char* LOC = "MEDSPLITTER::MESHCollectionDriver::readFileStruct()";
00126   BEGIN_OF_MED(LOC);
00127 
00128   const MEDMEM::MEDFILEBROWSER med_struct (_filename[0]);
00129   int nb_fields = med_struct.getNumberOfFields();
00130 
00131   MESSAGE_MED("found "<<nb_fields<<" fields in file");
00132   vector<string> names = med_struct.getFieldNames();
00133   for (int ifield = 0; ifield < nb_fields; ifield++)
00134   {
00135     MEDMEM::VEC_DT_IT_ dtit=med_struct.getFieldIteration(names[ifield]);
00136     for (unsigned i = 0; i < dtit.size(); i++)
00137     {
00138       field_names.push_back(names[ifield]);
00139       iternumber.push_back(dtit[i].dt);
00140       ordernumber.push_back(dtit[i].it);
00141       types.push_back( MED_REEL64 == med_struct.getFieldType( names[ifield] ));
00142     }
00143   }
00144   END_OF_MED(LOC);
00145 }
00146 
00148 int MESHCollectionDriver::getFieldType(const string& fieldname)
00149 {
00150   const char* LOC = "MEDSPLITTER::MESHCollectionDriver::getFieldType()";
00151   BEGIN_OF_MED(LOC);
00152   const MEDMEM::MEDFILEBROWSER med_struct (_filename[0]);
00153 
00154   int ret = ( MED_REEL64 == med_struct.getFieldType( fieldname ));
00155   END_OF_MED(LOC);
00156 
00157   return ret;
00158 }
00159 namespace
00160 {
00161 //================================================================================
00165 struct TJointData
00166 {
00167   char             _name[MED_NAME_SIZE+1];
00168   int              _nb_couples;
00169   med_2_3::med_int _distant_domain;
00170   med_2_3::med_geometry_type _geo_local, _geo_dist;
00171 };
00172 }
00173 
00174 //================================================================================
00182 //================================================================================
00183 
00184 void MESHCollectionDriver::readLoc2GlobCellConnect(int                 idomain,
00185                                                    const set<int>&     loc_domains,
00186                                                    ParaDomainSelector* domain_selector,
00187                                                    vector<int> &       loc2glob_corr)
00188 {
00189   using namespace med_2_3;
00190   med_2_3::med_err err;
00191 
00192   // find joints of domains touching idomain and loaded on other processors
00193 
00194   TJointData joint;
00195   list< TJointData > joints;
00196   int total_nb_couples = 0;
00197 
00198   MEDMEM::MESH* loc_mesh = _collection->getMesh()[idomain];
00199   char* meshname = (char*) _meshname[idomain].c_str();
00200   char* filename = (char*) _filename[idomain].c_str();
00201   //cout << "#" << domain_selector->rank() << ": mesh - " << meshname << endl;
00202 
00203   const int loc_mesh_dim = loc_mesh->getMeshDimension();
00204 
00205   const MED_EN::medGeometryElement* types = loc_mesh->getTypes(MED_EN::MED_CELL);
00206   int                             nbtypes = loc_mesh->getNumberOfTypes(MED_EN::MED_CELL);
00207   const list<MED_EN::medGeometryElement>& all_types = MED_EN::meshEntities[ MED_EN::MED_CELL ];
00208 
00209   med_idt fid = MEDfileOpen( filename, med_2_3::MED_ACC_RDONLY );
00210   int  njoint = MEDnSubdomainJoint(fid, meshname);
00211   for (int ijoint=0; ijoint<njoint; ijoint++)
00212   {
00213     char joint_description[MED_COMMENT_SIZE+1], name_distant[MED_NAME_SIZE+1];
00214     int nstep,nocstpncorr;
00215     err = med_2_3::MEDsubdomainJointInfo(fid,meshname, ijoint+1,
00216                                          joint._name, joint_description,
00217                                          &joint._distant_domain, name_distant,
00218                                          &nstep,&nocstpncorr);
00219     if ( err || loc_domains.count( joint._distant_domain ))
00220       continue;  // distant is on this proc
00221 
00222     for (int itype=0; itype<nbtypes;itype++)
00223     {
00224       joint._geo_local = (med_geometry_type) types[itype];
00225       list<MED_EN::medGeometryElement>::const_iterator dist_type = all_types.begin();
00226       for ( ; dist_type != all_types.end(); ++dist_type )
00227       {
00228         if ( !_collection->isDimensionOK( *dist_type, loc_mesh_dim ))
00229           continue;
00230         joint._geo_dist = (med_geometry_type) *dist_type;
00231         err = MEDsubdomainCorrespondenceSize(fid, meshname, joint._name,
00232                                              MED_NO_DT, MED_NO_DT,
00233                                              med_2_3::MED_CELL, joint._geo_local,
00234                                              med_2_3::MED_CELL, joint._geo_dist,
00235                                              & joint._nb_couples );
00236         if ( !err && joint._nb_couples > 0 )
00237         {
00238           joints.push_back( joint );
00239           total_nb_couples += joint._nb_couples;
00240         }
00241       }
00242     }
00243   }
00244 
00245   // read cell pairs of found joints and replace distant local ids with global ones
00246 
00247   loc2glob_corr.resize( 2 * total_nb_couples );
00248   if ( total_nb_couples > 0 )
00249   {
00250     int* cell_corresp = & loc2glob_corr[0];
00251 
00252     list< TJointData >::iterator jnt = joints.begin();
00253     for ( ; jnt != joints.end(); ++jnt )
00254     {
00255       // read cell couples
00256       err = MEDsubdomainCorrespondenceRd(fid, meshname, jnt->_name,
00257                                          MED_NO_DT,MED_NO_IT, 
00258                                          med_2_3::MED_CELL, jnt->_geo_local,
00259                                          med_2_3::MED_CELL, jnt->_geo_dist,
00260                                          cell_corresp);
00261       if ( err ) continue;
00262 
00263       // distant local ids -> global ids
00264       if ( int shift_to_global = domain_selector->getDomainShift( jnt->_distant_domain ))
00265         for ( int i = 0 ; i < jnt->_nb_couples; ++i )
00266           cell_corresp[ 2*i + 1 ] += shift_to_global;
00267 
00268       cell_corresp += 2 * jnt->_nb_couples;
00269     }
00270   }
00271   MEDfileClose( fid );
00272 }
00273 
00274 //================================================================================
00278 //================================================================================
00279 
00280 int MESHCollectionDriver::readMeshDimension() const
00281 {
00282   const char* LOC = "MESHCollectionDriver::readMeshDimension(): ";
00283   if ( _filename.empty() || _meshname.empty() )
00284     throw MED_EXCEPTION( MEDMEM::STRING(LOC) << "file name or mesh name not available");
00285 
00286   med_2_3::med_idt fid = med_2_3::MEDfileOpen(_filename[0].c_str(),med_2_3::MED_ACC_RDONLY);
00287   if ( fid < 0 )
00288     throw MED_EXCEPTION( MEDMEM::STRING(LOC) << "can't open file " << _filename[0]);
00289 
00290   med_2_3::med_int meshDimension, spaceDimension;
00291   med_2_3::med_mesh_type meshType;
00292   char commentp3[MED_COMMENT_SIZE+1];
00293   char dtunittp3[MED_LNAME_SIZE+1];
00294   med_2_3::med_sorting_type sttp3;
00295   int nstep;
00296   med_2_3::med_axis_type axtypp3;
00297   char *t1pp3=new char[10*MED_SNAME_SIZE+1]; // preview 10-dimensional space :)
00298   char *t2pp3=new char[10*MED_SNAME_SIZE+1];
00299   int err= med_2_3::MEDmeshInfoByName(fid,_meshname[0].c_str(),
00300                                       &spaceDimension, &meshDimension, &meshType,
00301                                       commentp3,dtunittp3,&sttp3,&nstep,&axtypp3,t1pp3,t2pp3);
00302   delete [] t1pp3;
00303   delete [] t2pp3;
00304 
00305   med_2_3::MEDfileClose(fid);
00306 
00307   if ( err )
00308     throw MED_EXCEPTION( MEDMEM::STRING(LOC) << "mesh name is invalid");
00309 
00310   return meshDimension;
00311 }
00312 
00313 void MESHCollectionDriver::readSubdomain(vector<int*>& cellglobal,
00314                                          vector<int*>& faceglobal,
00315                                          vector<int*>& nodeglobal, int idomain)
00316 {
00317   const char* LOC = "MEDSPLITTER::MESHCollectionDriver::readSubdomain()";
00318   BEGIN_OF_MED(LOC);
00319   char file[256];
00320   char meshname[MED_NAME_SIZE+1];
00321 
00322   strcpy(meshname,_meshname[idomain].c_str());
00323   strcpy(file,_filename[idomain].c_str());
00324   cout << "Reading "<<_meshname[idomain]<<" in "<<_filename[idomain]<<endl;
00325   MEDMEM::MESH* mesh = (_collection->getMesh())[idomain]=new MEDMEM::MESH;
00326   MED_MESH_RDONLY_DRIVER meshDrv(file,mesh);
00327   meshDrv.setMeshName( _meshname[idomain] );
00328   meshDrv.desactivateFacesComputation(); // else global face numbering becomes invalid
00329   meshDrv.open();
00330   meshDrv.read();
00331   meshDrv.close();
00332   cout <<"End of Read"<<endl;
00333   //reading MEDSPLITTER::CONNECTZONEs NODE/NODE and CELL/CELL
00334   med_2_3::med_idt fid = med_2_3::MEDfileOpen(file,med_2_3::MED_ACC_RDONLY);
00335   med_2_3::med_int njoint = med_2_3::MEDnSubdomainJoint(fid, meshname);
00336   for (int ijoint=0; ijoint<njoint; ijoint++)
00337   {
00338     int distant;
00339     char joint_description[MED_COMMENT_SIZE+1];
00340     char name[MED_NAME_SIZE+1];
00341     char name_distant[MED_NAME_SIZE+1];
00342     int nstep,nocstpncorr;
00343     int ncorr = med_2_3::MEDsubdomainJointInfo(fid,meshname, ijoint+1, name, 
00344                                                joint_description,
00345                                                &distant, name_distant,&nstep,&nocstpncorr);
00346     for (int ic=0; ic<ncorr; ic++)
00347     {
00348       med_2_3::med_entity_type cor_typent_local;
00349       med_2_3::med_geometry_type cor_typgeo_local;
00350       med_2_3::med_entity_type cor_typent_dist;
00351       med_2_3::med_geometry_type cor_typgeo_dist;
00352 
00353       int ncouples;
00354       med_2_3::MEDsubdomainCorrespondenceSizeInfo(fid, meshname, name, MED_NO_DT,MED_NO_IT, ic+1,
00355                                                   &cor_typent_local,  &cor_typgeo_local,
00356                                                   &cor_typent_dist, &cor_typgeo_dist,&ncouples
00357                                                   );
00358       int* node_corresp=new int[ncouples*2];
00359       if (cor_typent_local == med_2_3::MED_NODE && cor_typent_dist == med_2_3::MED_NODE)
00360       {
00361         med_2_3::MEDsubdomainCorrespondenceRd(fid, meshname, name,
00362                                               MED_NO_DT,MED_NO_IT, 
00363                                               cor_typent_local, cor_typgeo_local,
00364                                               cor_typent_dist,  cor_typgeo_dist,
00365                                               node_corresp);
00366       }
00367       //constructing the connect zone and adding it to the connect zone list
00368       MEDMEM::CONNECTZONE* cz = new MEDMEM::CONNECTZONE();
00369       cz->setName(string(name));
00370       cz->setDescription(joint_description);
00371       cz->setLocalDomainNumber(idomain);
00372       cz->setDistantDomainNumber(distant);
00373       cz->setLocalMesh((_collection->getMesh())[idomain]);
00374       cz->setDistantMesh((_collection->getMesh())[distant]);
00375       cz->setNodeCorresp(node_corresp,ncouples);
00376       (_collection->getCZ()).push_back(cz);
00377 
00378     }//loop on correspom_topology->nbDomain())ndances
00379   }//loop on joints 
00380 
00381   // Reading global numbering
00382   //
00383   int ncell=(_collection->getMesh())[idomain]->getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
00384   if (ncell>0)
00385   {
00386     int * array=new int[ncell];
00387     int offset=0;
00388     MESSAGE_MED("Reading cell global numbering for mesh "<< idomain);
00389     list<MED_EN::medGeometryElement>::const_iterator iter;
00390     char meshchar[MED_NAME_SIZE+1];
00391     strcpy(meshchar,(_collection->getMesh())[idomain]->getName().c_str());
00392     int nbtypes = (_collection->getMesh())[idomain]->getNumberOfTypes(MED_EN::MED_CELL);
00393     const MED_EN::medGeometryElement* types =(_collection->getMesh())[idomain]->getTypes(MED_EN::MED_CELL);
00394     for (int itype=0; itype<nbtypes;itype++)
00395     {
00396       MED_EN::medGeometryElement type=types[itype];
00397       if (!_collection->isDimensionOK(type,(_collection->getMesh())[idomain]->getMeshDimension()))
00398         continue;
00399       int ntype = (_collection->getMesh())[idomain]->getNumberOfElements(MED_EN::MED_CELL,type);
00400       if (ntype==0) continue;
00401       med_2_3::MEDmeshGlobalNumberRd(fid, meshname,
00402                                      MED_NO_DT, MED_NO_IT,
00403                                      med_2_3::MED_CELL, (med_2_3::med_geometry_type)type,
00404                                      array+offset );
00405       offset+=ntype;
00406     }
00407     cellglobal[idomain]=array;
00408   }
00409 
00410   MESSAGE_MED("Reading node global numbering");
00411   int nnode= (_collection->getMesh())[idomain]->getNumberOfNodes();
00412   {
00413     int* array=new int[nnode];
00414     med_2_3::MEDmeshGlobalNumberRd(fid,meshname,
00415                                    MED_NO_DT, MED_NO_IT,
00416                                    med_2_3::MED_NODE, MED_POINT1,
00417                                    array);
00418     nodeglobal[idomain]=array;
00419   }
00420 
00421   MESSAGE_MED("Reading face global numbering for mesh "<<idomain);
00422   MED_EN::medEntityMesh entity =
00423     (mesh->getMeshDimension()==3)?MED_EN::MED_FACE:MED_EN::MED_EDGE;
00424   int nbface=(_collection->getMesh())[idomain]->getNumberOfElements(entity,MED_EN::MED_ALL_ELEMENTS);
00425   if (nbface!=0)
00426   {
00427     int* array=new int[nbface];
00428     int offset=0;
00429     int nbtypes = mesh->getNumberOfTypes( entity );
00430     const MED_EN::medGeometryElement* types =mesh->getTypes(entity);
00431 
00432     for (int itype=0; itype< nbtypes; itype++)
00433     {
00434       MED_EN::medGeometryElement type = types[itype];
00435       if (!_collection->isDimensionOK(type,mesh->getMeshDimension()-1)) continue;
00436 
00437       int ntype = mesh->getNumberOfElements(entity,type);
00438       if (ntype==0) continue;
00439       med_2_3::MEDmeshGlobalNumberRd( fid, meshname,
00440                                       MED_NO_DT, MED_NO_IT,
00441                                       med_2_3::MED_CELL, (med_2_3::med_geometry_type)type,
00442                                       array+offset );
00443       offset+=ntype;
00444     }
00445     faceglobal[idomain]=array;
00446   }
00447   med_2_3::MEDfileClose(fid);
00448 
00449   END_OF_MED(LOC);
00450 }
00451 
00452 void MESHCollectionDriver::writeSubdomain(int idomain, int nbdomains, char* distfilename,
00453                                           ParaDomainSelector* /*domain_selector*/)
00454 {
00455   MESSAGE_MED(" Number of connect zones "<<(_collection->getCZ()).size());
00456 
00457   //writing connect zones in joints
00458 
00459   med_2_3::med_idt fid = med_2_3::MEDfileOpen(distfilename,med_2_3::MED_ACC_RDWR);
00460 
00461   int index_joint=0;
00462 
00463 
00464   for (unsigned icz=0; icz<(_collection->getCZ()).size(); icz++)
00465   {
00466     if ((_collection->getCZ())[icz]->getLocalDomainNumber()==idomain)
00467     {
00468       med_2_3::med_err error;
00469       int idistant=(_collection->getCZ())[icz]->getDistantDomainNumber();
00470       char joint_name[MED_NAME_SIZE+1];
00471       sprintf(joint_name,"joint_%i",idistant+1);
00472       char desc[MED_COMMENT_SIZE+1];
00473       sprintf(desc,"connect_zone_%d",icz+1);
00474 
00475       char distant_name[MED_NAME_SIZE+1];
00476       //sprintf(distant_name,"domain_%i",(_collection->getCZ())[icz]->getDistantDomainNumber());
00477 
00478       //        sprintf(distant_name,(_collection->getMesh())[idistant]->getName().c_str());
00479       sprintf(distant_name,"domain_%i",idistant);
00480       char mesh_name[MED_NAME_SIZE+1];
00481 
00482       strcpy (mesh_name, _collection->getMesh(idomain)->getName().c_str());
00483       SCRUTE_MED(_collection->getMesh(idomain)->getName());
00484       error = med_2_3::MEDsubdomainJointCr(fid,mesh_name, joint_name, desc, 
00485                                            idistant, distant_name);
00486       if (error==-1) cout << "erreur creation de joint "<<endl;
00487 
00489       //writing node/node correspondency
00491       int nbnodes=(_collection->getCZ())[icz]->getNodeNumber();
00492       int* node_corresp=const_cast<int*>((_collection->getCZ())[icz]->getNodeCorrespValue());
00493 
00494       /* Nodes are reordered so that the ordering on the local and the distant domain
00495          correspond. The chosen order is the natural ordering on the domain
00496          with lowest proc id*/
00497       if (_collection->getSubdomainBoundaryCreates())
00498       {
00499         if (idomain<idistant)
00500           jointSort(node_corresp, nbnodes, true);
00501         else
00502           jointSort(node_corresp, nbnodes, false);
00503       }
00504       error=
00505         med_2_3::MEDsubdomainCorrespondenceWr(fid, mesh_name, joint_name,
00506                                               MED_NO_DT, MED_NO_IT,
00507                                               med_2_3::MED_NODE, MED_POINT1,
00508                                               med_2_3::MED_NODE, MED_POINT1,
00509                                               nbnodes, node_corresp);
00510       if (error==-1) cout << "erreur creation de joint "<<endl;
00511 
00512       //writing cell/cell joint      
00513       writeElementJoint(MED_EN::MED_CELL, icz, idomain, idistant, mesh_name,joint_name,fid);
00514       //writing face/face joint
00515       if (_collection->getSubdomainBoundaryCreates())
00516       {
00517         MED_EN::medEntityMesh constituent_entity =
00518           (_collection->getMeshDimension()==3)?MED_EN::MED_FACE:MED_EN::MED_EDGE;
00519         writeElementJoint(constituent_entity, icz, idomain, idistant, mesh_name,joint_name,fid);                 
00520       }                   
00521       index_joint++;
00522     }
00523   }
00524 
00525   char meshchar[MED_NAME_SIZE+1];
00526   strcpy(meshchar,(_collection->getMesh())[idomain]->getName().c_str());
00527 
00528   // Writing cell global numbering
00529   // 
00530   {
00531     int ncell=_collection->getTopology()->getCellNumber(idomain);
00532     int* array=new int[ncell];
00533     _collection->getTopology()->getCellList(idomain,array);
00534     int offset=0;
00535 
00536     MED_EN::MESH_ENTITIES::const_iterator currentEntity;
00537     list<MED_EN::medGeometryElement>::const_iterator iter;
00538     currentEntity  = MED_EN::meshEntities.find(MED_EN::MED_CELL);
00539 
00540     int nbtypes = (_collection->getMesh())[idomain]->getNumberOfTypes(MED_EN::MED_CELL);
00541     const MED_EN::medGeometryElement* types =(_collection->getMesh())[idomain]->getTypes(MED_EN::MED_CELL);
00542     for (int itype=0; itype<nbtypes;itype++)
00543     {
00544       MED_EN::medGeometryElement type=types[itype];
00545       if (!_collection->isDimensionOK(type,(_collection->getMesh())[idomain]->getMeshDimension()))
00546         continue;
00547       int ntype = (_collection->getMesh())[idomain]->getNumberOfElements(MED_EN::MED_CELL,type);
00548       if (ntype==0) continue;
00549       med_2_3::MEDmeshGlobalNumberWr(fid,meshchar,
00550                                      MED_NO_DT, MED_NO_IT,
00551                                      med_2_3::MED_CELL, (med_2_3::med_geometry_type)type,
00552                                      ntype, array+offset);
00553       offset+=ntype;
00554     }
00555     delete[] array;
00556   }
00557 
00558   MED_EN::medEntityMesh constituent_entity;
00559   if (_collection->getMeshDimension()==3)
00560     constituent_entity=MED_EN::MED_FACE;
00561   else if (_collection->getMeshDimension()==2)
00562     constituent_entity=MED_EN::MED_EDGE;
00563   else throw MEDEXCEPTION("Wrong dimension");
00564 
00565 
00566   //writing face global numbering
00567   {
00568     int offset=0;
00569     int nface= _collection->getTopology()->getFaceNumber(idomain);
00570     int * array=new int[nface];
00571     _collection->getTopology()->getFaceList(idomain,array);
00572 
00573     int nbfacetypes = 0;
00574     const MED_EN::medGeometryElement* facetypes = 0;
00575     if ( _collection->getMesh()[idomain]->getNumberOfElements( constituent_entity, MED_ALL_ELEMENTS))
00576       {
00577         nbfacetypes = (_collection->getMesh())[idomain]->getNumberOfTypes(constituent_entity);
00578         if (nbfacetypes>0)
00579           facetypes =(_collection->getMesh())[idomain]->getTypes(constituent_entity);
00580       }
00581     for (int itype=0; itype<nbfacetypes;itype++)
00582       {
00583         MED_EN::medGeometryElement type=facetypes[itype];
00584         if (!_collection->isDimensionOK(type,_collection->getMeshDimension()-1)) continue;
00585 
00586         int ntype = (_collection->getMesh())[idomain]->getNumberOfElements(constituent_entity,type);
00587         if (ntype==0) continue;
00588         med_2_3::MEDmeshGlobalNumberWr(fid,meshchar,
00589                                        MED_NO_DT, MED_NO_IT,
00590                                        med_2_3::MED_CELL, (med_2_3::med_geometry_type)type,
00591                                        ntype, array+offset );
00592         offset+=ntype;
00593       }
00594     delete[] array;
00595   }
00596 
00597 
00598   //writing node global numbering
00599   {
00600     int nnode=_collection->getTopology()->getNodeNumber(idomain);
00601     int* array = new int[nnode];
00602     _collection->getTopology()->getNodeList(idomain,array);
00603 
00604     med_2_3::MEDmeshGlobalNumberWr(fid,meshchar,
00605                                    MED_NO_DT, MED_NO_IT,
00606                                    med_2_3::MED_NODE, MED_POINT1,
00607                                    nnode, array);
00608     delete[] array;
00609   }
00610 
00611   med_2_3::MEDfileClose(fid);
00612   MESSAGE_MED("End of writing");
00613 
00614 }
00615 
00616 void MESHCollectionDriver::writeElementJoint(medEntityMesh entity ,
00617                                              int icz, 
00618                                              int idomain, 
00619                                              int idistant, 
00620                                              char* mesh_name, 
00621                                              char* joint_name,  
00622                                              med_2_3::med_idt fid )
00623 {
00625   //writing cell/cell correspondency
00627   int nbcells=(_collection->getCZ())[icz]->getEntityCorrespNumber(entity,entity);
00628   const int* index = (_collection->getCZ())[icz]->getEntityCorrespIndex(entity,entity);
00629   const int* value = (_collection->getCZ())[icz]->getEntityCorrespValue(entity,entity);
00630   if ( nbcells==0 ) return; // e.g. domains have 1 common node
00631 
00632   map <pair <MED_EN::medGeometryElement, MED_EN::medGeometryElement> , vector<int> > cellmap;
00633   map <MED_EN::medGeometryElement, int> local_offset;
00634   map <MED_EN::medGeometryElement, int> distant_offset;
00635 
00636   //definition of the local offsets for the types present on local
00637   //and distant domains
00638   // for a mesh containing 2 triangles and 3 quads
00639   //local_offset[TRIA3]=0
00640   //local_offset[QUAD4]=2
00641 
00642   int nb_types_local=(_collection->getMesh())[idomain]-> getNumberOfTypes(entity);
00643   const MED_EN::medGeometryElement* local_types = (_collection->getMesh())[idomain]->getTypes(entity);
00644   const int* local_gni = (_collection->getMesh())[idomain]-> getGlobalNumberingIndex(entity);
00645   for (int i=0; i< nb_types_local; i++)
00646   {
00647     local_offset[local_types[i]]=local_gni[i]-1;
00648   }                                      
00649 
00650   int nb_types_distant=(_collection->getMesh())[idistant]-> getNumberOfTypes(entity);
00651   const MED_EN::medGeometryElement* distant_types = (_collection->getMesh())[idistant]->getTypes(entity);
00652   const int* distant_gni = (_collection->getMesh())[idistant]-> getGlobalNumberingIndex(entity);
00653   for (int i=0; i< nb_types_distant; i++)
00654   {
00655     distant_offset[distant_types[i]]=distant_gni[i]-1;
00656   } 
00657 
00658   if (nb_types_local==1 && nb_types_distant==1)
00659   {
00660     MED_EN::medGeometryElement local_type =  (_collection->getMesh())[idomain]->getElementType(entity,1);
00661     MED_EN::medGeometryElement distant_type = (_collection->getMesh())[idistant]->getElementType(entity,1);
00662     vector<int> corresp;
00663     for (int i=0; i<nbcells; i++)
00664       for (int icol = index[i]-1; icol<index[i+1]-1; icol++)
00665       {
00666         corresp.push_back(i+1);
00667         corresp.push_back(value[icol]);
00668       }
00669     int size_joint = corresp.size()/2;
00670     med_2_3::MEDsubdomainCorrespondenceWr(fid, mesh_name, joint_name,
00671                                           MED_NO_DT, MED_NO_IT,
00672                                           med_2_3::MED_CELL,(med_2_3::med_geometry_type)local_type,
00673                                           med_2_3::MED_CELL,(med_2_3::med_geometry_type)distant_type,
00674                                           size_joint,&corresp[0]);
00675   }
00676   else
00677   {
00678     //classifying all the cell/cell relationships into geomtype/geomtype relationships
00679     //there exists a vector for each geomtype/geomtype pair
00680     // the vectors are stored in cellmap, a std::map with a pair<geomtype,geomtype> key 
00681 
00682     for (int i=0; i<nbcells; i++)
00683       for (int icol = index[i]-1; icol<index[i+1]-1; icol++)
00684       {
00685         MED_EN::medGeometryElement local_type =  (_collection->getMesh())[idomain]->getElementType(entity,i+1);
00686         MED_EN::medGeometryElement distant_type = (_collection->getMesh())[idistant]->getElementType(entity,value[icol]);
00687 
00688         cellmap[make_pair(local_type, distant_type)].push_back(i+1-local_offset[local_type]);
00689         cellmap[make_pair(local_type, distant_type)].push_back(value[icol]-distant_offset[distant_type]);
00690 
00691       }
00692     map <pair <MED_EN::medGeometryElement, MED_EN::medGeometryElement> , vector<int> >::const_iterator iter;
00693 
00694     //going through all the (geom,geom) pairs and writing the joints
00695     for (iter= cellmap.begin(); iter != cellmap.end(); iter++)
00696     {
00697       int size= iter->second.size();
00698       int *corresp = new int[size];
00699       for (int ind=0; ind < size; ind++)
00700         corresp[ind]=(iter->second)[ind];
00701       med_2_3::med_geometry_type local_geo_elem=(med_2_3::med_geometry_type)iter->first.first;
00702       med_2_3::med_geometry_type distant_geo_elem=(med_2_3::med_geometry_type)iter->first.second;
00703       int size_joint=size/2;
00704       //med_2_3::med_err error =
00705       med_2_3::MEDsubdomainCorrespondenceWr(fid, mesh_name, joint_name,
00706                                             MED_NO_DT, MED_NO_IT,
00707                                             med_2_3::MED_CELL, local_geo_elem,
00708                                             med_2_3::MED_CELL, distant_geo_elem,
00709                                             size_joint, corresp);
00710       delete[] corresp;
00711     }
00712     // MED v 2.3.1 returns an error code when
00713     // writing a joint that is already present in the file.
00714     // Also, it returns an error code if a joint
00715     // concerns 3D elements.
00716     // Until these two items are not
00717     // changed, the following line must be commented out
00718 
00719     //if (error==-1) throw MEDEXCEPTION("Error filling joint");
00720 
00721 
00722   }
00723 }
00724 
00725 void MESHCollectionDriver::jointSort(int* elems, int nbelems, bool is_first)
00726 {
00727   //filling an ordered structure with the elem ids
00728   map <int,int> nodemap;
00729   if (is_first)
00730     for (int i=0; i<nbelems; i++) 
00731       nodemap.insert(make_pair(elems[2*i],elems[2*i+1]));
00732 
00733   else
00734     for (int i=0; i<nbelems; i++)
00735       nodemap.insert(make_pair(elems[2*i+1],elems[2*i]));
00736 
00737   int* ptr_elems=elems;      
00738 
00739   //filling the vector in appropriate order 
00740   for (map<int,int>::const_iterator iter=nodemap.begin(); iter!=nodemap.end(); iter++)
00741   {
00742     if (is_first)
00743     {
00744       *ptr_elems++=iter->first;
00745       *ptr_elems++=iter->second;
00746     }
00747     else
00748     {
00749       *ptr_elems++=iter->second;
00750       *ptr_elems++=iter->first;
00751     }
00752   }     
00753 }