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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2078 - (hide annotations)
Thu Nov 20 16:10:10 2008 UTC (10 years, 10 months ago) by phornby
File size: 32449 byte(s)
Two changes.

1. Move blocktimer from escript to esysUtils.
2. Make it possible to link to paso as a DLL or .so.

Should have no effect on 'nix's

In respect of 1., blocktimer had begun to spring up everywhere, so
for the moment I thought it best to move it to the only other library that
pops up all over the place.

In respect of 2., paso needed to be a DLL in order to use the windows intelc /fast
option, which does aggressive multi-file optimisations. Even in its current form, it either
vectorises or parallelises  hundreds more loops in the esys system than appear in the pragmas.

In achieving 2. I have not been too delicate in adding

PASO_DLL_API

declarations to the .h files in paso/src. Only toward the end of the process of
the conversion, when the number of linker errors dropped below 20, say, did I choosy about what
functions in a header I declared PASO_DLL_API. As a result, there are likely to be many routines
declared as external functions symbols that are in fact internal to the paso DLL. 
Why is this an issue?  It prevents the intelc compiler from getting aggressive on the paso module.
With pain there is sometimes gain. At least all the DLL rules in windows give good
(non-microsoft) compiler writers a chance to really shine.

So, if you should see a PASO_DLL_API on a function in a paso header file,
and think to yourself, "that function is only called in paso, why export it?", then feel free to
delete the PASO_DLL_API export declaration.

