1 // Copyright (C) 2007-2012 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
23 #ifndef DRIVERTOOLS_HXX
24 #define DRIVERTOOLS_HXX
28 #include "MEDMEM_define.hxx"
29 #include "MEDMEM_Exception.hxx"
30 #include "MEDMEM_FieldForward.hxx"
31 #include "MEDMEM_ArrayInterface.hxx"
32 #include "MEDMEM_Group.hxx"
33 #include "MEDMEM_GaussLocalization.hxx"
44 // Way of storing values with gauss points by CASTEM: full/no interlace.
45 // Remove code switched by this symbol as soon as the way is found out,
46 // from here and from gibi driver
47 #define CASTEM_FULL_INTERLACE
57 class _intermediateMED;
59 // ==============================================================================
60 struct MEDMEM_EXPORT _noeud
62 mutable int number; // negative if node is merged
63 std::vector<double> coord;
66 // ==============================================================================
67 typedef pair<int,int> _link; // a pair of node numbers
69 // ==============================================================================
70 struct MEDMEM_EXPORT _maille
72 typedef std::map<int,_noeud>::iterator TNoeud; //iter;
73 std::vector< TNoeud > sommets;
74 MED_EN::medGeometryElement geometricType;
75 mutable bool reverse; // to reverse sommets of a face
76 mutable int* sortedNodeIDs; // for comparison and merge
77 //mutable list<unsigned> groupes; // the GROUPs maille belongs to, used to create families
79 _maille(MED_EN::medGeometryElement type=MED_EN::MED_NONE, size_t nelem=0)
80 : geometricType(type),reverse(false),sortedNodeIDs(0),_ordre(0) { sommets.reserve(nelem); }
82 _maille(const _maille& ma);
83 void init() const { if ( sortedNodeIDs ) delete [] sortedNodeIDs; sortedNodeIDs = 0; }
84 ~_maille() { init(); }
86 int dimension() const // retourne la dimension de la maille
87 { return geometricType/100; }
89 int dimensionWithPoly() const // retourne la dimension de la maille
90 { return geometricType >= MED_EN::MED_POLYGON ? dimension()-2 : dimension(); }
92 const int* getSortedNodes() const; // creates if needed and return sortedNodeIDs
93 bool operator < (const _maille& ma) const;
94 MED_EN::medEntityMesh getEntity(const int meshDimension) const throw (MEDEXCEPTION);
95 _link link(int i) const;
97 //bool isMerged(int i) const { return sommets[i]->second.number < 0; }
98 int nodeNum(int i) const { return abs( sommets[i]->second.number ); }
99 int nodeID (int i) const { return sommets[i]->first; } // key in points map
101 unsigned ordre() const { return abs( _ordre ); }
102 bool isMerged() const { return _ordre < 0; }
103 void setMergedOrdre(unsigned o) const { _ordre = -o; }
104 void setOrdre(int o) const { _ordre = o; }
107 mutable int _ordre; // l'ordre est fixé après insertion dans le set, et ne change ni l'état, ni l'ordre -> mutable
108 _maille& operator=(const _maille& ma);
111 // ==============================================================================
112 struct MEDMEM_EXPORT _mailleIteratorCompare // pour ordonner le set d'iterateurs sur mailles
114 // this operator is used for ordering mailles thinin a group, which happens when
115 // right numbers are already given to _maille's, so we can use _maille::ordre
116 // for comparison instead of a heavy _maille::operator <
117 bool operator () (std::set<_maille>::iterator i1, std::set<_maille>::iterator i2)
120 return i1->ordre() < i2->ordre();
124 // ==============================================================================
125 struct MEDMEM_EXPORT _groupe
127 typedef std::set<_maille>::iterator TMaille;
128 typedef std::vector< TMaille >::const_iterator TMailleIter;
130 std::vector< TMaille > mailles; // iterateurs sur les mailles composant le groupe
131 std::vector<int> groupes; // indices des sous-groupes composant le groupe
132 std::map<unsigned,int> relocMap; // map _maille::ordre() -> index in GROUP, built by getGroups()
135 std::vector<std::string> refNames; /* names of groups referring this one;
136 refNames is resized according to nb of references
137 while reading a group (pile 1) and it is filled with
138 names while reading long names (pile 27); each named
139 reference is converted into a copy of the medGroup
142 const _maille& maille(int index) { return *mailles[index]; }
143 bool empty() const { return mailles.empty() && groupes.empty(); }
145 int size() const { return (mailles.size()>relocMap.size())?mailles.size():relocMap.size(); }
147 int size() const { return std::max( mailles.size(), relocMap.size() ); }
149 _groupe():medGroup(0) {}
152 // ==============================================================================
153 struct MEDMEM_EXPORT _fieldBase {
154 // a field contains several subcomponents each referring to its own support and
155 // having several named components
156 // ----------------------------------------------------------------------------
157 struct _sub_data // a subcomponent
158 // --------------------------------
160 int _supp_id; // group index within _intermediateMED::groupes
161 std::vector<std::string> _comp_names; // component names
162 std::vector<int> _nb_gauss; // nb values per element in a component
164 void setData( int nb_comp, int supp_id )
165 { _supp_id = supp_id - 1; _comp_names.resize(nb_comp); _nb_gauss.resize(nb_comp,1); }
166 int nbComponents() const { return _comp_names.size(); }
167 std::string & compName( int i_comp ) { return _comp_names[ i_comp ]; }
168 bool isValidNbGauss() const { return *std::max_element( _nb_gauss.begin(), _nb_gauss.end() ) ==
169 *std::min_element( _nb_gauss.begin(), _nb_gauss.end() ); }
171 int nbGauss() const { return (1>_nb_gauss[0])?1:_nb_gauss[0]; }
173 int nbGauss() const { return std::max( 1, _nb_gauss[0] ); }
175 bool hasGauss() const { return nbGauss() > 1; }
177 // ----------------------------------------------------------------------------
178 std::vector< _sub_data > _sub;
179 int _group_id; // group index within _intermediateMED::groupes
180 // if _group_id == -1 then each subcomponent makes a separate MEDMEM::FIELD, else all subcomponents
181 // are converted into one MEDMEM::FIELD. The latter is possible only if nb of components in all subs
182 // is the same and supports of subcomponents do not overlap
183 MED_EN::med_type_champ _type;
185 std::string _description;// field description
187 _fieldBase( MED_EN::med_type_champ theType, int nb_sub )
188 : _group_id(-1),_type(theType) { _sub.resize( nb_sub ); }
189 virtual std::list<std::pair< FIELD_*, int> > getField(std::vector<_groupe>& groupes) const = 0;
190 void getGroupIds( std::set<int> & ids, bool all ) const; // return ids of main and/or sub-groups
191 bool hasCommonSupport() const { return _group_id >= 0; } // true if there is one support for all subs
192 bool hasSameComponentsBySupport() const;
193 virtual void dump(std::ostream&) const;
194 virtual ~_fieldBase() {}
197 // ==============================================================================
198 template< class T > class _field: public _fieldBase
200 std::vector< std::vector< T > > comp_values;
202 _field< T > ( MED_EN::med_type_champ theType, int nb_sub, int total_nb_comp )
203 : _fieldBase( theType, nb_sub ) { comp_values.reserve( total_nb_comp ); }
204 std::vector< T >& addComponent( int nb_values ); // return a vector ready to fill in
205 std::list<std::pair< FIELD_*, int> > getField(std::vector<_groupe>& groupes) const;
206 virtual void dump(std::ostream&) const;
209 // ==============================================================================
212 * Iterator on set of _maille's of given dimension
215 class MEDMEM_EXPORT _maillageByDimIterator
218 // if (convertPoly) then treat poly as 2d and 3d, else as 4d and 5d (=medGeometryElement/100)
219 _maillageByDimIterator( const _intermediateMED & medi,
220 const int dim=-1, // dim=-1 - for all dimensions
221 const bool convertPoly = false )
222 { myImed = & medi; init( dim, convertPoly ); }
224 // if (convertPoly) then treat poly as 2d and 3d, else as 4d and 5d (=medGeometryElement/100)
225 void init(const int dim=-1, // dim=-1 - for all dimensions
226 const bool convertPoly = false );
228 //!< return next set of _maille's of required dimension
229 const std::set<_maille > * nextType() {
230 while ( myIt != myEnd )
231 if ( !myIt->second.empty() && ( myDim < 0 || dim(false) == myDim ))
232 return & (myIt++)->second;
237 //!< return dimension of mailles returned by the last or further next()
238 int dim(const bool last=true) const {
241 return myConvertPoly ?
242 it->second.begin()->dimensionWithPoly() :
243 it->second.begin()->dimension();
245 //!< return type of mailles returned by the last next()
246 MED_EN::medGeometryElement type() const { iterator it = myIt; return (--it)->first; }
248 //!< return number of mailles taking into account merged ones
249 int sizeWithoutMerged() const {
251 removed::const_iterator tNb = nbRemovedByType->find( (--it)->first );
252 return it->second.size() - ( tNb == nbRemovedByType->end() ? 0 : tNb->second );
255 typedef std::map<MED_EN::medGeometryElement, int > removed;
256 typedef std::map<MED_EN::medGeometryElement, std::set<_maille > > TMaillageByType;
257 typedef TMaillageByType::const_iterator iterator;
259 const _intermediateMED* myImed;
260 iterator myIt, myEnd;
261 int myDim, myConvertPoly;
262 const removed * nbRemovedByType;
265 // ==============================================================================
268 * Intermediate structure used by drivers to store data read from the other format file.
269 * The structure provides functions that transform the stored data to the MED format :
270 * getCoordinate(), getConnectivity(), getGroups().
271 * The elements inserted in maillage and points are automaticaly ordered.
272 * Renumbering are performed by getConnectivity & getGroups before writing the MED structures.
273 * Read the conception ducumentation for more details.
276 struct MEDMEM_EXPORT _intermediateMED
278 typedef std::map<MED_EN::medGeometryElement, std::set<_maille > > TMaillageByType;
279 typedef std::map<MED_EN::medGeometryElement, int > TNbByType;
280 typedef std::map< const _maille*, std::vector<int> > TPolyherdalNbFaceNodes;
282 TNbByType nbRemovedByType; // nb mailles removed by merge
283 std::vector<_groupe> groupes;
284 //std::vector<GROUP *> medGroupes;
285 std::map< int, _noeud > points;
286 std::list< _fieldBase* > fields;
287 bool hasMixedCells; // true if there are groups with mixed entity types
288 TPolyherdalNbFaceNodes polyherdalNbFaceNodes; // nb of nodes in faces for each polyhedron
290 inline _groupe::TMaille insert(const _maille& ma);
292 int getMeshDimension() const;
293 void mergeNodesAndElements(double tolerance); // optionally merge nodes and elements
294 CONNECTIVITY * getConnectivity(); // creates MED connectivity from the intermediate structure
295 COORDINATE * getCoordinate(const string & coordinateSystem="CARTESIAN"); // makes MED coordinate
296 // void getFamilies(std::vector<FAMILY *> & _famCell, std::vector<FAMILY *> & _famFace,
297 // std::vector<FAMILY *> & _famEdge, std::vector<FAMILY *> & _famNode,
299 void getGroups(std::vector<GROUP *> & _groupCell, std::vector<GROUP *> & _groupFace,
300 std::vector<GROUP *> & _groupEdge, std::vector<GROUP *> & _groupNode,
302 //GROUP * getGroup( int i );
304 void getFields(std::list< FIELD_* >& fields);
306 // used by previous functions to renumber points & mesh.
307 bool myGroupsTreated;
308 void treatGroupes(); // detect groupes of mixed dimension
309 void numerotationMaillage();
310 bool numerotationPoints(); // return true if renumeration done
311 int nbMerged(int geoType) const; //!< nb nodes removed by merge
314 { myNodesNumerated = myMaillesNumerated = myGroupsTreated = false; currentTypeMailles = 0; }
319 bool myNodesNumerated, myMaillesNumerated;
320 // mailles groupped by geom type; use insert() for filling in and
321 // _maillageByDimIterator for exploring it
322 //std::set<_maille> maillage;
323 TMaillageByType maillageByType;
324 TMaillageByType::value_type* currentTypeMailles; // for fast insertion
325 friend class _maillageByDimIterator;
327 //-----------------------------------------------------------------------
328 _groupe::TMaille _intermediateMED::insert(const _maille& ma)
330 if ( !currentTypeMailles || currentTypeMailles->first != ma.geometricType )
331 currentTypeMailles = & *maillageByType.insert
332 ( make_pair( ma.geometricType, std::set<_maille >())).first;
334 _groupe::TMaille res = currentTypeMailles->second.insert( ma ).first;
336 ((_maille&)ma).init(); // this method was added for the sake of this call which is needed to
337 // remove comparison key (sortedNodeIDs) from a temporary _maille ma
342 // ==============================================================================
344 std::ostream& operator << (std::ostream& , const _maille& );
345 std::ostream& operator << (std::ostream& , const _groupe& );
346 std::ostream& operator << (std::ostream& , const _noeud& );
347 std::ostream& operator << (std::ostream& , const _intermediateMED& );
348 std::ostream& operator << (std::ostream& , const _fieldBase* );
350 // ===========================================================
351 // field template implementation //
352 // ===========================================================
355 std::vector< T >& _field< T >::addComponent( int nb_values )
357 comp_values.push_back( std::vector< T >() );
358 std::vector< T >& res = comp_values.back();
359 res.resize( nb_values );
363 //=======================================================================
364 //function : getField
365 //purpose : return list of pairs< field, supporting_group_id >
366 //=======================================================================
369 std::list<std::pair< FIELD_*, int> > _field< T >::getField(std::vector<_groupe> & groupes) const
371 const char* LOC = "_field< T >::getField()";
373 std::list<std::pair< FIELD_*, int> > res;
377 vector<int> nbelgeoc(2,0), nbgaussgeo(2,0);
379 int i_comp_tot = 0, nb_fields = 0;
380 std::set<int> supp_id_set; // to create a new field when support repeats if hasCommonSupport()
381 std::vector< _sub_data >::const_iterator sub_data, sub_end = _sub.end();
385 if ( hasCommonSupport() ) // several subs are combined into one field
387 grp = & groupes[ _group_id ];
388 medGrp = grp->medGroup;
389 if ( !grp || grp->empty() || !medGrp || !medGrp->getNumberOfTypes())
392 // Make gauss array data
393 nbtypegeo = medGrp->getNumberOfTypes();
394 nbelgeoc .resize( nbtypegeo + 1, 0 );
395 nbgaussgeo.resize( nbtypegeo + 1, 0 );
396 const int * nbElemByType = medGrp->getNumberOfElements();
397 const MED_EN::medGeometryElement* types = medGrp->getTypes();
398 for (int iType = 0; iType < nbtypegeo; ++iType) {
400 nbelgeoc [ iType+1 ] = nbelgeoc[ iType ] + nbElemByType[ iType ];
401 // nb gauss by type; find a sub for a geo type
402 for ( sub_data = _sub.begin(); sub_data != sub_end; ++sub_data ) {
403 _groupe & sub_grp = groupes[ sub_data->_supp_id ];
404 if ( !sub_grp.empty() && sub_grp.mailles[0]->geometricType == types[ iType ])
407 ASSERT_MED( sub_data != sub_end );
408 nbgaussgeo[ iType+1 ] = sub_data->nbGauss();
411 typedef typename MEDMEM_ArrayInterface<T,FullInterlace,NoGauss>::Array TArrayNoGauss;
412 typedef typename MEDMEM_ArrayInterface<T,FullInterlace,Gauss>::Array TArrayGauss;
413 FIELD< T, FullInterlace > * f = 0;
414 TArrayNoGauss * arrayNoGauss = 0;
415 TArrayGauss * arrayGauss = 0;
417 // loop on subs of this field
419 for ( sub_data = _sub.begin(); sub_data != sub_end; ++sub_data, ++i_sub )
421 // nb values in a field
422 if ( !hasCommonSupport() ) {
423 grp = & groupes[ sub_data->_supp_id ];
424 medGrp = grp->medGroup;
426 int nb_val = grp->size();
428 // check validity of a sub_data
429 bool validSub = true;
431 INFOS_MED("Skip field <" << _name << ">: invalid supporting group "
432 << (hasCommonSupport() ? _group_id : sub_data->_supp_id )
433 << " of " << i_sub << "-th subcomponent" );
436 if ( !sub_data->isValidNbGauss() ) {
437 INFOS_MED("Skip field <" << _name << ">: different nb of gauss points in components ");
441 if ( hasCommonSupport() ) {
442 if ( !res.empty() ) {
444 f->removeReference();
449 i_comp_tot += sub_data->nbComponents();
455 if ( !f || !hasCommonSupport() || !supp_id_set.insert( sub_data->_supp_id ).second )
462 f = new FIELD< T, FullInterlace >();
464 f->setNumberOfComponents( sub_data->nbComponents() );
465 f->setComponentsNames( & sub_data->_comp_names[ 0 ] );
466 f->setNumberOfValues ( nb_val );
468 f->setDescription( _description );
469 vector<string> str( sub_data->nbComponents() );
470 f->setComponentsDescriptions( &str[0] );
471 f->setMEDComponentsUnits( &str[0] );
472 if ( !hasCommonSupport() && nb_fields > 1 )
474 f->setName( MEDMEM::STRING(_name) << "_Sub_" << nb_fields );
475 INFOS_MED("Warning: field |" <<_name<<"| is incompatible with MED format (has "
476 "sub-fields of different nature), so we map its sub-field #"<< nb_fields <<
477 " into a separate field |"<<f->getName() << "|");
479 res.push_back( make_pair( f , hasCommonSupport() ? _group_id : sub_data->_supp_id ));
480 MESSAGE_MED(" MAKE " << nb_fields << "-th field <" << _name << "> on group_id " << _group_id );
483 if ( !sub_data->hasGauss() ) {
484 arrayNoGauss = new TArrayNoGauss( sub_data->nbComponents(), nb_val );
485 f->setArray( arrayNoGauss );
488 if ( !hasCommonSupport() ) {
490 nbelgeoc [1] = nb_val;
491 nbgaussgeo[1] = sub_data->nbGauss();
493 arrayGauss = new TArrayGauss(sub_data->nbComponents(), nb_val,
494 nbtypegeo, & nbelgeoc[0], & nbgaussgeo[0]);
495 f->setArray( arrayGauss );
497 // PAL11040 "GIBI driver for Castem fields with Gauss point values"
498 const MED_EN::medGeometryElement* types = medGrp->getTypes();
499 for (int iGeom = 0; iGeom < nbtypegeo; ++iGeom) {
501 name << "Gauss_" << nbgaussgeo[iGeom+1] << "points_on" << types[iGeom] << "geom";
502 GAUSS_LOCALIZATION_* loc = GAUSS_LOCALIZATION_::makeDefaultLocalization
503 (name.str(), types[iGeom], nbgaussgeo[iGeom+1]);
504 f->setGaussLocalization( types[iGeom], loc );
511 // get nb elements in a group
512 _groupe & sub_grp = groupes[ sub_data->_supp_id ];
513 int nb_supp_elems = sub_grp.mailles.size();
514 int nb_gauss = sub_data->nbGauss();
515 MESSAGE_MED("insert sub data, group_id: " << sub_data->_supp_id <<
516 ", nb values: " << comp_values[ i_comp_tot ].size() <<
517 ", relocMap size: " << sub_grp.relocMap.size() <<
518 ", nb mailles: " << nb_supp_elems);
520 #ifdef CASTEM_FULL_INTERLACE
521 const int gauss_step = 1;
522 const int elem_step = nb_gauss;
524 const int gauss_step = nb_supp_elems;
525 const int elem_step = 1;
528 // loop on components of a sub
529 for ( int i_comp = 1; i_comp <= sub_data->nbComponents(); ++i_comp )
532 const std::vector< T > & values = comp_values[ i_comp_tot++ ];
533 bool oneValue = ( values.size() == 1 );
534 ASSERT_MED( oneValue || (int)values.size() == nb_supp_elems * nb_gauss );
535 for ( int k = 0; k < nb_supp_elems; ++k )
537 const T& val = oneValue ? values[ 0 ] : values[ k * elem_step ];
538 const _maille& ma = sub_grp.maille( k );
539 if ( medGrp->isOnAllElements() ) {
543 std::map<unsigned,int>::const_iterator ordre_i = grp->relocMap.find( ma.ordre() );
544 if ( ordre_i == grp->relocMap.end() )
545 throw MEDEXCEPTION (LOCALIZED(STRING(LOC) << ", cant find elem index. "
546 << k << "-th elem: " << ma));
549 if ( arrayNoGauss ) {
550 arrayNoGauss->setIJ( i, i_comp, val );
553 const T* pVal = & val;
554 for ( int iGauss = 1; iGauss <= nb_gauss; ++iGauss, pVal += gauss_step )
555 arrayGauss->setIJK( i, i_comp, iGauss, *pVal);
559 } // loop on subs of the field
565 // ==============================================================================
566 template <class T> void _field< T >::dump(std::ostream& os) const
567 // ==============================================================================
569 _fieldBase::dump(os);
571 for ( int i = 0 ; i < (int)comp_values.size(); ++i )
573 os << " " << i+1 << "-th component, nb values: " << comp_values[ i ].size() << endl;
577 } // namespace MEDMEM
579 #endif /* DRIVERTOOLS_HXX */