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

Diff of /trunk/dudley/src/CPPAdapter/MeshAdapterFactory.cpp

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

revision 3636 by caltinay, Thu Oct 28 00:50:41 2010 UTC revision 3637 by caltinay, Fri Oct 21 04:29:33 2011 UTC
# Line 11  Line 11 
11  *  *
12  *******************************************************/  *******************************************************/
13    
   
14  #include "MeshAdapterFactory.h"  #include "MeshAdapterFactory.h"
15  #include "DudleyError.h"  #include "DudleyError.h"
16  extern "C" {  extern "C" {
# Line 43  namespace dudley { Line 42  namespace dudley {
42      NcAtt *attr;      NcAtt *attr;
43      char error_msg[LenErrorMsg_MAX];      char error_msg[LenErrorMsg_MAX];
44      if (! (attr=dataFile->get_att(attr_name)) ) {      if (! (attr=dataFile->get_att(attr_name)) ) {
45        sprintf(error_msg,"Error retrieving integer attribute '%s' from NetCDF file '%s'", attr_name, fName);        sprintf(error_msg,"loadMesh: Error retrieving integer attribute '%s' from NetCDF file '%s'", attr_name, fName);
46        throw DataException(error_msg);        throw DataException(error_msg);
47      }      }
48      int temp = attr->as_int(0);      int temp = attr->as_int(0);
# Line 52  namespace dudley { Line 51  namespace dudley {
51    }    }
52  #endif  #endif
53    
54      inline void cleanupAndThrow(Dudley_Mesh* mesh, Esys_MPIInfo* info, string msg)
55      {
56          Dudley_Mesh_free(mesh);
57          Esys_MPIInfo_free(info);
58          string msgPrefix("loadMesh: NetCDF operation failed - ");
59          throw DataException(msgPrefix+msg);
60      }
61    
62  //   AbstractContinuousDomain* loadMesh(const std::string& fileName)  //   AbstractContinuousDomain* loadMesh(const std::string& fileName)
63    Domain_ptr loadMesh(const std::string& fileName)    Domain_ptr loadMesh(const std::string& fileName)
64    {    {
65  #ifdef USE_NETCDF  #ifdef USE_NETCDF
66      Esys_MPIInfo *mpi_info = Esys_MPIInfo_alloc( MPI_COMM_WORLD );      Esys_MPIInfo *mpi_info = Esys_MPIInfo_alloc( MPI_COMM_WORLD );
     AbstractContinuousDomain* temp;  
67      Dudley_Mesh *mesh_p=NULL;      Dudley_Mesh *mesh_p=NULL;
68      char error_msg[LenErrorMsg_MAX];      char error_msg[LenErrorMsg_MAX];
69    
# Line 77  namespace dudley { Line 83  namespace dudley {
83      // Create the NetCDF file.      // Create the NetCDF file.
84      NcFile dataFile(fName, NcFile::ReadOnly);      NcFile dataFile(fName, NcFile::ReadOnly);
85      if (!dataFile.is_valid()) {      if (!dataFile.is_valid()) {
86        sprintf(error_msg,"loadMesh: Opening file NetCDF %s for reading failed.", fName);        sprintf(error_msg,"loadMesh: Opening NetCDF file '%s' for reading failed.", fName);
87        Dudley_setError(IO_ERROR,error_msg);        Dudley_setError(IO_ERROR,error_msg);
88        Esys_MPIInfo_free( mpi_info );        Esys_MPIInfo_free( mpi_info );
89        throw DataException(error_msg);        throw DataException(error_msg);
90      }      }
91    
92      // Read NetCDF integer attributes      // Read NetCDF integer attributes
93      int mpi_size            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_size");      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");      int mpi_rank                        = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_rank");
95      int numDim              = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numDim");      int numDim                          = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numDim");
96      int numNodes            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numNodes");      int numNodes                        = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numNodes");
97      int num_Elements            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements");      int num_Elements                    = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements");
98      int num_FaceElements        = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements");      int num_FaceElements                = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements");
99      int num_Points          = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Points");      int num_Points                      = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Points");
100      int num_Elements_numNodes       = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements_numNodes");      int num_Elements_numNodes           = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements_numNodes");
101      int Elements_TypeId         = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Elements_TypeId");      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");      int num_FaceElements_numNodes       = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements_numNodes");
103      int FaceElements_TypeId     = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"FaceElements_TypeId");      int FaceElements_TypeId             = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"FaceElements_TypeId");
104      int Points_TypeId           = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Points_TypeId");      int Points_TypeId                   = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Points_TypeId");
105      int num_Tags            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Tags");      int num_Tags                        = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Tags");
106    
107      // Verify size and rank      // Verify size and rank
108      if (mpi_info->size != mpi_size) {      if (mpi_info->size != mpi_size) {
109        sprintf(error_msg, "Error loadMesh - The NetCDF file '%s' can only be read on %d CPUs instead of %d", fName, mpi_size, mpi_info->size);        sprintf(error_msg, "loadMesh: The NetCDF file '%s' can only be read on %d CPUs instead of %d", fName, mpi_size, mpi_info->size);
110        throw DataException(error_msg);        throw DataException(error_msg);
111      }      }
112      if (mpi_info->rank != mpi_rank) {      if (mpi_info->rank != mpi_rank) {
113        sprintf(error_msg, "Error loadMesh - The NetCDF file '%s' should be read on CPU #%d instead of %d", fName, mpi_rank, mpi_info->rank);        sprintf(error_msg, "loadMesh: The NetCDF file '%s' should be read on CPU #%d instead of %d", fName, mpi_rank, mpi_info->rank);
114        throw DataException(error_msg);        throw DataException(error_msg);
115      }      }
116    
117      // Read mesh name      // Read mesh name
118      if (! (attr=dataFile.get_att("Name")) ) {      if (! (attr=dataFile.get_att("Name")) ) {
119        sprintf(error_msg,"Error retrieving mesh name from NetCDF file '%s'", fName);        sprintf(error_msg,"loadMesh: Error retrieving mesh name from NetCDF file '%s'", fName);
120        throw DataException(error_msg);        throw DataException(error_msg);
121      }      }
122      char *name = attr->as_string(0);      char *name = attr->as_string(0);
123      delete attr;      delete attr;
124    
125      string msgPrefix("Error in loadMesh: NetCDF operation failed - ");      TMPMEMFREE(fName);
126    
127      /* allocate mesh */      /* allocate mesh */
128      mesh_p = Dudley_Mesh_alloc(name,numDim,mpi_info);      mesh_p = Dudley_Mesh_alloc(name,numDim,mpi_info);
# Line 126  namespace dudley { Line 132  namespace dudley {
132          Dudley_NodeFile_allocTable(mesh_p->Nodes, numNodes);          Dudley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
133          // Nodes_Id          // Nodes_Id
134          if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )          if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
135            throw DataException(msgPrefix+"get_var(Nodes_Id)");              cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_Id)");
136          if (! nc_var_temp->get(&mesh_p->Nodes->Id[0], numNodes) ) {          if (! nc_var_temp->get(&mesh_p->Nodes->Id[0], numNodes) )
137            TMPMEMFREE(mesh_p->Nodes->Id);              cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_Id)");
           throw DataException("get(Nodes_Id)");  
         }  
