Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ReadSms.cpp
Go to the documentation of this file.
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 
52 {
53  return new ReadSms( iface );
54 }
55 
56 ReadSms::ReadSms( Interface* impl ) : mdbImpl( impl ), globalId( 0 ), paramCoords( 0 ), geomDimension( 0 ), setId( 0 )
57 {
59 }
60 
62 {
63  if( readMeshIface )
64  {
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 
111 
112  ErrorCode result =
114  CHECK( "Failed to create param coords tag." );
115 
116  int negone = -1;
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, &regions[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 
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