/[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 3317 - (hide annotations)
Thu Oct 28 00:50:41 2010 UTC (8 years, 10 months ago) by caltinay
File size: 31640 byte(s)
Removed bogus mpi.h includes, replaced others by our esysUtils wrapper
and rearranged so that the wrapper always comes before netcdf which fixes
linking problems when disabling mpi on systems with netCDF 4.x.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26