/[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 2748 - (hide annotations)
Tue Nov 17 07:32:59 2009 UTC (9 years, 10 months ago) by gross
File size: 33657 byte(s)
Macro elements are implemented now. VTK writer for macro elements still needs testing.
1 ksteube 1312
2     /*******************************************************
3 ksteube 1811 *
4 jfenwick 2548 * Copyright (c) 2003-2009 by University of Queensland
5 ksteube 1811 * 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 gross 2748 Finley_ReferenceElementSet *refElements= Finley_ReferenceElementSet_alloc((ElementTypeId)Elements_TypeId,mesh_p->order, mesh_p->reduced_order);
196     if (Finley_noError()) {
197     mesh_p->Elements=Finley_ElementFile_alloc(refElements, mpi_info);
198     }
199     Finley_ReferenceElementSet_dealloc(refElements);
200 ksteube 1347 if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Elements, num_Elements);
201 gross 2748 if (Finley_noError()) {
202     mesh_p->Elements->minColor=0;
203     mesh_p->Elements->maxColor=num_Elements-1;
204     if (num_Elements>0) {
205 ksteube 1346 // Elements_Id
206     if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
207     throw DataException("Error - loadMesh:: unable to read Elements_Id from netCDF file: " + *fName);
208     if (! nc_var_temp->get(&mesh_p->Elements->Id[0], num_Elements) ) {
209 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Id);
210 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Elements_Id from NetCDF file: " + *fName);
211     }
212     // Elements_Tag
213     if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
214     throw DataException("Error - loadMesh:: unable to read Elements_Tag from netCDF file: " + *fName);
215     if (! nc_var_temp->get(&mesh_p->Elements->Tag[0], num_Elements) ) {
216 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Tag);
217 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Elements_Tag from NetCDF file: " + *fName);
218     }
219     // Elements_Owner
220     if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
221     throw DataException("Error - loadMesh:: unable to read Elements_Owner from netCDF file: " + *fName);
222     if (! nc_var_temp->get(&mesh_p->Elements->Owner[0], num_Elements) ) {
223 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Owner);
224 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Elements_Owner from NetCDF file: " + *fName);
225     }
226     // Elements_Color
227     if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
228     throw DataException("Error - loadMesh:: unable to read Elements_Color from netCDF file: " + *fName);
229     if (! nc_var_temp->get(&mesh_p->Elements->Color[0], num_Elements) ) {
230 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Color);
231 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Elements_Color from NetCDF file: " + *fName);
232     }
233     // Elements_Nodes
234 gross 2748 int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);
235 ksteube 1346 if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
236 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Nodes);
237 ksteube 1346 throw DataException("Error - loadMesh:: unable to read Elements_Nodes from netCDF file: " + *fName);
238     }
239     if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {
240 ksteube 1755 TMPMEMFREE(Elements_Nodes);
241 ksteube 1346 throw DataException("Error - load:: unable to recover Elements_Nodes from netCDF file: " + *fName);
242     }
243 gross 2748 // Copy temp array into mesh_p->Elements->Nodes
244     for (int i=0; i<num_Elements; i++) {
245     for (int j=0; j<num_Elements_numNodes; j++) {
246     mesh_p->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
247     = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
248     }
249     }
250     TMPMEMFREE(Elements_Nodes);
251    
252     } /* num_Elements>0 */
253     }
254 ksteube 1345 }
255    
256     /* get the face elements */
257 ksteube 1346 if (Finley_noError()) {
258 gross 2748 Finley_ReferenceElementSet *refFaceElements= Finley_ReferenceElementSet_alloc((ElementTypeId)FaceElements_TypeId ,mesh_p->order, mesh_p->reduced_order);
259     if (Finley_noError()) {
260     mesh_p->FaceElements=Finley_ElementFile_alloc(refFaceElements, mpi_info);
261     }
262     Finley_ReferenceElementSet_dealloc(refFaceElements);
263 ksteube 1347 if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->FaceElements, num_FaceElements);
264 gross 2748 if (Finley_noError()) {
265     mesh_p->FaceElements->minColor=0;
266     mesh_p->FaceElements->maxColor=num_FaceElements-1;
267     if (num_FaceElements>0) {
268 ksteube 1346 // FaceElements_Id
269     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
270     throw DataException("Error - loadMesh:: unable to read FaceElements_Id from netCDF file: " + *fName);
271     if (! nc_var_temp->get(&mesh_p->FaceElements->Id[0], num_FaceElements) ) {
272 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Id);
273 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover FaceElements_Id from NetCDF file: " + *fName);
274     }
275     // FaceElements_Tag
276     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
277     throw DataException("Error - loadMesh:: unable to read FaceElements_Tag from netCDF file: " + *fName);
278     if (! nc_var_temp->get(&mesh_p->FaceElements->Tag[0], num_FaceElements) ) {
279 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Tag);
280 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover FaceElements_Tag from NetCDF file: " + *fName);
281     }
282     // FaceElements_Owner
283     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
284     throw DataException("Error - loadMesh:: unable to read FaceElements_Owner from netCDF file: " + *fName);
285     if (! nc_var_temp->get(&mesh_p->FaceElements->Owner[0], num_FaceElements) ) {
286 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Owner);
287 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover FaceElements_Owner from NetCDF file: " + *fName);
288     }
289     // FaceElements_Color
290     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
291     throw DataException("Error - loadMesh:: unable to read FaceElements_Color from netCDF file: " + *fName);
292     if (! nc_var_temp->get(&mesh_p->FaceElements->Color[0], num_FaceElements) ) {
293 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Color);
294 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover FaceElements_Color from NetCDF file: " + *fName);
295     }
296     // FaceElements_Nodes
297 gross 2748 int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);
298 ksteube 1346 if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
299 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Nodes);
300 ksteube 1346 throw DataException("Error - loadMesh:: unable to read FaceElements_Nodes from netCDF file: " + *fName);
301     }
302     if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
303 ksteube 1755 TMPMEMFREE(FaceElements_Nodes);
304 ksteube 1346 throw DataException("Error - load:: unable to recover FaceElements_Nodes from netCDF file: " + *fName);
305     }
306 gross 2748 // Copy temp array into mesh_p->FaceElements->Nodes
307     for (int i=0; i<num_FaceElements; i++) {
308     for (int j=0; j<num_FaceElements_numNodes; j++) {
309     mesh_p->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
310     }
311     }
312     TMPMEMFREE(FaceElements_Nodes);
313     } /* num_FaceElements>0 */
314     }
315 ksteube 1346 }
316 ksteube 1345
317 ksteube 1346 /* get the Contact elements */
318     if (Finley_noError()) {
319 gross 2748 Finley_ReferenceElementSet *refContactElements= Finley_ReferenceElementSet_alloc((ElementTypeId)ContactElements_TypeId,mesh_p->order, mesh_p->reduced_order);
320     if (Finley_noError()) {
321     mesh_p->ContactElements=Finley_ElementFile_alloc(refContactElements, mpi_info);
322     }
323     Finley_ReferenceElementSet_dealloc(refContactElements);
324 ksteube 1347 if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->ContactElements, num_ContactElements);
325 gross 2748 if (Finley_noError()) {
326     mesh_p->ContactElements->minColor=0;
327     mesh_p->ContactElements->maxColor=num_ContactElements-1;
328     if (num_ContactElements>0) {
329 ksteube 1346 // ContactElements_Id
330     if (! ( nc_var_temp = dataFile.get_var("ContactElements_Id")) )
331     throw DataException("Error - loadMesh:: unable to read ContactElements_Id from netCDF file: " + *fName);
332     if (! nc_var_temp->get(&mesh_p->ContactElements->Id[0], num_ContactElements) ) {
333 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Id);
334 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover ContactElements_Id from NetCDF file: " + *fName);
335     }
336     // ContactElements_Tag
337     if (! ( nc_var_temp = dataFile.get_var("ContactElements_Tag")) )
338     throw DataException("Error - loadMesh:: unable to read ContactElements_Tag from netCDF file: " + *fName);
339     if (! nc_var_temp->get(&mesh_p->ContactElements->Tag[0], num_ContactElements) ) {
340 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Tag);
341 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover ContactElements_Tag from NetCDF file: " + *fName);
342     }
343     // ContactElements_Owner
344     if (! ( nc_var_temp = dataFile.get_var("ContactElements_Owner")) )
345     throw DataException("Error - loadMesh:: unable to read ContactElements_Owner from netCDF file: " + *fName);
346     if (! nc_var_temp->get(&mesh_p->ContactElements->Owner[0], num_ContactElements) ) {
347 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Owner);
348 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover ContactElements_Owner from NetCDF file: " + *fName);
349     }
350     // ContactElements_Color
351     if (! ( nc_var_temp = dataFile.get_var("ContactElements_Color")) )
352     throw DataException("Error - loadMesh:: unable to read ContactElements_Color from netCDF file: " + *fName);
353     if (! nc_var_temp->get(&mesh_p->ContactElements->Color[0], num_ContactElements) ) {
354 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Color);
355 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover ContactElements_Color from NetCDF file: " + *fName);
356     }
357     // ContactElements_Nodes
358 gross 2748 int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);
359 ksteube 1346 if (!(nc_var_temp = dataFile.get_var("ContactElements_Nodes"))) {
360 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Nodes);
361 ksteube 1346 throw DataException("Error - loadMesh:: unable to read ContactElements_Nodes from netCDF file: " + *fName);
362     }
363     if (! nc_var_temp->get(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes) ) {
364 ksteube 1755 TMPMEMFREE(ContactElements_Nodes);
365 ksteube 1346 throw DataException("Error - load:: unable to recover ContactElements_Nodes from netCDF file: " + *fName);
366     }
367 gross 2748 // Copy temp array into mesh_p->ContactElements->Nodes
368     for (int i=0; i<num_ContactElements; i++) {
369     for (int j=0; j<num_ContactElements_numNodes; j++) {
370     mesh_p->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]= ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)];
371     }
372     }
373     TMPMEMFREE(ContactElements_Nodes);
374     } /* num_ContactElements>0 */
375     }
376    
377 ksteube 1346 }
378 ksteube 1345
379 ksteube 1346 /* get the Points (nodal elements) */
380     if (Finley_noError()) {
381 gross 2748 Finley_ReferenceElementSet *refPoints= Finley_ReferenceElementSet_alloc((ElementTypeId)Points_TypeId,mesh_p->order, mesh_p->reduced_order);
382     if (Finley_noError()) {
383     mesh_p->Points=Finley_ElementFile_alloc(refPoints, mpi_info);
384     }
385     Finley_ReferenceElementSet_dealloc(refPoints);
386 ksteube 1347 if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Points, num_Points);
387 gross 2748 if (Finley_noError()) {
388     mesh_p->Points->minColor=0;
389     mesh_p->Points->maxColor=num_Points-1;
390     if (num_Points>0) {
391 ksteube 1346 // Points_Id
392     if (! ( nc_var_temp = dataFile.get_var("Points_Id")) )
393     throw DataException("Error - loadMesh:: unable to read Points_Id from netCDF file: " + *fName);
394     if (! nc_var_temp->get(&mesh_p->Points->Id[0], num_Points) ) {
395 ksteube 1755 TMPMEMFREE(mesh_p->Points->Id);
396 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Points_Id from NetCDF file: " + *fName);
397     }
398     // Points_Tag
399     if (! ( nc_var_temp = dataFile.get_var("Points_Tag")) )
400     throw DataException("Error - loadMesh:: unable to read Points_Tag from netCDF file: " + *fName);
401     if (! nc_var_temp->get(&mesh_p->Points->Tag[0], num_Points) ) {
402 ksteube 1755 TMPMEMFREE(mesh_p->Points->Tag);
403 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Points_Tag from NetCDF file: " + *fName);
404     }
405     // Points_Owner
406     if (! ( nc_var_temp = dataFile.get_var("Points_Owner")) )
407     throw DataException("Error - loadMesh:: unable to read Points_Owner from netCDF file: " + *fName);
408     if (! nc_var_temp->get(&mesh_p->Points->Owner[0], num_Points) ) {
409 ksteube 1755 TMPMEMFREE(mesh_p->Points->Owner);
410 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Points_Owner from NetCDF file: " + *fName);
411     }
412     // Points_Color
413     if (! ( nc_var_temp = dataFile.get_var("Points_Color")) )
414     throw DataException("Error - loadMesh:: unable to read Points_Color from netCDF file: " + *fName);
415     if (! nc_var_temp->get(&mesh_p->Points->Color[0], num_Points) ) {
416 ksteube 1755 TMPMEMFREE(mesh_p->Points->Color);
417 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Points_Color from NetCDF file: " + *fName);
418     }
419     // Points_Nodes
420     int *Points_Nodes = TMPMEMALLOC(num_Points,int);
421     if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
422 ksteube 1755 TMPMEMFREE(mesh_p->Points->Nodes);
423 ksteube 1346 throw DataException("Error - loadMesh:: unable to read Points_Nodes from netCDF file: " + *fName);
424     }
425     if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {
426 ksteube 1755 TMPMEMFREE(Points_Nodes);
427 ksteube 1346 throw DataException("Error - load:: unable to recover Points_Nodes from netCDF file: " + *fName);
428     }
429     // Copy temp array into mesh_p->Points->Nodes
430     for (int i=0; i<num_Points; i++) {
431     mesh_p->Nodes->Id[mesh_p->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
432     }
433     TMPMEMFREE(Points_Nodes);
434 gross 2748
435     } /* num_Points>0 */
436     }
437     }
438 ksteube 1345
439 ksteube 1347 /* get the tags */
440     if (Finley_noError()) {
441     if (num_Tags>0) {
442     // Temp storage to gather node IDs
443     int *Tags_keys = TMPMEMALLOC(num_Tags, int);
444     char name_temp[4096];
445     int i;
446 ksteube 1345
447 ksteube 1347 // Tags_keys
448     if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) )
449     throw DataException("Error - loadMesh:: unable to read Tags_keys from netCDF file: " + *fName);
450     if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {
451 ksteube 1755 TMPMEMFREE(Tags_keys);
452 ksteube 1347 throw DataException("Error - loadMesh:: unable to recover Tags_keys from NetCDF file: " + *fName);
453     }
454     for (i=0; i<num_Tags; i++) {
455     // Retrieve tag name
456     sprintf(name_temp, "Tags_name_%d", i);
457     if (! (attr=dataFile.get_att(name_temp)) ) {
458     sprintf(error_msg,"Error retrieving tag name from NetCDF file '%s'", fName);
459     throw DataException(error_msg);
460     }
461     char *name = attr->as_string(0);
462     delete attr;
463     Finley_Mesh_addTagMap(mesh_p, name, Tags_keys[i]);
464     }
465     }
466     }
467 ksteube 1990
468     if (Finley_noError()) Finley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent);
469     TMPMEMFREE(first_DofComponent);
470     TMPMEMFREE(first_NodeComponent);
471 ksteube 1347
472 ksteube 1345 } /* Finley_noError() after Finley_Mesh_alloc() */
473    
474 ksteube 1312 checkFinleyError();
475 ksteube 1345 temp=new MeshAdapter(mesh_p);
476    
477     if (! Finley_noError()) {
478     Finley_Mesh_free(mesh_p);
479     }
480    
481 ksteube 1312 /* win32 refactor */
482     TMPMEMFREE(fName);
483 ksteube 1345
484     blocktimer_increment("LoadMesh()", blocktimer_start);
485 jfenwick 1872 return temp->getPtr();
486 ksteube 1345 #else
487     throw DataException("Error - loadMesh: is not compiled with NetCFD. Please contact your installation manager.");
488     #endif /* USE_NETCDF */
489 ksteube 1312 }
490    
491 jfenwick 1872 Domain_ptr readMesh(const std::string& fileName,
492 gross 1059 int integrationOrder,
493     int reducedIntegrationOrder,
494 ksteube 1312 int optimize)
495 jgs 82 {
496     //
497     // create a copy of the filename to overcome the non-constness of call
498     // to Finley_Mesh_read
499 bcumming 751 Finley_Mesh* fMesh=0;
500 woo409 757 // Win32 refactor
501 phornby 1628 if( fileName.size() == 0 )
502     {
503     throw DataException("Null file name!");
504     }
505    
506     char *fName = TMPMEMALLOC(fileName.size()+1,char);
507    
508 jgs 82 strcpy(fName,fileName.c_str());
509 ksteube 1345 double blocktimer_start = blocktimer_time();
510 bcumming 751
511 gross 2748 fMesh=Finley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
512 ksteube 1360 checkFinleyError();
513     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
514    
515     /* win32 refactor */
516     TMPMEMFREE(fName);
517    
518     blocktimer_increment("ReadMesh()", blocktimer_start);
519 jfenwick 1872 return temp->getPtr();
520 ksteube 1360 }
521    
522 jfenwick 1872 Domain_ptr readGmsh(const std::string& fileName,
523 gross 934 int numDim,
524     int integrationOrder,
525     int reducedIntegrationOrder,
526 gross 2722 int optimize,
527     int useMacroElements)
528 gross 934 {
529     //
530     // create a copy of the filename to overcome the non-constness of call
531     // to Finley_Mesh_read
532     Finley_Mesh* fMesh=0;
533     // Win32 refactor
534 phornby 1628 if( fileName.size() == 0 )
535     {
536     throw DataException("Null file name!");
537     }
538    
539     char *fName = TMPMEMALLOC(fileName.size()+1,char);
540    
541 gross 934 strcpy(fName,fileName.c_str());
542 ksteube 1345 double blocktimer_start = blocktimer_time();
543 gross 934
544 gross 2722 fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE), (useMacroElements ? TRUE : FALSE));
545 gross 934 checkFinleyError();
546     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
547    
548     /* win32 refactor */
549     TMPMEMFREE(fName);
550    
551 ksteube 1345 blocktimer_increment("ReadGmsh()", blocktimer_start);
552 jfenwick 1872 return temp->getPtr();
553 gross 934 }
554    
555 jfenwick 1872 /* AbstractContinuousDomain* brick(int n0,int n1,int n2,int order,*/
556     Domain_ptr brick(int n0,int n1,int n2,int order,
557 jgs 82 double l0,double l1,double l2,
558     int periodic0,int periodic1,
559     int periodic2,
560     int integrationOrder,
561 gross 1059 int reducedIntegrationOrder,
562 ksteube 1312 int useElementsOnFace,
563     int useFullElementOrder,
564     int optimize)
565 jgs 82 {
566     int numElements[]={n0,n1,n2};
567     double length[]={l0,l1,l2};
568     int periodic[]={periodic0, periodic1, periodic2};
569    
570     //
571     // linearInterpolation
572 bcumming 751 Finley_Mesh* fMesh=NULL;
573    
574 jgs 82 if (order==1) {
575 gross 2722 fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
576 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
577 gross 2722 } else if (order==2) {
578     fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
579     useElementsOnFace,useFullElementOrder,FALSE, (optimize ? TRUE : FALSE)) ;
580     } else if (order==-1) {
581     fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
582     useElementsOnFace,useFullElementOrder,TRUE,(optimize ? TRUE : FALSE)) ;
583 jgs 82 } else {
584     stringstream temp;
585     temp << "Illegal interpolation order: " << order;
586     setFinleyError(VALUE_ERROR,temp.str().c_str());
587     }
588     //
589     // Convert any finley errors into a C++ exception
590     checkFinleyError();
591     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
592 jfenwick 1872 return temp->getPtr();
593 jgs 82 }
594 jfenwick 1872
595     /* AbstractContinuousDomain* rectangle(int n0,int n1,int order,*/
596     Domain_ptr rectangle(int n0,int n1,int order,
597 jgs 82 double l0, double l1,
598     int periodic0,int periodic1,
599     int integrationOrder,
600 gross 1059 int reducedIntegrationOrder,
601 ksteube 1312 int useElementsOnFace,
602     int useFullElementOrder,
603     int optimize)
604 jgs 82 {
605     int numElements[]={n0,n1};
606     double length[]={l0,l1};
607     int periodic[]={periodic0, periodic1};
608    
609 bcumming 751 Finley_Mesh* fMesh=0;
610 jgs 82 if (order==1) {
611 gross 2722 fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
612 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
613 gross 2722 } else if (order==2) {
614     fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
615     useElementsOnFace,useFullElementOrder,FALSE,(optimize ? TRUE : FALSE));
616     } else if (order==-1) {
617     fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
618     useElementsOnFace,useFullElementOrder,TRUE,(optimize ? TRUE : FALSE));
619     } else {
620 jgs 82 stringstream temp;
621     temp << "Illegal interpolation order: " << order;
622     setFinleyError(VALUE_ERROR,temp.str().c_str());
623     }
624     //
625     // Convert any finley errors into a C++ exception
626     checkFinleyError();
627     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
628 jfenwick 1872 return temp->getPtr();
629 jgs 82 }
630 gross 1062
631 jfenwick 1872 Domain_ptr meshMerge(const boost::python::list& meshList)
632 jgs 82 {
633 jgs 110 Finley_Mesh* fMesh=0;
634     //
635     // extract the meshes from meshList
636     int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
637 woo409 757 Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
638 jgs 110 for (int i=0;i<numMsh;++i) {
639     AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
640     const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
641     mshes[i]=finley_meshListMember->getFinley_Mesh();
642     }
643     //
644     // merge the meshes:
645     fMesh=Finley_Mesh_merge(numMsh,mshes);
646 bcumming 759 TMPMEMFREE(mshes);
647 jgs 110 //
648     // Convert any finley errors into a C++ exception
649     checkFinleyError();
650     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
651 woo409 757
652 jfenwick 1872 return temp->getPtr();
653 jgs 82 }
654 jfenwick 1872
655     Domain_ptr glueFaces(const boost::python::list& meshList,
656 ksteube 1312 double safety_factor,
657     double tolerance,
658     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    
665 jgs 110 //
666     // glue the faces:
667 jfenwick 1872 const MeshAdapter* merged_finley_meshes=dynamic_cast<const MeshAdapter*>(merged_meshes.get());
668 jgs 110 fMesh=merged_finley_meshes->getFinley_Mesh();
669 ksteube 1312 Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
670 bcumming 751
671 jgs 110 //
672     // Convert any finley errors into a C++ exception
673     checkFinleyError();
674 jfenwick 1872 return merged_meshes->getPtr();
675 jgs 82 }
676 jfenwick 1872 Domain_ptr joinFaces(const boost::python::list& meshList,
677 jgs 82 double safety_factor,
678 gross 1059 double tolerance,
679 ksteube 1312 int optimize)
680 jgs 82 {
681 jgs 110 Finley_Mesh* fMesh=0;
682     //
683     // merge the meshes:
684 jfenwick 1872 Domain_ptr merged_meshes=meshMerge(meshList);
685 jgs 110 //
686     // join the faces:
687 jfenwick 1872 const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes.get());
688 jgs 110 fMesh=merged_finley_meshes->getFinley_Mesh();
689 ksteube 1312 Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
690 jgs 110 //
691     // Convert any finley errors into a C++ exception
692     checkFinleyError();
693 jfenwick 1872 return merged_meshes->getPtr();
694 jgs 82 }
695    
696 gross 1062 // end of namespace
697    
698     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26