src/ParaMEDLoader/ParaMEDFileMesh.cxx

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