138          // Nodes_Tag          // Nodes_Tag
139          if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )          if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )
140            throw DataException("get_var(Nodes_Tag)");              cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_Tag)");
141          if (! nc_var_temp->get(&mesh_p->Nodes->Tag[0], numNodes) ) {          if (! nc_var_temp->get(&mesh_p->Nodes->Tag[0], numNodes) )
142            TMPMEMFREE(mesh_p->Nodes->Tag);              cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_Tag)");
           throw DataException("get(Nodes_Tag)");  
         }  
143          // Nodes_gDOF          // Nodes_gDOF
144          if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )          if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )
145            throw DataException("get_var(Nodes_gDOF)");              cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_gDOF)");
146          if (! nc_var_temp->get(&mesh_p->Nodes->globalDegreesOfFreedom[0], numNodes) ) {          if (! nc_var_temp->get(&mesh_p->Nodes->globalDegreesOfFreedom[0], numNodes) )
147            TMPMEMFREE(mesh_p->Nodes->globalDegreesOfFreedom);              cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_gDOF)");
           throw DataException("get(Nodes_gDOF)");  
         }  
148          // Nodes_gNI          // Nodes_gNI
149          if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )          if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )
150            throw DataException("get_var(Nodes_gNI)");              cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_gNI)");
151          if (! nc_var_temp->get(&mesh_p->Nodes->globalNodesIndex[0], numNodes) ) {          if (! nc_var_temp->get(&mesh_p->Nodes->globalNodesIndex[0], numNodes) )
152            TMPMEMFREE(mesh_p->Nodes->globalNodesIndex);              cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_gNI)");
           throw DataException("get(Nodes_gNI)");  
         }  
153          // Nodes_grDfI          // Nodes_grDfI
154          if (! ( nc_var_temp = dataFile.get_var("Nodes_grDfI")) )          if (! ( nc_var_temp = dataFile.get_var("Nodes_grDfI")) )
155            throw DataException("get_var(Nodes_grDfI)");              cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_grDfI)");
156          if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedDOFIndex[0], numNodes) ) {          if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedDOFIndex[0], numNodes) )
157            TMPMEMFREE(mesh_p->Nodes->globalReducedDOFIndex);              cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_grDfI)");
           throw DataException("get(Nodes_grDfI)");  
         }  
158          // Nodes_grNI          // Nodes_grNI
159          if (! ( nc_var_temp = dataFile.get_var("Nodes_grNI")) )          if (! ( nc_var_temp = dataFile.get_var("Nodes_grNI")) )
160            throw DataException("get_var(Nodes_grNI)");              cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_grNI)");
161          if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedNodesIndex[0], numNodes) ) {          if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedNodesIndex[0], numNodes) )
162            TMPMEMFREE(mesh_p->Nodes->globalReducedNodesIndex);              cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_grNI)");
           throw DataException("get(Nodes_grNI)");  
         }  
