src/FILTER/field2nodes.cxx

   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 #include<string>
  21 #include<deque>
  22
  23 #include "MEDMEM_Exception.hxx"
  24 #include "MEDMEM_define.hxx"
  25
  26 #include "MEDMEM_Med.hxx"
  27 #include "MEDMEM_Mesh.hxx"
  28 #include "MEDMEM_Family.hxx"
  29 #include "MEDMEM_Support.hxx"
  30 #include "MEDMEM_Field.hxx"
  31
  32 using namespace std;
  33 using namespace MEDMEM;
  34
  35 void usage(char * name)
  36 {
  37   cout << " ERROR ABOUT SYNTAX " << endl ;
  38   cout << "  " << name << " <input med file> <output med file> " << endl ;
  39   exit(-1);
  40 }
  41
  42 int main (int argc, char ** argv) {
  43
  44   string filenameIN ;
  45   string filenameOUT;
  46   set <int> listElements, supset;
  47   set <int>::iterator elemIt, ssIt;
  48
  49   if ( argc == 3 ) {
  50     filenameIN  = argv[1] ;
  51     filenameOUT = argv[2] ;
  52     cout << "-> reading all into the Med file " << filenameIN << " and writing all into the Ensight file " << filenameOUT <<  endl ;
  53
  54     MED myMed(MED_DRIVER,filenameIN) ;
  55
  56     cout << "-> Read all meshes "  ;
  57     int NumberOfMeshes = myMed.getNumberOfMeshes() ;
  58     cout << "( "<<NumberOfMeshes << " ) :" << endl ;
  59     FIELD<double> **volume = new FIELD<double>*[NumberOfMeshes];
  60     deque<string> MeshName = myMed.getMeshNames() ;
  61     for (int i=0; i<NumberOfMeshes; i++) {
  62       volume[i] = NULL;
  63       myMed.getMesh(MeshName[i])->read() ;
  64       cout << "-> Mesh "<<i+1<<", named "<<MeshName[i]<<" is read !" << endl;
  65       int id = myMed.getMesh(MeshName[i])->addDriver(MED_DRIVER,filenameOUT,MeshName[i],MED_EN::MED_ECRI);
  66       myMed.getMesh(MeshName[i])->write(id);
  67       cout << "-> Mesh "<<i+1<<", named "<<MeshName[i]<<" is written !" << endl;
  68     }
  69
  70     myMed.updateSupport() ;
  71
  72     cout << "-> Read all fields " ;
  73     int NumberOfFields = myMed.getNumberOfFields() ;
  74     cout << "( "<<NumberOfFields << " ) :" << endl;
  75     deque<string> FieldName = myMed.getFieldNames() ;
  76     for (int i=0; i<NumberOfFields; i++) {
  77       deque<DT_IT_> FieldIteration = myMed.getFieldIteration(FieldName[i]) ;
  78       cout << "-> Field "<<i+1<<", named "<<FieldName[i] << " :" << endl ;
  79       int NumberOfIteration = FieldIteration.size() ;
  80       cout << "    Number of iteration pair : "<< NumberOfIteration << endl;
  81       for (int j=0; j<NumberOfIteration; j++) {
  82         FIELD_ * myField = myMed.getField(FieldName[i],FieldIteration[j].dt,FieldIteration[j].it) ;
  83
  84         myField->read() ;
  85         cout << "    * Iteration "<<FieldIteration[j].dt<<" and  order number "<<FieldIteration[j].it<<" ) is read !" << endl;
  86
  87         int NumberOfComponents = myField->getNumberOfComponents();
  88         switch(myField->getSupport()->getEntity()){
  89         case MED_CELL:
  90           {
  91             cout << "*************************** CHAMP AUX CELLULES" << endl;
  92             MESH *mesh = myField->getSupport()->getMesh();
  93
  94             // create new support for new field
  95             SUPPORT *newSup;
  96             if( myField->getSupport()->isOnAllElements() )
  97               newSup = new SUPPORT(mesh,"Support",MED_NODE);
  98             else{
  99               int nbe = myField->getSupport()->getNumberOfElements(MED_ALL_ELEMENTS);
 100               const int *numb = myField->getSupport()->getNumber(MED_ALL_ELEMENTS);
 101               list<int> myList;
 102               for(int k=0;k<nbe;k++){
 103                 myList.push_back(numb[k]);
 104                 supset.insert(numb[k]);
 105               }
 106               newSup = mesh->buildSupportOnNodeFromElementList(myList,MED_CELL);
 107             }
 108
 109             // read number of nodes
 110             int NumberOfNodes = newSup->getNumberOfElements(MED_ALL_ELEMENTS);
 111             // calculate reverse connectivity to have the list of elements which contains node i
 112             const int *revC = myField->getSupport()->getMesh()->getReverseConnectivity(MED_NODAL,MED_CELL);
 113             const int *indC = myField->getSupport()->getMesh()->getReverseConnectivityIndex(MED_NODAL,MED_CELL);
 114             int ivol;
 115             // calculate volume field on mesh
 116             for(int k=0;k<NumberOfMeshes;k++)
 117               if(strcmp(mesh->getName().c_str(),MeshName[k].c_str())==0)
 118                 ivol = k;
 119             if( volume[ivol] == NULL )
 120               volume[ivol] = myField->getSupport()->getMesh()->getVolume(myField->getSupport());
 121             if (dynamic_cast<MEDMEM::FIELD<double>*>(myField)){
 122               FIELD<double> *myDField = (MEDMEM::FIELD<double>*)myField;
 123               FIELD<double> *newDField = new FIELD<double>(newSup,NumberOfComponents);
 124               newDField->setName(myField->getName());
 125               newDField->setIterationNumber(FieldIteration[j].dt);
 126               newDField->setOrderNumber(FieldIteration[j].it);
 127               newDField->setTime(myDField->getTime());
 128               double *val = new double[NumberOfComponents];
 129               if( newSup->isOnAllElements() ){
 130                 for (int k=1; k<NumberOfNodes+1; k++){
 131                   // listElements contains elements which contains a node of element i
 132                   listElements.clear();
 133                   for(int l=indC[k-1];l<indC[k];l++){
 134                     // c element contains node k
 135                     int c=revC[l-1];
 136                     listElements.insert(c);
 137                   }
 138
 139                   // calculate field value on node k
 140                   double sigmaV = 0.;
 141                   for(int l=0;l<NumberOfComponents;l++)
 142                     val[l] = 0.;
 143                   for(elemIt=listElements.begin();elemIt!=listElements.end();elemIt++){
 144                     int elem = *elemIt;
 145                     double vol = volume[ivol]->getValueIJ(elem,1);
 146                     if( vol != 0. ){
 147                       sigmaV += 1./vol;
 148                       for(int l=1;l<=NumberOfComponents;l++)
 149                         val[l-1] += myDField->getValueIJ(elem,l)/vol;
 150                     }
 151                   }
 152                   for(int l=1;l<=NumberOfComponents;l++)
 153                     newDField->setValueIJ(k,l,val[l-1]/sigmaV);
 154                 }
 155               }
 156               else{
 157                 const int *numb = newSup->getNumber(MED_ALL_ELEMENTS);
 158                 for (int k=1; k<=NumberOfNodes; k++){
 159                   int kk = numb[k-1];
 160                   // listElements contains elements which contains a node of element i
 161                   listElements.clear();
 162                   for(int l=indC[kk-1];l<indC[kk];l++){
 163                     // c element contains node k
 164                     int c=revC[l-1];
 165                     // add c element only if it is in partial support
 166                     ssIt = supset.find(c);
 167                     if(ssIt!=supset.end())
 168                       listElements.insert(c);
 169                   }
 170
 171                   // calculate field value on node k
 172                   double sigmaV = 0.;
 173                   for(int l=0;l<NumberOfComponents;l++)
 174                     val[l] = 0.;
 175                   for(elemIt=listElements.begin();elemIt!=listElements.end();elemIt++){
 176                     int elem = *elemIt;
 177                     double vol = volume[ivol]->getValueIJ(elem,1);
 178                     if( vol != 0. ){
 179                       sigmaV += 1./vol;
 180                       for(int l=1;l<=NumberOfComponents;l++)
 181                         val[l-1] += myDField->getValueIJ(elem,l)/vol;
 182                     }
 183                   }
 184                   for(int l=1;l<=NumberOfComponents;l++)
 185                     newDField->setValueIJ(k,l,val[l-1]/sigmaV);
 186                 }
 187               }
 188               delete [] val;
 189               int id = newDField->addDriver(MED_DRIVER,filenameOUT,FieldName[i],MED_EN::MED_ECRI);
 190               newDField->write(id);
 191               delete newSup;
 192               delete newDField;
 193               cout << "    * Iteration "<<FieldIteration[j].dt<<" and  order number "<<FieldIteration[j].it<<" ) is written !" << endl;
 194             }
 195             else{
 196               FIELD<int> *myIField = (MEDMEM::FIELD<int>*)myField;
 197               FIELD<int> *newIField = new FIELD<int>(newSup,NumberOfComponents);
 198               newIField->setName(myField->getName());
 199               newIField->setIterationNumber(FieldIteration[j].dt);
 200               newIField->setOrderNumber(FieldIteration[j].it);
 201               newIField->setTime(myIField->getTime());
 202               double *val = new double[NumberOfComponents];
 203               if( newSup->isOnAllElements() ){
 204                 for (int k=1; k<NumberOfNodes+1; k++){
 205                   // listElements contains elements which contains a node of element i
 206                   listElements.clear();
 207                   for(int l=indC[k-1];l<indC[k];l++){
 208                     // c element contains node i
 209                     int c=revC[l-1];
 210                     listElements.insert(c);
 211                   }
 212
 213                   // calculate field value on node k
 214                   double sigmaV = 0.;
 215                   for(int l=0;l<NumberOfComponents;l++)
 216                     val[l] = 0.;
 217                   for(elemIt=listElements.begin();elemIt!=listElements.end();elemIt++){
 218                     int elem = *elemIt;
 219                     double vol = volume[ivol]->getValueIJ(elem,1);
 220                     if( vol != 0. ){
 221                       sigmaV += 1./vol;
 222                       for(int l=1;l<=NumberOfComponents;l++)
 223                         val[l-1] += ((double)myIField->getValueIJ(elem,l))/vol;
 224                     }
 225                   }
 226                   for(int l=1;l<=NumberOfComponents;l++)
 227                     newIField->setValueIJ(k,l,(int)(val[l-1]/sigmaV));
 228                 }
 229               }
 230               else{
 231                 const int *numb = newSup->getNumber(MED_ALL_ELEMENTS);
 232                 for (int k=1; k<=NumberOfNodes; k++){
 233                   int kk = numb[k-1];
 234                   // listElements contains elements which contains a node of element i
 235                   listElements.clear();
 236                   for(int l=indC[kk-1];l<indC[kk];l++){
 237                     // c element contains node k
 238                     int c=revC[l-1];
 239                     // add c element only if it is in partial support
 240                     ssIt = supset.find(c);
 241                     if(ssIt!=supset.end())
 242                       listElements.insert(c);
 243                   }
 244
 245                   // calculate field value on node k
 246                   double sigmaV = 0.;
 247                   for(int l=0;l<NumberOfComponents;l++)
 248                     val[l] = 0.;
 249                   for(elemIt=listElements.begin();elemIt!=listElements.end();elemIt++){
 250                     int elem = *elemIt;
 251                     double vol = volume[ivol]->getValueIJ(elem,1);
 252                     if( vol != 0. ){
 253                       sigmaV += 1./vol;
 254                       for(int l=1;l<=NumberOfComponents;l++)
 255                         val[l-1] += myIField->getValueIJ(elem,l)/vol;
 256                     }
 257                   }
 258                   for(int l=1;l<=NumberOfComponents;l++)
 259                     newIField->setValueIJ(k,l,(int)(val[l-1]/sigmaV));
 260                 }
 261               }
 262               delete [] val;
 263               int id = newIField->addDriver(MED_DRIVER,filenameOUT,FieldName[i],MED_EN::MED_ECRI);
 264               newIField->write(id);
 265               delete newSup;
 266               delete newIField;
 267               cout << "    * Iteration "<<FieldIteration[j].dt<<" and  order number "<<FieldIteration[j].it<<" ) is written !" << endl;
 268             }
 269           }
 270           break;
 271         case MED_FACE:
 272           cout << "*************************** CHAMP AUX FACES" << endl;
 273           break;
 274         case MED_EDGE:
 275           cout << "*************************** CHAMP AUX ARETES" << endl;
 276           break;
 277         case MED_NODE:
 278           cout << "*************************** CHAMP AUX NOEUDS" << endl;
 279           int id = myField->addDriver(MED_DRIVER,filenameOUT,FieldName[i],MED_EN::MED_ECRI);
 280           myField->write(id);
 281           cout << "    * Iteration "<<FieldIteration[j].dt<<" and  order number "<<FieldIteration[j].it<<" ) is written !" << endl;
 282           break;
 283         }
 284
 285       }
 286     }
 287
 288   }
 289   else usage(argv[0]);
 290 }