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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1345 - (show annotations)
Wed Nov 14 07:53:34 2007 UTC (12 years, 4 months ago) by ksteube
Original Path: trunk/finley/src/CPPAdapter/MeshAdapterFactory.cpp
File size: 14836 byte(s)
Created LoadMesh to read a mesh from a distributed NetCDF file.
Can read nodes but not elements yet.

1
2 /* $Id$ */
3
4 /*******************************************************
5 *
6 * Copyright 2003-2007 by ACceSS MNRF
7 * Copyright 2007 by University of Queensland
8 *
9 * http://esscc.uq.edu.au
10 * Primary Business: Queensland, Australia
11 * Licensed under the Open Software License version 3.0
12 * http://www.opensource.org/licenses/osl-3.0.php
13 *
14 *******************************************************/
15
16 #ifdef PASO_MPI
17 #include <mpi.h>
18 #endif
19 #ifdef USE_NETCDF
20 #include <netcdfcpp.h>
21 #endif
22 #include "MeshAdapterFactory.h"
23 #include "FinleyError.h"
24 extern "C" {
25 #include "escript/blocktimer.h"
26 }
27
28 #include <boost/python/extract.hpp>
29
30 #include <sstream>
31
32 using namespace std;
33 using namespace escript;
34
35 namespace finley {
36
37 #ifdef USE_NETCDF
38 // A convenience method to retrieve an integer attribute from a NetCDF file
39 int NetCDF_Get_Int_Attribute(NcFile *dataFile, char *fName, char *attr_name) {
40 NcAtt *attr;
41 char error_msg[LenErrorMsg_MAX];
42 if (! (attr=dataFile->get_att(attr_name)) ) {
43 sprintf(error_msg,"Error retrieving integer attribute '%s' from NetCDF file '%s'", attr_name, fName);
44 throw DataException(error_msg);
45 }
46 int temp = attr->as_int(0);
47 delete attr;
48 return(temp);
49 }
50 #endif
51
52 AbstractContinuousDomain* loadMesh(const std::string& fileName)
53 {
54 #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 // create a copy of the filename to overcome the non-constness of call
60 // to Finley_Mesh_read
61 // Win32 refactor
62 char *fName = ((fileName.size()+1)>0) ? TMPMEMALLOC((fileName.size()+1),char) : (char*)NULL;
63 strcpy(fName,fileName.c_str());
64
65 printf("ksteube finley::loadMesh %s\n", fName);
66
67 double blocktimer_start = blocktimer_time();
68 Finley_resetError();
69
70 // Open NetCDF file for reading
71 NcAtt *attr;
72 NcVar *nc_var_temp;
73 // netCDF error handler
74 NcError err(NcError::silent_nonfatal);
75 // Create the NetCDF file.
76 NcFile dataFile(fName, NcFile::ReadOnly);
77 if (!dataFile.is_valid()) {
78 sprintf(error_msg,"loadMesh: Opening file NetCDF %s for reading failed.", fName);
79 Finley_setError(IO_ERROR,error_msg);
80 Paso_MPIInfo_free( mpi_info );
81 throw DataException("Error - loadMesh:: Could not read NetCDF file.");
82 }
83
84 // Read NetCDF integer attributes
85 int mpi_size = NetCDF_Get_Int_Attribute(&dataFile, fName, "mpi_size");
86 int mpi_rank = NetCDF_Get_Int_Attribute(&dataFile, fName, "mpi_rank");
87 int numDim = NetCDF_Get_Int_Attribute(&dataFile, fName, "numDim");
88 int order = NetCDF_Get_Int_Attribute(&dataFile, fName, "order");
89 int reduced_order = NetCDF_Get_Int_Attribute(&dataFile, fName, "reduced_order");
90 int numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, "numNodes");
91 int num_Elements = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_Elements");
92 int num_FaceElements = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_FaceElements");
93 int num_ContactElements = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_ContactElements");
94 int num_Points = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_Points");
95
96 // Read mesh name
97 if (! (attr=dataFile.get_att("Name")) ) {
98 sprintf(error_msg,"Error retrieving mesh name from NetCDF file '%s'", fName);
99 throw DataException(error_msg);
100 }
101 char *name = attr->as_string(0);
102 delete attr;
103
104 /* allocate mesh */
105 mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
106 if (Finley_noError()) {
107
108 /* read nodes */
109 Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
110 // Nodes_Id
111 if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
112 throw DataException("Error - loadMesh:: unable to read Nodes_Id from netCDF file: " + *fName);
113 if (! nc_var_temp->get(&mesh_p->Nodes->Id[0], numNodes) ) {
114 free(&mesh_p->Nodes->Id);
115 throw DataException("Error - loadMesh:: unable to recover Nodes_Id from NetCDF file: " + *fName);
116 }
117 // printf("ksteube Nodes_Id: "); for (int i=0; i<numNodes; i++) { printf(" %d", mesh_p->Nodes->Id[i]); } printf("\n");
118 // Nodes_Tag
119 if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )
120 throw DataException("Error - loadMesh:: unable to read Nodes_Tag from netCDF file: " + *fName);
121 if (! nc_var_temp->get(&mesh_p->Nodes->Tag[0], numNodes) ) {
122 free(&mesh_p->Nodes->Tag);
123 throw DataException("Error - loadMesh:: unable to recover Nodes_Tag from NetCDF file: " + *fName);
124 }
125 // Nodes_gDOF
126 if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )
127 throw DataException("Error - loadMesh:: unable to read Nodes_gDOF from netCDF file: " + *fName);
128 if (! nc_var_temp->get(&mesh_p->Nodes->globalDegreesOfFreedom[0], numNodes) ) {
129 free(&mesh_p->Nodes->globalDegreesOfFreedom);
130 throw DataException("Error - loadMesh:: unable to recover Nodes_gDOF from NetCDF file: " + *fName);
131 }
132 // Nodes_gNI
133 if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )
134 throw DataException("Error - loadMesh:: unable to read Nodes_gNI from netCDF file: " + *fName);
135 if (! nc_var_temp->get(&mesh_p->Nodes->globalNodesIndex[0], numNodes) ) {
136 free(&mesh_p->Nodes->globalNodesIndex);
137 throw DataException("Error - loadMesh:: unable to recover Nodes_gNI from NetCDF file: " + *fName);
138 }
139 // Nodes_grDfI
140 if (! ( nc_var_temp = dataFile.get_var("Nodes_grDfI")) )
141 throw DataException("Error - loadMesh:: unable to read Nodes_grDfI from netCDF file: " + *fName);
142 if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedDOFIndex[0], numNodes) ) {
143 free(&mesh_p->Nodes->globalReducedDOFIndex);
144 throw DataException("Error - loadMesh:: unable to recover Nodes_grDfI from NetCDF file: " + *fName);
145 }
146 // Nodes_grNI
147 if (! ( nc_var_temp = dataFile.get_var("Nodes_grNI")) )
148 throw DataException("Error - loadMesh:: unable to read Nodes_grNI from netCDF file: " + *fName);
149 if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedNodesIndex[0], numNodes) ) {
150 free(&mesh_p->Nodes->globalReducedNodesIndex);
151 throw DataException("Error - loadMesh:: unable to recover Nodes_grNI from NetCDF file: " + *fName);
152 }
153 // Nodes_Coordinates
154 if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates"))) {
155 free(&mesh_p->Nodes->Coordinates);
156 throw DataException("Error - loadMesh:: unable to read Nodes_Coordinates from netCDF file: " + *fName);
157 }
158 if (! nc_var_temp->get(&(mesh_p->Nodes->Coordinates[0]), numNodes, numDim) ) {
159 free(&mesh_p->Nodes->Coordinates);
160 throw DataException("Error - load:: unable to recover Nodes_Coordinates from netCDF file: " + *fName);
161 }
162
163 #if 0 /* Not yet...finish the above first */
164 /* read elements */
165 if (Finley_noError()) {
166 mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
167 if (Finley_noError()) {
168 Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
169 mesh_p->Elements->minColor=0;
170 mesh_p->Elements->maxColor=numEle-1;
171 if (Finley_noError()) {
172 }
173 }
174 }
175 #endif
176
177 /* get the face elements */
178
179 /* get the Contact face element */
180
181 /* get the nodal element */
182
183 /* get the name tags */
184
185 } /* Finley_noError() after Finley_Mesh_alloc() */
186
187 #if 0 /* Not yet...finish the above first */
188 if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
189 if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
190 #endif
191
192 checkFinleyError();
193 temp=new MeshAdapter(mesh_p);
194
195 if (! Finley_noError()) {
196 Finley_Mesh_free(mesh_p);
197 }
198
199 /* win32 refactor */
200 TMPMEMFREE(fName);
201
202 blocktimer_increment("LoadMesh()", blocktimer_start);
203 return temp;
204 #else
205 throw DataException("Error - loadMesh: is not compiled with NetCFD. Please contact your installation manager.");
206 #endif /* USE_NETCDF */
207 }
208
209 AbstractContinuousDomain* readMesh(const std::string& fileName,
210 int integrationOrder,
211 int reducedIntegrationOrder,
212 int optimize)
213 {
214 //
215 // create a copy of the filename to overcome the non-constness of call
216 // to Finley_Mesh_read
217 Finley_Mesh* fMesh=0;
218 // Win32 refactor
219 char *fName = ((fileName.size()+1)>0) ? TMPMEMALLOC((fileName.size()+1),char) : (char*)NULL;
220 strcpy(fName,fileName.c_str());
221 double blocktimer_start = blocktimer_time();
222
223 fMesh=Finley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
224 checkFinleyError();
225 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
226
227 /* win32 refactor */
228 TMPMEMFREE(fName);
229
230 blocktimer_increment("ReadMesh()", blocktimer_start);
231 return temp;
232 }
233
234 AbstractContinuousDomain* readGmsh(const std::string& fileName,
235 int numDim,
236 int integrationOrder,
237 int reducedIntegrationOrder,
238 int optimize)
239 {
240 //
241 // create a copy of the filename to overcome the non-constness of call
242 // to Finley_Mesh_read
243 Finley_Mesh* fMesh=0;
244 // Win32 refactor
245 char *fName = ((fileName.size()+1)>0) ? TMPMEMALLOC((fileName.size()+1),char) : (char*)NULL;
246 strcpy(fName,fileName.c_str());
247 double blocktimer_start = blocktimer_time();
248
249 fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
250 checkFinleyError();
251 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
252
253 /* win32 refactor */
254 TMPMEMFREE(fName);
255
256 blocktimer_increment("ReadGmsh()", blocktimer_start);
257 return temp;
258 }
259
260 AbstractContinuousDomain* brick(int n0,int n1,int n2,int order,
261 double l0,double l1,double l2,
262 int periodic0,int periodic1,
263 int periodic2,
264 int integrationOrder,
265 int reducedIntegrationOrder,
266 int useElementsOnFace,
267 int useFullElementOrder,
268 int optimize)
269 {
270 int numElements[]={n0,n1,n2};
271 double length[]={l0,l1,l2};
272 int periodic[]={periodic0, periodic1, periodic2};
273
274 //
275 // linearInterpolation
276 Finley_Mesh* fMesh=NULL;
277
278 if (order==1) {
279 fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
280 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
281 }
282 else if (order==2) {
283 fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
284 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
285 } else {
286 stringstream temp;
287 temp << "Illegal interpolation order: " << order;
288 setFinleyError(VALUE_ERROR,temp.str().c_str());
289 }
290 //
291 // Convert any finley errors into a C++ exception
292 checkFinleyError();
293 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
294 return temp;
295 }
296 AbstractContinuousDomain* rectangle(int n0,int n1,int order,
297 double l0, double l1,
298 int periodic0,int periodic1,
299 int integrationOrder,
300 int reducedIntegrationOrder,
301 int useElementsOnFace,
302 int useFullElementOrder,
303 int optimize)
304 {
305 int numElements[]={n0,n1};
306 double length[]={l0,l1};
307 int periodic[]={periodic0, periodic1};
308
309 Finley_Mesh* fMesh=0;
310 if (order==1) {
311 fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
313 }
314 else if (order==2) {
315 fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
316 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
317 }
318 else {
319 stringstream temp;
320 temp << "Illegal interpolation order: " << order;
321 setFinleyError(VALUE_ERROR,temp.str().c_str());
322 }
323 //
324 // Convert any finley errors into a C++ exception
325 checkFinleyError();
326 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
327 return temp;
328 }
329
330 AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)
331 {
332 Finley_Mesh* fMesh=0;
333 //
334 // extract the meshes from meshList
335 int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
336 Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
337 for (int i=0;i<numMsh;++i) {
338 AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
339 const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
340 mshes[i]=finley_meshListMember->getFinley_Mesh();
341 }
342 //
343 // merge the meshes:
344 fMesh=Finley_Mesh_merge(numMsh,mshes);
345 TMPMEMFREE(mshes);
346 //
347 // Convert any finley errors into a C++ exception
348 checkFinleyError();
349 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
350
351 return temp;
352 }
353 AbstractContinuousDomain* glueFaces(const boost::python::list& meshList,
354 double safety_factor,
355 double tolerance,
356 int optimize)
357 {
358 Finley_Mesh* fMesh=0;
359 //
360 // merge the meshes:
361 AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
362 //
363 // glue the faces:
364 const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
365 fMesh=merged_finley_meshes->getFinley_Mesh();
366 Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
367
368 //
369 // Convert any finley errors into a C++ exception
370 checkFinleyError();
371 return merged_meshes;
372 }
373 AbstractContinuousDomain* joinFaces(const boost::python::list& meshList,
374 double safety_factor,
375 double tolerance,
376 int optimize)
377 {
378 Finley_Mesh* fMesh=0;
379 //
380 // merge the meshes:
381 AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
382 //
383 // join the faces:
384 const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
385 fMesh=merged_finley_meshes->getFinley_Mesh();
386 Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
387 //
388 // Convert any finley errors into a C++ exception
389 checkFinleyError();
390 return merged_meshes;
391 }
392
393 // end of namespace
394
395 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26