← Back to branch summary

~vcs-imports/escript-finley/trunk

« back to all changes in this revision

Viewing changes to dudley/src/CPPAdapter/MeshAdapterFactory.cpp

  • Committer: jfenwick
  • Date: 2010-10-11 01:48:14 UTC
  • Revision ID: svn-v4:77569008-7704-0410-b7a0-a92fef0b09fd:trunk:3259
Merging dudley and scons updates from branches

Show diffs side-by-side

added added

removed removed

Lines of Context:
dudley/src/CPPAdapter/MeshAdapterFactory.cpp
 
1
 
 
2
/*******************************************************
 
3
*
 
4
* Copyright (c) 2003-2010 by University of Queensland
 
5
* Earth Systems Science Computational Center (ESSCC)
 
6
* http://www.uq.edu.au/esscc
 
7
*
 
8
* Primary Business: Queensland, Australia
 
9
* Licensed under the Open Software License version 3.0
 
10
* http://www.opensource.org/licenses/osl-3.0.php
 
11
*
 
12
*******************************************************/
 
13
 
 
14
 
 
15
#ifdef ESYS_MPI
 
16
#include <mpi.h>
 
17
#endif
 
18
#ifdef USE_NETCDF
 
19
#include <netcdfcpp.h>
 
20
#endif
 
21
#include "MeshAdapterFactory.h"
 
22
#include "DudleyError.h"
 
23
extern "C" {
 
24
#include "esysUtils/blocktimer.h"
 
25
#include "dudley/Dudley.h"
 
26
#include "dudley/Mesh.h"
 
27
#include "dudley/TriangularMesh.h"
 
28
 
 
29
}
 
30
 
 
31
#include <boost/python/extract.hpp>
 
32
 
 
33
#include <sstream>
 
34
 
 
35
using namespace std;
 
36
using namespace escript;
 
37
 
 
38
namespace dudley {
 
39
 
 
40
#ifdef USE_NETCDF
 
41
  // A convenience method to retrieve an integer attribute from a NetCDF file
 
42
  int NetCDF_Get_Int_Attribute(NcFile *dataFile, char *fName, char *attr_name) {
 
43
    NcAtt *attr;
 
44
    char error_msg[LenErrorMsg_MAX];
 
45
    if (! (attr=dataFile->get_att(attr_name)) ) {
 
46
      sprintf(error_msg,"Error retrieving integer attribute '%s' from NetCDF file '%s'", attr_name, fName);
 
47
      throw DataException(error_msg);
 
48
    }
 
49
    int temp = attr->as_int(0);
 
50
    delete attr;
 
51
    return(temp);
 
52
  }
 
53
#endif
 
54
 
 
55
//   AbstractContinuousDomain* loadMesh(const std::string& fileName)
 
56
  Domain_ptr loadMesh(const std::string& fileName)
 
57
  {
 
58
#ifdef USE_NETCDF
 
59
    Esys_MPIInfo *mpi_info = Esys_MPIInfo_alloc( MPI_COMM_WORLD );
 
60
    AbstractContinuousDomain* temp;
 
61
    Dudley_Mesh *mesh_p=NULL;
 
62
    char error_msg[LenErrorMsg_MAX];
 
63
 
 
64
    char *fName = Esys_MPI_appendRankToFileName(fileName.c_str(),
 
65
                                                mpi_info->size,
 
66
                                                mpi_info->rank);
 
67
 
 
68
    double blocktimer_start = blocktimer_time();
 
69
    Dudley_resetError();
 
70
    int *first_DofComponent, *first_NodeComponent;
 
71
 
 
72
    // Open NetCDF file for reading
 
73
    NcAtt *attr;
 
74
    NcVar *nc_var_temp;
 
75
    // netCDF error handler
 
76
    NcError err(NcError::silent_nonfatal);
 
77
    // Create the NetCDF file.
 
78
    NcFile dataFile(fName, NcFile::ReadOnly);
 
79
    if (!dataFile.is_valid()) {
 
80
      sprintf(error_msg,"loadMesh: Opening file NetCDF %s for reading failed.", fName);
 
81
      Dudley_setError(IO_ERROR,error_msg);
 
82
      Esys_MPIInfo_free( mpi_info );
 
83
      throw DataException(error_msg);
 
84
    }
 
85
 
 
86
    // Read NetCDF integer attributes
 
87
    int mpi_size                        = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_size");
 
88
    int mpi_rank                        = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_rank");
 
89
    int numDim                          = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numDim");
 
90
    int numNodes                        = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numNodes");
 
91
    int num_Elements                    = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements");
 
92
    int num_FaceElements                = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements");
 
93
    int num_Points                      = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Points");
 
94
    int num_Elements_numNodes           = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements_numNodes");
 
95
    int Elements_TypeId                 = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Elements_TypeId");
 
96
    int num_FaceElements_numNodes       = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements_numNodes");
 
97
    int FaceElements_TypeId             = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"FaceElements_TypeId");
 
98
    int Points_TypeId                   = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Points_TypeId");
 
99
    int num_Tags                        = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Tags");
 
100
 
 
101
    // Verify size and rank
 
102
    if (mpi_info->size != mpi_size) {
 
103
      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);
 
104
      throw DataException(error_msg);
 
105
    }
 
