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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26