Mesh Oriented datABase  (version 5.6.0)
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 
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 );
154  MB_CHK_ERR( result );
155  }
156 
157  EntityHandle this_gent, new_handle;
158  std::vector< EntityHandle > gentities[4];
159  int gent_id, dum_int;
160  int gent_type, num_connections;
161 
162  for( int i = 0; i < nvertices; i++ )
163  {
164  n = fscanf( file_ptr, "%d", &gent_id );
165  CHECKN( 1 );
166  if( !gent_id ) continue;
167 
168  n = fscanf( file_ptr, "%d %d %lf %lf %lf", &gent_type, &num_connections, coord_arrays[0] + i,
169  coord_arrays[1] + i, coord_arrays[2] + i );
170  CHECKN( 5 );
171 
172  result = get_set( gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag );
173  MB_CHK_ERR( result );
174 
175  new_handle = vstart + i;
176  result = mdbImpl->add_entities( this_gent, &new_handle, 1 );
177  CHECK( "Adding vertex to geom set failed." );
178 
179  switch( gent_type )
180  {
181  case 1:
182  n = fscanf( file_ptr, "%le", dum_params );
183  CHECKN( 1 );
184  result = mdbImpl->tag_set_data( paramCoords, &new_handle, 1, dum_params );
185  CHECK( "Failed to set param coords tag for vertex." );
186  break;
187  case 2:
188  n = fscanf( file_ptr, "%le %le %d", dum_params, dum_params + 1, &dum_int );
189  CHECKN( 3 );
190  dum_params[2] = dum_int;
191  result = mdbImpl->tag_set_data( paramCoords, &new_handle, 1, dum_params );
192  CHECK( "Failed to set param coords tag for vertex." );
193  break;
194  default:
195  break;
196  }
197  } // End of reading vertices
198 
199  // *******************************
200  // Read Edges
201  // *******************************
202 
203  int vert1, vert2, num_pts;
204  std::vector< EntityHandle > everts( 2 );
205  EntityHandle estart, *connect;
206  result = readMeshIface->get_element_connect( nedges, 2, MBEDGE, 1, estart, connect );
207  CHECK( "Failed to create array of edges." );
208 
209  if( file_id_tag )
210  {
211  result = add_entities( estart, nedges, file_id_tag );
212  if( MB_SUCCESS != result ) return result;
213  }
214 
215  for( int i = 0; i < nedges; i++ )
216  {
217  n = fscanf( file_ptr, "%d", &gent_id );
218  CHECKN( 1 );
219  if( !gent_id ) continue;
220 
221  n = fscanf( file_ptr, "%d %d %d %d %d", &gent_type, &vert1, &vert2, &num_connections, &num_pts );
222  CHECKN( 5 );
223  if( vert1 < 1 || vert1 > nvertices ) return MB_FILE_WRITE_ERROR;
224  if( vert2 < 1 || vert2 > nvertices ) return MB_FILE_WRITE_ERROR;
225 
226  connect[0] = vstart + vert1 - 1;
227  connect[1] = vstart + vert2 - 1;
228  if( num_pts > 1 && !warned )
229  {
230  std::cout << "Warning: num_points > 1 not supported; choosing last one." << std::endl;
231  warned = true;
232  }
233 
234  result = get_set( gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag );
235  CHECK( "Problem getting geom set for edge." );
236 
237  new_handle = estart + i;
238  result = mdbImpl->add_entities( this_gent, &new_handle, 1 );
239  CHECK( "Failed to add edge to geom set." );
240 
241  connect += 2;
242 
243  for( int j = 0; j < num_pts; j++ )
244  {
245  switch( gent_type )
246  {
247  case 1:
248  n = fscanf( file_ptr, "%le", dum_params );
249  CHECKN( 1 );
250  result = mdbImpl->tag_set_data( paramCoords, &new_handle, 1, dum_params );
251  CHECK( "Failed to set param coords tag for edge." );
252  break;
253  case 2:
254  n = fscanf( file_ptr, "%le %le %d", dum_params, dum_params + 1, &dum_int );
255  CHECKN( 3 );
256  dum_params[2] = dum_int;
257  result = mdbImpl->tag_set_data( paramCoords, &new_handle, 1, dum_params );
258  CHECK( "Failed to set param coords tag for edge." );
259  break;
260  default:
261  break;
262  }
263  }
264  } // End of reading edges
265 
266  // *******************************
267  // Read Faces
268  // *******************************
269  std::vector< EntityHandle > bound_ents, bound_verts, new_faces;
270  int bound_id;
271  Range shverts;
272  new_faces.resize( nfaces );
273  int num_bounding;
274 
275  for( int i = 0; i < nfaces; i++ )
276  {
277  n = fscanf( file_ptr, "%d", &gent_id );
278  CHECKN( 1 );
279  if( !gent_id ) continue;
280 
281  n = fscanf( file_ptr, "%d %d", &gent_type, &num_bounding );
282  CHECKN( 2 );
283 
284  result = get_set( gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag );
285  CHECK( "Problem getting geom set for face." );
286 
287  bound_ents.resize( num_bounding + 1 );
288  bound_verts.resize( num_bounding );
289  for( int j = 0; j < num_bounding; j++ )
290  {
291  n = fscanf( file_ptr, "%d ", &bound_id );
292  CHECKN( 1 );
293  if( 0 > bound_id ) bound_id = abs( bound_id );
294  assert( 0 < bound_id && bound_id <= nedges );
295  if( bound_id < 1 || bound_id > nedges ) return MB_FILE_WRITE_ERROR;
296  bound_ents[j] = estart + abs( bound_id ) - 1;
297  }
298 
299  // Convert edge-based model to vertex-based one
300  for( int j = 0; j < num_bounding; j++ )
301  {
302  if( j == num_bounding - 1 ) bound_ents[j + 1] = bound_ents[0];
303  result = mdbImpl->get_adjacencies( &bound_ents[j], 2, 0, false, shverts );
304  CHECK( "Failed to get vertices bounding edge." );
305  assert( shverts.size() == 1 );
306  bound_verts[j] = *shverts.begin();
307  shverts.clear();
308  }
309 
310  result = mdbImpl->create_element( (EntityType)( MBTRI + num_bounding - 3 ), &bound_verts[0], bound_verts.size(),
311  new_faces[i] );
312  CHECK( "Failed to create edge." );
313 
314  result = mdbImpl->add_entities( this_gent, &new_faces[i], 1 );
315  CHECK( "Failed to add edge to geom set." );
316 
317  int num_read = fscanf( file_ptr, "%d", &num_pts );
318  if( !num_pts || !num_read ) continue;
319 
320  for( int j = 0; j < num_pts; j++ )
321  {
322  switch( gent_type )
323  {
324  case 1:
325  n = fscanf( file_ptr, "%le", dum_params );
326  CHECKN( 1 );
327  result = mdbImpl->tag_set_data( paramCoords, &new_faces[i], 1, dum_params );
328  CHECK( "Failed to set param coords tag for face." );
329  break;
330  case 2:
331  n = fscanf( file_ptr, "%le %le %d", dum_params, dum_params + 1, &dum_int );
332  CHECKN( 3 );
333  dum_params[2] = dum_int;
334  result = mdbImpl->tag_set_data( paramCoords, &new_faces[i], 1, dum_params );
335  CHECK( "Failed to set param coords tag for face." );
336  break;
337  default:
338  break;
339  }
340  }
341  } // End of reading faces
342 
343  if( file_id_tag )
344  {
345  result = readMeshIface->assign_ids( *file_id_tag, &new_faces[0], new_faces.size(), 1 );
346  if( MB_SUCCESS != result ) return result;
347  }
348 
349  // *******************************
350  // Read Regions
351  // *******************************
352  int sense[MAX_SUB_ENTITIES];
353  bound_verts.resize( MAX_SUB_ENTITIES );
354 
355  std::vector< EntityHandle > regions;
356  if( file_id_tag ) regions.resize( nregions );
357  for( int i = 0; i < nregions; i++ )
358  {
359  n = fscanf( file_ptr, "%d", &gent_id );
360  CHECKN( 1 );
361  if( !gent_id ) continue;
362  result = get_set( gentities, 3, gent_id, geomDimension, this_gent, file_id_tag );
363  CHECK( "Couldn't get geom set for region." );
364  n = fscanf( file_ptr, "%d", &num_bounding );
365  CHECKN( 1 );
366  bound_ents.resize( num_bounding );
367  for( int j = 0; j < num_bounding; j++ )
368  {
369  n = fscanf( file_ptr, "%d ", &bound_id );
370  CHECKN( 1 );
371  assert( abs( bound_id ) < (int)new_faces.size() + 1 && bound_id );
372  if( !bound_id || abs( bound_id ) > nfaces ) return MB_FILE_WRITE_ERROR;
373  sense[j] = ( bound_id < 0 ) ? -1 : 1;
374  bound_ents[j] = new_faces[abs( bound_id ) - 1];
375  }
376 
377  EntityType etype;
378  result = readMeshIface->get_ordered_vertices( &bound_ents[0], sense, num_bounding, 3, &bound_verts[0], etype );
379  CHECK( "Failed in get_ordered_vertices." );
380 
381  // Make the element
382  result = mdbImpl->create_element( etype, &bound_verts[0], CN::VerticesPerEntity( etype ), new_handle );
383  CHECK( "Failed to create region." );
384 
385  result = mdbImpl->add_entities( this_gent, &new_handle, 1 );
386  CHECK( "Failed to add region to geom set." );
387 
388  if( file_id_tag ) regions[i] = new_handle;
389 
390  n = fscanf( file_ptr, "%d ", &dum_int );
391  CHECKN( 1 );
392  } // End of reading regions
393 
394  if( file_id_tag )
395  {
396  result = readMeshIface->assign_ids( *file_id_tag, &regions[0], regions.size(), 1 );
397  if( MB_SUCCESS != result ) return result;
398  }
399 
400  return MB_SUCCESS;
401 }
402 
403 ErrorCode ReadSms::get_set( std::vector< EntityHandle >* sets,
404  int set_dim,
405  int set_id,
406  Tag dim_tag,
407  EntityHandle& this_set,
408  const Tag* file_id_tag )
409 {
410  ErrorCode result = MB_SUCCESS;
411 
412  if( set_dim < 0 || set_dim > 3 ) return MB_FILE_WRITE_ERROR;
413 
414  if( (int)sets[set_dim].size() <= set_id || !sets[set_dim][set_id] )
415  {
416  if( (int)sets[set_dim].size() <= set_id ) sets[set_dim].resize( set_id + 1, 0 );
417 
418  if( !sets[set_dim][set_id] )
419  {
420  result = mdbImpl->create_meshset( MESHSET_SET, sets[set_dim][set_id] );
421  if( MB_SUCCESS != result ) return result;
422  result = mdbImpl->tag_set_data( globalId, &sets[set_dim][set_id], 1, &set_id );
423  if( MB_SUCCESS != result ) return result;
424  result = mdbImpl->tag_set_data( dim_tag, &sets[set_dim][set_id], 1, &set_dim );
425  if( MB_SUCCESS != result ) return result;
426 
427  if( file_id_tag )
428  {
429  result = mdbImpl->tag_set_data( *file_id_tag, &sets[set_dim][set_id], 1, &setId );
430  ++setId;
431  }
432  }
433  }
434 
435  this_set = sets[set_dim][set_id];
436 
437  return result;
438 }
439 
440 ErrorCode ReadSms::read_parallel_info( FILE* file_ptr )
441 {
442  // ErrorCode result;
443 
444  // Read partition info
445  int nparts, part_id, num_ifaces, num_corner_ents;
446  int num_read = fscanf( file_ptr, "%d %d %d %d", &nparts, &part_id, &num_ifaces, &num_corner_ents );
447  if( !num_read ) return MB_FAILURE;
448 
449  // Read interfaces
450  int iface_id, iface_dim, iface_own, num_iface_corners;
451  // EntityHandle this_iface;
452  std::vector< int >* iface_corners = NULL;
453  for( int i = 0; i < num_ifaces; i++ )
454  {
455  num_read = fscanf( file_ptr, "%d %d %d %d", &iface_id, &iface_dim, &iface_own, &num_iface_corners );
456  if( !num_read ) return MB_FAILURE;
457 
458  // result = get_set(sets, iface_dim, iface_id, dim_tag, iface_own, this_iface);
459  // CHECK("Failed to make iface set.");
460 
461  // Read the corner ids and store them on the set for now
462  iface_corners = new std::vector< int >( num_iface_corners );
463  for( int j = 0; j < num_iface_corners; j++ )
464  {
465  num_read = fscanf( file_ptr, "%d", &( *iface_corners )[j] );
466  if( !num_read )
467  {
468  delete iface_corners;
469  return MB_FAILURE;
470  }
471  }
472 
473  // result = tag_set_data(ifaceCornerTag, &this_iface, 1,
474  //&iface_corners);
475  // CHECK("Failed to set iface corner tag.");
476 
477  delete iface_corners;
478  iface_corners = NULL;
479  }
480 
481  // Interface data has been read
482  return MB_SUCCESS;
483 }
484 
485 ErrorCode ReadSms::add_entities( EntityHandle start, EntityHandle count, const Tag* file_id_tag )
486 {
487  if( !count || !file_id_tag ) return MB_FAILURE;
488 
489  Range range;
490  range.insert( start, start + count - 1 );
491  return readMeshIface->assign_ids( *file_id_tag, range, 1 );
492 }
493 
494 } // namespace moab