/[escript]/branches/trilinos_from_5897/dudley/src/CPPAdapter/MeshAdapterFactory.cpp
ViewVC logotype

Diff of /branches/trilinos_from_5897/dudley/src/CPPAdapter/MeshAdapterFactory.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 5962 by caltinay, Fri Feb 5 03:37:49 2016 UTC revision 5963 by caltinay, Mon Feb 22 06:59:27 2016 UTC
# Line 17  Line 17 
17  #define ESNEEDPYTHON  #define ESNEEDPYTHON
18  #include "esysUtils/first.h"  #include "esysUtils/first.h"
19    
   
20  #include "MeshAdapterFactory.h"  #include "MeshAdapterFactory.h"
21  #include "DudleyError.h"  #include "DudleyError.h"
 #include "esysUtils/blocktimer.h"  
22  #include "dudley/Dudley.h"  #include "dudley/Dudley.h"
23  #include "dudley/Mesh.h"  #include "dudley/Mesh.h"
24  #include "dudley/TriangularMesh.h"  #include "dudley/TriangularMesh.h"
 #ifdef ESYS_MPI  
 #include "esysUtils/Esys_MPI.h"  
 #endif  
25    
26    #include "esysUtils/Esys_MPI.h"
27  #include "escript/SubWorld.h"  #include "escript/SubWorld.h"
28    
29  #ifdef USE_NETCDF  #ifdef USE_NETCDF
# Line 45  using namespace escript; Line 41  using namespace escript;
41  namespace dudley {  namespace dudley {
42    
43  #ifdef USE_NETCDF  #ifdef USE_NETCDF
44    // A convenience method to retrieve an integer attribute from a NetCDF file  // A convenience method to retrieve an integer attribute from a NetCDF file
45    int NetCDF_Get_Int_Attribute(NcFile *dataFile, const std::string &fName, char *attr_name)  template<typename T>
46    {  T ncReadAtt(NcFile* dataFile, const string& fName, const string& attrName)
47      NcAtt *attr;  {
48      char error_msg[LenErrorMsg_MAX];      NcAtt* attr = dataFile->get_att(attrName.c_str());
49      if (! (attr=dataFile->get_att(attr_name)) ) {      if (!attr) {
50        sprintf(error_msg,"loadMesh: Error retrieving integer attribute '%s' from NetCDF file '%s'", attr_name, fName.c_str());          stringstream msg;
51        throw DataException(error_msg);          msg << "loadMesh: Error retrieving integer attribute '" << attrName
52                << "' from NetCDF file '" << fName << "'";
53            throw DudleyAdapterException(msg.str());
54      }      }
55      int temp = attr->as_int(0);      T value = (sizeof(T) > 4 ? attr->as_long(0) : attr->as_int(0));
56      delete attr;      delete attr;
57      return(temp);      return value;
58    }  }
59  #endif  #endif
60    
61    inline void cleanupAndThrow(Dudley_Mesh* mesh, string msg)  inline void cleanupAndThrow(Dudley_Mesh* mesh, string msg)
62    {  {
63        Dudley_Mesh_free(mesh);      Dudley_Mesh_free(mesh);
64        string msgPrefix("loadMesh: NetCDF operation failed - ");      string msgPrefix("loadMesh: NetCDF operation failed - ");
65        throw DataException(msgPrefix+msg);      throw DudleyAdapterException(msgPrefix+msg);
66    }  }
67    
68  //   AbstractContinuousDomain* loadMesh(const std::string& fileName)  Domain_ptr loadMesh(const std::string& fileName)
69    Domain_ptr loadMesh(const std::string& fileName)  {
   {  
70  #ifdef USE_NETCDF  #ifdef USE_NETCDF
71      esysUtils::JMPI mpi_info = esysUtils::makeInfo( MPI_COMM_WORLD );      esysUtils::JMPI mpi_info = esysUtils::makeInfo( MPI_COMM_WORLD );
72      Dudley_Mesh *mesh_p=NULL;      Dudley_Mesh *mesh_p=NULL;
73      char error_msg[LenErrorMsg_MAX];      const string fName(esysUtils::appendRankToFileName(fileName,
74                                            mpi_info->size, mpi_info->rank));
75    
     std::string fName(esysUtils::appendRankToFileName(fileName, mpi_info->size,  
                                                       mpi_info->rank));  
   
     double blocktimer_start = blocktimer_time();  
76      Dudley_resetError();      Dudley_resetError();
77      int *first_DofComponent, *first_NodeComponent;      int *first_DofComponent, *first_NodeComponent;
78    
# Line 90  namespace dudley { Line 84  namespace dudley {
84      // Create the NetCDF file.      // Create the NetCDF file.
85      NcFile dataFile(fName.c_str(), NcFile::ReadOnly);      NcFile dataFile(fName.c_str(), NcFile::ReadOnly);
86      if (!dataFile.is_valid()) {      if (!dataFile.is_valid()) {
87        sprintf(error_msg,"loadMesh: Opening NetCDF file '%s' for reading failed.", fName.c_str());          stringstream msg;
88        Dudley_setError(IO_ERROR,error_msg);          msg << "loadMesh: Opening NetCDF file '" << fName << "' for reading failed.";
89        throw DataException(error_msg);          throw DudleyAdapterException(msg.str());
90      }      }
91    
92      // Read NetCDF integer attributes      // Read NetCDF integer attributes
93      int mpi_size                        = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_size");  
94      int mpi_rank                        = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_rank");      // index_size was only introduced with 64-bit index support so fall back
95      int numDim                          = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numDim");      // to 32 bits if not found.
96      int numNodes                        = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numNodes");      int index_size;
97      int num_Elements                    = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements");      try {
98      int num_FaceElements                = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements");          index_size = ncReadAtt<int>(&dataFile, fName, "index_size");
99      int num_Points                      = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Points");      } catch (DudleyAdapterException& e) {
100      int num_Elements_numNodes           = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements_numNodes");          index_size = 4;
101      int Elements_TypeId                 = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Elements_TypeId");      }
102      int num_FaceElements_numNodes       = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements_numNodes");      // technically we could cast if reading 32-bit data on 64-bit escript
103      int FaceElements_TypeId             = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"FaceElements_TypeId");      // but cost-benefit analysis clearly favours this implementation for now
104      int Points_TypeId                   = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Points_TypeId");      if (sizeof(index_t) != index_size) {
105      int num_Tags                        = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Tags");          throw DudleyAdapterException("loadMesh: size of index types at runtime differ from dump file");
106        }
107    
108        int mpi_size = ncReadAtt<int>(&dataFile, fName, "mpi_size");
109        int mpi_rank = ncReadAtt<int>(&dataFile, fName, "mpi_rank");
110        int numDim = ncReadAtt<int>(&dataFile, fName, "numDim");
111        dim_t numNodes = ncReadAtt<dim_t>(&dataFile, fName, "numNodes");
112        dim_t num_Elements = ncReadAtt<dim_t>(&dataFile, fName, "num_Elements");
113        dim_t num_FaceElements = ncReadAtt<dim_t>(&dataFile, fName, "num_FaceElements");
114        dim_t num_Points = ncReadAtt<dim_t>(&dataFile, fName, "num_Points");
115        int num_Elements_numNodes = ncReadAtt<int>(&dataFile, fName, "num_Elements_numNodes");
116        int Elements_TypeId = ncReadAtt<int>(&dataFile, fName, "Elements_TypeId");
117        int num_FaceElements_numNodes = ncReadAtt<int>(&dataFile, fName, "num_FaceElements_numNodes");
118        int FaceElements_TypeId = ncReadAtt<int>(&dataFile, fName, "FaceElements_TypeId");
119        int Points_TypeId = ncReadAtt<int>(&dataFile, fName, "Points_TypeId");
120        int num_Tags = ncReadAtt<int>(&dataFile, fName, "num_Tags");
121    
122      // Verify size and rank      // Verify size and rank
123      if (mpi_info->size != mpi_size) {      if (mpi_info->size != mpi_size) {
124        sprintf(error_msg, "loadMesh: The NetCDF file '%s' can only be read on %d CPUs instead of %d", fName.c_str(), mpi_size, mpi_info->size);          stringstream msg;
125        throw DataException(error_msg);          msg << "loadMesh: The NetCDF file '" << fName
126                << "' can only be read on " << mpi_size
127                << " CPUs. Currently running: " << mpi_info->size;
128            throw DudleyAdapterException(msg.str());
129      }      }
130      if (mpi_info->rank != mpi_rank) {      if (mpi_info->rank != mpi_rank) {
131        sprintf(error_msg, "loadMesh: The NetCDF file '%s' should be read on CPU #%d instead of %d", fName.c_str(), mpi_rank, mpi_info->rank);          stringstream msg;
132        throw DataException(error_msg);          msg << "loadMesh: The NetCDF file '" << fName
133                << "' should be read on CPU #" << mpi_rank
134                << " and NOT on #" << mpi_info->rank;
135            throw DudleyAdapterException(msg.str());
136      }      }
137    
138      // Read mesh name      // Read mesh name
139      if (! (attr=dataFile.get_att("Name")) ) {      if (! (attr=dataFile.get_att("Name")) ) {
140        sprintf(error_msg,"loadMesh: Error retrieving mesh name from NetCDF file '%s'", fName.c_str());          stringstream msg;
141        throw DataException(error_msg);          msg << "loadMesh: Error retrieving mesh name from NetCDF file '"
142                << fName << "'";
143            throw DudleyAdapterException(msg.str());
144      }      }
145      boost::scoped_array<char> name(attr->as_string(0));      boost::scoped_array<char> name(attr->as_string(0));
146      delete attr;      delete attr;
147    
148      /* allocate mesh */      // allocate mesh
149      mesh_p = Dudley_Mesh_alloc(name.get(), numDim, mpi_info);      mesh_p = Dudley_Mesh_alloc(name.get(), numDim, mpi_info);
150      if (Dudley_noError()) {      if (Dudley_noError()) {
151            // read nodes
         /* read nodes */  
152          Dudley_NodeFile_allocTable(mesh_p->Nodes, numNodes);          Dudley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
153          // Nodes_Id          // Nodes_Id
154          if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )          if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
# Line 347  namespace dudley { Line 363  namespace dudley {
363                sprintf(name_temp, "Tags_name_%d", i);                sprintf(name_temp, "Tags_name_%d", i);
364                if (! (attr=dataFile.get_att(name_temp)) ) {                if (! (attr=dataFile.get_att(name_temp)) ) {
365                    delete[] Tags_keys;                    delete[] Tags_keys;
366                    sprintf(error_msg,"get_att(%s)", name_temp);                    stringstream msg;
367                    cleanupAndThrow(mesh_p, error_msg);                    msg << "get_att(" << name_temp << ")";
368                      cleanupAndThrow(mesh_p, msg.str());
369                }                }
370                boost::scoped_array<char> name(attr->as_string(0));                boost::scoped_array<char> name(attr->as_string(0));
371                delete attr;                delete attr;
# Line 396  namespace dudley { Line 413  namespace dudley {
413          Dudley_Mesh_free(mesh_p);          Dudley_Mesh_free(mesh_p);
414      }      }
415    
     blocktimer_increment("LoadMesh()", blocktimer_start);  
416      return dom->getPtr();      return dom->getPtr();
417  #else  #else
418      throw DataException("loadMesh: not compiled with NetCDF. Please contact your installation manager.");      throw DataException("loadMesh: not compiled with NetCDF. Please contact your installation manager.");
# Line 421  namespace dudley { Line 437  namespace dudley {
437      char *fName = new char[fileName.size()+1];      char *fName = new char[fileName.size()+1];
438                    
439      strcpy(fName,fileName.c_str());      strcpy(fName,fileName.c_str());
     double blocktimer_start = blocktimer_time();  
440    
441      fMesh=Dudley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));      fMesh=Dudley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
442      checkDudleyError();      checkDudleyError();
# Line 429  namespace dudley { Line 444  namespace dudley {
444            
445      delete[] fName;      delete[] fName;
446            
     blocktimer_increment("ReadMesh()", blocktimer_start);  
447      return temp->getPtr();      return temp->getPtr();
448    }    }
449    
# Line 453  namespace dudley { Line 467  namespace dudley {
467      char *fName = new char[fileName.size()+1];      char *fName = new char[fileName.size()+1];
468                    
469      strcpy(fName,fileName.c_str());      strcpy(fName,fileName.c_str());
     double blocktimer_start = blocktimer_time();  
470    
471      fMesh=Dudley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE), (useMacroElements ? TRUE : FALSE));      fMesh=Dudley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE), (useMacroElements ? TRUE : FALSE));
472      checkDudleyError();      checkDudleyError();
# Line 461  namespace dudley { Line 474  namespace dudley {
474            
475      delete[] fName;      delete[] fName;
476            
     blocktimer_increment("ReadGmsh()", blocktimer_start);  
477      return temp->getPtr();      return temp->getPtr();
478    }    }
479    

Legend:
Removed from v.5962  
changed lines
  Added in v.5963

  ViewVC Help
Powered by ViewVC 1.1.26