163          // Nodes_Coordinates          // Nodes_Coordinates
164          if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates"))) {          if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates")))
165            TMPMEMFREE(mesh_p->Nodes->Coordinates);              cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_Coordinates)");
166            throw DataException("get_var(Nodes_Coordinates)");          if (! nc_var_temp->get(&(mesh_p->Nodes->Coordinates[0]), numNodes, numDim) )
167          }              cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_Coordinates)");
         if (! nc_var_temp->get(&(mesh_p->Nodes->Coordinates[0]), numNodes, numDim) ) {  
           TMPMEMFREE(mesh_p->Nodes->Coordinates);  
           throw DataException("get(Nodes_Coordinates)");  
         }  
         // Nodes_DofDistribution  
         first_DofComponent = TMPMEMALLOC(mpi_size+1,index_t);  
         if (! ( nc_var_temp = dataFile.get_var("Nodes_DofDistribution")) )  
           throw DataException("get_var(Nodes_DofDistribution)");  
         if (! nc_var_temp->get(&first_DofComponent[0], mpi_size+1) ) {  
           throw DataException("get(Nodes_DofDistribution)");  
         }  
   
         // Nodes_NodeDistribution  
         first_NodeComponent = TMPMEMALLOC(mpi_size+1,index_t);  
         if (! ( nc_var_temp = dataFile.get_var("Nodes_NodeDistribution")) )  
           throw DataException("get_var(Nodes_NodeDistribution)");  
         if (! nc_var_temp->get(&first_NodeComponent[0], mpi_size+1) ) {  
           throw DataException("get(Nodes_NodeDistribution)");  
         }  
168    
169          /* read elements */          /* read elements */
170          if (Dudley_noError()) {          if (Dudley_noError()) {
171            if (Dudley_noError())  {              mesh_p->Elements=Dudley_ElementFile_alloc((Dudley_ElementTypeId)Elements_TypeId, mpi_info);
172                mesh_p->Elements=Dudley_ElementFile_alloc((Dudley_ElementTypeId)Elements_TypeId, mpi_info);              if (Dudley_noError())
173            }                  Dudley_ElementFile_allocTable(mesh_p->Elements, num_Elements);
174            if (Dudley_noError()) Dudley_ElementFile_allocTable(mesh_p->Elements, num_Elements);              if (Dudley_noError()) {
175            if (Dudley_noError()) {                  mesh_p->Elements->minColor=0;
176                mesh_p->Elements->minColor=0;                  mesh_p->Elements->maxColor=num_Elements-1;
177                mesh_p->Elements->maxColor=num_Elements-1;                  if (num_Elements>0) {
               if (num_Elements>0) {  
178                     // Elements_Id                     // Elements_Id
179                     if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )                     if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
180                       throw DataException("get_var(Elements_Id)");                         cleanupAndThrow(mesh_p, mpi_info, "get_var(Elements_Id)");
181                     if (! nc_var_temp->get(&mesh_p->Elements->Id[0], num_Elements) ) {                     if (! nc_var_temp->get(&mesh_p->Elements->Id[0], num_Elements) )
182                       TMPMEMFREE(mesh_p->Elements->Id);                         cleanupAndThrow(mesh_p, mpi_info, "get(Elements_Id)");
                      throw DataException("get(Elements_Id)");  
                    }  
183                     // Elements_Tag                     // Elements_Tag
184                     if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )                     if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
185                       throw DataException("get_var(Elements_Tag)");                         cleanupAndThrow(mesh_p, mpi_info, "get_var(Elements_Tag)");
186                     if (! nc_var_temp->get(&mesh_p->Elements->Tag[0], num_Elements) ) {                     if (! nc_var_temp->get(&mesh_p->Elements->Tag[0], num_Elements) )
187                       TMPMEMFREE(mesh_p->Elements->Tag);                         cleanupAndThrow(mesh_p, mpi_info, "get(Elements_Tag)");
                      throw DataException("get(Elements_Tag)");  
                    }  
188                     // Elements_Owner                     // Elements_Owner
189                     if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )                     if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
190                       throw DataException("get_var(Elements_Owner)");                         cleanupAndThrow(mesh_p, mpi_info, "get_var(Elements_Owner)");
191                     if (! nc_var_temp->get(&mesh_p->Elements->Owner[0], num_Elements) ) {                     if (! nc_var_temp->get(&mesh_p->Elements->Owner[0], num_Elements) )
192                       TMPMEMFREE(mesh_p->Elements->Owner);                         cleanupAndThrow(mesh_p, mpi_info, "get(Elements_Owner)");
                      throw DataException("get(Elements_Owner)");  
                    }  
193                     // Elements_Color                     // Elements_Color
194                     if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )                     if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
195                       throw DataException("get_var(Elements_Color)");                         cleanupAndThrow(mesh_p, mpi_info, "get_var(Elements_Color)");
196                     if (! nc_var_temp->get(&mesh_p->Elements->Color[0], num_Elements) ) {                     if (! nc_var_temp->get(&mesh_p->Elements->Color[0], num_Elements) )
197                       TMPMEMFREE(mesh_p->Elements->Color);                         cleanupAndThrow(mesh_p, mpi_info, "get(Elements_Color)");
                      throw DataException("get(Elements_Color)");  
                    }  
198                     // Elements_Nodes                     // Elements_Nodes
199                     int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);                     int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);
200                     if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {                     if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
201                       TMPMEMFREE(mesh_p->Elements->Nodes);                         TMPMEMFREE(Elements_Nodes);
202                       throw DataException("get_var(Elements_Nodes)");                         cleanupAndThrow(mesh_p, mpi_info, "get_var(Elements_Nodes)");
203                     }                     }
204                     if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {                     if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {
205                       TMPMEMFREE(Elements_Nodes);                         TMPMEMFREE(Elements_Nodes);
206                       throw DataException("get(Elements_Nodes)");                         cleanupAndThrow(mesh_p, mpi_info, "get(Elements_Nodes)");
207                     }                     }
208    
209                     // Copy temp array into mesh_p->Elements->Nodes                     // Copy temp array into mesh_p->Elements->Nodes
210                     for (int i=0; i<num_Elements; i++) {                     for (int i=0; i<num_Elements; i++) {
211                         for (int j=0; j<num_Elements_numNodes; j++) {                         for (int j=0; j<num_Elements_numNodes; j++) {
212                             mesh_p->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]                             mesh_p->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
213                             = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];                                  = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
214                         }                         }
215                     }                     }
216                     TMPMEMFREE(Elements_Nodes);                     TMPMEMFREE(Elements_Nodes);
217                                    } /* num_Elements>0 */
218                } /* num_Elements>0 */              }
219            }          }
     }  