106
    if (mpi_info->rank != mpi_rank) {
 
107
      sprintf(error_msg, "Error loadMesh - The NetCDF file '%s' should be read on CPU #%d instead of %d", fName, mpi_rank, mpi_info->rank);
 
108
      throw DataException(error_msg);
 
109
    }
 
110
 
 
111
    // Read mesh name
 
112
    if (! (attr=dataFile.get_att("Name")) ) {
 
113
      sprintf(error_msg,"Error retrieving mesh name from NetCDF file '%s'", fName);
 
114
      throw DataException(error_msg);
 
115
    }
 
116
    char *name = attr->as_string(0);
 
117
    delete attr;
 
118
 
 
119
    string msgPrefix("Error in loadMesh: NetCDF operation failed - ");
 
120
 
 
121
    /* allocate mesh */
 
122
    mesh_p = Dudley_Mesh_alloc(name,numDim,mpi_info);
 
123
    if (Dudley_noError()) {
 
124
 
 
125
        /* read nodes */
 
126
        Dudley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
 
127
        // Nodes_Id
 
128
        if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
 
129
          throw DataException(msgPrefix+"get_var(Nodes_Id)");
 
130
        if (! nc_var_temp->get(&mesh_p->Nodes->Id[0], numNodes) ) {
 
131
          TMPMEMFREE(mesh_p->Nodes->Id);
 
132
          throw DataException("get(Nodes_Id)");
 
133
        }
 
134
        // Nodes_Tag
 
135
        if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )
 
136
          throw DataException("get_var(Nodes_Tag)");
 
137
        if (! nc_var_temp->get(&mesh_p->Nodes->Tag[0], numNodes) ) {
 
138
          TMPMEMFREE(mesh_p->Nodes->Tag);
 
139
          throw DataException("get(Nodes_Tag)");
 
140
        }
 
141
        // Nodes_gDOF
 
142
        if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )
 
143
          throw DataException("get_var(Nodes_gDOF)");
 
144
        if (! nc_var_temp->get(&mesh_p->Nodes->globalDegreesOfFreedom[0], numNodes) ) {
 
145
          TMPMEMFREE(mesh_p->Nodes->globalDegreesOfFreedom);
 
146
          throw DataException("get(Nodes_gDOF)");
 
147
        }
 
148
        // Nodes_gNI
 
149
        if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )
 
150
          throw DataException("get_var(Nodes_gNI)");
 
151
        if (! nc_var_temp->get(&mesh_p->Nodes->globalNodesIndex[0], numNodes) ) {
 
152
          TMPMEMFREE(mesh_p->Nodes->globalNodesIndex);
 
153
          throw DataException("get(Nodes_gNI)");
 
154
        }
 
155
        // Nodes_grDfI
 
156
        if (! ( nc_var_temp = dataFile.get_var("Nodes_grDfI")) )
 
157
          throw DataException("get_var(Nodes_grDfI)");
 
