 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 /**
2 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3 * storing and accessing finite element mesh data.
4 *
5 * Copyright 2004 Sandia Corporation. Under the terms of Contract
6 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7 * retains certain rights in this software.
8 *
9 * This library is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public
11 * License as published by the Free Software Foundation; either
12 * version 2.1 of the License, or (at your option) any later version.
13 *
14 */
15
16 /**
17 * \class ReadSms
18 * \brief Sms (http://www.geuz.org/sms) file reader
19 * \author Jason Kraftcheck
20 */
21
22 #include "ReadSms.hpp"
23 #include "FileTokenizer.hpp" // For file tokenizer
24 #include "Internals.hpp"
25 #include "moab/Interface.hpp"
26 #include "moab/ReadUtilIface.hpp"
27 #include "moab/Range.hpp"
28 #include "MBTagConventions.hpp"
29 #include "MBParallelConventions.h"
30 #include "moab/CN.hpp"
31
32 #include <cerrno>
33 #include <cstring>
34 #include <map>
35 #include <set>
36 #include <iostream>
37
38 #define CHECK( a ) \
39 if( MB_SUCCESS != result ) \
40 { \
41 std::cerr << ( a ) << std::endl; \
42 return result; \
43 }
44
45 #define CHECKN( a ) \
46 if( n != ( a ) ) return MB_FILE_WRITE_ERROR
47
48 namespace moab
49 {
50
51 ReaderIface* ReadSms::factory( Interface* iface )
52 {
53 return new ReadSms( iface );
54 }
55
56 ReadSms::ReadSms( Interface* impl ) : mdbImpl( impl ), globalId( 0 ), paramCoords( 0 ), geomDimension( 0 ), setId( 0 )
57 {
58 mdbImpl->query_interface( readMeshIface );
59 }
60
61 ReadSms::~ReadSms()
62 {
63 if( readMeshIface )
64 {
65 mdbImpl->release_interface( readMeshIface );
66 readMeshIface = 0;
67 }
68 }
69
70 ErrorCode ReadSms::read_tag_values( const char* /* file_name */,
71 const char* /* tag_name */,
72 const FileOptions& /* opts */,
73 std::vector< int >& /* tag_values_out */,
74 const SubsetList* /* subset_list */ )
75 {
76 return MB_NOT_IMPLEMENTED;
77 }
78
79 ErrorCode ReadSms::load_file( const char* filename,
80 const EntityHandle* /* file_set */,
81 const FileOptions& /* opts */,
82 const ReaderIface::SubsetList* subset_list,
83 const Tag* file_id_tag )
84 {
85 if( subset_list )
86 {
87 MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for Sms" );
88 }
89
90 setId = 1;
91
92 // Open file
93 FILE* file_ptr = fopen( filename, "r" );
94 if( !file_ptr )
95 {
96 MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, filename << ": " << strerror( errno ) );
97 }
98
99 const ErrorCode result = load_file_impl( file_ptr, file_id_tag );
100 fclose( file_ptr );
101
102 return result;
103 }
104
105 ErrorCode ReadSms::load_file_impl( FILE* file_ptr, const Tag* file_id_tag )
106 {
107 bool warned = false;
108 double dum_params[] = { 0.0, 0.0, 0.0 };
109
110 globalId = mdbImpl->globalId_tag();
111
112 ErrorCode result =
113 mdbImpl->tag_get_handle( "PARAMETER_COORDS", 3, MB_TYPE_DOUBLE, paramCoords, MB_TAG_DENSE | MB_TAG_CREAT );
114 CHECK( "Failed to create param coords tag." );
115
116 int negone = -1;
117 result = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomDimension,
118 MB_TAG_SPARSE | MB_TAG_CREAT, &negone );
119 CHECK( "Failed to create geom dim tag." );
120
121 int n;
122 char line[256], all_line[1024];
123 int file_type;
124
125 if( fgets( all_line, sizeof( all_line ), file_ptr ) == NULL )
126 {
127 return MB_FAILURE;
128 }
129 if( sscanf( all_line, "%s %d", line, &file_type ) != 2 )
130 {
131 return MB_FAILURE;
132 }
133
134 if( 3 == file_type )
135 {
136 result = read_parallel_info( file_ptr );
137 CHECK( "Failed to read parallel info." );
138 }
139
140 int nregions, nfaces, nedges, nvertices, npoints;
141 n = fscanf( file_ptr, "%d %d %d %d %d", &nregions, &nfaces, &nedges, &nvertices, &npoints );
142 CHECKN( 5 );
143 if( nregions < 0 || nfaces < 0 || nedges < 0 || nvertices < 0 || npoints < 0 ) return MB_FILE_WRITE_ERROR;
144
145 // Create the vertices
146 std::vector< double* > coord_arrays;
147 EntityHandle vstart = 0;
148 result = readMeshIface->get_node_coords( 3, nvertices, MB_START_ID, vstart, coord_arrays );
149 CHECK( "Failed to get node arrays." );
150
151 if( file_id_tag )
152 {
153 result = add_entities( vstart, nvertices, file_id_tag );MB_CHK_ERR( result );
154 }
155
156 EntityHandle this_gent, new_handle;
157 std::vector< EntityHandle > gentities[4];
158 int gent_id, dum_int;
159 int gent_type, num_connections;
160
161 for( int i = 0; i < nvertices; i++ )
162 {
163 n = fscanf( file_ptr, "%d", &gent_id );
164 CHECKN( 1 );
165 if( !gent_id ) continue;
166
167 n = fscanf( file_ptr, "%d %d %lf %lf %lf", &gent_type, &num_connections, coord_arrays[0] + i,
168 coord_arrays[1] + i, coord_arrays[2] + i );
169 CHECKN( 5 );
170
171 result = get_set( gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag );MB_CHK_ERR( result );
172
173 new_handle = vstart + i;
174 result = mdbImpl->add_entities( this_gent, &new_handle, 1 );
175 CHECK( "Adding vertex to geom set failed." );
176
177 switch( gent_type )
178 {
179 case 1:
180 n = fscanf( file_ptr, "%le", dum_params );
181 CHECKN( 1 );
182 result = mdbImpl->tag_set_data( paramCoords, &new_handle, 1, dum_params );
183 CHECK( "Failed to set param coords tag for vertex." );
184 break;
185 case 2:
186 n = fscanf( file_ptr, "%le %le %d", dum_params, dum_params + 1, &dum_int );
187 CHECKN( 3 );
188 dum_params[2] = dum_int;
189 result = mdbImpl->tag_set_data( paramCoords, &new_handle, 1, dum_params );
190 CHECK( "Failed to set param coords tag for vertex." );
191 break;
192 default:
193 break;
194 }
195 } // End of reading vertices
196
197 // *******************************
198 // Read Edges
199 // *******************************
200
201 int vert1, vert2, num_pts;
202 std::vector< EntityHandle > everts( 2 );
203 EntityHandle estart, *connect;
204 result = readMeshIface->get_element_connect( nedges, 2, MBEDGE, 1, estart, connect );
205 CHECK( "Failed to create array of edges." );
206
207 if( file_id_tag )
208 {
209 result = add_entities( estart, nedges, file_id_tag );
210 if( MB_SUCCESS != result ) return result;
211 }
212
213 for( int i = 0; i < nedges; i++ )
214 {
215 n = fscanf( file_ptr, "%d", &gent_id );
216 CHECKN( 1 );
217 if( !gent_id ) continue;
218
219 n = fscanf( file_ptr, "%d %d %d %d %d", &gent_type, &vert1, &vert2, &num_connections, &num_pts );
220 CHECKN( 5 );
221 if( vert1 < 1 || vert1 > nvertices ) return MB_FILE_WRITE_ERROR;
222 if( vert2 < 1 || vert2 > nvertices ) return MB_FILE_WRITE_ERROR;
223
224 connect[0] = vstart + vert1 - 1;
225 connect[1] = vstart + vert2 - 1;
226 if( num_pts > 1 && !warned )
227 {
228 std::cout << "Warning: num_points > 1 not supported; choosing last one." << std::endl;
229 warned = true;
230 }
231
232 result = get_set( gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag );
233 CHECK( "Problem getting geom set for edge." );
234
235 new_handle = estart + i;
236 result = mdbImpl->add_entities( this_gent, &new_handle, 1 );
237 CHECK( "Failed to add edge to geom set." );
238
239 connect += 2;
240
241 for( int j = 0; j < num_pts; j++ )
242 {
243 switch( gent_type )
244 {
245 case 1:
246 n = fscanf( file_ptr, "%le", dum_params );
247 CHECKN( 1 );
248 result = mdbImpl->tag_set_data( paramCoords, &new_handle, 1, dum_params );
249 CHECK( "Failed to set param coords tag for edge." );
250 break;
251 case 2:
252 n = fscanf( file_ptr, "%le %le %d", dum_params, dum_params + 1, &dum_int );
253 CHECKN( 3 );
254 dum_params[2] = dum_int;
255 result = mdbImpl->tag_set_data( paramCoords, &new_handle, 1, dum_params );
256 CHECK( "Failed to set param coords tag for edge." );
257 break;
258 default:
259 break;
260 }
261 }
262 } // End of reading edges
263
264 // *******************************
265 // Read Faces
266 // *******************************
267 std::vector< EntityHandle > bound_ents, bound_verts, new_faces;
268 int bound_id;
269 Range shverts;
270 new_faces.resize( nfaces );
271 int num_bounding;
272
273 for( int i = 0; i < nfaces; i++ )
274 {
275 n = fscanf( file_ptr, "%d", &gent_id );
276 CHECKN( 1 );
277 if( !gent_id ) continue;
278
279 n = fscanf( file_ptr, "%d %d", &gent_type, &num_bounding );
280 CHECKN( 2 );
281
282 result = get_set( gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag );
283 CHECK( "Problem getting geom set for face." );
284
285 bound_ents.resize( num_bounding + 1 );
286 bound_verts.resize( num_bounding );
287 for( int j = 0; j < num_bounding; j++ )
288 {
289 n = fscanf( file_ptr, "%d ", &bound_id );
290 CHECKN( 1 );
291 if( 0 > bound_id ) bound_id = abs( bound_id );
292 assert( 0 < bound_id && bound_id <= nedges );
293 if( bound_id < 1 || bound_id > nedges ) return MB_FILE_WRITE_ERROR;
294 bound_ents[j] = estart + abs( bound_id ) - 1;
295 }
296
297 // Convert edge-based model to vertex-based one
298 for( int j = 0; j < num_bounding; j++ )
299 {
300 if( j == num_bounding - 1 ) bound_ents[j + 1] = bound_ents[0];
301 result = mdbImpl->get_adjacencies( &bound_ents[j], 2, 0, false, shverts );
302 CHECK( "Failed to get vertices bounding edge." );
303 assert( shverts.size() == 1 );
304 bound_verts[j] = *shverts.begin();
305 shverts.clear();
306 }
307
308 result = mdbImpl->create_element( (EntityType)( MBTRI + num_bounding - 3 ), &bound_verts[0], bound_verts.size(),
309 new_faces[i] );
310 CHECK( "Failed to create edge." );
311
312 result = mdbImpl->add_entities( this_gent, &new_faces[i], 1 );
313 CHECK( "Failed to add edge to geom set." );
314
315 int num_read = fscanf( file_ptr, "%d", &num_pts );
316 if( !num_pts || !num_read ) continue;
317
318 for( int j = 0; j < num_pts; j++ )
319 {
320 switch( gent_type )
321 {
322 case 1:
323 n = fscanf( file_ptr, "%le", dum_params );
324 CHECKN( 1 );
325 result = mdbImpl->tag_set_data( paramCoords, &new_faces[i], 1, dum_params );
326 CHECK( "Failed to set param coords tag for face." );
327 break;
328 case 2:
329 n = fscanf( file_ptr, "%le %le %d", dum_params, dum_params + 1, &dum_int );
330 CHECKN( 3 );
331 dum_params[2] = dum_int;
332 result = mdbImpl->tag_set_data( paramCoords, &new_faces[i], 1, dum_params );
333 CHECK( "Failed to set param coords tag for face." );
334 break;
335 default:
336 break;
337 }
338 }
339 } // End of reading faces
340
341 if( file_id_tag )
342 {
343 result = readMeshIface->assign_ids( *file_id_tag, &new_faces[0], new_faces.size(), 1 );
344 if( MB_SUCCESS != result ) return result;
345 }
346
347 // *******************************
348 // Read Regions
349 // *******************************
350 int sense[MAX_SUB_ENTITIES];
351 bound_verts.resize( MAX_SUB_ENTITIES );
352
353 std::vector< EntityHandle > regions;
354 if( file_id_tag ) regions.resize( nregions );
355 for( int i = 0; i < nregions; i++ )
356 {
357 n = fscanf( file_ptr, "%d", &gent_id );
358 CHECKN( 1 );
359 if( !gent_id ) continue;
360 result = get_set( gentities, 3, gent_id, geomDimension, this_gent, file_id_tag );
361 CHECK( "Couldn't get geom set for region." );
362 n = fscanf( file_ptr, "%d", &num_bounding );
363 CHECKN( 1 );
364 bound_ents.resize( num_bounding );
365 for( int j = 0; j < num_bounding; j++ )
366 {
367 n = fscanf( file_ptr, "%d ", &bound_id );
368 CHECKN( 1 );
369 assert( abs( bound_id ) < (int)new_faces.size() + 1 && bound_id );
370 if( !bound_id || abs( bound_id ) > nfaces ) return MB_FILE_WRITE_ERROR;
371 sense[j] = ( bound_id < 0 ) ? -1 : 1;
372 bound_ents[j] = new_faces[abs( bound_id ) - 1];
373 }
374
375 EntityType etype;
376 result = readMeshIface->get_ordered_vertices( &bound_ents[0], sense, num_bounding, 3, &bound_verts[0], etype );
377 CHECK( "Failed in get_ordered_vertices." );
378
379 // Make the element
380 result = mdbImpl->create_element( etype, &bound_verts[0], CN::VerticesPerEntity( etype ), new_handle );
381 CHECK( "Failed to create region." );
382
383 result = mdbImpl->add_entities( this_gent, &new_handle, 1 );
384 CHECK( "Failed to add region to geom set." );
385
386 if( file_id_tag ) regions[i] = new_handle;
387
388 n = fscanf( file_ptr, "%d ", &dum_int );
389 CHECKN( 1 );
390 } // End of reading regions
391
392 if( file_id_tag )
393 {
394 result = readMeshIface->assign_ids( *file_id_tag, ®ions[0], regions.size(), 1 );
395 if( MB_SUCCESS != result ) return result;
396 }
397
398 return MB_SUCCESS;
399 }
400
401 ErrorCode ReadSms::get_set( std::vector< EntityHandle >* sets,
402 int set_dim,
403 int set_id,
404 Tag dim_tag,
405 EntityHandle& this_set,
406 const Tag* file_id_tag )
407 {
408 ErrorCode result = MB_SUCCESS;
409
410 if( set_dim < 0 || set_dim > 3 ) return MB_FILE_WRITE_ERROR;
411
412 if( (int)sets[set_dim].size() <= set_id || !sets[set_dim][set_id] )
413 {
414 if( (int)sets[set_dim].size() <= set_id ) sets[set_dim].resize( set_id + 1, 0 );
415
416 if( !sets[set_dim][set_id] )
417 {
418 result = mdbImpl->create_meshset( MESHSET_SET, sets[set_dim][set_id] );
419 if( MB_SUCCESS != result ) return result;
420 result = mdbImpl->tag_set_data( globalId, &sets[set_dim][set_id], 1, &set_id );
421 if( MB_SUCCESS != result ) return result;
422 result = mdbImpl->tag_set_data( dim_tag, &sets[set_dim][set_id], 1, &set_dim );
423 if( MB_SUCCESS != result ) return result;
424
425 if( file_id_tag )
426 {
427 result = mdbImpl->tag_set_data( *file_id_tag, &sets[set_dim][set_id], 1, &setId );
428 ++setId;
429 }
430 }
431 }
432
433 this_set = sets[set_dim][set_id];
434
435 return result;
436 }
437
438 ErrorCode ReadSms::read_parallel_info( FILE* file_ptr )
439 {
440 // ErrorCode result;
441
442 // Read partition info
443 int nparts, part_id, num_ifaces, num_corner_ents;
444 int num_read = fscanf( file_ptr, "%d %d %d %d", &nparts, &part_id, &num_ifaces, &num_corner_ents );
445 if( !num_read ) return MB_FAILURE;
446
447 // Read interfaces
448 int iface_id, iface_dim, iface_own, num_iface_corners;
449 // EntityHandle this_iface;
450 std::vector< int >* iface_corners = NULL;
451 for( int i = 0; i < num_ifaces; i++ )
452 {
453 num_read = fscanf( file_ptr, "%d %d %d %d", &iface_id, &iface_dim, &iface_own, &num_iface_corners );
454 if( !num_read ) return MB_FAILURE;
455
456 // result = get_set(sets, iface_dim, iface_id, dim_tag, iface_own, this_iface);
457 // CHECK("Failed to make iface set.");
458
459 // Read the corner ids and store them on the set for now
460 iface_corners = new std::vector< int >( num_iface_corners );
461 for( int j = 0; j < num_iface_corners; j++ )
462 {
463 num_read = fscanf( file_ptr, "%d", &( *iface_corners )[j] );
464 if( !num_read )
465 {
466 delete iface_corners;
467 return MB_FAILURE;
468 }
469 }
470
471 // result = tag_set_data(ifaceCornerTag, &this_iface, 1,
472 //&iface_corners);
473 // CHECK("Failed to set iface corner tag.");
474
475 delete iface_corners;
476 iface_corners = NULL;
477 }
478
479 // Interface data has been read
480 return MB_SUCCESS;
481 }
482
483 ErrorCode ReadSms::add_entities( EntityHandle start, EntityHandle count, const Tag* file_id_tag )
484 {
485 if( !count || !file_id_tag ) return MB_FAILURE;
486
487 Range range;
488 range.insert( start, start + count - 1 );
489 return readMeshIface->assign_ids( *file_id_tag, range, 1 );
490 }
491
492 } // namespace moab
1.9.1.