220    
221          /* get the face elements */          /* get the face elements */
222          if (Dudley_noError()) {          if (Dudley_noError()) {
223            if (Dudley_noError())  {              mesh_p->FaceElements=Dudley_ElementFile_alloc((Dudley_ElementTypeId)FaceElements_TypeId, mpi_info);
224                mesh_p->FaceElements=Dudley_ElementFile_alloc((Dudley_ElementTypeId)FaceElements_TypeId, mpi_info);              if (Dudley_noError())
225            }                  Dudley_ElementFile_allocTable(mesh_p->FaceElements, num_FaceElements);
226            if (Dudley_noError()) Dudley_ElementFile_allocTable(mesh_p->FaceElements, num_FaceElements);              if (Dudley_noError()) {
227            if (Dudley_noError()) {                  mesh_p->FaceElements->minColor=0;
228                mesh_p->FaceElements->minColor=0;                  mesh_p->FaceElements->maxColor=num_FaceElements-1;
229                mesh_p->FaceElements->maxColor=num_FaceElements-1;                  if (num_FaceElements>0) {
               if (num_FaceElements>0) {  
230                     // FaceElements_Id                     // FaceElements_Id
231                     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )                     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
232                       throw DataException("get_var(FaceElements_Id)");                         cleanupAndThrow(mesh_p, mpi_info, "get_var(FaceElements_Id)");
233                     if (! nc_var_temp->get(&mesh_p->FaceElements->Id[0], num_FaceElements) ) {                     if (! nc_var_temp->get(&mesh_p->FaceElements->Id[0], num_FaceElements) )
234                       TMPMEMFREE(mesh_p->FaceElements->Id);                         cleanupAndThrow(mesh_p, mpi_info, "get(FaceElements_Id)");
                      throw DataException("get(FaceElements_Id)");  
                    }  
235                     // FaceElements_Tag                     // FaceElements_Tag
236                     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )                     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
237                       throw DataException("get_var(FaceElements_Tag)");                         cleanupAndThrow(mesh_p, mpi_info, "get_var(FaceElements_Tag)");
238                     if (! nc_var_temp->get(&mesh_p->FaceElements->Tag[0], num_FaceElements) ) {                     if (! nc_var_temp->get(&mesh_p->FaceElements->Tag[0], num_FaceElements) )
239                       TMPMEMFREE(mesh_p->FaceElements->Tag);                         cleanupAndThrow(mesh_p, mpi_info, "get(FaceElements_Tag)");
                      throw DataException("get(FaceElements_Tag)");  
                    }  
240                     // FaceElements_Owner                     // FaceElements_Owner
241                     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )                     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
242                       throw DataException("get_var(FaceElements_Owner)");                         cleanupAndThrow(mesh_p, mpi_info, "get_var(FaceElements_Owner)");
243                     if (! nc_var_temp->get(&mesh_p->FaceElements->Owner[0], num_FaceElements) ) {                     if (! nc_var_temp->get(&mesh_p->FaceElements->Owner[0], num_FaceElements) )
244                       TMPMEMFREE(mesh_p->FaceElements->Owner);                         cleanupAndThrow(mesh_p, mpi_info, "get(FaceElements_Owner)");
                      throw DataException("get(FaceElements_Owner)");  
                    }  
245                     // FaceElements_Color                     // FaceElements_Color
246                     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )                     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
247                       throw DataException("get_var(FaceElements_Color)");                         cleanupAndThrow(mesh_p, mpi_info, "get_var(FaceElements_Color)");
248                     if (! nc_var_temp->get(&mesh_p->FaceElements->Color[0], num_FaceElements) ) {                     if (! nc_var_temp->get(&mesh_p->FaceElements->Color[0], num_FaceElements) )
249                       TMPMEMFREE(mesh_p->FaceElements->Color);                         cleanupAndThrow(mesh_p, mpi_info, "get(FaceElements_Color)");
                      throw DataException("get(FaceElements_Color)");  
                    }  
250                     // FaceElements_Nodes                     // FaceElements_Nodes
251                     int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);                     int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);
252                     if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {                     if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
253                       TMPMEMFREE(mesh_p->FaceElements->Nodes);                         TMPMEMFREE(FaceElements_Nodes);
254                       throw DataException("get_var(FaceElements_Nodes)");                         cleanupAndThrow(mesh_p, mpi_info, "get_var(FaceElements_Nodes)");
255                     }                     }
256                     if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {                     if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
257                       TMPMEMFREE(FaceElements_Nodes);                         TMPMEMFREE(FaceElements_Nodes);
258                       throw DataException("get(FaceElements_Nodes)");                         cleanupAndThrow(mesh_p, mpi_info, "get(FaceElements_Nodes)");
259                     }                     }
260                     // Copy temp array into mesh_p->FaceElements->Nodes                     // Copy temp array into mesh_p->FaceElements->Nodes
261                     for (int i=0; i<num_FaceElements; i++) {                     for (int i=0; i<num_FaceElements; i++) {
262                          for (int j=0; j<num_FaceElements_numNodes; j++) {                         for (int j=0; j<num_FaceElements_numNodes; j++) {
263                              mesh_p->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];                             mesh_p->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
264                          }                         }
265                      }                     }
266                      TMPMEMFREE(FaceElements_Nodes);                     TMPMEMFREE(FaceElements_Nodes);
267                } /* num_FaceElements>0 */                  } /* num_FaceElements>0 */
268            }              }
269      }          }
270    
271          /* get the Points (nodal elements) */          /* get the Points (nodal elements) */
272          if (Dudley_noError()) {          if (Dudley_noError()) {
273            if (Dudley_noError())  {              mesh_p->Points=Dudley_ElementFile_alloc((Dudley_ElementTypeId)Points_TypeId, mpi_info);
274                mesh_p->Points=Dudley_ElementFile_alloc((Dudley_ElementTypeId)Points_TypeId, mpi_info);              if (Dudley_noError())
275            }                  Dudley_ElementFile_allocTable(mesh_p->Points, num_Points);
276            if (Dudley_noError()) Dudley_ElementFile_allocTable(mesh_p->Points, num_Points);              if (Dudley_noError()) {
277            if (Dudley_noError()) {                  mesh_p->Points->minColor=0;
278                mesh_p->Points->minColor=0;                  mesh_p->Points->maxColor=num_Points-1;
279                mesh_p->Points->maxColor=num_Points-1;                  if (num_Points>0) {
280                if (num_Points>0) {                     // Points_Id
281                 // Points_Id                     if (! ( nc_var_temp = dataFile.get_var("Points_Id")))
282                     if (! ( nc_var_temp = dataFile.get_var("Points_Id")) )                         cleanupAndThrow(mesh_p, mpi_info, "get_var(Points_Id)");
283                       throw DataException("get_var(Points_Id)");                     if (! nc_var_temp->get(&mesh_p->Points->Id[0], num_Points))
284                     if (! nc_var_temp->get(&mesh_p->Points->Id[0], num_Points) ) {                         cleanupAndThrow(mesh_p, mpi_info, "get(Points_Id)");
285                       TMPMEMFREE(mesh_p->Points->Id);                     // Points_Tag
286                       throw DataException("get(Points_Id)");                     if (! ( nc_var_temp = dataFile.get_var("Points_Tag")))
287                     }                         cleanupAndThrow(mesh_p, mpi_info, "get_var(Points_Tag)");
288                 // Points_Tag                     if (! nc_var_temp->get(&mesh_p->Points->Tag[0], num_Points))
289                     if (! ( nc_var_temp = dataFile.get_var("Points_Tag")) )                         cleanupAndThrow(mesh_p, mpi_info, "get(Points_Tag)");
290                       throw DataException("get_var(Points_Tag)");                     // Points_Owner
291                     if (! nc_var_temp->get(&mesh_p->Points->Tag[0], num_Points) ) {                     if (! ( nc_var_temp = dataFile.get_var("Points_Owner")))
292                       TMPMEMFREE(mesh_p->Points->Tag);                         cleanupAndThrow(mesh_p, mpi_info, "get_var(Points_Owner)");
293                       throw DataException("get(Points_Tag)");                     if (!nc_var_temp->get(&mesh_p->Points->Owner[0], num_Points))
294                     }                         cleanupAndThrow(mesh_p, mpi_info, "get(Points_Owner)");
295                 // Points_Owner                     // Points_Color
296                     if (! ( nc_var_temp = dataFile.get_var("Points_Owner")) )                     if (! ( nc_var_temp = dataFile.get_var("Points_Color")))
297                       throw DataException("get_var(Points_Owner)");                         cleanupAndThrow(mesh_p, mpi_info, "get_var(Points_Color)");
298                     if (! nc_var_temp->get(&mesh_p->Points->Owner[0], num_Points) ) {                     if (!nc_var_temp->get(&mesh_p->Points->Color[0], num_Points))
299                       TMPMEMFREE(mesh_p->Points->Owner);                         cleanupAndThrow(mesh_p, mpi_info, "get(Points_Color)");
300                       throw DataException("get(Points_Owner)");                     // Points_Nodes
301                     }                     int *Points_Nodes = TMPMEMALLOC(num_Points,int);
                // Points_Color  
                    if (! ( nc_var_temp = dataFile.get_var("Points_Color")) )  
                      throw DataException("get_var(Points_Color)");  
                    if (! nc_var_temp->get(&mesh_p->Points->Color[0], num_Points) ) {  
                      TMPMEMFREE(mesh_p->Points->Color);  
                      throw DataException("get(Points_Color)");  
                    }  
                // Points_Nodes  
            int *Points_Nodes = TMPMEMALLOC(num_Points,int);  
302                     if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {                     if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
303                       TMPMEMFREE(mesh_p->Points->Nodes);                         TMPMEMFREE(Points_Nodes);
304                       throw DataException("get_var(Points_Nodes)");                         cleanupAndThrow(mesh_p, mpi_info, "get_var(Points_Nodes)");
305                     }                     }
306                     if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {                     if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {
307                       TMPMEMFREE(Points_Nodes);                         TMPMEMFREE(Points_Nodes);
308                       throw DataException("get(Points_Nodes)");                         cleanupAndThrow(mesh_p, mpi_info, "get(Points_Nodes)");
309                       }
310                       // Copy temp array into mesh_p->Points->Nodes
311                       for (int i=0; i<num_Points; i++) {
312                           mesh_p->Nodes->Id[mesh_p->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
313                     }                     }
314             // Copy temp array into mesh_p->Points->Nodes                     TMPMEMFREE(Points_Nodes);
315             for (int i=0; i<num_Points; i++) {                  } /* num_Points>0 */
316               mesh_p->Nodes->Id[mesh_p->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];              }
317             }          }
            TMPMEMFREE(Points_Nodes);  
           
               } /* num_Points>0 */  
           }  
         }  
