src/FILTER/field2nodes.cxx

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