158
        if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedDOFIndex[0], numNodes) ) {
 
159
          TMPMEMFREE(mesh_p->Nodes->globalReducedDOFIndex);
 
160
          throw DataException("get(Nodes_grDfI)");
 
161
        }
 
162
        // Nodes_grNI
 
163
        if (! ( nc_var_temp = dataFile.get_var("Nodes_grNI")) )
 
164
          throw DataException("get_var(Nodes_grNI)");
 
165
        if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedNodesIndex[0], numNodes) ) {
 
166
          TMPMEMFREE(mesh_p->Nodes->globalReducedNodesIndex);
 
167
          throw DataException("get(Nodes_grNI)");
 
168
        }
 
169
        // Nodes_Coordinates
 
170
        if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates"))) {
 
171
          TMPMEMFREE(mesh_p->Nodes->Coordinates);
 
172
          throw DataException("get_var(Nodes_Coordinates)");
 
173
        }
 
174
        if (! nc_var_temp->get(&(mesh_p->Nodes->Coordinates[0]), numNodes, numDim) ) {
 
175
          TMPMEMFREE(mesh_p->Nodes->Coordinates);
 
176
          throw DataException("get(Nodes_Coordinates)");
 
177
        }
 
178
        // Nodes_DofDistribution
 
179
        first_DofComponent = TMPMEMALLOC(mpi_size+1,index_t);
 
180
        if (! ( nc_var_temp = dataFile.get_var("Nodes_DofDistribution")) )
 
181
          throw DataException("get_var(Nodes_DofDistribution)");
 
182
        if (! nc_var_temp->get(&first_DofComponent[0], mpi_size+1) ) {
 
183
          throw DataException("get(Nodes_DofDistribution)");
 
184
        }
 
185
 
 
186
        // Nodes_NodeDistribution
 
187
        first_NodeComponent = TMPMEMALLOC(mpi_size+1,index_t);
 
188
        if (! ( nc_var_temp = dataFile.get_var("Nodes_NodeDistribution")) )
 
189
          throw DataException("get_var(Nodes_NodeDistribution)");
 
190
        if (! nc_var_temp->get(&first_NodeComponent[0], mpi_size+1) ) {
 
191
          throw DataException("get(Nodes_NodeDistribution)");
 
192
        }
 
193
 
 
194
        /* read elements */
 
195
        if (Dudley_noError()) {
 
196
                  if (Dudley_noError())  {
 
197
                          mesh_p->Elements=Dudley_ElementFile_alloc((Dudley_ElementTypeId)Elements_TypeId, mpi_info);
 
198
                  }
 
199
          if (Dudley_noError()) Dudley_ElementFile_allocTable(mesh_p->Elements, num_Elements);
 
200
                  if (Dudley_noError()) {
 
201
                          mesh_p->Elements->minColor=0;
 
202
                          mesh_p->Elements->maxColor=num_Elements-1;
 
203
                          if (num_Elements>0) {
 
204
                   // Elements_Id
 
205
                   if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
 
206
                     throw DataException("get_var(Elements_Id)");
 
207
                   if (! nc_var_temp->get(&mesh_p->Elements->Id[0], num_Elements) ) {
 
208
                     TMPMEMFREE(mesh_p->Elements->Id);
 
209
                     throw DataException("get(Elements_Id)");
 
210
                   }
 
211
                   // Elements_Tag
 
212
                   if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
 
213
                     throw DataException("get_var(Elements_Tag)");
 
214
                   if (! nc_var_temp->get(&mesh_p->Elements->Tag[0], num_Elements) ) {
 
215
                     TMPMEMFREE(mesh_p->Elements->Tag);
 
216
                     throw DataException("get(Elements_Tag)");
 
217
                   }
 
218
                   // Elements_Owner
 
219
                   if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
 
220
                     throw DataException("get_var(Elements_Owner)");
 
221
                   if (! nc_var_temp->get(&mesh_p->Elements->Owner[0], num_Elements) ) {
 
222
                     TMPMEMFREE(mesh_p->Elements->Owner);
 
223
                     throw DataException("get(Elements_Owner)");
 
224
                   }
 
225
                   // Elements_Color
 
226
                   if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
 
227
                     throw DataException("get_var(Elements_Color)");
 
228
                   if (! nc_var_temp->get(&mesh_p->Elements->Color[0], num_Elements) ) {
 
229
                     TMPMEMFREE(mesh_p->Elements->Color);
 
230
                     throw DataException("get(Elements_Color)");
 
231
                   }
 
232
                   // Elements_Nodes
 
233
                   int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);
 
234
                   if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
 
235
                     TMPMEMFREE(mesh_p->Elements->Nodes);
 
236
                     throw DataException("get_var(Elements_Nodes)");
 
237
                   }
 
238
                   if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {
 
239
                     TMPMEMFREE(Elements_Nodes);
 
240
                     throw DataException("get(Elements_Nodes)");
 
241
                   }
 
242
                   // Copy temp array into mesh_p->Elements->Nodes
 
243
                                   for (int i=0; i<num_Elements; i++) {
 
244
                                           for (int j=0; j<num_Elements_numNodes; j++) {
 
245
                                                   mesh_p->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
 
246
                                                   = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
 
247
                                           }
 
248
                                   }
 
249
                                   TMPMEMFREE(Elements_Nodes);
 
250
                                
 
251
                          } /* num_Elements>0 */
 