318    
319          /* get the tags */          /* get the tags */
320          if (Dudley_noError()) {          if (Dudley_noError()) {
# Line 375  namespace dudley { Line 322  namespace dudley {
322              // Temp storage to gather node IDs              // Temp storage to gather node IDs
323              int *Tags_keys = TMPMEMALLOC(num_Tags, int);              int *Tags_keys = TMPMEMALLOC(num_Tags, int);
324              char name_temp[4096];              char name_temp[4096];
325          int i;              int i;
326    
327          // Tags_keys              // Tags_keys
328              if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) )              if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) ) {
329                throw DataException("get_var(Tags_keys)");                  TMPMEMFREE(Tags_keys);
330                    cleanupAndThrow(mesh_p, mpi_info, "get_var(Tags_keys)");
331                }
332              if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {              if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {
333                TMPMEMFREE(Tags_keys);                  TMPMEMFREE(Tags_keys);
334                throw DataException("get(Tags_keys)");                  cleanupAndThrow(mesh_p, mpi_info, "get(Tags_keys)");
335              }              }
336          for (i=0; i<num_Tags; i++) {              for (i=0; i<num_Tags; i++) {
337                // Retrieve tag name                // Retrieve tag name
338                sprintf(name_temp, "Tags_name_%d", i);                sprintf(name_temp, "Tags_name_%d", i);
339                if (! (attr=dataFile.get_att(name_temp)) ) {                if (! (attr=dataFile.get_att(name_temp)) ) {
340                  sprintf(error_msg,"Error retrieving tag name from NetCDF file '%s'", fName);                    TMPMEMFREE(Tags_keys);
341                  throw DataException(error_msg);                    sprintf(error_msg,"get_att(%s)", name_temp);
342                      cleanupAndThrow(mesh_p, mpi_info, error_msg);
343                }                }
344                char *name = attr->as_string(0);                char *name = attr->as_string(0);
345                delete attr;                delete attr;
346                Dudley_Mesh_addTagMap(mesh_p, name, Tags_keys[i]);                Dudley_Mesh_addTagMap(mesh_p, name, Tags_keys[i]);
347          }              }
348        }              TMPMEMFREE(Tags_keys);
349      }            }
350            }
351        
352          if (Dudley_noError()) Dudley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent);          if (Dudley_noError()) {
353          TMPMEMFREE(first_DofComponent);              // Nodes_DofDistribution
354          TMPMEMFREE(first_NodeComponent);              first_DofComponent = TMPMEMALLOC(mpi_size+1,index_t);
355                if (! ( nc_var_temp = dataFile.get_var("Nodes_DofDistribution")) ) {
356                    TMPMEMFREE(first_DofComponent);
357                    cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_DofDistribution)");
358                }
359                if (! nc_var_temp->get(&first_DofComponent[0], mpi_size+1) ) {
360                    TMPMEMFREE(first_DofComponent);
361                    cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_DofDistribution)");
362                }
363    
364                // Nodes_NodeDistribution
365                first_NodeComponent = TMPMEMALLOC(mpi_size+1,index_t);
366                if (! ( nc_var_temp = dataFile.get_var("Nodes_NodeDistribution")) ) {
367                    TMPMEMFREE(first_DofComponent);
368                    TMPMEMFREE(first_NodeComponent);
369                    cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_NodeDistribution)");
370                }
371                if (! nc_var_temp->get(&first_NodeComponent[0], mpi_size+1) ) {
372                    TMPMEMFREE(first_DofComponent);
373                    TMPMEMFREE(first_NodeComponent);
374                    cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_NodeDistribution)");
375                }
376                Dudley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent);
377                TMPMEMFREE(first_DofComponent);
378                TMPMEMFREE(first_NodeComponent);
379            }
380    
381      } /* Dudley_noError() after Dudley_Mesh_alloc() */      } /* Dudley_noError() after Dudley_Mesh_alloc() */
382    
383      checkDudleyError();      checkDudleyError();
384      temp=new MeshAdapter(mesh_p);      AbstractContinuousDomain* dom=new MeshAdapter(mesh_p);
385    
386      if (! Dudley_noError()) {      if (! Dudley_noError()) {
387        Dudley_Mesh_free(mesh_p);          Dudley_Mesh_free(mesh_p);
388      }      }
389    
     /* win32 refactor */  
     TMPMEMFREE(fName);  
   
