src/ParaMEDLoader/ParaMEDFileMesh.cxx

   1 // Copyright (C) 2007-2023  CEA, EDF
   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, or (at your option) any later version.
   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 // Author : Anthony Geay (EDF R&D)
  20
  21 #include "ParaMEDFileMesh.hxx"
  22 #include "MCAuto.hxx"
  23 #include "MEDFileMesh.hxx"
  24 #include "MEDFileMeshLL.hxx"
  25 #include "MEDLoader.hxx"
  26 #include "MEDFileField1TS.hxx"
  27 #include "MEDFileUtilities.hxx"
  28 #include "MEDFileEntities.hxx"
  29 #include <iostream>
  30 #include <fstream>
  31
  32
  33 // From MEDLOader.cxx TU
  34 extern med_geometry_type typmai3[INTERP_KERNEL::NORM_MAXTYPE];
  35
  36 using namespace MEDCoupling;
  37
  38 void getSingleGeometricType(const std::string& fileName, const std::string& mName, INTERP_KERNEL::NormalizedCellType& geoType)
  39 {
  40   int meshDim, spaceDim;
  41    mcIdType numberOfNodes;
  42    std::vector< std::vector< std::pair<INTERP_KERNEL::NormalizedCellType,int> > > typesDistrib(GetUMeshGlobalInfo(fileName,mName,meshDim,spaceDim,numberOfNodes));
  43    std::size_t numberOfTypesMaxDimension = typesDistrib[0].size();
  44    if(numberOfTypesMaxDimension != 1)
  45      throw INTERP_KERNEL::Exception("ParaMEDFileMesh : only mesh with single geometrical type are supported with given distribution !");
  46    geoType = typesDistrib[0][0].first;
  47 }
  48
  49 void checkDistribution(const MPI_Comm& com, mcIdType totalNumberOfElements, const std::vector<mcIdType>& distrib)
  50 {
  51   mcIdType nbEltsInDistribLoc = distrib.size();
  52   mcIdType nbEltsInDistribTot = -1;
  53 #ifdef HAVE_MPI
  54   MPI_Allreduce(&nbEltsInDistribLoc, &nbEltsInDistribTot, 1, MPI_LONG, MPI_SUM, com);
  55 #else
  56       throw INTERP_KERNEL::Exception("not(HAVE_MPI) incompatible with MPI_World_Size>1");
  57 #endif
  58   if(nbEltsInDistribTot != totalNumberOfElements)
  59     {
  60       if(nbEltsInDistribTot > totalNumberOfElements)
  61         throw INTERP_KERNEL::Exception("ParaMEDFileMesh : Some of your partitions overlap each other ! Each element in your distribution vector must appear only once ! ");
  62       else
  63         throw INTERP_KERNEL::Exception("ParaMEDFileMesh : The distribution does not cover the whole mesh ! Each element of the mesh must appear once in your distribution vector ");
  64     }
  65 }
  66
  67
  68 MEDFileMesh *ParaMEDFileMesh::New(int iPart, int nbOfParts, const std::string& fileName, const std::string& mName, int dt, int it, MEDFileMeshReadSelector *mrs)
  69 {
  70   MEDFileUtilities::CheckFileForRead(fileName);
  71   MEDCoupling::MEDCouplingMeshType meshType;
  72   MEDFileUtilities::AutoFid fid(MEDfileOpen(fileName.c_str(),MED_ACC_RDONLY));
  73   int dummy0,dummy1;
  74   std::string dummy2;
  75   MEDCoupling::MEDCouplingAxisType dummy3;
  76   MEDFileMeshL2::GetMeshIdFromName(fid,mName,meshType,dummy3,dummy0,dummy1,dummy2);
  77   switch(meshType)
  78   {
  79     case UNSTRUCTURED:
  80       {
  81         return ParaMEDFileUMesh::New(iPart,nbOfParts,fileName,mName,dt,it,mrs);
  82       }
  83     default:
  84       throw INTERP_KERNEL::Exception("ParaMEDFileMesh::New : only unstructured mesh supported for the moment !");
  85   }
  86 }
  87
  88 MEDFileMesh *ParaMEDFileMesh::ParaNew(int iPart, int nbOfParts, const MPI_Comm& com, const MPI_Info& nfo, const std::string& fileName, const std::string& mName, int dt, int it, MEDFileMeshReadSelector *mrs)
  89 {
  90   MEDFileUtilities::CheckFileForRead(fileName);
  91   MEDCoupling::MEDCouplingMeshType meshType;
  92   MEDFileUtilities::AutoFid fid(MEDfileOpen(fileName.c_str(),MED_ACC_RDONLY));
  93   int dummy0,dummy1;
  94   std::string dummy2;
  95   MEDCoupling::MEDCouplingAxisType dummy3;
  96   MEDFileMeshL2::GetMeshIdFromName(fid,mName,meshType,dummy3,dummy0,dummy1,dummy2);
  97   switch(meshType)
  98   {
  99     case UNSTRUCTURED:
 100       {
 101         return ParaMEDFileUMesh::ParaNew(iPart,nbOfParts,com,nfo,fileName,mName,dt,it,mrs);
 102       }
 103     default:
 104       throw INTERP_KERNEL::Exception("ParaMEDFileMesh::ParaNew : only unstructured mesh supported for the moment !");
 105   }
 106 }
 107
 108 MEDFileUMesh *ParaMEDFileUMesh::New(int iPart, int nbOfParts, const std::string& fileName, const std::string& mName, int dt, int it, MEDFileMeshReadSelector *mrs)
 109 {
 110   MEDFileUtilities::CheckFileForRead(fileName);
 111   MEDFileUtilities::AutoFid fid(MEDfileOpen(fileName.c_str(),MED_ACC_RDONLY));
 112   return ParaMEDFileUMesh::NewPrivate(fid,iPart,nbOfParts,fileName,mName,dt,it,mrs);
 113 }
 114
 115 /*!
 116  * Opens the given file in parallel so that each processor can load a specific part of the mesh \a mName.
 117  * Each processor will load the cells contained in the vector \a distrib (only nodes lying on those cells will be loaded),
 118  * in order to read the entire mesh in parallel (memory consumption is thus distributed among all the processes).
 119  * \param [in] distrib - map defining for each geometric type, the corresponding vector of cells we want to load with c-type indexing (starting from zero).
 120  * Each vector in this map has an independant numerotation, which means on one processor, vectors of different types may contain the same numbers, they will not refer to the same cells
 121  * (the i-th cell of a type A does not correspond to the i-th cell of type B)
 122  * However they have to differ from one processor to another, as to ensure that:
 123  * 1) each processor only loads a unique part of the mesh
 124  * 2) the combined distribution vectors cover the entire mesh
 125  * \param [in] com - group of MPI processes that will read the file
 126  * \param [in] nfo- MPI info object (used to manage MPI routines)
 127  * \param [in] filename - name of the file we want to read
 128  * \param [in] mName - name of the mesh we want to read
 129  * \param [in] dt - order at which to read the mesh
 130  * \param [in] it - iteration at which to read the mesh
 131  * \param [in] mrs - object used to read additional low-level information
 132  * \return MEDFileUMesh* - a new instance of MEDFileUMesh. The
 133  *         caller is to delete this mesh using decrRef() as it is no more needed.
 134  * \throw exception if the mesh contains multiple types of cells
 135  * \throw exception if the partition of the mesh cells defined by \a distrib does not cover the whole mesh
 136  */
 137 MEDFileUMesh *ParaMEDFileUMesh::ParaNew(const std::map<INTERP_KERNEL::NormalizedCellType,std::vector<mcIdType>> &distrib, const MPI_Comm& com, const MPI_Info& nfo, const std::string& fileName, const std::string& mName, int dt, int it, MEDFileMeshReadSelector *mrs)
 138 {
 139   MEDFileUtilities::CheckFileForRead(fileName);
 140 #ifdef HDF5_IS_PARALLEL
 141   MEDFileUtilities::AutoFid fid(MEDparFileOpen(fileName.c_str(),MED_ACC_RDONLY,com,nfo));
 142 #else
 143   MEDFileUtilities::AutoFid fid(MEDfileOpen(fileName.c_str(),MED_ACC_RDONLY));
 144 #endif
 145   return ParaMEDFileUMesh::NewPrivate(fid,com,distrib,fileName,mName,dt,it,mrs);
 146 }
 147
 148
 149 /*!
 150  * Opens the given file in parallel so that each processor can load its part of the mesh \a mName.
 151  * The mesh will be equally and linearly distributed among all processes:
 152  * the list of cells will be divided into \a nbOfParts slices and only slice \a iPart (cells and nodes lying on those cells) will be loaded by the current processor.
 153  * The entire mesh is thus read in parallel and memory consumption is divided among the group of processes.
 154  * \param [in] iPart - part of the mesh that will be loaded
 155  * \param [in] nbOfParts - total number of parts in which to divide the mesh
 156  * \param [in] com - group of MPI processes that will read the file
 157  * \param [in] nfo- MPI info object (used to manage MPI routines)
 158  * \param [in] filename - name of the file we want to read
 159  * \param [in] mName - name of the mesh we want to read
 160  * \param [in] dt - Time order at which to read the mesh
 161  * \param [in] it - Time iteration at which to read the mesh
 162  * \param [in] mrs - object used to read additional low-level information
 163  * \return MEDFileUMesh* - a new instance of MEDFileUMesh. The
 164  *         caller is to delete this mesh using decrRef() as it is no more needed.
 165  */
 166 MEDFileUMesh *ParaMEDFileUMesh::ParaNew(int iPart, int nbOfParts, const MPI_Comm& com, const MPI_Info& nfo, const std::string& fileName, const std::string& mName, int dt, int it, MEDFileMeshReadSelector *mrs)
 167 {
 168   MEDFileUtilities::CheckFileForRead(fileName);
 169 #ifdef HDF5_IS_PARALLEL
 170   MEDFileUtilities::AutoFid fid(MEDparFileOpen(fileName.c_str(),MED_ACC_RDONLY,com,nfo)); // MPI_COMM_WORLD, MPI_INFO_NULL
 171 #else
 172   MEDFileUtilities::AutoFid fid(MEDfileOpen(fileName.c_str(),MED_ACC_RDONLY));
 173 #endif
 174   return ParaMEDFileUMesh::NewPrivate(fid,iPart,nbOfParts,fileName,mName,dt,it,mrs);
 175 }
 176
 177 /*!
 178  * Loads mesh \a mName in parallel using a custom partition of the mesh cells among the processes.
 179  * See ParaMEDFileUMesh::ParaNew for detailed description.
 180  */
 181 MEDFileUMesh *ParaMEDFileUMesh::NewPrivate(med_idt fid, const MPI_Comm& com, const std::map<INTERP_KERNEL::NormalizedCellType,std::vector<mcIdType>>& distrib, const std::string& fileName, const std::string& mName, int dt, int it, MEDFileMeshReadSelector *mrs)
 182 {
 183   MCAuto<MEDFileUMesh> ret;
 184   for(std::map<INTERP_KERNEL::NormalizedCellType,std::vector<mcIdType>>::const_iterator iter=distrib.begin(); iter!= distrib.end(); iter++)
 185     {
 186         med_geometry_type geoMedType(typmai3[iter->first /*current geometric type*/]);
 187         med_bool changement,transformation;
 188         med_int totalNumberOfElements(MEDmeshnEntity(fid,mName.c_str(),dt,it,MED_CELL,geoMedType,MED_CONNECTIVITY,MED_NODAL,&changement,&transformation));
 189         checkDistribution(com,totalNumberOfElements,iter->second /*distrib over this geometric type*/);
 190     }
 191   ret=MEDFileUMesh::LoadPartOfFromUserDistrib(fid,mName,distrib,dt,it,mrs);
 192   return ret.retn();
 193 }
 194
 195 /*!
 196  * Loads mesh \a mName in parallel using a slice partition of the mesh cells among the processes
 197  * See ParaMEDFileUMesh::ParaNew for detailed description.
 198  */
 199 MEDFileUMesh *ParaMEDFileUMesh::NewPrivate(med_idt fid, int iPart, int nbOfParts, const std::string& fileName, const std::string& mName, int dt, int it, MEDFileMeshReadSelector *mrs)
 200 {
 201   MCAuto<MEDFileUMesh> ret;
 202   int meshDim, spaceDim;
 203   mcIdType numberOfNodes;
 204   std::vector< std::vector< std::pair<INTERP_KERNEL::NormalizedCellType,int> > > typesDistrib(GetUMeshGlobalInfo(fileName,mName,meshDim,spaceDim,numberOfNodes));
 205   std::vector<INTERP_KERNEL::NormalizedCellType> types;
 206   std::vector<mcIdType> distrib;
 207   for(std::vector< std::vector< std::pair<INTERP_KERNEL::NormalizedCellType,int> > >::const_iterator it0=typesDistrib.begin();it0!=typesDistrib.end();it0++)
 208     for(std::vector< std::pair<INTERP_KERNEL::NormalizedCellType,int> >::const_iterator it1=(*it0).begin();it1!=(*it0).end();it1++)
 209       {
 210         types.push_back((*it1).first);
 211         mcIdType tmp[3];
 212         DataArray::GetSlice(0,(*it1).second,1,iPart,nbOfParts,tmp[0],tmp[1]);
 213         tmp[2]=1;
 214         distrib.insert(distrib.end(),tmp,tmp+3);
 215       }
 216   ret=MEDFileUMesh::LoadPartOf(fid,mName,types,distrib,dt,it,mrs);
 217   return ret.retn();
 218 }
 219
 220 /*!
 221  * Loads field \a fName laying on mesh \a mName from the filename \a fileName in parallel:
 222  * each processor will load their portion of the field (ie the portion laying on the cells in the vector \a distrib given in the parameters).
 223  * WARNING : this will only load the array of values of the field, additionnal information of the local field such as the number of its tuples might be incorrect
 224  * \param [in] com - group of MPI processes that will read the file
 225  * \param [in] nfo- MPI info object (used to manage MPI routines)
 226  * \param [in] fileName - name of the file containing the field
 227  * \param [in] fName - name of the field we want to load
 228  * \param [in] mName - name of the mesh on which the field is defined
 229  * \param [in] distrib - vector of cells on which we want to load the field \a fName (with c-type indexing, so starting from zero).
 230  * \param [in] dt - Time order at which to read the field
 231  * \param [in] it - Time iteration at which to read the field
 232  * \return MEDFileField1TS* - a new instance of MEDFileField1TS. The
 233  *         caller is to delete it using decrRef() as it is no more needed.
 234  * \throw exception if the field is not of type FLOAT64
 235  * \throw exception if the mesh contains more than one geometric type
 236  * \throw exception if the given distribution does not cover the entire mesh on which the field is defined
 237  */
 238 MEDFileField1TS *ParaMEDFileField1TS::ParaNew(const MPI_Comm& com, const MPI_Info& nfo, const std::string& fileName, const std::string& fName, const std::string& mName, const std::vector<mcIdType>& distrib, TypeOfField loc, int dt, int it)
 239 {
 240   MEDFileUtilities::CheckFileForRead(fileName);
 241 #ifdef HDF5_IS_PARALLEL
 242   MEDFileUtilities::AutoFid fid(MEDparFileOpen(fileName.c_str(),MED_ACC_RDONLY,com,nfo));
 243 #else
 244   MEDFileUtilities::AutoFid fid(MEDfileOpen(fileName.c_str(),MED_ACC_RDONLY));
 245 #endif
 246   return ParaMEDFileField1TS::NewPrivate(fid,com,fName,mName,distrib,loc,dt,it);
 247 }
 248
 249 /*!
 250  * Loads field \a fName in parallel using a custom partition of the mesh cells on which the field is defined among the processes.
 251  * See ParaMEDFileField1TS::ParaNew for detailed description.
 252  */
 253 MEDFileField1TS *ParaMEDFileField1TS::NewPrivate(med_idt fid, const MPI_Comm& com, const std::string& fName, const std::string& mName, const std::vector<mcIdType>& distrib, TypeOfField loc, int dt, int it)
 254 {
 255   INTERP_KERNEL::NormalizedCellType geoType;
 256   getSingleGeometricType(MEDFileWritable::FileNameFromFID(fid),mName,geoType);
 257   if(loc==ON_CELLS) //if distribution is on nodes, no fast way to check it (as a node can be shared by multiple processors)
 258     {
 259       med_geometry_type geoMedType(typmai3[geoType]);
 260       med_bool changement,transformation;
 261       med_int totalNumberOfElements(MEDmeshnEntity(fid,mName.c_str(),dt,it,MED_CELL,geoMedType,MED_FAMILY_NUMBER,MED_NODAL,&changement,&transformation));
 262       checkDistribution(com,totalNumberOfElements,distrib);
 263     }
 264   std::vector<std::pair<TypeOfField,INTERP_KERNEL::NormalizedCellType>> tmp={ {loc, geoType} };
 265   INTERP_KERNEL::AutoCppPtr<MEDFileEntities> entities(MEDFileEntities::BuildFrom(&tmp));
 266   MCAuto<MEDFileAnyTypeField1TS> partFile(MEDFileAnyTypeField1TS::NewAdv(fid,fName,dt,it,entities,distrib));
 267
 268   MCAuto<MEDFileField1TS> ret(MEDCoupling::DynamicCast<MEDFileAnyTypeField1TS,MEDFileField1TS>(partFile));
 269   if(ret.isNotNull())
 270     return ret.retn();
 271   else
 272     throw INTERP_KERNEL::Exception("ParaMEDFileField1TS::ParaNew : only FLOAT64 field supported for the moment !");
 273 }
 274
 275
 276 MEDFileMeshes *ParaMEDFileMeshes::New(int iPart, int nbOfParts, const std::string& fileName)
 277 {
 278   std::vector<std::string> ms(GetMeshNames(fileName));
 279   MCAuto<MEDFileMeshes> ret(MEDFileMeshes::New());
 280   for(std::vector<std::string>::const_iterator it=ms.begin();it!=ms.end();it++)
 281     {
 282       MCAuto<MEDFileMesh> mesh(ParaMEDFileMesh::New(iPart,nbOfParts,fileName,(*it)));
 283       ret->pushMesh(mesh);
 284     }
 285   return ret.retn();
 286 }
 287
 288 MEDFileMeshes *ParaMEDFileMeshes::ParaNew(int iPart, int nbOfParts, const MPI_Comm& com, const MPI_Info& nfo, const std::string& fileName)
 289 {
 290   std::vector<std::string> ms(GetMeshNames(fileName));
 291   MCAuto<MEDFileMeshes> ret(MEDFileMeshes::New());
 292   for(std::vector<std::string>::const_iterator it=ms.begin();it!=ms.end();it++)
 293     {
 294       MCAuto<MEDFileMesh> mesh(ParaMEDFileMesh::ParaNew(iPart,nbOfParts,com,nfo,fileName,(*it)));
 295       ret->pushMesh(mesh);
 296     }
 297   return ret.retn();
 298 }