252
                  }
 
253
        }
 
254
 
 
255
        /* get the face elements */
 
256
        if (Dudley_noError()) {
 
257
                  if (Dudley_noError())  {
 
258
                          mesh_p->FaceElements=Dudley_ElementFile_alloc((Dudley_ElementTypeId)FaceElements_TypeId, mpi_info);
 
259
                  }
 
260
          if (Dudley_noError()) Dudley_ElementFile_allocTable(mesh_p->FaceElements, num_FaceElements);
 
261
                  if (Dudley_noError()) {
 
262
                          mesh_p->FaceElements->minColor=0;
 
263
                          mesh_p->FaceElements->maxColor=num_FaceElements-1;
 
264
                          if (num_FaceElements>0) {
 
265
                   // FaceElements_Id
 
266
                   if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
 
267
                     throw DataException("get_var(FaceElements_Id)");
 
268
                   if (! nc_var_temp->get(&mesh_p->FaceElements->Id[0], num_FaceElements) ) {
 
269
                     TMPMEMFREE(mesh_p->FaceElements->Id);
 
270
                     throw DataException("get(FaceElements_Id)");
 
271
                   }
 
272
                   // FaceElements_Tag
 
273
                   if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
 
274
                     throw DataException("get_var(FaceElements_Tag)");
 
275
                   if (! nc_var_temp->get(&mesh_p->FaceElements->Tag[0], num_FaceElements) ) {
 
276
                     TMPMEMFREE(mesh_p->FaceElements->Tag);
 
277
                     throw DataException("get(FaceElements_Tag)");
 
278
                   }
 
279
                   // FaceElements_Owner
 
280
                   if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
 
281
                     throw DataException("get_var(FaceElements_Owner)");
 
282
                   if (! nc_var_temp->get(&mesh_p->FaceElements->Owner[0], num_FaceElements) ) {
 
283
                     TMPMEMFREE(mesh_p->FaceElements->Owner);
 
284
                     throw DataException("get(FaceElements_Owner)");
 
285
                   }
 
286
                   // FaceElements_Color
 
287
                   if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
 
288
                     throw DataException("get_var(FaceElements_Color)");
 
289
                   if (! nc_var_temp->get(&mesh_p->FaceElements->Color[0], num_FaceElements) ) {
 
290
                     TMPMEMFREE(mesh_p->FaceElements->Color);
 
291
                     throw DataException("get(FaceElements_Color)");
 
292
                   }
 
293
                   // FaceElements_Nodes
 
294
                   int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);
 
295
                   if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
 
296
                     TMPMEMFREE(mesh_p->FaceElements->Nodes);
 
297
                     throw DataException("get_var(FaceElements_Nodes)");
 
298
                   }
 
299
                   if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
 
300
                     TMPMEMFREE(FaceElements_Nodes);
 
301
                     throw DataException("get(FaceElements_Nodes)");
 
302
                   }
 