390      blocktimer_increment("LoadMesh()", blocktimer_start);      blocktimer_increment("LoadMesh()", blocktimer_start);
391      return temp->getPtr();      return dom->getPtr();
392  #else  #else
393      throw DataException("Error - loadMesh: is not compiled with NetCDF. Please contact your installation manager.");      throw DataException("loadMesh: not compiled with NetCDF. Please contact your installation manager.");
394  #endif /* USE_NETCDF */  #endif /* USE_NETCDF */
395    }    }
396    
397    Domain_ptr readMesh(const std::string& fileName,    Domain_ptr readMesh(const std::string& fileName,
398                       int integrationOrder,                        int integrationOrder,
399                                       int reducedIntegrationOrder,                        int reducedIntegrationOrder,
400                                       int optimize)                        int optimize)
401    {    {
402      //      //
403      // create a copy of the filename to overcome the non-constness of call      // create a copy of the filename to overcome the non-constness of call
# Line 437  namespace dudley { Line 410  namespace dudley {
410      }      }
411    
412      char *fName = TMPMEMALLOC(fileName.size()+1,char);      char *fName = TMPMEMALLOC(fileName.size()+1,char);
413                
414      strcpy(fName,fileName.c_str());      strcpy(fName,fileName.c_str());
415      double blocktimer_start = blocktimer_time();      double blocktimer_start = blocktimer_time();
416    
# Line 457  namespace dudley { Line 430  namespace dudley {
430                                       int integrationOrder,                                       int integrationOrder,
431                                       int reducedIntegrationOrder,                                       int reducedIntegrationOrder,
432                                       int optimize,                                       int optimize,
433                               int useMacroElements)                                       int useMacroElements)
434    {    {
435      //      //
436      // create a copy of the filename to overcome the non-constness of call      // create a copy of the filename to overcome the non-constness of call
# Line 470  namespace dudley { Line 443  namespace dudley {
443      }      }
444    
445      char *fName = TMPMEMALLOC(fileName.size()+1,char);      char *fName = TMPMEMALLOC(fileName.size()+1,char);
446                
447      strcpy(fName,fileName.c_str());      strcpy(fName,fileName.c_str());
448      double blocktimer_start = blocktimer_time();      double blocktimer_start = blocktimer_time();
449    
# Line 485  namespace dudley { Line 458  namespace dudley {
458      return temp->getPtr();      return temp->getPtr();
459    }    }
460    
   
