/[escript]/branches/more_shared_ptrs_from_1812/finley/src/CPPAdapter/MeshAdapterFactory.cpp
ViewVC logotype

Annotation of /branches/more_shared_ptrs_from_1812/finley/src/CPPAdapter/MeshAdapterFactory.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1813 - (hide annotations)
Fri Sep 26 00:58:05 2008 UTC (11 years, 6 months ago) by jfenwick
File size: 32376 byte(s)
Branching to experiment with proxies and shared pointers.

1 ksteube 1312
2     /*******************************************************
3 ksteube 1811 *
4     * Copyright (c) 2003-2008 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 ksteube 1312
14 ksteube 1811
15 ksteube 817 #ifdef PASO_MPI
16     #include <mpi.h>
17     #endif
18 ksteube 1345 #ifdef USE_NETCDF
19     #include <netcdfcpp.h>
20     #endif
21 jgs 203 #include "MeshAdapterFactory.h"
22 jgs 480 #include "FinleyError.h"
23 ksteube 1345 extern "C" {
24     #include "escript/blocktimer.h"
25     }
26 jgs 82
27 jgs 480 #include <boost/python/extract.hpp>
28    
29     #include <sstream>
30    
31 jgs 82 using namespace std;
32     using namespace escript;
33    
34     namespace finley {
35    
36 ksteube 1345 #ifdef USE_NETCDF
37     // A convenience method to retrieve an integer attribute from a NetCDF file
38     int NetCDF_Get_Int_Attribute(NcFile *dataFile, char *fName, char *attr_name) {
39     NcAtt *attr;
40     char error_msg[LenErrorMsg_MAX];
41     if (! (attr=dataFile->get_att(attr_name)) ) {
42     sprintf(error_msg,"Error retrieving integer attribute '%s' from NetCDF file '%s'", attr_name, fName);
43     throw DataException(error_msg);
44     }
45     int temp = attr->as_int(0);
46     delete attr;
47     return(temp);
48     }
49     #endif
50    
51     AbstractContinuousDomain* loadMesh(const std::string& fileName)
52 ksteube 1312 {
53 ksteube 1345 #ifdef USE_NETCDF
54 ksteube 1347 bool optimize=FALSE; // Don't optimize since this would cause problems with Data().dump()
55 ksteube 1345 Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
56     AbstractContinuousDomain* temp;
57     Finley_Mesh *mesh_p=NULL;
58     char error_msg[LenErrorMsg_MAX];
59 ksteube 1312
60 phornby 1628 char *fName = Paso_MPI_appendRankToFileName(fileName.c_str(),
61     mpi_info->size,
62     mpi_info->rank);
63    
64 ksteube 1345 double blocktimer_start = blocktimer_time();
65     Finley_resetError();
66 ksteube 1755 int *first_DofComponent, *first_NodeComponent;
67 ksteube 1345
68     // Open NetCDF file for reading
69     NcAtt *attr;
70     NcVar *nc_var_temp;
71     // netCDF error handler
72     NcError err(NcError::silent_nonfatal);
73     // Create the NetCDF file.
74     NcFile dataFile(fName, NcFile::ReadOnly);
75     if (!dataFile.is_valid()) {
76     sprintf(error_msg,"loadMesh: Opening file NetCDF %s for reading failed.", fName);
77     Finley_setError(IO_ERROR,error_msg);
78     Paso_MPIInfo_free( mpi_info );
79 ksteube 1347 throw DataException(error_msg);
80 ksteube 1345 }
81    
82     // Read NetCDF integer attributes
83 phornby 1634 int mpi_size = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_size");
84     int mpi_rank = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_rank");
85     int numDim = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numDim");
86     int order = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"order");
87     int reduced_order = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"reduced_order");
88     int numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numNodes");
89     int num_Elements = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements");
90     int num_FaceElements = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements");
91     int num_ContactElements = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_ContactElements");
92     int num_Points = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Points");
93     int num_Elements_numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements_numNodes");
94     int Elements_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Elements_TypeId");
95     int num_FaceElements_numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements_numNodes");
96     int FaceElements_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"FaceElements_TypeId");
97     int num_ContactElements_numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_ContactElements_numNodes");
98     int ContactElements_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"ContactElements_TypeId");
99     int Points_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Points_TypeId");
100     int num_Tags = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Tags");
101 ksteube 1345
102 ksteube 1347 // Verify size and rank
103     if (mpi_info->size != mpi_size) {
104     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);
105     throw DataException(error_msg);
106     }
107     if (mpi_info->rank != mpi_rank) {
108     sprintf(error_msg, "Error loadMesh - The NetCDF file '%s' should be read on CPU #%d instead of %d", fName, mpi_rank, mpi_info->rank);
109     throw DataException(error_msg);
110     }
111    
112 ksteube 1345 // Read mesh name
113     if (! (attr=dataFile.get_att("Name")) ) {
114     sprintf(error_msg,"Error retrieving mesh name from NetCDF file '%s'", fName);
115     throw DataException(error_msg);
116     }
117     char *name = attr->as_string(0);
118     delete attr;
119    
120     /* allocate mesh */
121     mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
122     if (Finley_noError()) {
123    
124     /* read nodes */
125     Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
126     // Nodes_Id
127     if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
128     throw DataException("Error - loadMesh:: unable to read Nodes_Id from netCDF file: " + *fName);
129     if (! nc_var_temp->get(&mesh_p->Nodes->Id[0], numNodes) ) {
130 ksteube 1755 TMPMEMFREE(mesh_p->Nodes->Id);
131 ksteube 1345 throw DataException("Error - loadMesh:: unable to recover Nodes_Id from NetCDF file: " + *fName);
132     }
133     // Nodes_Tag
134     if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )
135     throw DataException("Error - loadMesh:: unable to read Nodes_Tag from netCDF file: " + *fName);
136     if (! nc_var_temp->get(&mesh_p->Nodes->Tag[0], numNodes) ) {
137 ksteube 1755 TMPMEMFREE(mesh_p->Nodes->Tag);
138 ksteube 1345 throw DataException("Error - loadMesh:: unable to recover Nodes_Tag from NetCDF file: " + *fName);
139     }
140     // Nodes_gDOF
141     if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )
142     throw DataException("Error - loadMesh:: unable to read Nodes_gDOF from netCDF file: " + *fName);
143     if (! nc_var_temp->get(&mesh_p->Nodes->globalDegreesOfFreedom[0], numNodes) ) {
144 ksteube 1755 TMPMEMFREE(mesh_p->Nodes->globalDegreesOfFreedom);
145 ksteube 1345 throw DataException("Error - loadMesh:: unable to recover Nodes_gDOF from NetCDF file: " + *fName);
146     }
147     // Nodes_gNI
148     if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )
149     throw DataException("Error - loadMesh:: unable to read Nodes_gNI from netCDF file: " + *fName);
150     if (! nc_var_temp->get(&mesh_p->Nodes->globalNodesIndex[0], numNodes) ) {
151 ksteube 1755 TMPMEMFREE(mesh_p->Nodes->globalNodesIndex);
152 ksteube 1345 throw DataException("Error - loadMesh:: unable to recover Nodes_gNI from NetCDF file: " + *fName);
153     }
154     // Nodes_grDfI
155     if (! ( nc_var_temp = dataFile.get_var("Nodes_grDfI")) )
156     throw DataException("Error - loadMesh:: unable to read Nodes_grDfI from netCDF file: " + *fName);
157     if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedDOFIndex[0], numNodes) ) {
158 ksteube 1755 TMPMEMFREE(mesh_p->Nodes->globalReducedDOFIndex);
159 ksteube 1345 throw DataException("Error - loadMesh:: unable to recover Nodes_grDfI from NetCDF file: " + *fName);
160     }
161     // Nodes_grNI
162     if (! ( nc_var_temp = dataFile.get_var("Nodes_grNI")) )
163     throw DataException("Error - loadMesh:: unable to read Nodes_grNI from netCDF file: " + *fName);
164     if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedNodesIndex[0], numNodes) ) {
165 ksteube 1755 TMPMEMFREE(mesh_p->Nodes->globalReducedNodesIndex);
166 ksteube 1345 throw DataException("Error - loadMesh:: unable to recover Nodes_grNI from NetCDF file: " + *fName);
167     }
168     // Nodes_Coordinates
169     if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates"))) {
170 ksteube 1755 TMPMEMFREE(mesh_p->Nodes->Coordinates);
171 ksteube 1345 throw DataException("Error - loadMesh:: unable to read Nodes_Coordinates from netCDF file: " + *fName);
172     }
173     if (! nc_var_temp->get(&(mesh_p->Nodes->Coordinates[0]), numNodes, numDim) ) {
174 ksteube 1755 TMPMEMFREE(mesh_p->Nodes->Coordinates);
175 ksteube 1345 throw DataException("Error - load:: unable to recover Nodes_Coordinates from netCDF file: " + *fName);
176     }
177 ksteube 1347 // Nodes_DofDistribution
178 ksteube 1755 first_DofComponent = TMPMEMALLOC(mpi_size+1,index_t);
179 ksteube 1347 if (! ( nc_var_temp = dataFile.get_var("Nodes_DofDistribution")) )
180     throw DataException("Error - loadMesh:: unable to read Nodes_DofDistribution from netCDF file: " + *fName);
181 ksteube 1755 if (! nc_var_temp->get(&first_DofComponent[0], mpi_size+1) ) {
182 ksteube 1347 throw DataException("Error - loadMesh:: unable to recover Nodes_DofDistribution from NetCDF file: " + *fName);
183     }
184 ksteube 1345
185 ksteube 1755 // Nodes_NodeDistribution
186     first_NodeComponent = TMPMEMALLOC(mpi_size+1,index_t);
187     if (! ( nc_var_temp = dataFile.get_var("Nodes_NodeDistribution")) )
188     throw DataException("Error - loadMesh:: unable to read Nodes_NodeDistribution from netCDF file: " + *fName);
189     if (! nc_var_temp->get(&first_NodeComponent[0], mpi_size+1) ) {
190     throw DataException("Error - loadMesh:: unable to recover Nodes_NodeDistribution from NetCDF file: " + *fName);
191     }
192    
193 ksteube 1345 /* read elements */
194     if (Finley_noError()) {
195 ksteube 1347 mesh_p->Elements=Finley_ElementFile_alloc((ElementTypeId)Elements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
196     if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Elements, num_Elements);
197     mesh_p->Elements->minColor=0;
198     mesh_p->Elements->maxColor=num_Elements-1;
199 ksteube 1346 if (num_Elements>0) {
200 ksteube 1345 if (Finley_noError()) {
201 ksteube 1346 // Elements_Id
202     if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
203     throw DataException("Error - loadMesh:: unable to read Elements_Id from netCDF file: " + *fName);
204     if (! nc_var_temp->get(&mesh_p->Elements->Id[0], num_Elements) ) {
205 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Id);
206 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Elements_Id from NetCDF file: " + *fName);
207     }
208     // Elements_Tag
209     if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
210     throw DataException("Error - loadMesh:: unable to read Elements_Tag from netCDF file: " + *fName);
211     if (! nc_var_temp->get(&mesh_p->Elements->Tag[0], num_Elements) ) {
212 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Tag);
213 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Elements_Tag from NetCDF file: " + *fName);
214     }
215     // Elements_Owner
216     if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
217     throw DataException("Error - loadMesh:: unable to read Elements_Owner from netCDF file: " + *fName);
218     if (! nc_var_temp->get(&mesh_p->Elements->Owner[0], num_Elements) ) {
219 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Owner);
220 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Elements_Owner from NetCDF file: " + *fName);
221     }
222     // Elements_Color
223     if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
224     throw DataException("Error - loadMesh:: unable to read Elements_Color from netCDF file: " + *fName);
225     if (! nc_var_temp->get(&mesh_p->Elements->Color[0], num_Elements) ) {
226 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Color);
227 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Elements_Color from NetCDF file: " + *fName);
228     }
229     // Elements_Nodes
230     int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);
231     if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
232 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Nodes);
233 ksteube 1346 throw DataException("Error - loadMesh:: unable to read Elements_Nodes from netCDF file: " + *fName);
234     }
235     if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {
236 ksteube 1755 TMPMEMFREE(Elements_Nodes);
237 ksteube 1346 throw DataException("Error - load:: unable to recover Elements_Nodes from netCDF file: " + *fName);
238     }
239     // Copy temp array into mesh_p->Elements->Nodes
240     for (int i=0; i<num_Elements; i++) {
241     for (int j=0; j<num_Elements_numNodes; j++) {
242     mesh_p->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
243     = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
244     }
245     }
246     TMPMEMFREE(Elements_Nodes);
247 ksteube 1345 }
248 ksteube 1346 } /* num_Elements>0 */
249 ksteube 1345 }
250    
251     /* get the face elements */
252 ksteube 1346 if (Finley_noError()) {
253 ksteube 1347 mesh_p->FaceElements=Finley_ElementFile_alloc((ElementTypeId)FaceElements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
254     if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->FaceElements, num_FaceElements);
255     mesh_p->FaceElements->minColor=0;
256     mesh_p->FaceElements->maxColor=num_FaceElements-1;
257 ksteube 1346 if (num_FaceElements>0) {
258     if (Finley_noError()) {
259     // FaceElements_Id
260     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
261     throw DataException("Error - loadMesh:: unable to read FaceElements_Id from netCDF file: " + *fName);
262     if (! nc_var_temp->get(&mesh_p->FaceElements->Id[0], num_FaceElements) ) {
263 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Id);
264 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover FaceElements_Id from NetCDF file: " + *fName);
265     }
266     // FaceElements_Tag
267     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
268     throw DataException("Error - loadMesh:: unable to read FaceElements_Tag from netCDF file: " + *fName);
269     if (! nc_var_temp->get(&mesh_p->FaceElements->Tag[0], num_FaceElements) ) {
270 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Tag);
271 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover FaceElements_Tag from NetCDF file: " + *fName);
272     }
273     // FaceElements_Owner
274     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
275     throw DataException("Error - loadMesh:: unable to read FaceElements_Owner from netCDF file: " + *fName);
276     if (! nc_var_temp->get(&mesh_p->FaceElements->Owner[0], num_FaceElements) ) {
277 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Owner);
278 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover FaceElements_Owner from NetCDF file: " + *fName);
279     }
280     // FaceElements_Color
281     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
282     throw DataException("Error - loadMesh:: unable to read FaceElements_Color from netCDF file: " + *fName);
283     if (! nc_var_temp->get(&mesh_p->FaceElements->Color[0], num_FaceElements) ) {
284 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Color);
285 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover FaceElements_Color from NetCDF file: " + *fName);
286     }
287     // FaceElements_Nodes
288     int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);
289     if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
290 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Nodes);
291 ksteube 1346 throw DataException("Error - loadMesh:: unable to read FaceElements_Nodes from netCDF file: " + *fName);
292     }
293     if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
294 ksteube 1755 TMPMEMFREE(FaceElements_Nodes);
295 ksteube 1346 throw DataException("Error - load:: unable to recover FaceElements_Nodes from netCDF file: " + *fName);
296     }
297     // Copy temp array into mesh_p->FaceElements->Nodes
298     for (int i=0; i<num_FaceElements; i++) {
299     for (int j=0; j<num_FaceElements_numNodes; j++) {
300     mesh_p->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)]
301     = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
302     }
303     }
304     TMPMEMFREE(FaceElements_Nodes);
305     }
306     } /* num_FaceElements>0 */
307     }
308 ksteube 1345
309 ksteube 1346 /* get the Contact elements */
310     if (Finley_noError()) {
311 ksteube 1347 mesh_p->ContactElements=Finley_ElementFile_alloc((ElementTypeId)ContactElements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
312     if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->ContactElements, num_ContactElements);
313     mesh_p->ContactElements->minColor=0;
314     mesh_p->ContactElements->maxColor=num_ContactElements-1;
315 ksteube 1346 if (num_ContactElements>0) {
316     if (Finley_noError()) {
317     // ContactElements_Id
318     if (! ( nc_var_temp = dataFile.get_var("ContactElements_Id")) )
319     throw DataException("Error - loadMesh:: unable to read ContactElements_Id from netCDF file: " + *fName);
320     if (! nc_var_temp->get(&mesh_p->ContactElements->Id[0], num_ContactElements) ) {
321 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Id);
322 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover ContactElements_Id from NetCDF file: " + *fName);
323     }
324     // ContactElements_Tag
325     if (! ( nc_var_temp = dataFile.get_var("ContactElements_Tag")) )
326     throw DataException("Error - loadMesh:: unable to read ContactElements_Tag from netCDF file: " + *fName);
327     if (! nc_var_temp->get(&mesh_p->ContactElements->Tag[0], num_ContactElements) ) {
328 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Tag);
329 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover ContactElements_Tag from NetCDF file: " + *fName);
330     }
331     // ContactElements_Owner
332     if (! ( nc_var_temp = dataFile.get_var("ContactElements_Owner")) )
333     throw DataException("Error - loadMesh:: unable to read ContactElements_Owner from netCDF file: " + *fName);
334     if (! nc_var_temp->get(&mesh_p->ContactElements->Owner[0], num_ContactElements) ) {
335 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Owner);
336 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover ContactElements_Owner from NetCDF file: " + *fName);
337     }
338     // ContactElements_Color
339     if (! ( nc_var_temp = dataFile.get_var("ContactElements_Color")) )
340     throw DataException("Error - loadMesh:: unable to read ContactElements_Color from netCDF file: " + *fName);
341     if (! nc_var_temp->get(&mesh_p->ContactElements->Color[0], num_ContactElements) ) {
342 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Color);
343 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover ContactElements_Color from NetCDF file: " + *fName);
344     }
345     // ContactElements_Nodes
346     int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);
347     if (!(nc_var_temp = dataFile.get_var("ContactElements_Nodes"))) {
348 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Nodes);
349 ksteube 1346 throw DataException("Error - loadMesh:: unable to read ContactElements_Nodes from netCDF file: " + *fName);
350     }
351     if (! nc_var_temp->get(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes) ) {
352 ksteube 1755 TMPMEMFREE(ContactElements_Nodes);
353 ksteube 1346 throw DataException("Error - load:: unable to recover ContactElements_Nodes from netCDF file: " + *fName);
354     }
355     // Copy temp array into mesh_p->ContactElements->Nodes
356     for (int i=0; i<num_ContactElements; i++) {
357     for (int j=0; j<num_ContactElements_numNodes; j++) {
358     mesh_p->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]
359     = ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)];
360     }
361     }
362     TMPMEMFREE(ContactElements_Nodes);
363     }
364     } /* num_ContactElements>0 */
365     }
366 ksteube 1345
367 ksteube 1346 /* get the Points (nodal elements) */
368     if (Finley_noError()) {
369 ksteube 1347 mesh_p->Points=Finley_ElementFile_alloc((ElementTypeId)Points_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
370     if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Points, num_Points);
371     mesh_p->Points->minColor=0;
372     mesh_p->Points->maxColor=num_Points-1;
373 ksteube 1346 if (num_Points>0) {
374     if (Finley_noError()) {
375     // Points_Id
376     if (! ( nc_var_temp = dataFile.get_var("Points_Id")) )
377     throw DataException("Error - loadMesh:: unable to read Points_Id from netCDF file: " + *fName);
378     if (! nc_var_temp->get(&mesh_p->Points->Id[0], num_Points) ) {
379 ksteube 1755 TMPMEMFREE(mesh_p->Points->Id);
380 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Points_Id from NetCDF file: " + *fName);
381     }
382     // Points_Tag
383     if (! ( nc_var_temp = dataFile.get_var("Points_Tag")) )
384     throw DataException("Error - loadMesh:: unable to read Points_Tag from netCDF file: " + *fName);
385     if (! nc_var_temp->get(&mesh_p->Points->Tag[0], num_Points) ) {
386 ksteube 1755 TMPMEMFREE(mesh_p->Points->Tag);
387 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Points_Tag from NetCDF file: " + *fName);
388     }
389     // Points_Owner
390     if (! ( nc_var_temp = dataFile.get_var("Points_Owner")) )
391     throw DataException("Error - loadMesh:: unable to read Points_Owner from netCDF file: " + *fName);
392     if (! nc_var_temp->get(&mesh_p->Points->Owner[0], num_Points) ) {
393 ksteube 1755 TMPMEMFREE(mesh_p->Points->Owner);
394 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Points_Owner from NetCDF file: " + *fName);
395     }
396     // Points_Color
397     if (! ( nc_var_temp = dataFile.get_var("Points_Color")) )
398     throw DataException("Error - loadMesh:: unable to read Points_Color from netCDF file: " + *fName);
399     if (! nc_var_temp->get(&mesh_p->Points->Color[0], num_Points) ) {
400 ksteube 1755 TMPMEMFREE(mesh_p->Points->Color);
401 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Points_Color from NetCDF file: " + *fName);
402     }
403     // Points_Nodes
404     int *Points_Nodes = TMPMEMALLOC(num_Points,int);
405     if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
406 ksteube 1755 TMPMEMFREE(mesh_p->Points->Nodes);
407 ksteube 1346 throw DataException("Error - loadMesh:: unable to read Points_Nodes from netCDF file: " + *fName);
408     }
409     if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {
410 ksteube 1755 TMPMEMFREE(Points_Nodes);
411 ksteube 1346 throw DataException("Error - load:: unable to recover Points_Nodes from netCDF file: " + *fName);
412     }
413     // Copy temp array into mesh_p->Points->Nodes
414     for (int i=0; i<num_Points; i++) {
415     mesh_p->Nodes->Id[mesh_p->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
416     }
417     TMPMEMFREE(Points_Nodes);
418     }
419     } /* num_Points>0 */
420     }
421 ksteube 1345
422 ksteube 1347 /* get the tags */
423     if (Finley_noError()) {
424     if (num_Tags>0) {
425     // Temp storage to gather node IDs
426     int *Tags_keys = TMPMEMALLOC(num_Tags, int);
427     char name_temp[4096];
428     int i;
429 ksteube 1345
430 ksteube 1347 // Tags_keys
431     if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) )
432     throw DataException("Error - loadMesh:: unable to read Tags_keys from netCDF file: " + *fName);
433     if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {
434 ksteube 1755 TMPMEMFREE(Tags_keys);
435 ksteube 1347 throw DataException("Error - loadMesh:: unable to recover Tags_keys from NetCDF file: " + *fName);
436     }
437     for (i=0; i<num_Tags; i++) {
438     // Retrieve tag name
439     sprintf(name_temp, "Tags_name_%d", i);
440     if (! (attr=dataFile.get_att(name_temp)) ) {
441     sprintf(error_msg,"Error retrieving tag name from NetCDF file '%s'", fName);
442     throw DataException(error_msg);
443     }
444     char *name = attr->as_string(0);
445     delete attr;
446     Finley_Mesh_addTagMap(mesh_p, name, Tags_keys[i]);
447     }
448     }
449     }
450    
451 ksteube 1345 } /* Finley_noError() after Finley_Mesh_alloc() */
452    
453 ksteube 1755 if (Finley_noError()) Finley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent);
454     TMPMEMFREE(first_DofComponent);
455     TMPMEMFREE(first_NodeComponent);
456 ksteube 1345
457 ksteube 1312 checkFinleyError();
458 ksteube 1345 temp=new MeshAdapter(mesh_p);
459    
460     if (! Finley_noError()) {
461     Finley_Mesh_free(mesh_p);
462     }
463    
464 ksteube 1312 /* win32 refactor */
465     TMPMEMFREE(fName);
466 ksteube 1345
467     blocktimer_increment("LoadMesh()", blocktimer_start);
468 ksteube 1312 return temp;
469 ksteube 1345 #else
470     throw DataException("Error - loadMesh: is not compiled with NetCFD. Please contact your installation manager.");
471     #endif /* USE_NETCDF */
472 ksteube 1312 }
473    
474 jgs 82 AbstractContinuousDomain* readMesh(const std::string& fileName,
475 gross 1059 int integrationOrder,
476     int reducedIntegrationOrder,
477 ksteube 1312 int optimize)
478 jgs 82 {
479     //
480     // create a copy of the filename to overcome the non-constness of call
481     // to Finley_Mesh_read
482 bcumming 751 Finley_Mesh* fMesh=0;
483 woo409 757 // Win32 refactor
484 phornby 1628 if( fileName.size() == 0 )
485     {
486     throw DataException("Null file name!");
487     }
488    
489     char *fName = TMPMEMALLOC(fileName.size()+1,char);
490    
491 jgs 82 strcpy(fName,fileName.c_str());
492 ksteube 1345 double blocktimer_start = blocktimer_time();
493 bcumming 751
494 ksteube 1360 fMesh=Finley_Mesh_read_MPI(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
495     checkFinleyError();
496     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
497    
498     /* win32 refactor */
499     TMPMEMFREE(fName);
500    
501     blocktimer_increment("ReadMesh()", blocktimer_start);
502     return temp;
503     }
504    
505 gross 934 AbstractContinuousDomain* readGmsh(const std::string& fileName,
506     int numDim,
507     int integrationOrder,
508     int reducedIntegrationOrder,
509 ksteube 1312 int optimize)
510 gross 934 {
511     //
512     // create a copy of the filename to overcome the non-constness of call
513     // to Finley_Mesh_read
514     Finley_Mesh* fMesh=0;
515     // Win32 refactor
516 phornby 1628 if( fileName.size() == 0 )
517     {
518     throw DataException("Null file name!");
519     }
520    
521     char *fName = TMPMEMALLOC(fileName.size()+1,char);
522    
523 gross 934 strcpy(fName,fileName.c_str());
524 ksteube 1345 double blocktimer_start = blocktimer_time();
525 gross 934
526 ksteube 1312 fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
527 gross 934 checkFinleyError();
528     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
529    
530     /* win32 refactor */
531     TMPMEMFREE(fName);
532    
533 ksteube 1345 blocktimer_increment("ReadGmsh()", blocktimer_start);
534 gross 934 return temp;
535     }
536    
537 jgs 82 AbstractContinuousDomain* brick(int n0,int n1,int n2,int order,
538     double l0,double l1,double l2,
539     int periodic0,int periodic1,
540     int periodic2,
541     int integrationOrder,
542 gross 1059 int reducedIntegrationOrder,
543 ksteube 1312 int useElementsOnFace,
544     int useFullElementOrder,
545     int optimize)
546 jgs 82 {
547     int numElements[]={n0,n1,n2};
548     double length[]={l0,l1,l2};
549     int periodic[]={periodic0, periodic1, periodic2};
550    
551     //
552     // linearInterpolation
553 bcumming 751 Finley_Mesh* fMesh=NULL;
554    
555 jgs 82 if (order==1) {
556 gross 1062 fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
557 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
558 bcumming 751 }
559     else if (order==2) {
560 gross 1062 fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
561 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
562 jgs 82 } else {
563     stringstream temp;
564     temp << "Illegal interpolation order: " << order;
565     setFinleyError(VALUE_ERROR,temp.str().c_str());
566     }
567     //
568     // Convert any finley errors into a C++ exception
569     checkFinleyError();
570     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
571     return temp;
572     }
573     AbstractContinuousDomain* rectangle(int n0,int n1,int order,
574     double l0, double l1,
575     int periodic0,int periodic1,
576     int integrationOrder,
577 gross 1059 int reducedIntegrationOrder,
578 ksteube 1312 int useElementsOnFace,
579     int useFullElementOrder,
580     int optimize)
581 jgs 82 {
582     int numElements[]={n0,n1};
583     double length[]={l0,l1};
584     int periodic[]={periodic0, periodic1};
585    
586 bcumming 751 Finley_Mesh* fMesh=0;
587 jgs 82 if (order==1) {
588 gross 1062 fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
589 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
590 bcumming 751 }
591     else if (order==2) {
592 gross 1062 fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
593 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
594 bcumming 751 }
595     else {
596 jgs 82 stringstream temp;
597     temp << "Illegal interpolation order: " << order;
598     setFinleyError(VALUE_ERROR,temp.str().c_str());
599     }
600     //
601     // Convert any finley errors into a C++ exception
602     checkFinleyError();
603     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
604     return temp;
605     }
606 gross 1062
607 jgs 110 AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)
608 jgs 82 {
609 jgs 110 Finley_Mesh* fMesh=0;
610     //
611     // extract the meshes from meshList
612     int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
613 woo409 757 Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
614 jgs 110 for (int i=0;i<numMsh;++i) {
615     AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
616     const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
617     mshes[i]=finley_meshListMember->getFinley_Mesh();
618     }
619     //
620     // merge the meshes:
621     fMesh=Finley_Mesh_merge(numMsh,mshes);
622 bcumming 759 TMPMEMFREE(mshes);
623 jgs 110 //
624     // Convert any finley errors into a C++ exception
625     checkFinleyError();
626     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
627 woo409 757
628 jgs 82 return temp;
629     }
630     AbstractContinuousDomain* glueFaces(const boost::python::list& meshList,
631 ksteube 1312 double safety_factor,
632     double tolerance,
633     int optimize)
634 jgs 82 {
635 jgs 110 Finley_Mesh* fMesh=0;
636     //
637     // merge the meshes:
638     AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
639     //
640     // glue the faces:
641     const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
642     fMesh=merged_finley_meshes->getFinley_Mesh();
643 ksteube 1312 Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
644 bcumming 751
645 jgs 110 //
646     // Convert any finley errors into a C++ exception
647     checkFinleyError();
648     return merged_meshes;
649 jgs 82 }
650     AbstractContinuousDomain* joinFaces(const boost::python::list& meshList,
651     double safety_factor,
652 gross 1059 double tolerance,
653 ksteube 1312 int optimize)
654 jgs 82 {
655 jgs 110 Finley_Mesh* fMesh=0;
656     //
657     // merge the meshes:
658     AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
659     //
660     // join the faces:
661     const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
662     fMesh=merged_finley_meshes->getFinley_Mesh();
663 ksteube 1312 Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
664 jgs 110 //
665     // Convert any finley errors into a C++ exception
666     checkFinleyError();
667     return merged_meshes;
668 jgs 82 }
669    
670 gross 1062 // end of namespace
671    
672     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26