303
                   // Copy temp array into mesh_p->FaceElements->Nodes
 
304
                   for (int i=0; i<num_FaceElements; i++) {
 
305
                                                for (int j=0; j<num_FaceElements_numNodes; j++) {
 
306
                                                        mesh_p->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
 
307
                                                }
 
308
                                        }
 
309
                                        TMPMEMFREE(FaceElements_Nodes);
 
310
                          } /* num_FaceElements>0 */
 
311
                  }
 
312
        }
 
313
 
 
314
        /* get the Points (nodal elements) */
 
315
        if (Dudley_noError()) {
 
316
                  if (Dudley_noError())  {
 
317
                          mesh_p->Points=Dudley_ElementFile_alloc((Dudley_ElementTypeId)Points_TypeId, mpi_info);
 
318
                  }
 
319
          if (Dudley_noError()) Dudley_ElementFile_allocTable(mesh_p->Points, num_Points);
 
320
                  if (Dudley_noError()) {
 
321
                          mesh_p->Points->minColor=0;
 
322
                          mesh_p->Points->maxColor=num_Points-1;
 
323
                          if (num_Points>0) {
 
324
                   // Points_Id
 
325
                   if (! ( nc_var_temp = dataFile.get_var("Points_Id")) )
 
326
                     throw DataException("get_var(Points_Id)");
 
327
                   if (! nc_var_temp->get(&mesh_p->Points->Id[0], num_Points) ) {
 
328
                     TMPMEMFREE(mesh_p->Points->Id);
 
329
                     throw DataException("get(Points_Id)");
 
330
                   }
 
331
                   // Points_Tag
 
332
                   if (! ( nc_var_temp = dataFile.get_var("Points_Tag")) )
 
333
                     throw DataException("get_var(Points_Tag)");
 
334
                   if (! nc_var_temp->get(&mesh_p->Points->Tag[0], num_Points) ) {
 
335
                     TMPMEMFREE(mesh_p->Points->Tag);
 
336
                     throw DataException("get(Points_Tag)");
 
337
                   }
 
338
                   // Points_Owner
 
339
                   if (! ( nc_var_temp = dataFile.get_var("Points_Owner")) )
 
340
                     throw DataException("get_var(Points_Owner)");
 
341
                   if (! nc_var_temp->get(&mesh_p->Points->Owner[0], num_Points) ) {
 
342
                     TMPMEMFREE(mesh_p->Points->Owner);
 
343
                     throw DataException("get(Points_Owner)");
 
344
                   }
 
345
                   // Points_Color
 
346
                   if (! ( nc_var_temp = dataFile.get_var("Points_Color")) )
 
347
                     throw DataException("get_var(Points_Color)");
 
348
                   if (! nc_var_temp->get(&mesh_p->Points->Color[0], num_Points) ) {
 
349
                     TMPMEMFREE(mesh_p->Points->Color);
 
350
                     throw DataException("get(Points_Color)");
 
351
                   }
 
352
                   // Points_Nodes
 
353
                   int *Points_Nodes = TMPMEMALLOC(num_Points,int);
 
354
                   if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
 
355
                     TMPMEMFREE(mesh_p->Points->Nodes);
 
356
                     throw DataException("get_var(Points_Nodes)");
 
357
                   }
 
358
                   if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {
 
359
                     TMPMEMFREE(Points_Nodes);
 
360
                     throw DataException("get(Points_Nodes)");
 
361
                   }
 
362
                   // Copy temp array into mesh_p->Points->Nodes
 
363
                   for (int i=0; i<num_Points; i++) {
 
364
                     mesh_p->Nodes->Id[mesh_p->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
 
365
                   }
 
366
                   TMPMEMFREE(Points_Nodes);
 
367
             
 
368
                          } /* num_Points>0 */
 
369
                  }
 
370
                }
 
371
 
 
372
        /* get the tags */
 