461    Domain_ptr brick(int n0,int n1,int n2,int order,    Domain_ptr brick(int n0,int n1,int n2,int order,
462                      double l0,double l1,double l2,                     double l0,double l1,double l2,
463                      int periodic0,int periodic1,                     int periodic0,int periodic1,
464                      int periodic2,                     int periodic2,
465                      int integrationOrder,                     int integrationOrder,
466                      int reducedIntegrationOrder,                     int reducedIntegrationOrder,
467                      int useElementsOnFace,                     int useElementsOnFace,
468                      int useFullElementOrder,                     int useFullElementOrder,
469                      int optimize)                     int optimize)
470    {    {
   
   
471      int numElements[]={n0,n1,n2};      int numElements[]={n0,n1,n2};
472      double length[]={l0,l1,l2};      double length[]={l0,l1,l2};
473    
474      if (periodic0 || periodic1) // we don't support periodic boundary conditions      if (periodic0 || periodic1) // we don't support periodic boundary conditions
475      {      {
476      throw DudleyAdapterException("Dudley does not support periodic boundary conditions.");          throw DudleyAdapterException("Dudley does not support periodic boundary conditions.");
477      }      }
478      else if (integrationOrder>3 || reducedIntegrationOrder>1)      else if (integrationOrder>3 || reducedIntegrationOrder>1)
479      {      {
480      throw DudleyAdapterException("Dudley does not support the requested integrationOrders.");          throw DudleyAdapterException("Dudley does not support the requested integrationOrders.");
481      }      }
482      else if (useElementsOnFace || useFullElementOrder)      else if (useElementsOnFace || useFullElementOrder)
483      {      {
484      throw DudleyAdapterException("Dudley does not support useElementsOnFace or useFullElementOrder.");          throw DudleyAdapterException("Dudley does not support useElementsOnFace or useFullElementOrder.");
485      }      }
486      if (order>1)      if (order>1)
487      {      {
488      throw DudleyAdapterException("Dudley does not support element order greater than 1.");          throw DudleyAdapterException("Dudley does not support element order greater than 1.");
489      }      }
490    
491      //      //
492      // linearInterpolation      // linearInterpolation
493      Dudley_Mesh* fMesh=NULL;      Dudley_Mesh* fMesh=NULL;
494    
495      fMesh=Dudley_TriangularMesh_Tet4(numElements,length,integrationOrder,reducedIntegrationOrder,      fMesh=Dudley_TriangularMesh_Tet4(numElements, length, integrationOrder,
496                      (optimize ? TRUE : FALSE)) ;                          reducedIntegrationOrder, (optimize ? TRUE : FALSE));
497    
498      //      //
499      // Convert any dudley errors into a C++ exception      // Convert any dudley errors into a C++ exception
# Line 531  namespace dudley { Line 502  namespace dudley {
502      return temp->getPtr();      return temp->getPtr();
503    }    }
504    
505    Domain_ptr  rectangle(int n0,int n1,int order,    Domain_ptr rectangle(int n0,int n1,int order,
506                          double l0, double l1,                         double l0, double l1,
507                          int periodic0,int periodic1,                         int periodic0,int periodic1,
508                          int integrationOrder,                         int integrationOrder,
509                          int reducedIntegrationOrder,                         int reducedIntegrationOrder,
510                          int useElementsOnFace,                         int useElementsOnFace,
511                          int useFullElementOrder,                         int useFullElementOrder,
512                          int optimize)                         int optimize)
513    {    {
514      int numElements[]={n0,n1};      int numElements[]={n0,n1};
515      double length[]={l0,l1};      double length[]={l0,l1};
516    
517      if (periodic0 || periodic1) // we don't support periodic boundary conditions      if (periodic0 || periodic1) // we don't support periodic boundary conditions
518      {      {
519      throw DudleyAdapterException("Dudley does not support periodic boundary conditions.");          throw DudleyAdapterException("Dudley does not support periodic boundary conditions.");
520      }      }
521      else if (integrationOrder>3 || reducedIntegrationOrder>1)      else if (integrationOrder>3 || reducedIntegrationOrder>1)
522      {      {
523      throw DudleyAdapterException("Dudley does not support the requested integrationOrders.");          throw DudleyAdapterException("Dudley does not support the requested integrationOrders.");
524      }      }
525      else if (useElementsOnFace || useFullElementOrder)      else if (useElementsOnFace || useFullElementOrder)
526      {      {
527      throw DudleyAdapterException("Dudley does not support useElementsOnFace or useFullElementOrder.");          throw DudleyAdapterException("Dudley does not support useElementsOnFace or useFullElementOrder.");
528      }      }
529    
530      if (order>1)      if (order>1)
531      {      {
532      throw DudleyAdapterException("Dudley does not support element order greater than 1.");          throw DudleyAdapterException("Dudley does not support element order greater than 1.");
533      }      }
534      Dudley_Mesh* fMesh=Dudley_TriangularMesh_Tri3(numElements, length,integrationOrder,reducedIntegrationOrder,      Dudley_Mesh* fMesh=Dudley_TriangularMesh_Tri3(numElements, length,
535                      (optimize ? TRUE : FALSE));            integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
536      //      //
537      // Convert any DUDLEY errors into a C++ exception      // Convert any dudley errors into a C++ exception
538      checkDudleyError();      checkDudleyError();
539      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
540    
541      return temp->getPtr();      return temp->getPtr();
542    }    }
543    
544    // end of namespace    // end of namespace
545    
546  }  }
547    

Legend:
Removed from v.3636  
changed lines
  Added in v.3637

  ViewVC Help
Powered by ViewVC 1.1.26