1 // Copyright (C) 2007-2012 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
27 #include <libxml/tree.h>
28 #include <libxml/parser.h>
29 #include <libxml/xpath.h>
30 #include <libxml/xpathInternals.h>
38 #include "MEDMEM_Utilities.hxx"
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"
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"
60 using namespace MEDSPLITTER;
63 //#include "MEDSPLITTER_MESHCollectionDriver.H"
66 MESHCollectionDriver::MESHCollectionDriver(MESHCollection* collection):_collection(collection)
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
76 int MESHCollectionDriver::readSeq(char* filename, char* meshname)
78 const char* LOC = "MEDSPLITTER::MESHCollectionDriver::readSeq()";
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
92 (_collection->getMesh()).push_back(new_mesh);
94 _collection->setName(meshname);
95 (_collection->getCZ()).clear();
96 vector<int*> cellglobal,nodeglobal,faceglobal;
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);
113 * Reads the file structure to determine the list
114 * of all the available fields
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
123 void MESHCollectionDriver::readFileStruct(vector <string>& field_names,vector<int>& iternumber,vector <int>& ordernumber, vector <int>& types)
125 const char* LOC = "MEDSPLITTER::MESHCollectionDriver::readFileStruct()";
128 const MEDMEM::MEDFILEBROWSER med_struct (_filename[0]);
129 int nb_fields = med_struct.getNumberOfFields();
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++)
135 MEDMEM::VEC_DT_IT_ dtit=med_struct.getFieldIteration(names[ifield]);
136 for (unsigned i = 0; i < dtit.size(); i++)
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] ));
147 //!retrieves the type of a field for a given fieldname
148 int MESHCollectionDriver::getFieldType(const string& fieldname)
150 const char* LOC = "MEDSPLITTER::MESHCollectionDriver::getFieldType()";
152 const MEDMEM::MEDFILEBROWSER med_struct (_filename[0]);
154 int ret = ( MED_REEL64 == med_struct.getFieldType( fieldname ));
161 //================================================================================
163 * \brief Structure used in the method below
167 char _name[MED_NAME_SIZE+1];
169 med_2_3::med_int _distant_domain;
170 med_2_3::med_geometry_type _geo_local, _geo_dist;
174 //================================================================================
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
182 //================================================================================
184 void MESHCollectionDriver::readLoc2GlobCellConnect(int idomain,
185 const set<int>& loc_domains,
186 ParaDomainSelector* domain_selector,
187 vector<int> & loc2glob_corr)
189 using namespace med_2_3;
190 med_2_3::med_err err;
192 // find joints of domains touching idomain and loaded on other processors
195 list< TJointData > joints;
196 int total_nb_couples = 0;
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;
203 const int loc_mesh_dim = loc_mesh->getMeshDimension();
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 ];
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++)
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
222 for (int itype=0; itype<nbtypes;itype++)
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 )
228 if ( !_collection->isDimensionOK( *dist_type, loc_mesh_dim ))
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 )
238 joints.push_back( joint );
239 total_nb_couples += joint._nb_couples;
245 // read cell pairs of found joints and replace distant local ids with global ones
247 loc2glob_corr.resize( 2 * total_nb_couples );
248 if ( total_nb_couples > 0 )
250 int* cell_corresp = & loc2glob_corr[0];
252 list< TJointData >::iterator jnt = joints.begin();
253 for ( ; jnt != joints.end(); ++jnt )
256 err = MEDsubdomainCorrespondenceRd(fid, meshname, jnt->_name,
258 med_2_3::MED_CELL, jnt->_geo_local,
259 med_2_3::MED_CELL, jnt->_geo_dist,
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;
268 cell_corresp += 2 * jnt->_nb_couples;
274 //================================================================================
276 * \brief Return mesh dimension from the distributed med file having been read
278 //================================================================================
280 int MESHCollectionDriver::readMeshDimension() const
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");
286 med_2_3::med_idt fid = med_2_3::MEDfileOpen(_filename[0].c_str(),med_2_3::MED_ACC_RDONLY);
288 throw MED_EXCEPTION( MEDMEM::STRING(LOC) << "can't open file " << _filename[0]);
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;
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);
305 med_2_3::MEDfileClose(fid);
308 throw MED_EXCEPTION( MEDMEM::STRING(LOC) << "mesh name is invalid");
310 return meshDimension;
313 void MESHCollectionDriver::readSubdomain(vector<int*>& cellglobal,
314 vector<int*>& faceglobal,
315 vector<int*>& nodeglobal, int idomain)
317 const char* LOC = "MEDSPLITTER::MESHCollectionDriver::readSubdomain()";
320 char meshname[MED_NAME_SIZE+1];
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
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++)
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,
345 &distant, name_distant,&nstep,&nocstpncorr);
346 for (int ic=0; ic<ncorr; ic++)
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;
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
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)
361 med_2_3::MEDsubdomainCorrespondenceRd(fid, meshname, name,
363 cor_typent_local, cor_typgeo_local,
364 cor_typent_dist, cor_typgeo_dist,
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);
378 }//loop on correspom_topology->nbDomain())ndances
381 // Reading global numbering
383 int ncell=(_collection->getMesh())[idomain]->getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
386 int * array=new int[ncell];
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++)
396 MED_EN::medGeometryElement type=types[itype];
397 if (!_collection->isDimensionOK(type,(_collection->getMesh())[idomain]->getMeshDimension()))
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,
407 cellglobal[idomain]=array;
410 MESSAGE_MED("Reading node global numbering");
411 int nnode= (_collection->getMesh())[idomain]->getNumberOfNodes();
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,
418 nodeglobal[idomain]=array;
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);
427 int* array=new int[nbface];
429 int nbtypes = mesh->getNumberOfTypes( entity );
430 const MED_EN::medGeometryElement* types =mesh->getTypes(entity);
432 for (int itype=0; itype< nbtypes; itype++)
434 MED_EN::medGeometryElement type = types[itype];
435 if (!_collection->isDimensionOK(type,mesh->getMeshDimension()-1)) continue;
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,
445 faceglobal[idomain]=array;
447 med_2_3::MEDfileClose(fid);
452 void MESHCollectionDriver::writeSubdomain(int idomain, int nbdomains, char* distfilename,
453 ParaDomainSelector* /*domain_selector*/)
455 MESSAGE_MED(" Number of connect zones "<<(_collection->getCZ()).size());
457 //writing connect zones in joints
459 med_2_3::med_idt fid = med_2_3::MEDfileOpen(distfilename,med_2_3::MED_ACC_RDWR);
464 for (unsigned icz=0; icz<(_collection->getCZ()).size(); icz++)
466 if ((_collection->getCZ())[icz]->getLocalDomainNumber()==idomain)
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);
475 char distant_name[MED_NAME_SIZE+1];
476 //sprintf(distant_name,"domain_%i",(_collection->getCZ())[icz]->getDistantDomainNumber());
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];
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;
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());
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())
499 if (idomain<idistant)
500 jointSort(node_corresp, nbnodes, true);
502 jointSort(node_corresp, nbnodes, false);
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;
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())
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);
525 char meshchar[MED_NAME_SIZE+1];
526 strcpy(meshchar,(_collection->getMesh())[idomain]->getName().c_str());
528 // Writing cell global numbering
531 int ncell=_collection->getTopology()->getCellNumber(idomain);
532 int* array=new int[ncell];
533 _collection->getTopology()->getCellList(idomain,array);
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);
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++)
544 MED_EN::medGeometryElement type=types[itype];
545 if (!_collection->isDimensionOK(type,(_collection->getMesh())[idomain]->getMeshDimension()))
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);
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");
566 //writing face global numbering
569 int nface= _collection->getTopology()->getFaceNumber(idomain);
570 int * array=new int[nface];
571 _collection->getTopology()->getFaceList(idomain,array);
574 const MED_EN::medGeometryElement* facetypes = 0;
575 if ( _collection->getMesh()[idomain]->getNumberOfElements( constituent_entity, MED_ALL_ELEMENTS))
577 nbfacetypes = (_collection->getMesh())[idomain]->getNumberOfTypes(constituent_entity);
579 facetypes =(_collection->getMesh())[idomain]->getTypes(constituent_entity);
581 for (int itype=0; itype<nbfacetypes;itype++)
583 MED_EN::medGeometryElement type=facetypes[itype];
584 if (!_collection->isDimensionOK(type,_collection->getMeshDimension()-1)) continue;
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 );
598 //writing node global numbering
600 int nnode=_collection->getTopology()->getNodeNumber(idomain);
601 int* array = new int[nnode];
602 _collection->getTopology()->getNodeList(idomain,array);
604 med_2_3::MEDmeshGlobalNumberWr(fid,meshchar,
605 MED_NO_DT, MED_NO_IT,
606 med_2_3::MED_NODE, MED_POINT1,
611 med_2_3::MEDfileClose(fid);
612 MESSAGE_MED("End of writing");
616 void MESHCollectionDriver::writeElementJoint(medEntityMesh entity ,
622 med_2_3::med_idt fid )
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
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;
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
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++)
647 local_offset[local_types[i]]=local_gni[i]-1;
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++)
655 distant_offset[distant_types[i]]=distant_gni[i]-1;
658 if (nb_types_local==1 && nb_types_distant==1)
660 MED_EN::medGeometryElement local_type = (_collection->getMesh())[idomain]->getElementType(entity,1);
661 MED_EN::medGeometryElement distant_type = (_collection->getMesh())[idistant]->getElementType(entity,1);
663 for (int i=0; i<nbcells; i++)
664 for (int icol = index[i]-1; icol<index[i+1]-1; icol++)
666 corresp.push_back(i+1);
667 corresp.push_back(value[icol]);
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]);
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
682 for (int i=0; i<nbcells; i++)
683 for (int icol = index[i]-1; icol<index[i+1]-1; icol++)
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]);
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]);
692 map <pair <MED_EN::medGeometryElement, MED_EN::medGeometryElement> , vector<int> >::const_iterator iter;
694 //going through all the (geom,geom) pairs and writing the joints
695 for (iter= cellmap.begin(); iter != cellmap.end(); iter++)
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);
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
719 //if (error==-1) throw MEDEXCEPTION("Error filling joint");
725 void MESHCollectionDriver::jointSort(int* elems, int nbelems, bool is_first)
727 //filling an ordered structure with the elem ids
728 map <int,int> nodemap;
730 for (int i=0; i<nbelems; i++)
731 nodemap.insert(make_pair(elems[2*i],elems[2*i+1]));
734 for (int i=0; i<nbelems; i++)
735 nodemap.insert(make_pair(elems[2*i+1],elems[2*i]));
737 int* ptr_elems=elems;
739 //filling the vector in appropriate order
740 for (map<int,int>::const_iterator iter=nodemap.begin(); iter!=nodemap.end(); iter++)
744 *ptr_elems++=iter->first;
745 *ptr_elems++=iter->second;
749 *ptr_elems++=iter->second;
750 *ptr_elems++=iter->first;