373
        if (Dudley_noError()) {
 
374
          if (num_Tags>0) {
 
375
            // Temp storage to gather node IDs
 
376
            int *Tags_keys = TMPMEMALLOC(num_Tags, int);
 
377
            char name_temp[4096];
 
378
            int i;
 
379
 
 
380
            // Tags_keys
 
381
            if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) )
 
382
              throw DataException("get_var(Tags_keys)");
 
383
            if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {
 
384
              TMPMEMFREE(Tags_keys);
 
385
              throw DataException("get(Tags_keys)");
 
386
            }
 
387
            for (i=0; i<num_Tags; i++) {
 
388
              // Retrieve tag name
 
389
              sprintf(name_temp, "Tags_name_%d", i);
 
390
              if (! (attr=dataFile.get_att(name_temp)) ) {
 
391
                sprintf(error_msg,"Error retrieving tag name from NetCDF file '%s'", fName);
 
392
                throw DataException(error_msg);
 
393
              }
 
394
              char *name = attr->as_string(0);
 
395
              delete attr;
 
396
              Dudley_Mesh_addTagMap(mesh_p, name, Tags_keys[i]);
 
397
            }
 
398
          }
 
399
        }
 
400
   
 
401
        if (Dudley_noError()) Dudley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent);
 
402
        TMPMEMFREE(first_DofComponent);
 
403
        TMPMEMFREE(first_NodeComponent);
 
404
 
 
405
    } /* Dudley_noError() after Dudley_Mesh_alloc() */
 
406
 
 
407
    checkDudleyError();
 
408
    temp=new MeshAdapter(mesh_p);
 
409
 
 
410
    if (! Dudley_noError()) {
 
411
      Dudley_Mesh_free(mesh_p);
 
412
    }
 
413
 
 
414
    /* win32 refactor */
 
415
    TMPMEMFREE(fName);
 
416
 
 
417
    blocktimer_increment("LoadMesh()", blocktimer_start);
 
418
    return temp->getPtr();
 
419
#else
 
420
    throw DataException("Error - loadMesh: is not compiled with NetCDF. Please contact your installation manager.");
 
421
#endif /* USE_NETCDF */
 
422
  }
 
423
 
 
424
  Domain_ptr readMesh(const std::string& fileName,
 
425
                                     int integrationOrder,
 
426
                                     int reducedIntegrationOrder,
 
427
                                     int optimize)
 
428
  {
 
429
    //
 
430
    // create a copy of the filename to overcome the non-constness of call
 
431
    // to Dudley_Mesh_read
 
432
    Dudley_Mesh* fMesh=0;
 
433
    // Win32 refactor
 
434
    if( fileName.size() == 0 )
 
435
    {
 
436
       throw DataException("Null file name!");
 
437
    }
 
438
 
 
439
    char *fName = TMPMEMALLOC(fileName.size()+1,char);
 
440
        
 
441
    strcpy(fName,fileName.c_str());
 
442
    double blocktimer_start = blocktimer_time();
 
443
 
 
444
    fMesh=Dudley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
 
445
    checkDudleyError();
 
446
    AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
 
447
    
 
448
    /* win32 refactor */
 
449
    TMPMEMFREE(fName);
 
450
    
 
451
    blocktimer_increment("ReadMesh()", blocktimer_start);
 
452
    return temp->getPtr();
 
453
  }
 
454
 
 
455
  Domain_ptr readGmsh(const std::string& fileName,
 
456
                                     int numDim,
 
457
                                     int integrationOrder,
 
458
                                     int reducedIntegrationOrder,
 
459
                                     int optimize,
 
460
                                     int useMacroElements)
 
461
  {
 
462
    //
 
463
    // create a copy of the filename to overcome the non-constness of call
 
464
    // to Dudley_Mesh_read
 
465
    Dudley_Mesh* fMesh=0;
 
466
    // Win32 refactor
 
467
    if( fileName.size() == 0 )
 
468
    {
 
469
       throw DataException("Null file name!");
 
470
    }
 
471
 
 
472
    char *fName = TMPMEMALLOC(fileName.size()+1,char);
 
473
        
 
474
    strcpy(fName,fileName.c_str());
 
475
    double blocktimer_start = blocktimer_time();
 
476
 
 
477
    fMesh=Dudley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE), (useMacroElements ? TRUE : FALSE));
 
