src/MEDSPLITTER/MEDSPLITTER_MESHCollection.H

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