src/MEDSPLITTER/MEDSPLITTER_MESHCollection.H

   1 // Copyright (C) 2007-2012  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 #ifndef __MEDSPLITTER_MESHCOLLECTION_H__
  20 #define __MEDSPLITTER_MESHCOLLECTION_H__
  21
  22 /*! Projects an integer or double field onto a new mesh
  23  *
  24  * \param old_collection initial mesh collection supporting the field
  25  * \param fieldname name of the field to be cast
  26  * \param itnumber time iteration number of the field to be cast
  27  * \param ordernumber inner iteration number of the field to be cast
  28  * \param type_of_template dummy argument passed for separate template instanciation
  29  */
  30
  31 template <class T>
  32 void MESHCollection::castFields (const MESHCollection& old_collection,
  33                                  const string& fieldname, int itnumber, int ordernumber)
  34 {
  35   char fieldchar[80];
  36   strcpy(fieldchar,fieldname.c_str());
  37
  38   int nb_old_domains=old_collection._topology->nbDomain();
  39   int nb_new_domains=_topology->nbDomain();
  40
  41   vector <list<int> > element_array  (_topology->nbDomain());
  42
  43   vector <MEDMEM::FIELD<T> *> old_fields ;
  44
  45   //cout << "MEDSPLITTER -  reading fields from old collection"<<endl;
  46   old_collection.getDriver()->readFields(old_fields,fieldchar, itnumber, ordernumber);
  47
  48   //if (dynamic_cast<MEDMEM::FIELD<int>*> (old_fields[0])==0)
  49   //  old_collection.getDriver()->readFieldsDouble(old_fields,fieldchar, itnumber, ordernumber);
  50   //else
  51   //  old_collection.getDriver()->readFieldsInt(old_fields,fieldchar, itnumber, ordernumber);
  52   //cout <<"MEDSPLITTER - end of read"<<endl;
  53
  54   vector <const MEDMEM::SUPPORT*> old_supports(nb_old_domains);
  55   vector <MEDMEM::SUPPORT*> new_supports(nb_new_domains);
  56   vector <MEDMEM::FIELD<T> *> new_fields(nb_new_domains);
  57
  58   for (int iold = 0; iold < nb_old_domains; iold++)
  59   {
  60     old_supports[iold]=old_fields[iold]->getSupport();
  61   }
  62   for (int inew = 0; inew < nb_new_domains; inew++)
  63   {
  64     new_supports[inew]=new MEDMEM::SUPPORT();
  65   }
  66
  67   //cout << "MEDSPLITTER - casting supports"<<endl;
  68   castSupport(old_collection,old_supports,new_supports);
  69   //cout << "MEDSPLITTER - end of cast"<<endl;
  70
  71   int nb_components = old_fields[0]->getNumberOfComponents();
  72   const string* components_names = old_fields[0]->getComponentsNames();
  73   const string* components_description = old_fields[0]->getComponentsDescriptions();
  74   const string* components_units = old_fields[0]->getMEDComponentsUnits();
  75   if (itnumber != old_fields[0]->getIterationNumber()) {cout << "PB with iteration number"<<endl;exit (1);}
  76   int iteration_number=old_fields[0]->getIterationNumber();
  77   int order_number=old_fields[0]->getOrderNumber();
  78   double time=old_fields[0]->getTime();
  79   bool has_gauss_pts = old_fields[0]->getGaussPresence();
  80   //bool has_gauss_pts=true;
  81
  82   // Creating the fields objects
  83   // Two different procedures are used whether the field contains Gauss points
  84   // or not
  85
  86   vector <MEDMEM::MEDMEM_Array<T, FullInterlaceGaussPolicy>*> medarray (nb_new_domains);
  87
  88   for (int inew=0; inew < nb_new_domains; inew++)
  89   {
  90     if (!has_gauss_pts) {
  91       new_fields[inew] = new MEDMEM::FIELD<T>(new_supports[inew],nb_components);
  92       new_fields[inew]->setName(fieldname);
  93       new_fields[inew]->setComponentsNames(components_names);
  94       new_fields[inew]->setComponentsDescriptions(components_description);
  95       new_fields[inew]->setMEDComponentsUnits(components_units);
  96       new_fields[inew]->setIterationNumber(iteration_number);
  97       new_fields[inew]->setOrderNumber(order_number);
  98       new_fields[inew]->setTime(time);
  99     }
 100     if (has_gauss_pts)
 101     {
 102       new_fields[inew]=new MEDMEM::FIELD<T>();
 103       //copying the structures describing the field
 104       new_fields[inew]->setNumberOfComponents(nb_components);
 105       new_fields[inew]->setSupport(new_supports[inew]);
 106       new_fields[inew]->setName(fieldname);
 107       new_fields[inew]->setComponentsNames(components_names);
 108       new_fields[inew]->setComponentsDescriptions(components_description);
 109       new_fields[inew]->setMEDComponentsUnits(components_units);
 110       new_fields[inew]->setIterationNumber(iteration_number);
 111       new_fields[inew]->setOrderNumber(order_number);
 112       new_fields[inew]->setTime(time);
 113
 114       //counters for the gauss points
 115       //nbtypegeo is the number of geometric types on the field
 116       //nbelgeoc is the count of element for each type
 117       //nbgaussgeo is the number of gauss points for each type
 118       int nbtypegeo=0;
 119       int* nbelgeoc=0;
 120       int* nbgaussgeo=0;
 121
 122       vector<int> gauss_pts_number;
 123       vector<int> nb_elem;
 124
 125       // the GaussLocalization structures are browsed in the old field
 126       // and copied to the new one
 127       // the nbtypegeo counter is incremented so that
 128       // it contains the number of types for which a gauss localization is defined
 129       MED_EN::MESH_ENTITIES::const_iterator currentEntity;
 130       std::list<MED_EN::medGeometryElement>::const_iterator iter;
 131       currentEntity  = MED_EN::meshEntities.find(MED_EN::MED_CELL);
 132       for (iter = (*currentEntity).second.begin();iter != (*currentEntity).second.end(); iter++)
 133       {
 134         // getGaussLocalization throws an exception
 135         // if the GaussLocalization is not available
 136         try
 137         {
 138           MEDMEM::GAUSS_LOCALIZATION<FullInterlace>
 139             gaussloc(old_fields[0]->getGaussLocalization(*iter));
 140           new_fields[inew]->setGaussLocalization(*iter,gaussloc);
 141
 142           nbtypegeo++;
 143           gauss_pts_number.push_back(gaussloc.getNbGauss());
 144
 145           nb_elem.push_back(new_supports[inew]->getNumberOfElements(*iter));
 146         }
 147         catch(...)
 148         {
 149           continue;
 150         }
 151       }
 152
 153       // les tableaux nbelgeoc commencent a 1
 154       nbelgeoc = new int [nbtypegeo+1];
 155       nbgaussgeo= new int [nbtypegeo+1];
 156       int size=0;
 157       nbelgeoc[0]=0;
 158       nbgaussgeo[0]=-1;
 159       for (int i=1; i<=nbtypegeo;i++)
 160       {
 161         size+=nb_elem[i-1];
 162         nbelgeoc [i]=nb_elem[i-1]+nbelgeoc[i-1];
 163         nbgaussgeo[i]=gauss_pts_number[i-1];
 164       }
 165
 166       //a MEDMEM_Array structure is created to contain
 167       //the data on the Gauss points
 168
 169       medarray[inew]=new MEDMEM::MEDMEM_Array<T,FullInterlaceGaussPolicy>
 170         (new_fields[inew]->getNumberOfComponents(), size, nbtypegeo,
 171          static_cast<const int* const>(nbelgeoc),
 172          static_cast<const int* const> (nbgaussgeo));
 173       new_fields[inew]->setArray(medarray[inew]);
 174       //                        delete[] nbelgeoc;
 175       //                        delete[] nbgaussgeo;
 176     }
 177   }
 178
 179   for (int idomain=0; idomain < old_collection._topology->nbDomain(); idomain++)
 180   {
 181     //retrieves the group igroup on domain idomain
 182
 183     const MEDMEM::SUPPORT* support = old_supports[idomain];
 184     int nbelem = support->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
 185     if (nbelem==0) continue;
 186     int* list_of_elems=0;
 187
 188     if (support->isOnAllElements())
 189     {
 190       list_of_elems = new int[nbelem];
 191       for (int i=0; i<nbelem;i++)
 192         list_of_elems[i]=i+1;
 193     }
 194     else
 195       list_of_elems = const_cast<int*> (support->getNumber(MED_EN::MED_ALL_ELEMENTS));
 196
 197     int* array=new int[nbelem];
 198     int* ip=0;
 199     int* local=0;
 200     int* initial_local=0;
 201     int* full_array=0;
 202
 203     int size=0;
 204     MED_EN::medEntityMesh entity = support->getEntity();
 205     switch (entity)
 206     {
 207     case MED_EN::MED_CELL :
 208       ip=new int[nbelem];
 209       local=new int[nbelem];
 210       initial_local=list_of_elems;
 211       old_collection.getTopology()->convertCellToGlobal(idomain,list_of_elems,nbelem,array);
 212       _topology->convertGlobalCellList(array,nbelem,local,ip);
 213       size=nbelem;
 214       break;
 215     case MED_EN::MED_FACE :
 216       old_collection.getTopology()->convertFaceToGlobal(idomain,list_of_elems,nbelem,array);
 217       _topology->convertGlobalFaceListWithTwins(array,nbelem,local,ip,full_array,size);
 218       initial_local=new int[size];
 219       old_collection.getTopology()->convertGlobalFaceList(full_array,size,initial_local,idomain);
 220       delete[] full_array;
 221       break;
 222     case MED_EN::MED_NODE :
 223       old_collection.getTopology()->convertNodeToGlobal(idomain,list_of_elems,nbelem,array);
 224       _topology->convertGlobalNodeListWithTwins(array,nbelem,local,ip,full_array,size);
 225       initial_local=new int[size];
 226       old_collection.getTopology()->convertGlobalNodeList(full_array,size,initial_local,idomain);
 227       delete[] full_array;
 228       break;
 229     }
 230
 231     if (!has_gauss_pts) {
 232       for (int i=0; i<size; i++)
 233       {
 234         for (int j=0; j<nb_components;j++)
 235         {
 236           T value = old_fields[idomain]->getValueIJ(initial_local[i],j+1);
 237
 238           new_fields[ip[i]]->setValueIJ(local[i],j+1, value);
 239         }
 240       }
 241     }
 242     else
 243     {
 244       for (int i=0; i<size; i++)
 245       {
 246         MED_EN::medGeometryElement type =
 247           old_collection._mesh[idomain]->getElementType(entity,initial_local[i]);
 248         int nb_gauss_points=old_fields[idomain]->getNumberOfGaussPoints(type);
 249         for (int j=0; j<nb_components;j++)
 250           for (int k=0;k<nb_gauss_points; k++)
 251           {
 252             T value = old_fields[idomain]->getValueIJK(initial_local[i],j+1,k+1);
 253             medarray[ip[i]]->setIJK(local[i],j+1,k+1, value);
 254           }
 255       }
 256     }
 257     delete[]array;
 258     delete[]ip;
 259     delete[]local;
 260     if (support->isOnAllElements())
 261     {
 262       delete[] list_of_elems;
 263       list_of_elems=0;
 264     }
 265     if (entity==MED_EN::MED_FACE || entity==MED_EN::MED_NODE) delete[] initial_local;
 266   }
 267
 268   retrieveDriver()->writeFields(new_fields,fieldchar);
 269   //if (dynamic_cast<MEDMEM::FIELD<int>*>(new_fields[0])==0)
 270   //  retrieveDriver()->writeFieldsDouble(new_fields,fieldchar);
 271   //else
 272   //  retrieveDriver()->writeFieldsInt(new_fields,fieldchar);
 273   for (int i=0; i<nb_new_domains; i++)
 274   {
 275     new_fields[i]->removeReference();
 276     new_supports[i]->removeReference();
 277     //delete medarray[i];
 278   }
 279   for (unsigned i=0; i<old_fields.size(); i++)
 280   {
 281     cout << "old field deletion" <<endl;
 282     old_fields[i]->removeReference();
 283     //old_supports[i]->removeReference();//as we did not addReference() neither created it
 284   }
 285 }
 286
 287 #endif /*MEDSPLITTER_MESHCOLLECTION_H_*/