478
    checkDudleyError();
 
479
    AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
 
480
    
 
481
    /* win32 refactor */
 
482
    TMPMEMFREE(fName);
 
483
    
 
484
    blocktimer_increment("ReadGmsh()", blocktimer_start);
 
485
    return temp->getPtr();
 
486
  }
 
487
 
 
488
 
 
489
  Domain_ptr brick(int n0,int n1,int n2,int order,
 
490
                    double l0,double l1,double l2,
 
491
                    int periodic0,int periodic1,
 
492
                    int periodic2,
 
493
                    int integrationOrder,
 
494
                    int reducedIntegrationOrder,
 
495
                    int useElementsOnFace,
 
496
                    int useFullElementOrder,
 
497
                    int optimize)
 
498
  {
 
499
 
 
500
 
 
501
    int numElements[]={n0,n1,n2};
 
502
    double length[]={l0,l1,l2};
 
503
 
 
504
    if (periodic0 || periodic1) // we don't support periodic boundary conditions
 
505
    {
 
506
        throw DudleyAdapterException("Dudley does not support periodic boundary conditions.");
 
507
    }
 
508
    else if (integrationOrder>3 || reducedIntegrationOrder>1)
 
509
    {
 
510
        throw DudleyAdapterException("Dudley does not support the requested integrationOrders.");
 
511
    }
 
512
    else if (useElementsOnFace || useFullElementOrder)
 
513
    {
 
514
        throw DudleyAdapterException("Dudley does not support useElementsOnFace or useFullElementOrder.");
 
515
    }
 
516
    if (order>1)
 
517
    {
 
518
        throw DudleyAdapterException("Dudley does not support element order greater than 1.");
 
519
    }
 
520
    //
 
521
    // linearInterpolation
 
522
    Dudley_Mesh* fMesh=NULL;
 
523
 
 
524
    fMesh=Dudley_TriangularMesh_Tet4(numElements,length,integrationOrder,reducedIntegrationOrder,
 
525
                                        (optimize ? TRUE : FALSE)) ;
 
526
 
 
527
    //
 
528
    // Convert any dudley errors into a C++ exception
 
529
    checkDudleyError();
 
530
    AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
 
531
    return temp->getPtr();
 
532
  }
 
533
 
 
534
  Domain_ptr  rectangle(int n0,int n1,int order,
 
535
                        double l0, double l1,
 
536
                        int periodic0,int periodic1,
 
537
                        int integrationOrder,
 
538
                        int reducedIntegrationOrder,
 
539
                        int useElementsOnFace,
 
540
                        int useFullElementOrder,
 
541
                        int optimize)
 
542
  {
 
543
    int numElements[]={n0,n1};
 
544
    double length[]={l0,l1};
 
545
 
 
546
    if (periodic0 || periodic1) // we don't support periodic boundary conditions
 
547
    {
 
548
        throw DudleyAdapterException("Dudley does not support periodic boundary conditions.");
 
549
    }
 
550
    else if (integrationOrder>3 || reducedIntegrationOrder>1)
 
551
    {
 
552
        throw DudleyAdapterException("Dudley does not support the requested integrationOrders.");
 
553
    }
 
554
    else if (useElementsOnFace || useFullElementOrder)
 
555
    {
 
556
        throw DudleyAdapterException("Dudley does not support useElementsOnFace or useFullElementOrder.");
 
557
    }
 
558
 
 
559
    if (order>1)
 
560
    {
 
561
        throw DudleyAdapterException("Dudley does not support element order greater than 1.");
 
562
    }
 
563
    Dudley_Mesh* fMesh=Dudley_TriangularMesh_Tri3(numElements, length,integrationOrder,reducedIntegrationOrder,
 
564
                                        (optimize ? TRUE : FALSE));
 
565
    //
 
566
    // Convert any DUDLEY errors into a C++ exception
 
567
    checkDudleyError();
 
568
    AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
 
569
    return temp->getPtr();
 
570
  }
 
571
 
 
572
  // end of namespace
 
573
 
 
574
}