/[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 3259 - (hide annotations)
Mon Oct 11 01:48:14 2010 UTC (8 years, 11 months ago) by jfenwick
File size: 31626 byte(s)
Merging dudley and scons updates from branches

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26