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