Back to index

salome-med  6.5.0
MEDMEM_DriverTools.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00004 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 
00023 #include "MEDMEM_DriverTools.hxx"
00024 #include "MEDMEM_STRING.hxx"
00025 #include "MEDMEM_Exception.hxx"
00026 #include "MEDMEM_Mesh.hxx"
00027 #include "MEDMEM_Group.hxx"
00028 #include "MEDMEM_Field.hxx"
00029 
00030 #include <iomanip>
00031 #include <algorithm>
00032 
00033 using namespace std;
00034 using namespace MED_EN;
00035 
00036 #define DUMP_LINES_LIMIT 20
00037 
00038 namespace MEDMEM {
00039 
00040 // avoid coping sortedNodeIDs
00041 _maille::_maille(const _maille& ma)
00042   : sommets(ma.sommets), geometricType(ma.geometricType), 
00043     reverse(ma.reverse), sortedNodeIDs(0), _ordre(ma._ordre)
00044 {
00045 }
00046 
00047 // Cet opérateur permet d'ordonner les mailles dans un set suivant l'ordre requis par MED
00048 bool _maille::operator < (const _maille& ma) const
00049 {
00050   // si le type géométrique differe, la comparaison est basée dessus
00051   // sinon on se base sur une comparaison des numéros de sommets
00052   // we compare _mailles of only same geometricType due to maillageByType usage
00053   size_t l=sommets.size();
00054   if ( dimension() > 3 ) { // poly
00055     size_t l2 = ma.sommets.size();
00056     if ( l != l2 )
00057       return l < l2;
00058   }
00059   const int* v1 = getSortedNodes();
00060   const int* v2 = ma.getSortedNodes();
00061   for ( const int* vEnd = v1 + l; v1 < vEnd; ++v1, ++v2 )
00062     if(*v1 != *v2)
00063       return *v1 < *v2;
00064   return false; // cas d'égalité
00065 
00066 //   if(geometricType==ma.geometricType)
00067 //   {
00068 //     // construction de deux vecteur temporaire contenant les numeros de sommets
00069 //     // pour faire le tri et les comparaisons
00070 //     size_t l=sommets.size();
00071 //     if ( dimension() > 3 ) { // poly
00072 //       size_t l2 = ma.sommets.size();
00073 //       if ( l != l2 )
00074 //         return l < l2;
00075 //     }
00076 //     std::vector<int> v1(l);
00077 //     std::vector<int> v2(l);
00078 //     for (unsigned int i=0; i!=l; ++i)
00079 //     {
00080 //       v1[i]=nodeNum(i);
00081 //       v2[i]=ma.nodeNum(i);
00082 //     }
00083 //     std::sort(v1.begin(), v1.end());
00084 //     std::sort(v2.begin(), v2.end());
00085 //     for(std::vector<int>::const_iterator i1=v1.begin(), i2=v2.begin(); i1!=v1.end(); ++i1, ++i2)
00086 //       if(*i1 != *i2)
00087 //         return *i1 < *i2;
00088 //     return false; // cas d'égalité
00089 //   }
00090 //   else
00091 //     return geometricType<ma.geometricType;
00092 }
00093 
00094 // creates if needed and return sortedNodeIDs
00095 const int* _maille::getSortedNodes() const
00096 {
00097   if ( !sortedNodeIDs )
00098   {
00099     size_t l=sommets.size();
00100     sortedNodeIDs = new int[ l ];
00101 
00102     for (size_t i=0; i!=l; ++i)
00103       sortedNodeIDs[i]=nodeID(i);
00104     std::sort( sortedNodeIDs, sortedNodeIDs + l );
00105   }
00106   return sortedNodeIDs;
00107 }
00108 
00109 _link _maille::link(int i) const
00110 {
00111   ASSERT_MED ( i >= 0 && i < (int)sommets.size() );
00112   int i2 = ( i + 1 == (int)sommets.size() ) ? 0 : i + 1;
00113   if ( reverse )
00114     return make_pair( sommets[i2]->first, sommets[i]->first );
00115   else
00116     return make_pair( sommets[i]->first, sommets[i2]->first );
00117 }
00118 
00119 // retourne l'entité d'une maille en fonction de la dimension du maillage.
00120 MED_EN::medEntityMesh _maille::getEntity(const int meshDimension) const throw (MEDEXCEPTION)
00121 {
00122   const char * LOC = "_maille::getEntity(const int meshDimension)" ;
00123 //   BEGIN_OF_MED(LOC);
00124 
00125   int mailleDimension = this->dimensionWithPoly();
00126   medEntityMesh entity;
00127   if (mailleDimension == meshDimension)
00128     entity = MED_CELL;
00129   else
00130     switch (mailleDimension)
00131       {
00132       case 0 :
00133         entity = MED_NODE;
00134         break;
00135       case 1 :
00136         entity = MED_EDGE;
00137         break;
00138       case 2 :
00139         entity = MED_FACE;
00140         break;
00141       default :
00142         throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "Impossible de determiner l'entite de la maille."));
00143       }
00144 return entity;
00145 
00146 //END_OF_MED(LOC);
00147 }
00148 
00149 void _maillageByDimIterator::init(const int dim, const bool convertPoly )
00150 {
00151   myIt = myImed->maillageByType.begin();
00152   myEnd = myImed->maillageByType.end();
00153   myDim = dim;
00154   myConvertPoly = convertPoly;
00155   nbRemovedByType = & myImed->nbRemovedByType;
00156 }
00157 
00158 std::ostream& operator << (std::ostream& os, const _maille& ma)
00159 {
00160     os << "maille " << ma.ordre() << " (" << ma.geometricType << ") : < ";
00161     os << ma.nodeNum(0);
00162     for( unsigned i=1; i!=ma.sommets.size(); ++i)
00163         os << ", " << ma.nodeNum( i );
00164     os << " > sortedNodeIDs: ";
00165     if ( ma.sortedNodeIDs ) {
00166       os << "< ";
00167       for( unsigned i=0; i!=ma.sommets.size(); ++i)
00168         os << ( i ? ", " : "" ) << ma.sortedNodeIDs[ i ];
00169       os << " >";
00170     }
00171     else {
00172       os << "NULL";
00173     }
00174     if ( ma.isMerged() )
00175       os << " MERGED ";
00176     return os;
00177 }
00178 
00179 std::ostream& operator << (std::ostream& os, const _groupe& gr)
00180 {
00181     os << "--- Groupe " << gr.nom << " --- " << std::endl ;
00182     os << " -> liste des sous-groupes : ";
00183     for( std::vector<int>::const_iterator i=gr.groupes.begin(); i!=gr.groupes.end(); ++i)
00184             os << *i << " ";
00185 
00186     os << std::endl << " -> liste des "<< gr.mailles.size() << " mailles : " << std::endl;
00187 
00188     _groupe::TMailleIter i1=gr.mailles.begin();
00189     int l;
00190     for(l = 0; l < DUMP_LINES_LIMIT && i1!=gr.mailles.end(); i1++, l++)
00191             os << setw(3) << l+1 << " " << *(*i1) << std::endl;
00192 
00193     if ( l == DUMP_LINES_LIMIT )
00194       os << "   ... skip " << gr.mailles.size() - l << " mailles" << endl;
00195 
00196     os << " relocMap, size=" << gr.relocMap.size() << endl;
00197     map<unsigned,int>::const_iterator it = gr.relocMap.begin();
00198     for ( l = 0; l < DUMP_LINES_LIMIT && it != gr.relocMap.end(); ++it, ++l )
00199       os << " (" << it->first << "," << it->second << ")";
00200     if ( gr.relocMap.size() > 0 )
00201       os << endl;
00202     return os;
00203 }
00204 
00205 std::ostream& operator << (std::ostream& os, const _noeud& no)
00206 {
00207     os << "noeud " << no.number << " : < ";
00208     std::vector<double>::const_iterator i=no.coord.begin();
00209     os << *i++ ;
00210     for( ; i!=no.coord.end(); ++i)
00211         os << ", " << *i;
00212     os << " >";
00213     return os;
00214 }
00215 
00216 void MEDMEM::_fieldBase::dump(std::ostream& os) const
00217 {
00218   os << "field " << "<" << _name << ">" << endl <<
00219     "  nb sub: " << _sub.size() << endl <<
00220     "  group index: " << _group_id << endl <<
00221     "  type: " << _type << endl;
00222   os << "  subcomponents:" << endl;
00223   vector< _sub_data >::const_iterator sub_data = _sub.begin();
00224   for ( ; sub_data != _sub.end(); ++sub_data ) {
00225     os << "    group index: " << sub_data->_supp_id <<
00226       ", " << sub_data->nbComponents() << " comp names: ";
00227     for ( int i_comp = 0; i_comp < sub_data->nbComponents(); ++i_comp )
00228       os << " |" << sub_data->_comp_names[ i_comp ] << "|";
00229     os << endl;
00230   }
00231 }
00232 
00233 std::ostream& operator << (std::ostream& os, const _fieldBase * f)
00234 {
00235   f->dump( os );
00236   return os;
00237 }
00238 
00239 std::ostream& operator << (std::ostream& os, const _intermediateMED& mi)
00240 {
00241   int l;
00242   _maillageByDimIterator maIt( mi );
00243   while ( const set<_maille >* maillage = maIt.nextType() )
00244   {
00245     os << "Set des " << maillage->size()
00246        << " mailles of type " << maIt.type() << ": "<< std::endl;
00247     std::set<_maille>::const_iterator i=maillage->begin();
00248     for( l = 0; l < DUMP_LINES_LIMIT && i!=maillage->end(); ++i, ++l)
00249       os << setw(3) << l+1 << " " <<*i << std::endl;
00250     if ( l == DUMP_LINES_LIMIT )
00251       os << "   ... skip " << maillage->size() - l << " mailles" << endl;
00252   }
00253   os << std::endl << "Vector des " << mi.groupes.size() << " groupes : " << std::endl;
00254   for (unsigned int k=0; k!=mi.groupes.size(); k++)
00255     os << std::setw(3) << k << " " << mi.groupes[k] << std::endl;
00256 
00257   os << std::endl << "map des " << mi.points.size() << " noeuds : " << std::endl;
00258   std::map<int,_noeud>::const_iterator j=mi.points.begin();
00259   for( l = 0; l < DUMP_LINES_LIMIT && j!=mi.points.end(); ++j, ++l)
00260     os << j->second << std::endl;
00261   if ( l == DUMP_LINES_LIMIT )
00262     os << "   ... skip " << mi.points.size() - l << " noeuds" << endl;
00263 
00264   os << endl << mi.fields.size() << " fields:" << endl;
00265   std::list<_fieldBase* >::const_iterator fIt = mi.fields.begin();
00266   for ( l = 0; fIt != mi.fields.end(); ++fIt, ++l )
00267     os << " - " << l+1 << " " << *fIt << endl;
00268 
00269   return os;
00270 }
00271 
00272 
00273 //=======================================================================
00274 //function : treatGroupes
00275 //purpose  : detect groupes of mixed dimension and erase groupes that
00276 //           won't be converted
00277 //=======================================================================
00278 
00279 void _intermediateMED::treatGroupes()
00280 {
00281   const char * LOC = "_intermediateMED::treatGroupes() : ";
00282   BEGIN_OF_MED(LOC);
00283 
00284   if ( myGroupsTreated )
00285     return;
00286   myGroupsTreated = true;
00287 
00288   // --------------------
00289   // erase useless group
00290   // --------------------
00291 
00292   // decrease hierarchical depth of subgroups
00293   vector<int>::iterator j;
00294   for (unsigned int i=0; i!=this->groupes.size(); ++i)
00295   {
00296     _groupe& grp = groupes[i];
00297     //INFOS_MED( i << " " << grp.nom );
00298     j = grp.groupes.begin();
00299     while( j!=grp.groupes.end() ) {
00300       int grpInd = *j-1;
00301       if ( grpInd < 0 || grpInd >= (int)groupes.size() ) {
00302         throw MEDEXCEPTION(LOCALIZED(STRING(LOC) << "Bad subgroup index: " << grpInd <<
00303                                      ", in " << i << " groupe.nom=" << grp.nom));
00304       }
00305       _groupe & sub_grp = groupes[ grpInd ];
00306       if ( !sub_grp.groupes.empty() ) {
00307         MESSAGE_MED("High hierarchical depth of subgroups in group " << i );
00308         *j = sub_grp.groupes[0]; // replace j with its 1st subgroup
00309         // push back the rest subs
00310         for ( int k = 1; k < (int)sub_grp.groupes.size(); ++k )
00311           grp.groupes.push_back( sub_grp.groupes[ k ]);
00312         // vector maybe is reallocated: restart iterator
00313         j = grp.groupes.begin();
00314       }
00315       else
00316         j++;
00317     }
00318     // remove empty sub-groupes
00319     j = grp.groupes.begin();
00320     while ( j!=grp.groupes.end() ) {
00321       if ( groupes[*j-1].empty() ) {
00322         grp.groupes.erase( j );
00323         j = grp.groupes.begin();
00324       }
00325       else
00326         j++;
00327     }
00328   }
00329   // get indices of groups that are field support -
00330   // do not erase them and their subgroups
00331   std::set<int> groups2convert;
00332   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
00333   for ( ; fIt != fields.end(); ++fIt )
00334     (*fIt)->getGroupIds( groups2convert, true );
00335 
00336   // keep named groups and their subgroups
00337   for (unsigned int i=0; i!=this->groupes.size(); ++i)
00338   {
00339     _groupe& grp = groupes[i];
00340     if ( grp.empty() || grp.nom.empty() )
00341       continue;
00342     groups2convert.insert( i );
00343     for( j = grp.groupes.begin(); j!=grp.groupes.end(); ++j)
00344       groups2convert.insert( *j-1 );
00345   }
00346   // erase groups that are not in groups2convert
00347   for (unsigned int i=0; i!=this->groupes.size(); ++i)
00348   {
00349     if ( groups2convert.find( i ) == groups2convert.end() ) {
00350       _groupe& grp = groupes[i];
00351       grp.mailles.clear();
00352       grp.groupes.clear();
00353       MESSAGE_MED( "Erase " << i << "-th group " << grp.nom );
00354     }
00355   }
00356 
00357   // ---------------------------------------------------
00358   // define if there are groups with mixed entity types
00359   // ---------------------------------------------------
00360 
00361   hasMixedCells = false;
00362   for (unsigned int i=0; i!=this->groupes.size(); ++i)
00363   {
00364     _groupe& grp = groupes[i];
00365     if ( grp.groupes.empty() )
00366       continue;
00367 
00368     // check if sub-groups have different dimension
00369     j = grp.groupes.begin();
00370     int dim = groupes[*j-1].mailles[0]->dimension();
00371     for( j++; j!=grp.groupes.end(); ++j) {
00372       int dim2 = groupes[*j-1].mailles[0]->dimension();
00373       if ( dim != dim2 ) {
00374         if ( dim == 0 || dim2 == 0 || dim + dim2 == 9 ) {
00375           // cant create a group of nodes plus anything else,
00376           // neither a group of polygones + polyhedrons
00377           grp.mailles.clear();
00378           grp.groupes.clear();
00379           MESSAGE_MED( "Erase mixed dim group with nodes:" << i << "-th group " << grp.nom );
00380           break;
00381         }
00382         int d1 = dim  < 4 ? dim  : dim  - 2; // care of poly
00383         int d2 = dim2 < 4 ? dim2 : dim2 - 2;
00384         if ( d1 != d2 ) {
00385           hasMixedCells = true;
00386           MESSAGE_MED( "Mixed dim group: " << i << "-th " << grp.nom <<
00387                    "  dim1 = " << dim << " dim2 = " << dim2 );
00388         }
00389       }
00390     }
00391   }
00392 
00393 //   if ( hasMixedCells )
00394 //     INFOS_MED( "There will be groups of mixed dimention" );
00395   END_OF_MED(LOC);
00396 }
00397 
00398 void _intermediateMED::numerotationMaillage()
00399 {
00400   const char* LOC = "_intermediateMED::numerotationMaillage() : ";
00401   BEGIN_OF_MED(LOC);
00402   if ( myMaillesNumerated )
00403     return;
00404   myMaillesNumerated = true;
00405 
00406   treatGroupes();
00407 
00408   set<_maille>::const_iterator i, iEnd;
00409 
00410   // numerotation mailles of type MED_POINT1 by node number
00411   bool hasPointMailles = false;
00412   if ( const set<_maille> * points = _maillageByDimIterator( *this, 0 ).nextType() ) {
00413     hasPointMailles = true;
00414     numerotationPoints();
00415     for ( i = points->begin(), iEnd = points->end(); i != iEnd; ++i )
00416       i->setOrdre( i->sommets[0]->second.number ); // is merged if point is merged
00417   }
00418 
00419   // loop on entities
00420   for ( int dim = 1; dim <= 3; ++dim )
00421   {
00422     int iterDim = hasMixedCells ? -1 : dim;
00423     const bool skipFirstType = ( hasPointMailles && hasMixedCells );
00424 
00425     // check if any numeration is present
00426     bool hasNumeration = true;
00427     _maillageByDimIterator entityMailles( *this, iterDim, true );
00428     if ( skipFirstType ) entityMailles.nextType();
00429     while ( const set<_maille> * typeMailles = entityMailles.nextType() ) {
00430       if ( typeMailles->begin()->ordre() == 0 || typeMailles->rbegin()->ordre() == 0 ) {
00431         hasNumeration = false;
00432         break;
00433       }
00434     }
00435     // check if re-numeration is needed
00436     bool ok = false, renumEntity = false;
00437     if ( hasNumeration )
00438     {
00439       ok = true;
00440       _maillageByDimIterator entityMailles( *this, iterDim, true );
00441       if ( skipFirstType ) entityMailles.nextType();
00442 
00443       int prevNbElems = 0;
00444       while ( const set<_maille> * typeMailles = entityMailles.nextType() )
00445       {
00446         unsigned minOrdre = INT_MAX, maxOrdre = 0;
00447         for ( i = typeMailles->begin(), iEnd = typeMailles->end(); i!=iEnd; ++i) {
00448           if ( i->ordre() < minOrdre ) minOrdre = i->ordre();
00449           if ( i->ordre() > maxOrdre ) maxOrdre = i->ordre();
00450         }
00451         unsigned typeSize = entityMailles.sizeWithoutMerged();
00452         if ( typeSize != maxOrdre - minOrdre + 1 )
00453           ok = false;
00454         if ( prevNbElems != 0 ) {
00455           if ( minOrdre == 1 )
00456             renumEntity = true;
00457           else if ( prevNbElems+1 != (int)minOrdre )
00458             ok = false;
00459         }
00460         prevNbElems += typeSize;
00461       }
00462 
00463       if ( ok && renumEntity ) // each geom type was numerated separately
00464       {
00465         entityMailles.init( iterDim, true );
00466         if ( skipFirstType ) entityMailles.nextType();
00467         prevNbElems = entityMailles.nextType()->size(); // no need to renumber the first type
00468         while ( const set<_maille> * typeMailles = entityMailles.nextType() ) {
00469           for ( i = typeMailles->begin(), iEnd = typeMailles->end(); i!=iEnd; ++i)
00470             i->setOrdre( i->ordre() + prevNbElems );
00471           prevNbElems += typeMailles->size();
00472         }
00473       }
00474     }
00475     if ( !ok )
00476     {
00477       int i_maille=0;
00478       entityMailles.init( iterDim, true );
00479       if ( skipFirstType ) entityMailles.nextType();
00480       while ( const set<_maille> * typeMailles = entityMailles.nextType() )
00481         for ( i = typeMailles->begin(), iEnd = typeMailles->end(); i!=iEnd; ++i)
00482           i->setOrdre( ++i_maille );
00483     }
00484   }
00485   END_OF_MED(LOC);
00486 }
00487 
00488 bool _intermediateMED::numerotationPoints()
00489 {
00490   if ( !myNodesNumerated ) // is negative if numerated by merge
00491   {
00492     int i_noeud=0;
00493     for( map<int,_noeud>::iterator i=points.begin(); i!=points.end(); ++i)
00494       i->second.number = ++i_noeud ;
00495     myNodesNumerated = true;
00496     return true;
00497   }
00498   return false;
00499 }
00500 
00501 int _intermediateMED::nbMerged(int type) const 
00502 {
00503   TNbByType::const_iterator typeNb = nbRemovedByType.find( type );
00504   return ( typeNb == nbRemovedByType.end() ? 0 : typeNb->second );
00505 }
00506 
00507 
00513 COORDINATE * _intermediateMED::getCoordinate(const string & coordinateSystem)
00514 {
00515     const medModeSwitch mode=MED_FULL_INTERLACE;
00516     int spaceDimension=points.begin()->second.coord.size();
00517     int numberOfNodes=points.size() - nbMerged( MED_POINT1 );
00518 
00519     // creation du tableau des coordonnees en mode MED_FULL_INTERLACE
00520     double * coord = new double[spaceDimension*numberOfNodes];
00521     double * xyz = coord;
00522     for( std::map<int,_noeud>::const_iterator i=points.begin(); i!=points.end(); ++i )
00523       if ( i->second.number > 0 ) {
00524         std::copy(i->second.coord.begin(), i->second.coord.end(), xyz );
00525         xyz += spaceDimension;
00526       }
00527 
00528     // creation de l'objet COORDINATE
00529     COORDINATE * coordinate = new COORDINATE(spaceDimension, numberOfNodes, mode);
00530     coordinate->setCoordinates(mode,coord);
00531     delete [] coord;
00532     coordinate->setCoordinatesSystem(coordinateSystem);
00533     return coordinate;
00534 }
00535 
00536 
00542 CONNECTIVITY * _intermediateMED::getConnectivity()
00543 {
00544   const char * LOC = "_intermediateMED::getConnectivity() : ";
00545   BEGIN_OF_MED(LOC);
00546 
00547   int numberOfNodes=points.size() - nbMerged( MED_POINT1 ); // number of nodes in the mesh
00548   medEntityMesh entity;
00549   CONNECTIVITY *Connectivity = NULL, *Constituent = NULL;
00550 
00551   set<_maille>::const_iterator i, iEnd; // iterateurs sur les mailles
00552 
00553   // find out mesh dimension
00554   medGeometryElement meshDim;
00555   _maillageByDimIterator allMailles( *this, -1, true );
00556   while ( allMailles.nextType() )
00557     meshDim = allMailles.dim();
00558 
00559   // renumerote les points de 1 a n (pour le cas ou certains points ne sont pas presents dans le maillage d'origine)
00560   numerotationPoints();
00561 
00562   // loop on entities
00563   for ( int dim = 0; dim <= 3; ++dim )
00564   {
00565     // skip nodes and elements of <dimension_maillage - 2> or less dimension
00566     // Unfortunately, it is impossible because of MESH::createFamilies() that requires
00567     // presence of connectivity even for nodes!
00568     //int dimension_maillage_moin_2=maillage.rbegin()->dimension() - 2;
00569 
00570     // tableau de travail : nombre d'elements pour chaque type geometrique
00571     vector<int> count;
00572     count.reserve( maillageByType.size() );
00573     count.push_back( 1 );
00574 
00575     // tableau de travail : stockage des types geometriques pour UNE entite
00576     vector<medGeometryElement> types;
00577     types.reserve( maillageByType.size() );
00578 
00579     // iterator returning mailles of each type of an entity,
00580     // but if hasMixedCells, we iterate on all types at every dim, since
00581     // in this case we store POINT1 elems as MED_NODE and
00582     // elems of all the rest types as MED_CELL
00583     int iterDim = hasMixedCells ? -1 : dim;
00584     _maillageByDimIterator entityMailles( *this, iterDim, /*convertPoly=*/true );
00585 
00586     // count nb of types and nb mailles of each type
00587     int dimension=0;
00588     if ( dim == 0 ) {
00589       if ( entityMailles.nextType() && entityMailles.dim() == 0 )
00590       {
00591         count.push_back( count.back() + numberOfNodes );
00592         types.push_back( entityMailles.type() );
00593       }
00594     }
00595     else {
00596       while ( entityMailles.nextType() )
00597       {
00598         //if ( entityMailles.dim() > 3 ) break; // ignore poly
00599 
00600         dimension = entityMailles.dim();
00601         if ( dimension == 0 ) continue; // if hasMixedCells, iterator returns all types
00602 
00603         count.push_back( count.back() + entityMailles.sizeWithoutMerged() );
00604         types.push_back( entityMailles.type() );
00605       }
00606     }
00607     int numberOfTypes = types.size(); // nombre de types de l'entite
00608     if ( numberOfTypes == 0 )
00609       continue;
00610 
00611     if ( dimension == meshDim ) entity=MED_CELL;
00612     else if (dimension==2 )     entity=MED_FACE;
00613     else if (dimension==1 )     entity=MED_EDGE;
00614     else if (dimension==0 )     entity=MED_NODE;
00615 
00616     Connectivity = new CONNECTIVITY ( numberOfTypes, entity );
00617     Connectivity->setEntityDimension( dimension );
00618     Connectivity->setNumberOfNodes  ( numberOfNodes );
00619     Connectivity->setGeometricTypes ( &types[0], entity);
00620     Connectivity->setCount          ( &count[0], entity );
00621 
00622     int prevNbElems = 1; // in previous type
00623     for (int k=0; k!=numberOfTypes; ++k )
00624     {
00625       set<_maille> & typeMailles = maillageByType[ types[k] ];
00626       i = typeMailles.begin(), iEnd = typeMailles.end();
00627       int nbMailles = count[k+1]-count[k];
00628       int* connectivity = 0, *index = 0;
00629 
00630       switch ( types[k] )
00631       {
00632       case MED_POLYGON:
00633         {
00634           // put polygones in order of increasing number
00635           vector<const _maille*> orderedPoly( nbMailles );
00636           for ( ; i != iEnd; ++i )
00637             if ( !i->isMerged() )
00638               orderedPoly[ i->ordre() - prevNbElems ] = &(*i);
00639 
00640           // make index
00641           int* polyIndex = index = new int[ nbMailles + 1 ];
00642           vector<const _maille*>::iterator poly = orderedPoly.begin(), polyEnd = orderedPoly.end();
00643           for ( *polyIndex++ = 1; polyIndex < index+nbMailles+1; ++poly, ++polyIndex)
00644             *polyIndex = polyIndex[-1] + (*poly)->sommets.size();
00645 
00646           // make connectivity
00647           int nbNodes = polyIndex[-1];
00648           int* conn = connectivity = new int[ nbNodes ];
00649           for ( poly = orderedPoly.begin(); poly != polyEnd; ++poly) {
00650             for ( int j = 0, nbNodes = (*poly)->sommets.size(); j < nbNodes; ++j )
00651               *conn++ = (*poly)->nodeNum( j );
00652           }
00653           break;
00654         }
00655 
00656       case MED_POLYHEDRA:
00657         {
00658           if ( typeMailles.size() != polyherdalNbFaceNodes.size() )
00659             throw MEDEXCEPTION (LOCALIZED(STRING(LOC) << "Missing info on polyhedron faces"));
00660 
00661           typedef TPolyherdalNbFaceNodes::iterator TPolyFaNoIter;
00662           TPolyFaNoIter polyFaNo, polyFaNoEnd = polyherdalNbFaceNodes.end();
00663 
00664           // put poly's in order of increasing number and count size of connectivity
00665           vector<TPolyFaNoIter> orderedPolyFaNo( nbMailles );
00666           int connSize = 0;
00667           for ( polyFaNo = polyherdalNbFaceNodes.begin(); polyFaNo != polyFaNoEnd; ++polyFaNo )
00668             if ( !polyFaNo->first->isMerged() )
00669             {
00670               orderedPolyFaNo[ polyFaNo->first->ordre() - prevNbElems ] = polyFaNo;
00671               connSize += polyFaNo->first->sommets.size() + polyFaNo->second.size() - 1;
00672             }
00673           vector<TPolyFaNoIter>::iterator pfnIter, pfnEnd = orderedPolyFaNo.end();
00674 
00675           // make index and connectivity
00676           int* conn = connectivity = new int[ connSize ];
00677           int* ind  = index        = new int[ nbMailles+1 ];
00678           *ind++ = 1;
00679           for ( pfnIter = orderedPolyFaNo.begin(); pfnIter != pfnEnd; ++pfnIter)
00680           {
00681             const _maille * poly = (*pfnIter)->first;
00682             const vector<int> & nbFaceNodes = (*pfnIter)->second;
00683             int nbNodes = 0;
00684             for ( unsigned iFace = 0; iFace < nbFaceNodes.size(); ++iFace )
00685             {
00686               for ( int j = 0, nbFNodes = nbFaceNodes[iFace]; j < nbFNodes; ++j )
00687                 *conn++ = poly->nodeNum( nbNodes++ );
00688               *conn++ = -1;
00689             }
00690             conn--;
00691             *ind = ind[-1] + nbNodes;
00692             ++ind;
00693           }
00694           break;
00695         }
00696 
00697       default: // CLASSIC TYPES
00698 
00699         // copie des sommets dans connectivity et set dans Connectivity
00700         int nbSommetsParMaille = i->sommets.size();
00701         int nbSommets = nbMailles * nbSommetsParMaille;
00702         connectivity = new int[ nbSommets ];
00703         if ( entity==MED_NODE )
00704         {
00705           for (int l=0; l!=nbSommets; ++l)
00706             connectivity[l] = l+1;
00707         }
00708         else
00709         {
00710           for ( ; i != iEnd; ++i ) { // loop on elements of geom type
00711             if ( i->isMerged() )
00712               continue;
00713             int* mailleConn = connectivity + nbSommetsParMaille * ( i->ordre() - prevNbElems );
00714             if ( i->reverse )
00715               for ( int n=nbSommetsParMaille-1; n!=-1; --n)
00716                 *mailleConn++ = i->nodeNum( n );
00717             else
00718               for ( int n=0; n != nbSommetsParMaille; ++n)
00719                 *mailleConn++ = i->nodeNum( n );
00720           }
00721           // DO NOT ERASE, maillage will be used while fields construction
00722           //maillage.erase(j);    ; // dangereux, mais optimise la memoire consommee!
00723         }
00724       }
00725 
00726       Connectivity->setNodal (connectivity, entity, types[k], index);
00727       delete [] connectivity;
00728       delete [] index; index = 0;
00729 
00730       prevNbElems += nbMailles;
00731     }
00732 
00733     if ( Constituent )
00734       Connectivity->setConstituent (Constituent);
00735     // stocke Connectivity pour utilisation dans setConstituent lors de la boucle suivante
00736     Constituent = Connectivity;
00737 
00738     if ( entity == MED_CELL )
00739       break; // necessary if hasMixedCells
00740   }
00741 
00742   END_OF_MED(LOC);
00743   return Connectivity;
00744 }
00745 
00746 
00756 void _intermediateMED::getGroups(vector<GROUP *> & _groupCell,
00757                                  vector<GROUP *> & _groupFace,
00758                                  vector<GROUP *> & _groupEdge,
00759                                  vector<GROUP *> & _groupNode, MESH * _ptrMesh)
00760 {
00761   const char* LOC = "_intermediateMED::getGroups() : ";
00762   BEGIN_OF_MED(LOC);
00763 
00764   //medGroupes.resize( groupes.size(), 0 );
00765   if (maillageByType.size() == 0) {
00766     INFOS_MED( "Erreur : il n'y a plus de mailles");
00767     return;
00768   }
00769 
00770   // get indices of groups that are field support - do not skip them
00771   std::set<int> support_groups;
00772   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
00773   for ( ; fIt != fields.end(); ++fIt )
00774     (*fIt)->getGroupIds( support_groups, false );
00775 
00776   numerotationMaillage(); // Renumerotation des mailles par entite
00777 
00778   int dimension_maillage=getMeshDimension();
00779 
00780   for (size_t i=0; i!=this->groupes.size(); ++i)
00781   {
00782     _groupe& grp = groupes[i];
00783     if ( grp.medGroup /*medGroupes[ i ]*/ )
00784       continue; // grp already converted into a MED group
00785 
00786     bool isFieldSupport = ( support_groups.find( i ) != support_groups.end() );
00787     // si le groupe est vide, ou si le groupe n'est pas nommé : on passe au suivant
00788     if ( grp.empty() || ( grp.nom.empty() && !isFieldSupport )) {
00789       if ( !grp.nom.empty() ) {
00790         INFOS_MED("Skip group " << i << " " << grp.nom);
00791 //       } else {
00792 //         MESSAGE_MED("Skip group " << i << " <" << grp.nom << "> " << ( grp.empty() ? ": empty" : ""));
00793       }
00794       continue;
00795     }
00796 
00797     // Build a set of mailles: sort mailles by type and exclude maille doubling
00798     bool isSelfIntersect = false;
00799     typedef set< set<_maille>::iterator, _mailleIteratorCompare > TMailleSet;
00800     TMailleSet mailleSet;
00801     if( grp.groupes.size() ) {// le groupe i contient des sous-maillages
00802       int nb_elem = 0;
00803       for( vector<int>::iterator j=grp.groupes.begin(); j!=grp.groupes.end(); ++j)
00804       {
00805         nb_elem += groupes[*j-1].mailles.size();
00806         _groupe::TMailleIter maIt=groupes[*j-1].mailles.begin();
00807         for( ; maIt!=groupes[*j-1].mailles.end(); ++maIt) {
00808 //           TMailleSet::const_iterator ma_it = mailleSet.find( *maIt );
00809 //           if ( ma_it != mailleSet.end() ) {
00810 //             MESSAGE_MED("EQUAL ELEMS: " << *ma_it << " AND " << *maIt);
00811 //           }
00812 //           else
00813             mailleSet.insert( *maIt );
00814         }
00815       }
00816       if ( nb_elem != (int)mailleSet.size() ) { // Self intersecting compound group
00817         isSelfIntersect = true;
00818         INFOS_MED("Self intersecting group: " << i << " <" << grp.nom << ">"
00819               << ", mailleSet.size = " << mailleSet.size() << ", sum nb elems = " << nb_elem);
00820         for( vector<int>::iterator j=grp.groupes.begin(); j!=grp.groupes.end(); ++j)
00821           INFOS_MED(" in sub-group "<<  *j << " <" << groupes[*j-1].nom << "> "
00822                 << groupes[*j-1].mailles.size() << " mailles of type "
00823                 << groupes[*j-1].mailles[0]->geometricType);
00824         // if a group serve as a support to a field, then the _field is to be converted
00825         // into several MEDMEM::FIELDs, one per sub-group; if we already skipped some sub-groups,
00826         // we roll back the loop on groups to create MED groups for skipped ones
00827         if ( isFieldSupport ) {
00828           if ( grp.nom.empty() ) // clear group existing for the sake of field only
00829             grp.groupes.clear();
00830           std::set<int> sub_grps;
00831           for ( fIt = fields.begin(); fIt != fields.end(); ++fIt ) {
00832             _fieldBase * field = (*fIt);
00833             if ( field->_group_id == (int)i ) {
00834               field->_group_id = -1; // -> a field by support
00835               field->getGroupIds( sub_grps, false );
00836             }
00837           }
00838           if ( (int)i > *sub_grps.begin() ) { // roll back
00839             support_groups.erase( i );
00840             support_groups.insert( sub_grps.begin(), sub_grps.end() );
00841             i = *sub_grps.begin() - 1;
00842             continue;
00843           }
00844         }
00845       }
00846     }
00847     else {
00848       _groupe::TMailleIter maIt=grp.mailles.begin();
00849       for(; maIt!=grp.mailles.end(); ++maIt)
00850         mailleSet.insert( *maIt );
00851       if ( grp.mailles.size() != mailleSet.size() )
00852         INFOS_MED( "Self intersecting group: " << i << " <" << grp.nom << ">"
00853                << ", mailleSet.size = " << mailleSet.size()
00854                << ", nb elems = " << grp.mailles.size());
00855     }
00856     // MEDMEM does not allow constituents of <mesh_dim>-2 and less dimension
00857     // but do not skip node group
00858     int group_min_dim = (**mailleSet.begin()).dimensionWithPoly();
00859     int group_max_dim = (**(--mailleSet.end())).dimensionWithPoly();
00860     if ( group_max_dim != 0 && group_min_dim <= dimension_maillage - 2  ) {
00861       MESSAGE_MED("Skip group: " << i << " <" << grp.nom << "> - too small dimension: "
00862               << group_min_dim);
00863       continue;
00864     }
00865 
00866     // initialise groupe_entity a l'entite de la premiere maille du groupe
00867     medEntityMesh groupe_entity = (**mailleSet.rbegin()).getEntity(dimension_maillage);
00868     if ( hasMixedCells && group_min_dim > 0 )
00869       groupe_entity = MED_CELL;
00870 
00871     // total nb of elements in mesh of groupe_entity
00872     int totalNbElements = 0;
00873     if ( groupe_entity == MED_NODE )
00874       totalNbElements = points.size() - nbMerged( MED_POINT1 );
00875     else {
00876       int entityDim = hasMixedCells ? -1 : group_min_dim;
00877       _maillageByDimIterator allMailles( *this, entityDim, true );
00878       while ( allMailles.nextType() )
00879         if ( allMailles.dim() > 0 )
00880           totalNbElements += allMailles.sizeWithoutMerged();
00881     }
00882     const bool isOnAll = ((int) mailleSet.size() == totalNbElements );
00883 
00884     // if !isOnAll, build a map _maille::ordre() -> index in GROUP.getNumber(MED_ALL_ELEMENTS).
00885     // It is used while fields building.
00886     if ( !isOnAll || isSelfIntersect || isFieldSupport ) {
00887       TMailleSet::iterator maIt = mailleSet.begin();
00888       for ( int iMa = 0; maIt != mailleSet.end(); maIt++ )
00889         grp.relocMap.insert( make_pair( (*maIt)->ordre(), ++iMa ));
00890     }
00891     //Parcours des mailles (a partir de la deuxieme) pour compter les types geometriques
00892     int nb_geometric_types=1;
00893     TMailleSet::iterator j=mailleSet.begin();
00894     medGeometryElement geometrictype=(**j).geometricType;
00895     for ( ++j ; j!=mailleSet.end(); ++j )
00896     {
00897       //Compte nombre de types geometriques
00898       if ( (**j).geometricType != geometrictype ) // si on change de type geometrique
00899       {
00900         nb_geometric_types++;
00901         geometrictype=(**j).geometricType;
00902       }
00903     }
00904 
00905     MED_EN::medGeometryElement * tab_types_geometriques = new MED_EN::medGeometryElement[nb_geometric_types];
00906     int * tab_index_types_geometriques = new int[nb_geometric_types+1];
00907     int * tab_numeros_elements = new int[mailleSet.size()];
00908     int * tab_nombres_elements = new int[nb_geometric_types];
00909 
00910     //Remplit tableaux entree des methodes set
00911     int indice_mailles=0/*, maxOrdre = -1*/;
00912     j=mailleSet.begin();
00913     geometrictype=(**j).geometricType;
00914     tab_index_types_geometriques[0]=1;
00915     int indice_types_geometriques=1;
00916     tab_types_geometriques[0]=geometrictype;
00917     //parcours des mailles du groupe
00918     for (  ; j!=mailleSet.end(); ++j , ++indice_mailles)
00919     {
00920       const _maille& ma = **j;
00921       tab_numeros_elements[indice_mailles]= ma.ordre();
00922 //       if ( maxOrdre < tab_numeros_elements[indice_mailles] )
00923 //         maxOrdre = tab_numeros_elements[indice_mailles];
00924       if (ma.geometricType != geometrictype)
00925       {
00926         tab_index_types_geometriques[indice_types_geometriques]=indice_mailles+1;
00927         geometrictype=ma.geometricType;
00928         tab_types_geometriques[indice_types_geometriques]=geometrictype;
00929         ++indice_types_geometriques;
00930       }
00931     }
00932     tab_index_types_geometriques[indice_types_geometriques]=indice_mailles+1;
00933     for (int k=0; k != nb_geometric_types; ++k)
00934     {
00935       tab_nombres_elements[k] = tab_index_types_geometriques[k+1]-tab_index_types_geometriques[k];
00936     }
00937     //INFOS_MED( "MAX ORDRE in grp " << grp.nom << " entity " << groupe_entity << " : " << maxOrdre);
00938 
00939     //Determination type entite du groupe
00940     vector <GROUP *> * vect_group;
00941     switch ( groupe_entity )
00942     {
00943     case MED_CELL :
00944       vect_group= & _groupCell;
00945       break;
00946     case MED_FACE :
00947       vect_group= & _groupFace;
00948       break;
00949     case MED_EDGE :
00950       vect_group= & _groupEdge;
00951       break;
00952     case MED_NODE :
00953       vect_group= & _groupNode;
00954       break;
00955     }
00956 
00957     //Creation nouveau groupe MED
00958     GROUP * new_group = new GROUP();
00959     grp.medGroup = new_group;
00960     //Appel methodes set
00961     //new_group->setTotalNumberOfElements(mailleSet.size());
00962     new_group->setName(grp.nom);
00963     new_group->setMesh(_ptrMesh);
00964     if ( _ptrMesh )
00965       _ptrMesh->removeReference();
00966     new_group->setNumberOfGeometricType(nb_geometric_types);
00967     new_group->setGeometricType(tab_types_geometriques);
00968     new_group->setNumberOfElements(tab_nombres_elements);
00969     new_group->setNumber(tab_index_types_geometriques,tab_numeros_elements);
00970     new_group->setEntity(groupe_entity);
00971     new_group->setAll( isOnAll );
00972 
00973     vect_group->push_back(new_group);
00974 
00975     // Issue 0021311. Use case: a gibi group has references (recorded in pile 1)
00976     // and several names (pile 27) refer (pile 10) to this group.
00977     // We create a copy of this group per each named reference
00978     for ( unsigned iRef = 0 ; iRef < grp.refNames.size(); ++iRef )
00979       if ( !grp.refNames[ iRef ].empty() )
00980       {
00981         vect_group->push_back( new GROUP( *new_group ));
00982         vect_group->back()->setName( grp.refNames[ iRef ] );
00983       }
00984 
00985     delete [] tab_types_geometriques;
00986     delete [] tab_index_types_geometriques;
00987     delete [] tab_numeros_elements;
00988     delete [] tab_nombres_elements;
00989   }
00990 
00991   END_OF_MED(LOC);
00992 }
00993 
00994 //=======================================================================
00995 //function : getFamilies
00996 //purpose  : create families like MESH::createFamilies() but preserves
00997 //           the order of elements in GROUPs defined by constituent families
00998 //           order. Call it after getGroups()
00999 //=======================================================================
01000 
01001 // void _intermediateMED::getFamilies(std::vector<FAMILY *> & _famCell,
01002 //                                    std::vector<FAMILY *> & _famFace,
01003 //                                    std::vector<FAMILY *> & _famEdge,
01004 //                                    std::vector<FAMILY *> & _famNode, MESH * _ptrMesh)
01005 // {
01006 //   const char * LOC = "_intermediateMED::getFamilies() : ";
01007 //   BEGIN_OF_MED(LOC);
01008 
01009 //   int nbElemFam = 0, nbNodeFam = 0;
01010 //   std::map< GROUP*, vector< FAMILY * > > grpFamsMap;
01011 //   int dimension_maillage=maillage.rbegin()->dimension();
01012 
01013 //   std::set<_maille>::const_iterator i=maillage.begin(); // iterateurs sur les mailles
01014 //   std::set<_maille>::const_iterator j=maillage.begin();
01015 
01016 //   do
01017 //   {
01018 //     // make a family containing mailles shared by the same set of groups
01019 //     std::list<unsigned>&  grList = i->groupes;  // to define the family end
01020 //     int           dimension = i->dimension();        // to define the entity end
01021 //     medGeometryElement type = i->geometricType;
01022 //     medEntityMesh    entity = i->getEntity( dimension_maillage );
01023 
01024 //     std::vector<medGeometryElement> tab_types_geometriques;
01025 //     std::vector<int> tab_index_types_geometriques;
01026 //     std::vector<int> tab_nombres_elements;
01027 //     std::vector<int> tab_numeros_elements;
01028 
01029 //     int iMa = 1, nbtype = 0;
01030 //     tab_types_geometriques.push_back( type );
01031 //     tab_index_types_geometriques.push_back( iMa );
01032 
01033 //     // scan family cells and fill the tab that are needed by the create a MED FAMILY
01034 //     while (i != maillage.end() &&
01035 //            i->groupes == grList &&
01036 //            i->dimension() == dimension)
01037 //     {
01038 //       if (type != i->geometricType) // si changement de type geometrique
01039 //       {
01040 //         tab_index_types_geometriques.push_back(iMa);
01041 //         tab_nombres_elements.push_back(nbtype);
01042 //         nbtype=0;
01043 //         type=i->geometricType;
01044 //         tab_types_geometriques.push_back(type); // stocke le nouveau type geometrique rencontre
01045 //       }
01046 //       ++nbtype;
01047 //       ++iMa;
01048 //       ++i;
01049 //     }
01050 //     tab_index_types_geometriques.push_back(iMa);
01051 //     tab_nombres_elements.push_back(nbtype); // n'a pas été stocké dans la boucle
01052 
01053 //     tab_numeros_elements.resize( iMa - 1 );
01054 //     for ( iMa = 0; j != i; j++, iMa++ )
01055 //       tab_numeros_elements[ iMa ] = j->ordre;
01056 
01057 //     int id = ( entity == MED_NODE ? ++nbNodeFam : -(++nbElemFam) );
01058 
01059 //     ostringstream name;
01060 //     name << "FAM_" << id;
01061 
01062 //     // create a empty MED FAMILY and fill it with the tabs we constructed
01063 //     FAMILY* newFam = new FAMILY();
01064 //     newFam->setTotalNumberOfElements( iMa );
01065 //     newFam->setName( name.str() );
01066 //     newFam->setMesh( _ptrMesh );
01067 //     newFam->setNumberOfGeometricType( tab_types_geometriques.size() );
01068 //     newFam->setGeometricType( &tab_types_geometriques[0] ); // we know the tab is not empy
01069 //     newFam->setNumberOfElements( &tab_nombres_elements[0] );
01070 //     newFam->setNumber( &tab_index_types_geometriques[0], &tab_numeros_elements[0] );
01071 //     newFam->setEntity( entity );
01072 //     newFam->setAll( false );
01073 //     newFam->setIdentifier( id );
01074 //     newFam->setNumberOfGroups( grList.size() );
01075 
01076 //     // Update links between families and groups
01077 //     if ( ! grList.empty() )
01078 //     {
01079 //       std::string * groupNames = new string[ grList.size() ];
01080 //       std::list<unsigned>::iterator g = grList.begin();
01081 //       for ( int i = 0; g != grList.end(); ++g, ++i ) {
01082 //         GROUP * medGROUP = getGroup( *g );
01083 //         groupNames[ i ] = medGROUP->getName();
01084 //         grpFamsMap[ medGROUP ].push_back( newFam );
01085 //       }
01086 //       newFam->setGroupsNames(groupNames);
01087 //     }
01088 //     // store newFam
01089 //     std::vector<FAMILY*>* families = 0;
01090 //     switch ( entity )
01091 //     {
01092 //     case MED_CELL :
01093 //       families = & _famCell; break;
01094 //     case MED_FACE :
01095 //       families = & _famFace; break;
01096 //     case MED_EDGE :
01097 //       families = & _famEdge; break;
01098 //     case MED_NODE :
01099 //       families = & _famNode; break;
01100 //     }
01101 //     if ( families )
01102 //       families->push_back( newFam );
01103 
01104 //   } while ( i != maillage.end() );
01105 
01106 //   // update references in groups
01107 //   std::map< GROUP*, vector< FAMILY * > >::iterator gf = grpFamsMap.begin();
01108 //   for ( ; gf != grpFamsMap.end(); ++gf )
01109 //   {
01110 //     gf->first->setNumberOfFamilies( gf->second.size() );
01111 //     gf->first->setFamilies( gf->second );
01112 //   }
01113 //}
01114 
01115 //=======================================================================
01116 //function : getGroup
01117 //purpose  :
01118 //=======================================================================
01119 
01120 // GROUP * _intermediateMED::getGroup( int i )
01121 // {
01122 //   if ( i <(int) medGroupes.size() )
01123 //     return medGroupes[ i ];
01124 //   throw MEDEXCEPTION
01125 //     (LOCALIZED(STRING("_intermediateMED::getGroup(): WRONG GROUP INDEX: ")
01126 //                << medGroupes.size() << " <= " << i ));
01127 // }
01128 
01129 //=======================================================================
01130 //function : getFields
01131 //purpose  :
01132 //=======================================================================
01133 
01134 void _intermediateMED::getFields(std::list< FIELD_* >& theFields)
01135 {
01136   const char * LOC = "_intermediateMED::getFields() : ";
01137   BEGIN_OF_MED(LOC);
01138 
01139   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
01140   for ( ; fIt != fields.end(); fIt++ )
01141   {
01142     const _fieldBase* fb = *fIt;
01143     list<pair< FIELD_*, int> >  ff = fb->getField(groupes);
01144     list<pair< FIELD_*, int> >::iterator f_sup = ff.begin();
01145     for (int j = 1 ; f_sup != ff.end(); f_sup++, ++j )
01146     {
01147       FIELD_    * f = f_sup->first;
01148       SUPPORT * sup = groupes[ f_sup->second ].medGroup;
01149       if ( !sup )
01150         throw MEDEXCEPTION
01151           (LOCALIZED(STRING(LOC) <<"_intermediateMED::getFields(), NULL field support: "
01152                      << " group index: " << fb->_group_id));
01153       int nb_elems = sup->getNumberOfElements( MED_ALL_ELEMENTS );
01154       if ( nb_elems != f->getNumberOfValues() )
01155         throw MEDEXCEPTION
01156           (LOCALIZED(STRING("_intermediateMED::getFields(), field support size (")
01157                      << nb_elems  << ") != NumberOfValues (" << f->getNumberOfValues()));
01158       theFields.push_back( f );
01159       if ( sup->getName().empty() ) {
01160         ostringstream name;
01161         name << "GRP_" << f->getName() << "_" << j;
01162         sup->setName( name.str() );
01163       }
01164       f->setSupport( sup );
01165       //f->setIterationNumber( j );
01166       f->setOrderNumber( j );
01167     }
01168   }
01169   END_OF_MED(LOC);
01170 }
01171 
01172 //=======================================================================
01173 //function : ~_intermediateMED
01174 //purpose  : destructor
01175 //=======================================================================
01176 
01177 _intermediateMED::~_intermediateMED()
01178 {
01179   MESSAGE_MED( "~_intermediateMED()");
01180   std::list< _fieldBase* >::const_iterator fIt = fields.begin();
01181   for ( ; fIt != fields.end(); fIt++ )
01182     delete *fIt;
01183 }
01184 
01185 //=======================================================================
01186 //function : getGroupIds
01187 //purpose  : return ids of main and/or sub- groups
01188 //=======================================================================
01189 
01190 void _fieldBase::getGroupIds( std::set<int> & ids, bool all ) const
01191 {
01192   if ( hasCommonSupport() )
01193     ids.insert( _group_id );
01194   if ( all || !hasCommonSupport() ) {
01195     vector< _sub_data >::const_iterator sub = _sub.begin();
01196     for ( ; sub != _sub.end(); ++sub )
01197       ids.insert( sub->_supp_id );
01198   }
01199 }
01200 
01201 //=======================================================================
01202 //function : hasSameComponentsBySupport
01203 //purpose  :
01204 //=======================================================================
01205 
01206 bool _fieldBase::hasSameComponentsBySupport() const
01207 {
01208   vector< _sub_data >::const_iterator sub_data = _sub.begin();
01209   const _sub_data& first_sub_data = *sub_data;
01210   for ( ++sub_data ; sub_data != _sub.end(); ++sub_data )
01211   {
01212     if ( first_sub_data._comp_names != sub_data->_comp_names )
01213       return false; // diff names of components
01214 
01215     if ( first_sub_data._nb_gauss != sub_data->_nb_gauss )
01216       return false; // diff nb of gauss points
01217   }
01218   return true;
01219 }
01220 
01221 //=======================================================================
01222 //
01223 //function : mergeNodes
01224 //
01225 //=======================================================================
01226 namespace {
01227 
01228 struct __NOEUD
01229 {
01230   typedef map<int,_noeud>::iterator Inoeud;
01231   Inoeud i_noeud;
01232 
01233   __NOEUD() {}
01234   __NOEUD(Inoeud n): i_noeud(n) {}
01235   const double & operator[] (int i) const { return i_noeud->second.coord[i]; }
01236   double         operator[] (int i)       { return i_noeud->second.coord[i]; }
01237   int dim() const { return i_noeud->second.coord.size(); }
01238   int& num() { return i_noeud->second.number; }
01239   int id() const { return i_noeud->first; }
01240   bool isMerged() const { return i_noeud->second.number < 1; }
01241 };
01242 //-----------------------------------------------------------------------
01243 double DistanceL2(const __NOEUD &A,const __NOEUD &B)
01244 {
01245   if ( B.isMerged() ) return DBL_MAX;
01246   const double* cooA = &A[0], *cooB = &B[0], *cooEnd = cooA + A.dim();
01247   double dist, somme=0;
01248   while ( cooA < cooEnd ) {
01249     dist=((*cooA++) - (*cooB++));
01250     somme+=dist*dist;
01251   }
01252   return sqrt(somme);
01253 }
01254 //-----------------------------------------------------------------------
01255 struct __NUAGENOEUD
01256 {
01257   vector< __NOEUD > nodes;
01258   __NUAGENOEUD(_intermediateMED& imed);
01259   __NOEUD & operator [] (int i) { return nodes[i]; }
01260   int size() const { return nodes.size(); }
01261 };
01262 //-----------------------------------------------------------------------
01263 __NUAGENOEUD::__NUAGENOEUD(_intermediateMED& imed)
01264 {
01265   nodes.resize( imed.points.size() );
01266   map<int,_noeud>::iterator i_noeud = imed.points.begin();
01267   for( int i = 0; i_noeud!=imed.points.end(); ++i_noeud, ++i ) {
01268     i_noeud->second.number = i+1;
01269     nodes[ i ] = i_noeud;
01270   }
01271 }
01272 //-----------------------------------------------------------------------
01273 // mergeNodes
01274 template<int DIM> int mergeNodes(double            tolerance,
01275                                  _intermediateMED& imed,
01276                                  vector< int > &   /*newNodeIDs*/)
01277 {
01278   /*typedef dTree<__NOEUD,__NUAGENOEUD,DIM > DTree;
01279   __NUAGENOEUD aNUAGENOEUD( imed );
01280   DTree tree( &aNUAGENOEUD );
01281 
01282 //   int maxID = imed.points.rbegin()->first;
01283 //   newNodeIDs.resize( maxID + 1, 0 );
01284 
01285   int num = 0, nbRemoved = 0;
01286   int nbNodes = aNUAGENOEUD.size();
01287   for ( int i = 0; i < nbNodes; ++i )
01288   {
01289     __NOEUD & node = aNUAGENOEUD[i];
01290     int & number = node.num();
01291     if ( number < 0 )
01292       continue; // already merged
01293     number = ++num;
01294 
01295     list<int> closeNumbers;
01296     int nbCoinsident = tree.get_all_close( node, tolerance, closeNumbers );
01297     if ( nbCoinsident > 1 ) // not only node it-self found
01298     {
01299       nbRemoved += nbCoinsident-1; // one of coincident nodes remains
01300       int id = node.id();
01301       list<int>::iterator n = closeNumbers.begin(), nEnd = closeNumbers.end();
01302       while ( n != nEnd ) {
01303         __NOEUD & coincNode = aNUAGENOEUD[ *n++ ];
01304         int coincID = coincNode.id();
01305         if ( coincID != id ) {
01306           coincNode.num() = -number;
01307           //newNodeIDs[ coincID ] = id;
01308         }
01309       }
01310     }
01311   }
01312   return nbRemoved;*/
01313   return 0;
01314 }
01315 //-----------------------------------------------------------------------
01316 // wrapper of _maille used after merging nodes to find equal mailles
01317 struct _mailleWrap
01318 {
01319   const _maille* _ma;
01320 
01321   _mailleWrap(const _maille* ma=0): _ma(ma) { if (_ma) _ma->init(); }
01322   ~_mailleWrap()                            { if (_ma) _ma->init(); }
01323 
01324   bool operator < (const _mailleWrap& mw) const
01325   {
01326     size_t l=_ma->sommets.size();
01327     const int* v1 = getSortedNodeNums(_ma);
01328     const int* v2 = getSortedNodeNums(mw._ma);
01329     for ( const int* vEnd = v1 + l; v1 < vEnd; ++v1, ++v2 )
01330       if(*v1 != *v2)
01331         return *v1 < *v2;
01332     return false; // cas d'égalité
01333   }
01334   static const int* getSortedNodeNums(const _maille* ma)
01335   {
01336     if ( !ma->sortedNodeIDs ) {
01337       size_t l=ma->sommets.size();
01338       ma->sortedNodeIDs = new int[ l ];
01339       for (size_t i=0; i!=l; ++i)
01340         ma->sortedNodeIDs[i]=ma->nodeNum(i);
01341       std::sort( ma->sortedNodeIDs, ma->sortedNodeIDs + l );
01342     }
01343     return ma->sortedNodeIDs;
01344   }
01345 };
01346 
01347 } // namespace
01348 
01349 //=======================================================================
01350 //function : mergeNodesAndElements
01351 //purpose  : optionally merge nodes and elements
01352 //=======================================================================
01353 
01354 void _intermediateMED::mergeNodesAndElements(double tolerance)
01355 {
01356   //const char * LOC = "_intermediateMED::mergeNodesAndElements(): ";
01357   vector< int > newNodeIDs; // for i-th node id, id to replace with
01358 
01359   int nbRemovedNodes = 0;
01360   const int spaceDimension=points.begin()->second.coord.size();
01361   if ( spaceDimension == 3 )
01362     nbRemovedNodes = mergeNodes<3>( tolerance, *this, newNodeIDs );
01363   else if ( spaceDimension == 2 )
01364     nbRemovedNodes = mergeNodes<2>( tolerance, *this, newNodeIDs );
01365 
01366   myNodesNumerated = true;
01367 
01368   if ( nbRemovedNodes == 0 )
01369     return;
01370 
01371   // all commented code here was used to keep initial numeration but it was too slow
01372   //numerotationMaillage();
01373   treatGroupes();
01374 
01375   nbRemovedByType[ MED_NONE   ] = nbRemovedNodes;
01376   nbRemovedByType[ MED_POINT1 ] = nbRemovedNodes;
01377 
01378   bool hasPointMailles = false;
01379   _maillageByDimIterator entityMailles( *this, 0 );
01380   if ( const set<_maille> * points = entityMailles.nextType() ) {
01381     hasPointMailles = true;
01382     set<_maille>::const_iterator i, iEnd = points->end();
01383     for ( i = points->begin(); i != iEnd; ++i )
01384       i->setOrdre( i->sommets[0]->second.number ); // is merged if point is merged
01385   }
01386   const bool skipFirstType = ( hasPointMailles && hasMixedCells );
01387   // loop on entities
01388   for ( int dim = 1; dim <= 3; ++dim )
01389   {
01390     int iterDim = hasMixedCells ? -1 : dim;
01391     //int nbRemovedInEntity = 0;
01392 
01393     // count total nb of mailles in entity
01394 //     int nbMailles = 0;
01395 //     entityMailles.init( iterDim, true );
01396 //     if ( skipFirstType ) entityMailles.nextType(); // merged POINT1's are same as nodes
01397 //     while ( const set<_maille> * typeMailles = entityMailles.nextType() )
01398 //       nbMailles += typeMailles->size();
01399 
01400     // for each maille number, count shift due to removed mailles with lower numbers
01401     //vector<int> shift( nbMailles+1, 0 );
01402 
01403     // merge and numerate mailles
01404     int num = 0;
01405     entityMailles.init( iterDim, true );
01406     if ( skipFirstType ) entityMailles.nextType(); // merged POINT1's are same as nodes
01407     while ( const set<_maille> * typeMailles = entityMailles.nextType() )
01408     {
01409       int nbRemoved = 0;
01410       set<_mailleWrap> maNodeNumSorted;
01411       pair< set<_mailleWrap>::const_iterator, bool > mw_isunique;
01412 
01413       set<_maille>::const_iterator ma = typeMailles->begin(), maEnd = typeMailles->end();
01414       while ( ma != maEnd )
01415       {
01416         const _maille & m = *ma++;
01417         mw_isunique = maNodeNumSorted.insert( &m );
01418         if ( mw_isunique.second ) {
01419           m.setOrdre( ++num );
01420         }
01421         else {
01422           const _maille* equalMa = mw_isunique.first->_ma;
01423           //unsigned ordreToRemove = m.ordre();
01424           m.setMergedOrdre( equalMa->ordre() );
01425           nbRemoved++;
01426 //           while ( ordreToRemove <= nbMailles )
01427 //             shift[ ordreToRemove++ ]++;
01428         }
01429       }
01430 
01431       if ( nbRemoved ) {
01432         nbRemovedByType[ entityMailles.type() ] = nbRemoved;
01433         //nbRemovedInEntity += nbRemoved;
01434       }
01435 
01436       // update maille ordre
01437 //       if ( nbRemovedInEntity ) {
01438 //         for ( ma = typeMailles->begin(); ma != maEnd; ++ma ) {
01439 //           unsigned newOrdre = ma->ordre() - shift[ ma->ordre() ];
01440 //           if ( ma->isMerged() )
01441 //             ma->setMergedOrdre( newOrdre );
01442 //           else
01443 //             ma->setOrdre( newOrdre );
01444 //         }
01445 //       }
01446     }
01447   }
01448 
01449   myMaillesNumerated = true;
01450 
01451 }
01452 
01453 //=======================================================================
01454 //function : getMeshDimension
01455 //purpose  : count mesh dimension
01456 //=======================================================================
01457 
01458 int _intermediateMED::getMeshDimension() const
01459 {
01460   int dim = 0;
01461   _maillageByDimIterator allMailles( *this, -1, true );
01462   while ( allMailles.nextType() )
01463     dim = allMailles.dim();
01464   return dim;
01465 }
01466 
01467 }
01468