src/MEDMEM/MEDMEM_DriverTools.hxx

   1 // Copyright (C) 2005  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   2 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   3 //
   4 // This library is free software; you can redistribute it and/or
   5 // modify it under the terms of the GNU Lesser General Public
   6 // License as published by the Free Software Foundation; either
   7 // version 2.1 of the License.
   8 //
   9 // This library is distributed in the hope that it will be useful
  10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 // Lesser General Public License for more details.
  13 //
  14 // You should have received a copy of the GNU Lesser General Public
  15 // License along with this library; if not, write to the Free Software
  16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 //
  18 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19 //
  20 #ifndef DRIVERTOOLS_HXX
  21 #define DRIVERTOOLS_HXX
  22
  23 #include <MEDMEM.hxx>
  24
  25 #include "MEDMEM_define.hxx"
  26 #include "MEDMEM_Exception.hxx"
  27 #include "MEDMEM_FieldForward.hxx"
  28 #include "MEDMEM_ArrayInterface.hxx"
  29 #include "MEDMEM_Group.hxx"
  30 #include "MEDMEM_GaussLocalization.hxx"
  31
  32 #include <string>
  33 #include <vector>
  34 #include <set>
  35 #include <list>
  36 #include <map>
  37 #include <iostream>
  38 #include <iomanip>
  39
  40 // Way of storing values with gauss points by CASTEM: full/no interlace.
  41 // Remove code switched by this symbol as soon as the way is found out,
  42 // from here and from gibi driver
  43 #define CASTEM_FULL_INTERLACE
  44
  45
  46 namespace MEDMEM {
  47 class MESH;
  48 class CONNECTIVITY;
  49 class COORDINATE;
  50 class GROUP;
  51 class FAMILY;
  52 class FIELD_;
  53
  54 struct MEDMEM_EXPORT _noeud
  55 {
  56     mutable int number;
  57     std::vector<double> coord;
  58 };
  59
  60 typedef pair<int,int> _link; // a pair of node numbers
  61
  62 struct MEDMEM_EXPORT _maille
  63 {
  64     typedef std::map<int,_noeud>::iterator iter;
  65     MED_EN::medGeometryElement geometricType;
  66     std::vector< iter > sommets;
  67     mutable unsigned ordre; // l'ordre est fixé après insertion dans le set, et ne change ni l'état, ni l'ordre -> mutable
  68     mutable bool reverse; // to reverse sommets of a face
  69     mutable list<unsigned> groupes; // the GROUPs maille belongs to, used to create families
  70
  71     _maille() : geometricType(MED_EN::MED_NONE),ordre(0),reverse(false)
  72     {
  73     };
  74     _maille(MED_EN::medGeometryElement _geometricType, size_t nelem) : geometricType(_geometricType),ordre(0),reverse(false)
  75     {
  76         sommets.reserve(nelem);
  77     };
  78     int dimension() const // retourne la dimension de la maille
  79     {
  80         return geometricType/100;
  81     };
  82     bool operator < (const _maille& ma) const;
  83     MED_EN::medEntityMesh getEntity(const int meshDimension) const throw (MEDEXCEPTION);
  84    _link link(int i) const;
  85 };
  86
  87 struct MEDMEM_EXPORT _mailleIteratorCompare // pour ordonner le set d'iterateurs sur mailles
  88 {
  89     bool operator () (std::set<_maille>::iterator i1, std::set<_maille>::iterator i2)
  90     {
  91         return *i1<*i2;
  92     }
  93 };
  94
  95 struct MEDMEM_EXPORT _groupe
  96 {
  97     typedef std::vector< std::set<_maille>::iterator>::const_iterator mailleIter;
  98     std::string nom;
  99     std::vector< std::set<_maille>::iterator > mailles; // iterateurs sur les mailles composant le groupe
 100     std::vector<int> groupes; // indices des sous-groupes composant le groupe
 101     std::map<const _maille*,int> relocMap; // map _maille* -> index in MEDMEM::GROUP.getNumber(MED_ALL_ELEMENTS). It is built by _intermediateMED::getGroups()
 102     bool empty() const { return mailles.empty() && groupes.empty(); }
 103 };
 104
 105 struct MEDMEM_EXPORT _fieldBase {
 106   // a field contains several subcomponents each referring to its own support and
 107   // having several named components
 108   struct _sub_data // a subcomponent
 109   {
 110     int                      _supp_id;    // group index within _intermediateMED::groupes
 111     std::vector<std::string> _comp_names; // component names
 112     std::vector<int>         _nb_gauss;   // nb values per element in a component
 113
 114     void setData( int nb_comp, int supp_id )
 115     { _supp_id = supp_id - 1; _comp_names.resize(nb_comp); _nb_gauss.resize(nb_comp,1); }
 116     int nbComponents() const { return _comp_names.size(); }
 117     std::string & compName( int i_comp ) { return _comp_names[ i_comp ]; }
 118     bool isValidNbGauss() const { return *std::max_element( _nb_gauss.begin(), _nb_gauss.end() ) ==
 119                                     *std::min_element( _nb_gauss.begin(), _nb_gauss.end() ); }
 120     int nbGauss() const { return std::max( 1, _nb_gauss[0] ); }
 121     bool hasGauss() const { return nbGauss() > 1; }
 122   };
 123   std::vector< _sub_data > _sub;
 124   int                      _group_id; // group index within _intermediateMED::groupes
 125   // if _group_id == -1 then each subcomponent makes a separate MEDMEM::FIELD, else all subcomponents
 126   // are converted into one MEDMEM::FIELD. The latter is possible only if nb of components in all subs
 127   // is the same and supports of subcomponents do not overlap
 128   MED_EN::med_type_champ   _type;
 129   std::string              _name;
 130   std::string              _description;// field description
 131
 132   _fieldBase( MED_EN::med_type_champ theType, int nb_sub )
 133     : _group_id(-1),_type(theType) { _sub.resize( nb_sub ); }
 134   virtual std::list<std::pair< FIELD_*, int> > getField(std::vector<_groupe>& groupes,
 135                                                         std::vector<GROUP *>& medGroups) const = 0;
 136   void getGroupIds( std::set<int> & ids, bool all ) const; // return ids of main and/or sub-groups
 137   bool hasCommonSupport() const { return _group_id >= 0; } // true if there is one support for all subs
 138   bool hasSameComponentsBySupport() const;
 139   virtual void dump(std::ostream&) const;
 140   virtual ~_fieldBase() {}
 141 };
 142
 143 template< class T > class _field: public _fieldBase
 144 {
 145   std::vector< std::vector< T > > comp_values;
 146  public:
 147   _field< T > ( MED_EN::med_type_champ theType, int nb_sub, int total_nb_comp )
 148     : _fieldBase( theType, nb_sub ) { comp_values.reserve( total_nb_comp ); }
 149   std::vector< T >& addComponent( int nb_values ); // return a vector ready to fill in
 150   std::list<std::pair< FIELD_*, int> > getField(std::vector<_groupe>& groupes,
 151                                                 std::vector<GROUP *>& medGroups) const;
 152   virtual void dump(std::ostream&) const;
 153 };
 154
 155 /*!
 156  * \if developper
 157  * Intermediate structure used by drivers to store data read from the other format file.
 158  * The structure provides functions that transform the stored data to the MED format :
 159  * getCoordinate(), getConnectivity(), getGroups().
 160  * The elements inserted in maillage and points are automaticaly ordered.
 161  * Renumbering are performed by getConnectivity & getGroups before writing the MED structures.
 162  * Read the conception ducumentation for more details.
 163  * \endif
 164  */
 165 struct MEDMEM_EXPORT _intermediateMED
 166 {
 167   std::set<_maille>        maillage;
 168   std::vector<_groupe>     groupes;
 169   std::vector<GROUP *>     medGroupes;
 170   std::map< int, _noeud >  points;
 171   std::list< _fieldBase* > fields;
 172   bool hasMixedCells; // true if there are groups with mixed entity types
 173
 174   CONNECTIVITY * getConnectivity(); // set MED connectivity from the intermediate structure
 175   COORDINATE * getCoordinate(const string & coordinateSystem="CARTESIAN"); // set MED coordinate from the
 176   // intermediate structure
 177   void getFamilies(std::vector<FAMILY *> & _famCell, std::vector<FAMILY *> & _famFace,
 178                    std::vector<FAMILY *> & _famEdge, std::vector<FAMILY *> & _famNode, MESH * _ptrMesh);
 179   void getGroups(std::vector<GROUP *> & _groupCell, std::vector<GROUP *> & _groupFace,
 180                  std::vector<GROUP *> & _groupEdge, std::vector<GROUP *> & _groupNode, MESH * _ptrMesh);
 181   GROUP * getGroup( int i );
 182
 183   void getFields(std::list< FIELD_* >& fields);
 184
 185   // used by previous functions to renumber points & mesh.
 186   void treatGroupes(); // detect groupes of mixed dimention
 187   void numerotationMaillage();
 188   void numerotationPoints();
 189   bool myGroupsTreated;
 190   bool myPointsNumerated;
 191
 192   _intermediateMED() { myGroupsTreated = myPointsNumerated = false; }
 193     ~_intermediateMED();
 194 };
 195
 196 std::ostream& operator << (std::ostream& , const _maille& );
 197 std::ostream& operator << (std::ostream& , const _groupe& );
 198 std::ostream& operator << (std::ostream& , const _noeud& );
 199 std::ostream& operator << (std::ostream& , const _intermediateMED& );
 200 std::ostream& operator << (std::ostream& , const _fieldBase* );
 201
 202 // ===========================================================
 203 //                 field template implementation           //
 204 // ===========================================================
 205
 206 template <class T>
 207   std::vector< T >& _field< T >::addComponent( int nb_values )
 208 {
 209   comp_values.push_back( std::vector< T >() );
 210   std::vector< T >& res = comp_values.back();
 211   res.resize( nb_values );
 212   return res;
 213 }
 214
 215 //=======================================================================
 216 //function : getField
 217 //purpose  : return list of pairs< field, supporting_group_id >
 218 //=======================================================================
 219
 220 template <class T>
 221 std::list<std::pair< FIELD_*, int> > _field< T >::getField(std::vector<_groupe> & groupes,
 222                                                            std::vector<GROUP *> & medGroups) const
 223 {
 224   const char* LOC = "_field< T >::getField()";
 225
 226   std::list<std::pair< FIELD_*, int> > res;
 227
 228   // gauss array data
 229   int nbtypegeo = 0;
 230   vector<int> nbelgeoc(2,0), nbgaussgeo(2,0);
 231
 232   int i_comp_tot = 0, nb_fields = 0;
 233   std::set<int> supp_id_set; // to create a new field when support repeats if hasCommonSupport()
 234   std::vector< _sub_data >::const_iterator sub_data, sub_end = _sub.end();
 235
 236   _groupe*  grp = 0;
 237   GROUP* medGrp = 0;
 238   if ( hasCommonSupport() ) // several subs are combined into one field
 239   {
 240     grp    = & groupes[ _group_id ];
 241     medGrp = medGroups[ _group_id ];
 242     if ( !grp || grp->empty() || !medGrp || !medGrp->getNumberOfTypes())
 243       return res;
 244
 245     // Make gauss array data
 246     nbtypegeo = medGrp->getNumberOfTypes();
 247     nbelgeoc  .resize( nbtypegeo + 1, 0 );
 248     nbgaussgeo.resize( nbtypegeo + 1, 0 );
 249     const int * nbElemByType = medGrp->getNumberOfElements();
 250     sub_data = _sub.begin();
 251     for (int iType = 0; iType < nbtypegeo; ++iType) {
 252       nbelgeoc  [ iType+1 ] = nbelgeoc[ iType ] + nbElemByType[ iType ];
 253       nbgaussgeo[ iType+1 ] = sub_data->nbGauss();
 254       int nbElemInSubs = 0;
 255       while ( nbElemInSubs < nbElemByType[ iType ] && sub_data != sub_end ) {
 256         nbElemInSubs += groupes[ sub_data->_supp_id ].relocMap.size();
 257         ++sub_data;
 258       }
 259     }
 260   }
 261   typedef typename MEDMEM_ArrayInterface<T,FullInterlace,NoGauss>::Array TArrayNoGauss;
 262   typedef typename MEDMEM_ArrayInterface<T,FullInterlace,Gauss>::Array   TArrayGauss;
 263   FIELD< T, FullInterlace > * f = 0;
 264   TArrayNoGauss * arrayNoGauss = 0;
 265   TArrayGauss   * arrayGauss = 0;
 266
 267   // loop on subs of this field
 268   int i_sub = 1;
 269   for ( sub_data = _sub.begin(); sub_data != sub_end; ++sub_data, ++i_sub )
 270   {
 271     // nb values in a field
 272     if ( !hasCommonSupport() ) {
 273       grp    = & groupes[ sub_data->_supp_id ];
 274       medGrp = medGroups[ sub_data->_supp_id ];
 275     }
 276     int nb_val = grp->relocMap.size();
 277
 278     // check validity of a sub_data
 279     bool validSub = true;
 280     if ( !nb_val ) {
 281       INFOS("Skip field <" << _name << ">: invalid supporting group "
 282             << (hasCommonSupport() ? _group_id : sub_data->_supp_id )
 283             << " of " << i_sub << "-th subcomponent" );
 284       validSub = false;
 285     }
 286     if ( !sub_data->isValidNbGauss() ) {
 287       INFOS("Skip field <" << _name << ">: different nb of gauss points in components ");
 288       validSub = false;
 289     }
 290     if ( !validSub ) {
 291       if ( hasCommonSupport() ) {
 292         if ( !res.empty() ) {
 293           delete f;
 294           res.clear();
 295         }
 296         return res;
 297       }
 298       i_comp_tot += sub_data->nbComponents();
 299       continue;
 300     }
 301
 302     // Create a field
 303
 304     if ( !f || !hasCommonSupport() || !supp_id_set.insert( sub_data->_supp_id ).second )
 305     {
 306       ++nb_fields;
 307       supp_id_set.clear();
 308       arrayNoGauss = 0;
 309       arrayGauss = 0;
 310
 311       f = new FIELD< T, FullInterlace >();
 312
 313       f->setNumberOfComponents( sub_data->nbComponents() );
 314       f->setComponentsNames( & sub_data->_comp_names[ 0 ] );
 315       f->setNumberOfValues ( nb_val );
 316       f->setName( _name );
 317       f->setDescription( _description );
 318       vector<string> str( sub_data->nbComponents() );
 319       f->setComponentsDescriptions( &str[0] );
 320       f->setMEDComponentsUnits( &str[0] );
 321
 322       res.push_back( make_pair( f , hasCommonSupport() ? _group_id : sub_data->_supp_id ));
 323       MESSAGE(" MAKE " << nb_fields << "-th field <" << _name << "> on group_id " << _group_id );
 324
 325       // make an array
 326       if ( !sub_data->hasGauss() ) {
 327         arrayNoGauss = new TArrayNoGauss( sub_data->nbComponents(), nb_val );
 328         f->setArray( arrayNoGauss );
 329       }
 330       else {
 331         if ( !hasCommonSupport() ) {
 332           nbtypegeo = 1;
 333           nbelgeoc  [1] = nb_val;
 334           nbgaussgeo[1] = sub_data->nbGauss();
 335         }
 336         arrayGauss = new TArrayGauss(sub_data->nbComponents(), nb_val,
 337                                      nbtypegeo, & nbelgeoc[0], & nbgaussgeo[0]);
 338         f->setArray( arrayGauss );
 339
 340         // PAL11040 "GIBI driver for Castem fields with Gauss point values"
 341         const MED_EN::medGeometryElement* types = medGrp->getTypes();
 342         for (int iGeom = 0; iGeom < nbtypegeo; ++iGeom) {
 343           ostringstream name;
 344           name << "Gauss_" << nbgaussgeo[iGeom+1] << "points_on" << types[iGeom] << "geom";
 345           GAUSS_LOCALIZATION_* loc = GAUSS_LOCALIZATION_::makeDefaultLocalization
 346             (name.str(), types[iGeom], nbgaussgeo[iGeom+1]);
 347           f->setGaussLocalization( types[iGeom], loc );
 348         }
 349       }
 350     }
 351
 352     // Set values
 353
 354     // get nb elements in a group
 355     _groupe & sub_grp = groupes[ sub_data->_supp_id ];
 356     int nb_supp_elems = sub_grp.mailles.size();
 357     int nb_gauss      = sub_data->nbGauss();
 358     MESSAGE("insert sub data, group_id: " << sub_data->_supp_id <<
 359             ", nb values: "               << comp_values[ i_comp_tot ].size() <<
 360             ", relocMap size: "           << sub_grp.relocMap.size() <<
 361             ", nb mailles: "              << nb_supp_elems);
 362
 363 #ifdef CASTEM_FULL_INTERLACE
 364     const int gauss_step = 1;
 365     const int elem_step = nb_gauss;
 366 #else
 367     const int gauss_step = nb_supp_elems;
 368     const int elem_step = 1;
 369 #endif
 370     int i; // elem index
 371     // loop on components of a sub
 372     for ( int i_comp = 1; i_comp <= sub_data->nbComponents(); ++i_comp )
 373     {
 374       // store values
 375       const std::vector< T > & values = comp_values[ i_comp_tot++ ];
 376       bool oneValue = ( values.size() == 1 );
 377       ASSERT( oneValue || values.size() == nb_supp_elems * nb_gauss );
 378       for ( int k = 0; k < nb_supp_elems; ++k )
 379       {
 380         const T& val = oneValue ? values[ 0 ] : values[ k * elem_step ];
 381         const _maille* ma = &(*sub_grp.mailles[ k ]);
 382         if ( medGrp->isOnAllElements() ) {
 383           i = ma->ordre;
 384           if ( i > nb_val )
 385             throw MEDEXCEPTION (LOCALIZED(STRING(LOC) << ", wrong elem position. "
 386                                           << k << "-th elem: " << ma
 387                                           << ", pos (" << i << ") must be <= " << nb_val));
 388         }
 389         else {
 390           std::map<const _maille*,int>::const_iterator ma_i = grp->relocMap.find( ma );
 391           if ( ma_i == grp->relocMap.end() )
 392             throw MEDEXCEPTION (LOCALIZED(STRING(LOC) << ", cant find elem index. "
 393                                           << k << "-th elem: " << ma));
 394           i = ma_i->second;
 395         }
 396         if ( arrayNoGauss ) {
 397           arrayNoGauss->setIJ( i, i_comp, val );
 398         }
 399         else {
 400           const T* pVal = & val;
 401           for ( int iGauss = 1; iGauss <= nb_gauss; ++iGauss, pVal += gauss_step )
 402             arrayGauss->setIJK( i, i_comp, iGauss, *pVal);
 403         }
 404       }
 405     }
 406   } // loop on subs of the field
 407
 408   return res;
 409 }
 410
 411
 412 template <class T> void _field< T >::dump(std::ostream& os) const
 413 {
 414   _fieldBase::dump(os);
 415   os << endl;
 416   for ( int i = 0 ; i < (int)comp_values.size(); ++i )
 417   {
 418     os << "    " << i+1 << "-th component, nb values: " << comp_values[ i ].size() << endl;
 419   }
 420 }
 421
 422 }; // namespace MEDMEM
 423
 424 #endif /* DRIVERTOOLS_HXX */