src/MEDMEM/MEDMEM_AsciiFieldDriver.hxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   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.
  10 //
  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.
  15 //
  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
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #ifndef ASCII_FIELD_DRIVER_HXX
  24 #define ASCII_FIELD_DRIVER_HXX
  25
  26 #include "MEDMEM.hxx"
  27 #include "MEDMEM_GenDriver.hxx"
  28 #include "MEDMEM_Exception.hxx"
  29 #include "MEDMEM_Unit.hxx"
  30 #include "MEDMEM_Support.hxx"
  31 #include "MEDMEM_Mesh.hxx"
  32 #include "MEDMEM_ArrayInterface.hxx"
  33 #include "MEDMEM_ArrayConvert.hxx"
  34
  35 #include <list>
  36 #include <string>
  37 #include <ctype.h>
  38 #include <iomanip>
  39 #include <stdlib.h>
  40 #include <string.h>
  41 #include <fstream>
  42
  43 namespace MEDMEM {
  44
  45   const int    PRECISION_IN_ASCII_FILE = 10;
  46   const double PRECISION_IN_COMPARE    = 1e-10;
  47   const int    SPACE_BETWEEN_NBS       = 19;
  48
  49   template<int N,unsigned int CODE>
  50   void fill(double *a, const double *b)
  51   {
  52     a[N]=b[CODE & 0x3 ];
  53     fill<N-1, (CODE>>2) >(a,b);
  54   }
  55
  56   template<int N>
  57   bool compare(const double* a, const double* b)
  58   {
  59     double sign=b[N]<0 ? -1 : 1;
  60     if(a[N]<b[N]*(1-sign*PRECISION_IN_COMPARE))
  61       return true;
  62     if(a[N]>b[N]*(1+sign*PRECISION_IN_COMPARE))
  63       return false;
  64     return compare<N-1>(a,b);
  65   }
  66
  67   template<> MEDMEM_EXPORT
  68   void fill<-1,0x3>(double *a, const double *b);
  69
  70   template<> MEDMEM_EXPORT
  71   bool compare<-1>(const double *a, const double *b);
  72
  73   template <class T, int SPACEDIMENSION, unsigned int SORTSTRATEGY>
  74   class SDForSorting
  75   {
  76   private:
  77     double _coords[SPACEDIMENSION];
  78     T *_components;
  79     int _nbComponents;
  80   public:
  81     SDForSorting(const double *coords, const T* comp, int nbComponents);
  82     SDForSorting(const SDForSorting& other);
  83     ~SDForSorting();
  84     bool operator< (const SDForSorting<T,SPACEDIMENSION,SORTSTRATEGY>& other) const;
  85     void writeLine(ofstream& file) const;
  86   };
  87
  88   template <class T>
  89   class ASCII_FIELD_DRIVER : public GENDRIVER
  90   {
  91   private:
  92     MESH                   *_mesh;
  93     SUPPORT                *_support;
  94     mutable FIELD<T>       *_ptrField;
  95     std::string             _fileName;
  96     mutable ofstream        _file;
  97     unsigned int            _code;
  98     MED_EN::med_sort_direc  _direc;
  99     int                     _nbComponents;
 100     int                     _spaceDimension;
 101     //static int           _nbComponentsForCpyInfo;
 102
 103   public:
 104     template <class INTERLACING_TAG>
 105     ASCII_FIELD_DRIVER():GENDRIVER(ASCII_DRIVER),
 106                          _ptrField((FIELD<T>)MED_NULL),
 107                          _fileName("") {}
 108
 109     template <class INTERLACING_TAG>
 110     ASCII_FIELD_DRIVER(const string & fileName, FIELD<T,INTERLACING_TAG> * ptrField,
 111                        MED_EN::med_sort_direc direction=MED_EN::ASCENDING,
 112                        const char *priority="");
 113
 114
 115     ASCII_FIELD_DRIVER(const ASCII_FIELD_DRIVER<T>& other);
 116     void open() throw (MEDEXCEPTION);
 117     void close();
 118     void read ( void ) throw (MEDEXCEPTION);
 119     void write( void ) const throw (MEDEXCEPTION);
 120     GENDRIVER* copy() const;
 121   private:
 122     void buildIntroduction() const;
 123     template<int SPACEDIMENSION, unsigned int SORTSTRATEGY>
 124     void sortAndWrite() const;
 125     //template<int SPACEDIMENSION, unsigned int SORTSTRATEGY>//, std::list< SDForSorting<T,SPACEDIMENSION,SORTSTRATEGY > > >
 126     //static void copyInfo(const double *a,T *b);
 127     //static void copyInfo2(const double *,T *);
 128   };
 129 }
 130
 131
 132 namespace MEDMEM {
 133
 134   template <class T, int SPACEDIMENSION, unsigned int SORTSTRATEGY>
 135   SDForSorting<T,SPACEDIMENSION,SORTSTRATEGY>::SDForSorting(const double *coords, const T* comp, int nbComponents):_nbComponents(nbComponents)
 136   {
 137     fill<SPACEDIMENSION-1,SORTSTRATEGY>(_coords,coords);
 138     _components=new T[_nbComponents];
 139     memcpy(_components,comp,sizeof(T)*_nbComponents);
 140   }
 141
 142   template <class T, int SPACEDIMENSION, unsigned int SORTSTRATEGY>
 143   SDForSorting<T,SPACEDIMENSION,SORTSTRATEGY>::SDForSorting(const SDForSorting& other):_nbComponents(other._nbComponents)
 144   {
 145     memcpy(_coords,other._coords,sizeof(double)*SPACEDIMENSION);
 146     _components=new T[_nbComponents];
 147     memcpy(_components,other._components,sizeof(T)*_nbComponents);
 148   }
 149
 150   template <class T, int SPACEDIMENSION, unsigned int SORTSTRATEGY>
 151   SDForSorting<T,SPACEDIMENSION,SORTSTRATEGY>::~SDForSorting()
 152   {
 153     delete [] _components;
 154   }
 155
 156   template <class T, int SPACEDIMENSION, unsigned int SORTSTRATEGY>
 157   bool SDForSorting<T,SPACEDIMENSION,SORTSTRATEGY>::operator< (const SDForSorting<T,SPACEDIMENSION,SORTSTRATEGY>& other) const
 158   {
 159     return compare<SPACEDIMENSION-1>(_coords,other._coords);
 160   }
 161
 162   template <class T, int SPACEDIMENSION, unsigned int SORTSTRATEGY>
 163   void SDForSorting<T,SPACEDIMENSION,SORTSTRATEGY>::writeLine(ofstream& file) const
 164   {
 165     int i;
 166     double temp[SPACEDIMENSION];
 167     fill<SPACEDIMENSION-1,SORTSTRATEGY>(temp,_coords);
 168     for(i=0;i<SPACEDIMENSION;i++)
 169       file << setw(SPACE_BETWEEN_NBS) << temp[i];
 170     for(i=0;i<_nbComponents;i++)
 171       file << setw(SPACE_BETWEEN_NBS) << _components[i];
 172     file << endl;
 173   }
 174
 175   template <class T>
 176   template <class INTERLACING_TAG>
 177   ASCII_FIELD_DRIVER<T>::ASCII_FIELD_DRIVER(const string & fileName,
 178                                             FIELD<T,INTERLACING_TAG> * ptrField,
 179                                             MED_EN::med_sort_direc direction,
 180                                             const char *priority)
 181     :GENDRIVER(fileName, MED_EN::WRONLY, ASCII_DRIVER),
 182      _ptrField((FIELD<T>*)ptrField),
 183      _fileName(fileName),
 184      _direc(direction)
 185     {
 186       _nbComponents=_ptrField->getNumberOfComponents();
 187       if(_nbComponents<=0)
 188         throw MEDEXCEPTION("ASCII_FIELD_DRIVER : No components in FIELD<T>");
 189       _support=(SUPPORT *)_ptrField->getSupport();
 190       _mesh=(MESH *)_support->getMesh();
 191       _spaceDimension=_mesh->getSpaceDimension();
 192       _code=3;
 193       int i;
 194       if(priority[0]=='\0')
 195         for(i=_spaceDimension-1;i>=0;i--)
 196           {
 197             _code<<=2;
 198             _code+=i;
 199           }
 200       else
 201         {
 202           if(_spaceDimension != (int)strlen(priority))
 203             throw MEDEXCEPTION("ASCII_FIELD_DRIVER : Coordinate priority invalid with spaceDim");
 204           for(i=_spaceDimension-1;i>=0;i--)
 205             {
 206               char c=toupper(priority[i]);
 207               if(int(c-'X')>(_spaceDimension-1) || int(c-'X')<0)
 208                 throw MEDEXCEPTION("ASCII_FIELD_DRIVER : Invalid priority definition");
 209               _code<<=2;
 210               _code+=c-'X';
 211             }
 212         }
 213     }
 214
 215
 216   template <class T>
 217   ASCII_FIELD_DRIVER<T>::ASCII_FIELD_DRIVER(const ASCII_FIELD_DRIVER<T>& other):
 218     GENDRIVER(ASCII_DRIVER),
 219     _mesh(other._mesh),
 220     _support(other._support),
 221     _ptrField(other._ptrField),
 222     _fileName(other._fileName),
 223     _code(other._code),
 224     _direc(other._direc),
 225     _nbComponents(other._nbComponents),
 226     _spaceDimension(other._spaceDimension)
 227   {
 228   }
 229
 230   template <class T>
 231   void ASCII_FIELD_DRIVER<T>::open() throw (MEDEXCEPTION)
 232   {
 233                 if (_file.is_open())
 234                         throw MEDEXCEPTION("ASCII_FIELD_DRIVER::open() : file is already open !");
 235     _file.open(_fileName.c_str(),ofstream::out | ofstream::app);
 236     // for MEDMEMTest_AsciiFieldDriver.cxx:208 :
 237     // must throw because the file is opened
 238     //CPPUNIT_ASSERT_MED_THROW(aDriver1->setFileName("anyfile2"), MEDEXCEPTION);
 239     _status = _file.is_open() ? MED_OPENED : MED_INVALID;
 240   }
 241
 242   template <class T>
 243   void ASCII_FIELD_DRIVER<T>::close()
 244   {
 245     _file.close();
 246     _status = MED_CLOSED;
 247   }
 248
 249   template <class T>
 250   void ASCII_FIELD_DRIVER<T>::read ( void ) throw (MEDEXCEPTION)
 251   {
 252     throw MEDEXCEPTION("ASCII_FIELD_DRIVER::read : Can't read with a WRONLY driver !");
 253   }
 254
 255   template <class T>
 256   GENDRIVER* ASCII_FIELD_DRIVER<T>::copy() const
 257   {
 258     return new ASCII_FIELD_DRIVER<T>(*this);
 259   }
 260
 261   template <class T>
 262   void ASCII_FIELD_DRIVER<T>::write( void ) const throw (MEDEXCEPTION)
 263   {
 264                 if (!_file.is_open())
 265                         throw MEDEXCEPTION("ASCII_FIELD_DRIVER::write : can't write a file that was not opened !");
 266
 267     buildIntroduction();
 268     switch(_spaceDimension)
 269       {
 270       case 2:
 271         {
 272           switch(_code)
 273             {
 274             case 52: //XY
 275               {
 276                 sortAndWrite<2,52>();
 277                 break;
 278               }
 279             case 49: //YX
 280               {
 281                 sortAndWrite<2,49>();
 282                 break;
 283               }
 284             default:
 285               MEDEXCEPTION("ASCII_FIELD_DRIVER : Invalid priority definition");
 286             }
 287           break;
 288         }
 289         case 3:
 290         {
 291           switch(_code)
 292             {
 293             case 228: //XYZ
 294               {
 295                 sortAndWrite<3,228>();
 296                 break;
 297               }
 298             case 216: //XZY
 299               {
 300                 sortAndWrite<3,216>();
 301                 break;
 302               }
 303             case 225://YXZ
 304               {
 305                 sortAndWrite<3,225>();
 306                 break;
 307               }
 308             case 201://YZX
 309               {
 310                 sortAndWrite<3,201>();
 311                 break;
 312               }
 313             case 210://ZXY
 314               {
 315                 sortAndWrite<3,210>();
 316                 break;
 317               }
 318             case 198://ZYX
 319               {
 320                 sortAndWrite<3,198>();
 321                 break;
 322               }
 323              default:
 324               MEDEXCEPTION("ASCII_FIELD_DRIVER : Invalid priority definition");
 325             }
 326           break;
 327           }
 328       default:
 329         MEDEXCEPTION("ASCII_FIELD_DRIVER : Invalid space dimension must be 2 or 3");
 330         }
 331   }
 332
 333   template <class T>
 334   void ASCII_FIELD_DRIVER<T>::buildIntroduction() const
 335   {
 336
 337     int i;
 338     _file  << setiosflags(ios::scientific);
 339     _file << "#TITLE: table " << _ptrField->getName() << " TIME: " << _ptrField->getTime() << " IT: " << _ptrField->getIterationNumber() << endl;
 340     _file << "#COLUMN_TITLES: ";
 341     for(i=0;i<_spaceDimension;i++)
 342       _file << char('X'+i) << " | ";
 343     const std::string *compoNames=_ptrField->getComponentsNames();
 344     for(i=0;i<_nbComponents;i++)
 345       {
 346         if(!compoNames)
 347           _file << compoNames[i];
 348         else
 349           _file << "None";
 350         if(i<_nbComponents-1)
 351           _file << " | ";
 352       }
 353     _file << endl;
 354     _file << "#COLUMN_UNITS: ";
 355     compoNames=_mesh->getCoordinateptr()->getCoordinatesUnits();
 356     for(i=0;i<_spaceDimension;i++)
 357       {
 358         if(!compoNames)
 359           _file << compoNames[i];
 360         else
 361           _file << "None";
 362         _file << " | ";
 363       }
 364     const UNIT *compoUnits=_ptrField->getComponentsUnits();
 365     for(i=0;i<_nbComponents;i++)
 366       {
 367         if(!compoUnits)
 368           _file << compoUnits[i].getName();
 369         else
 370           _file << "None";
 371         if(i<_nbComponents-1)
 372           _file << " | ";
 373       }
 374     _file << endl;
 375   }
 376 }
 377
 378 #include "MEDMEM_Field.hxx"
 379
 380 namespace MEDMEM
 381 {
 382   template <class T>
 383   template<int SPACEDIMENSION, unsigned int SORTSTRATEGY>
 384   void ASCII_FIELD_DRIVER<T>::sortAndWrite() const
 385   {
 386     typedef typename MEDMEM_ArrayInterface<double,NoInterlace,NoGauss>::Array    ArrayDoubleNo;
 387     typedef typename MEDMEM_ArrayInterface<double,FullInterlace,NoGauss>::Array  ArrayDoubleFull;
 388     typedef typename MEDMEM_ArrayInterface<T,NoInterlace,NoGauss>::Array         ArrayNo;
 389     typedef typename MEDMEM_ArrayInterface<T,NoInterlaceByType,NoGauss>::Array   ArrayNoByType;
 390     typedef typename MEDMEM_ArrayInterface<T,FullInterlace,NoGauss>::Array       ArrayFull;
 391
 392     int i,j;
 393     int numberOfValues=_ptrField->getNumberOfValues();
 394     std::list< SDForSorting<T,SPACEDIMENSION,SORTSTRATEGY > > li;
 395     const double * coord;
 396     FIELD<double,FullInterlace> * barycenterField=0;
 397     ArrayDoubleNo * baryArrayTmp = NULL;
 398     double * xyz[SPACEDIMENSION];
 399     bool deallocateXyz=false;
 400
 401     if(_support->getEntity()==MED_EN::MED_NODE) {
 402       if (_support->isOnAllElements()) {
 403
 404         coord=_mesh->getCoordinates(MED_EN::MED_NO_INTERLACE);
 405         for(i=0; i<SPACEDIMENSION; i++)
 406           xyz[i]=(double *)coord+i*numberOfValues;
 407
 408       } else {
 409
 410         coord = _mesh->getCoordinates(MED_EN::MED_FULL_INTERLACE);
 411         const int * nodesNumber=_support->getNumber(MED_EN::MED_ALL_ELEMENTS);
 412         for(i=0; i<SPACEDIMENSION; i++)
 413           xyz[i]=new double[numberOfValues];
 414             deallocateXyz=true;
 415             for(i=0;i<numberOfValues;i++) {
 416               for(j=0;j<SPACEDIMENSION;j++)
 417                 xyz[j][i]=coord[(nodesNumber[i]-1)*SPACEDIMENSION+j];
 418             }
 419       }
 420     } else {
 421
 422       barycenterField = _mesh->getBarycenter(_support);
 423       baryArrayTmp = ArrayConvert
 424         ( *( static_cast<ArrayDoubleFull*>(barycenterField->getArray()) ) );
 425       coord = baryArrayTmp->getPtr();
 426       for(i=0; i<SPACEDIMENSION; i++)
 427         xyz[i]=(double *)(coord+i*numberOfValues);
 428     }
 429
 430     const T * valsToSet;
 431     ArrayFull * tmpArray = NULL;
 432     if ( _ptrField->getInterlacingType() == MED_EN::MED_FULL_INTERLACE )
 433       valsToSet= _ptrField->getValue();
 434     else if ( _ptrField->getInterlacingType() == MED_EN::MED_NO_INTERLACE_BY_TYPE ) {
 435       tmpArray = ArrayConvert
 436         ( *( static_cast<ArrayNoByType*>(_ptrField->getArray()) ) );
 437       valsToSet= tmpArray->getPtr();
 438     }
 439     else {
 440       tmpArray = ArrayConvert
 441         ( *( static_cast<ArrayNo*>(_ptrField->getArray()) ) );
 442       valsToSet= tmpArray->getPtr();
 443     }
 444     double temp[SPACEDIMENSION];
 445     for(i=0;i<numberOfValues;i++) {
 446       for(j=0;j<SPACEDIMENSION;j++)
 447         temp[j]=xyz[j][i];
 448       li.push_back(SDForSorting<T,SPACEDIMENSION,SORTSTRATEGY>(temp,valsToSet+i*_nbComponents,_nbComponents));
 449     }
 450
 451     if (barycenterField) barycenterField->removeReference();
 452     if (baryArrayTmp)    delete baryArrayTmp;
 453     if (tmpArray)        delete tmpArray;
 454
 455     if(deallocateXyz)
 456       for(j=0;j<SPACEDIMENSION;j++)
 457         delete [] xyz[j];
 458
 459     li.sort();
 460     _file << setprecision(PRECISION_IN_ASCII_FILE);
 461     if(_direc==MED_EN::ASCENDING) {
 462       typename std::list< SDForSorting<T,SPACEDIMENSION,SORTSTRATEGY > >::iterator iter;
 463       for(iter=li.begin();iter!=li.end();iter++)
 464         (*iter).writeLine(_file);
 465       _file << endl;
 466     } else if (_direc==MED_EN::DESCENDING) {
 467
 468       typename std::list< SDForSorting<T,SPACEDIMENSION,SORTSTRATEGY > >::reverse_iterator iter;
 469       for(iter=li.rbegin();iter!=li.rend();iter++)
 470         (*iter).writeLine(_file);
 471       _file << endl;
 472     } else
 473       MEDEXCEPTION("ASCII_FIELD_DRIVER : Invalid sort direction");
 474   }
 475
 476   //{
 477     //_nbComponentsForCpyInfo=_nbComponents;
 478     //_ptrField->fillFromAnalytic <TEST<T>::copyInfo3> ();
 479     //_ptrField->fillFromAnalytic<  ASCII_FIELD_DRIVER<T>::copyInfo<SPACEDIMENSION,SORTSTRATEGY> > ();
 480     //li.sort();
 481     //typename std::list< SDForSorting<T,SPACEDIMENSION,SORTSTRATEGY > >::iterator iter;
 482     //_file << setiosflags(ios::scientific) << setprecision(PRECISION_IN_ASCII_FILE);
 483     //for(iter=li.begin();iter!=li.end();iter++)
 484     //  {
 485     //(*iter).writeLine(_file);
 486     // }
 487
 488 //   template <class T>
 489 //   template<int SPACEDIMENSION, unsigned int SORTSTRATEGY>//, std::list< SDForSorting<T,SPACEDIMENSION,SORTSTRATEGY> > lis>
 490 //   void ASCII_FIELD_DRIVER<T>::copyInfo(const double *a,T *b)
 491 //   {
 492 //     //lis.push_back(SDForSorting<T,SPACEDIMENSION,SORTSTRATEGY>(a,b,_nbComponentsForCpyInfo));
 493 //   }
 494
 495 //   template <class T>
 496 //   int ASCII_FIELD_DRIVER<T>::_nbComponentsForCpyInfo=0;
 497 }
 498
 499 #endif