Here's hoping for no breakage.....
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 phornby 2078 #include "esysUtils/blocktimer.h"
25 ksteube 1345 }
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 jfenwick 1872 // AbstractContinuousDomain* loadMesh(const std::string& fileName)
52     Domain_ptr loadMesh(const std::string& fileName)
53 ksteube 1312 {
54 ksteube 1345 #ifdef USE_NETCDF
55     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 ksteube 1990
451     if (Finley_noError()) Finley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent);
452     TMPMEMFREE(first_DofComponent);
453     TMPMEMFREE(first_NodeComponent);
454 ksteube 1347
455 ksteube 1345 } /* Finley_noError() after Finley_Mesh_alloc() */
456    
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 jfenwick 1872 return temp->getPtr();
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 jfenwick 1872 Domain_ptr 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 jfenwick 1872 return temp->getPtr();
503 ksteube 1360 }
504    
505 jfenwick 1872 Domain_ptr readGmsh(const std::string& fileName,
506 gross 934 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 jfenwick 1872 return temp->getPtr();
535 gross 934 }
536    
537 jfenwick 1872 /* AbstractContinuousDomain* brick(int n0,int n1,int n2,int order,*/
538     Domain_ptr brick(int n0,int n1,int n2,int order,
539 jgs 82 double l0,double l1,double l2,
540     int periodic0,int periodic1,
541     int periodic2,
542     int integrationOrder,
543 gross 1059 int reducedIntegrationOrder,
544 ksteube 1312 int useElementsOnFace,
545     int useFullElementOrder,
546     int optimize)
547 jgs 82 {
548     int numElements[]={n0,n1,n2};
549     double length[]={l0,l1,l2};
550     int periodic[]={periodic0, periodic1, periodic2};
551    
552     //
553     // linearInterpolation
554 bcumming 751 Finley_Mesh* fMesh=NULL;
555    
556 jgs 82 if (order==1) {
557 gross 1062 fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
558 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
559 bcumming 751 }
560     else if (order==2) {
561 gross 1062 fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
562 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
563 jgs 82 } else {
564     stringstream temp;
565     temp << "Illegal interpolation order: " << order;
566     setFinleyError(VALUE_ERROR,temp.str().c_str());
567     }
568     //
569     // Convert any finley errors into a C++ exception
570     checkFinleyError();
571     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
572 jfenwick 1872 return temp->getPtr();
573 jgs 82 }
574 jfenwick 1872
575     /* AbstractContinuousDomain* rectangle(int n0,int n1,int order,*/
576     Domain_ptr rectangle(int n0,int n1,int order,
577 jgs 82 double l0, double l1,
578     int periodic0,int periodic1,
579     int integrationOrder,
580 gross 1059 int reducedIntegrationOrder,
581 ksteube 1312 int useElementsOnFace,
582     int useFullElementOrder,
583     int optimize)
584 jgs 82 {
585     int numElements[]={n0,n1};
586     double length[]={l0,l1};
587     int periodic[]={periodic0, periodic1};
588    
589 bcumming 751 Finley_Mesh* fMesh=0;
590 jgs 82 if (order==1) {
591 gross 1062 fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
592 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
593 bcumming 751 }
594     else if (order==2) {
595 gross 1062 fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
596 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
597 bcumming 751 }
598     else {
599 jgs 82 stringstream temp;
600     temp << "Illegal interpolation order: " << order;
601     setFinleyError(VALUE_ERROR,temp.str().c_str());
602     }
603     //
604     // Convert any finley errors into a C++ exception
605     checkFinleyError();
606     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
607 jfenwick 1872 return temp->getPtr();
608 jgs 82 }
609 gross 1062
610 jfenwick 1872 Domain_ptr meshMerge(const boost::python::list& meshList)
611 jgs 82 {
612 jgs 110 Finley_Mesh* fMesh=0;
613     //
614     // extract the meshes from meshList
615     int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
616 woo409 757 Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
617 jgs 110 for (int i=0;i<numMsh;++i) {
618     AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
619     const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
620     mshes[i]=finley_meshListMember->getFinley_Mesh();
621     }
622     //
623     // merge the meshes:
624     fMesh=Finley_Mesh_merge(numMsh,mshes);
625 bcumming 759 TMPMEMFREE(mshes);
626 jgs 110 //
627     // Convert any finley errors into a C++ exception
628     checkFinleyError();
629     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
630 woo409 757
631 jfenwick 1872 return temp->getPtr();
632 jgs 82 }
633 jfenwick 1872
634     Domain_ptr glueFaces(const boost::python::list& meshList,
635 ksteube 1312 double safety_factor,
636     double tolerance,
637     int optimize)
638 jgs 82 {
639 jgs 110 Finley_Mesh* fMesh=0;
640     //
641     // merge the meshes:
642 jfenwick 1872 Domain_ptr merged_meshes=meshMerge(meshList);
643    
644 jgs 110 //
645     // glue the faces:
646 jfenwick 1872 const MeshAdapter* merged_finley_meshes=dynamic_cast<const MeshAdapter*>(merged_meshes.get());
647 jgs 110 fMesh=merged_finley_meshes->getFinley_Mesh();
648 ksteube 1312 Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
649 bcumming 751
650 jgs 110 //
651     // Convert any finley errors into a C++ exception
652     checkFinleyError();
653 jfenwick 1872 return merged_meshes->getPtr();
654 jgs 82 }
655 jfenwick 1872 Domain_ptr joinFaces(const boost::python::list& meshList,
656 jgs 82 double safety_factor,
657 gross 1059 double tolerance,
658 ksteube 1312 int optimize)
659 jgs 82 {
660 jgs 110 Finley_Mesh* fMesh=0;
661     //
662     // merge the meshes:
663 jfenwick 1872 Domain_ptr merged_meshes=meshMerge(meshList);
664 jgs 110 //
665     // join the faces:
666 jfenwick 1872 const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes.get());
667 jgs 110 fMesh=merged_finley_meshes->getFinley_Mesh();
668 ksteube 1312 Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
669 jgs 110 //
670     // Convert any finley errors into a C++ exception
671     checkFinleyError();
672 jfenwick 1872 return merged_meshes->getPtr();
673 jgs 82 }
674    
675 gross 1062 // end of namespace
676    
677     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26