Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ReadIDEAS.cpp
Go to the documentation of this file.
1 #include <iostream>
2 #include <fstream>
3 #include <vector>
4 #include <cstdlib>
5 #include <sstream>
6 #include <cassert>
7 
8 #include "ReadIDEAS.hpp"
9 #include "moab/Interface.hpp"
10 #include "Internals.hpp"
11 #include "moab/ReadUtilIface.hpp"
12 #include "FileTokenizer.hpp"
13 #include "MBTagConventions.hpp"
14 #include "moab/Range.hpp"
15 #include "moab/CN.hpp"
16 
17 namespace moab
18 {
19 
21 {
22  return new ReadIDEAS( iface );
23 }
24 
26 {
28 }
29 
30 ErrorCode ReadIDEAS::read_tag_values( const char* /* file_name */,
31  const char* /* tag_name */,
32  const FileOptions& /* opts */,
33  std::vector< int >& /* tag_values_out */,
34  const SubsetList* /* subset_list */ )
35 {
36  return MB_NOT_IMPLEMENTED;
37 }
38 
39 ErrorCode ReadIDEAS::load_file( const char* fname,
40  const EntityHandle*,
41  const FileOptions& /*options*/,
42  const ReaderIface::SubsetList* subset_list,
43  const Tag* file_id_tag )
44 {
45  if( subset_list )
46  {
47  MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for IDEAS" );
48  }
49 
50  file.open( fname );
51  if( !file.good() )
52  {
53  MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "Failed to open file: " << fname );
54  }
55 
56  ErrorCode rval;
57 
58  char line[10000];
59  file.getline( line, 10000 );
60  char* liter = line;
61  while( *liter && isspace( *liter ) )
62  ++liter;
63  if( *liter != '-' ) return MB_FAILURE;
64  ++liter;
65  if( *liter != '1' ) return MB_FAILURE;
66  while( *++liter )
67  if( !isspace( *liter ) ) return MB_FAILURE;
68 
69  EntityHandle first_vertex = 0;
70  while( !file.eof() )
71  {
72  file.getline( line, 10000 );
73  unsigned int header_id = (unsigned int)strtol( line, NULL, 10 );
74 
75  // Create vertices
76  if( DOUBLE_PRECISION_NODES0 == header_id || DOUBLE_PRECISION_NODES1 == header_id )
77  {
78  if( first_vertex ) // multiple vertex blocks?
79  return MB_FAILURE;
80  rval = create_vertices( first_vertex, file_id_tag );MB_CHK_SET_ERR( rval, "Failed to read vertices" );
81  }
82  // Create elements
83  else if( ELEMENTS0 == header_id || ELEMENTS1 == header_id || ELEMENTS2 == header_id )
84  {
85  if( !first_vertex ) // Need to read vertices first
86  return MB_FAILURE;
87  rval = create_elements( first_vertex, file_id_tag );MB_CHK_SET_ERR( rval, "Failed to read elements" );
88  }
89  // Skip everything else
90  else
91  {
92  rval = skip_header();
93  if( MB_SUCCESS != rval ) return MB_FAILURE;
94  }
95  }
96 
97  file.close();
98  return MB_SUCCESS;
99 }
100 
102 {
103  // Go until finding a pair of -1 lines
104  char* ctmp;
105  char line[10000];
106  std::string s;
107 
108  int end_of_block = 0;
109 
110  long int il;
111 
112  while( file.getline( line, 10000 ) )
113  {
114  il = std::strtol( line, &ctmp, 10 );
115  if( il == -1 )
116  {
117  s = ctmp;
118  if( s.empty() ) end_of_block++;
119  }
120  else
121  end_of_block = 0;
122 
123  if( end_of_block >= 2 ) return MB_SUCCESS;
124  }
125 
126  return MB_SUCCESS;
127 }
128 
129 ErrorCode ReadIDEAS::create_vertices( EntityHandle& first_vertex, const Tag* file_id_tag )
130 {
131  // Read two lines: first has some data, second has coordinates
132  char line1[10000], line2[10000];
133  int il1, il2;
134  char *ctmp1, *ctmp2;
135  std::string s1, s2;
136 
137  ErrorCode rval;
138 
139  std::streampos top_of_block = file.tellg();
140  unsigned int num_verts = 0;
141 
142  for( ;; )
143  {
144  // Read both lines
145  if( !file.getline( line1, 10000 ) ) return MB_FAILURE;
146  if( !file.getline( line2, 10000 ) ) return MB_FAILURE;
147 
148  // Check if we are at the end of the block
149  il1 = std::strtol( line1, &ctmp1, 10 );
150  il2 = std::strtol( line2, &ctmp2, 10 );
151  if( ( il1 == -1 ) && ( il2 == -1 ) )
152  {
153  s1 = ctmp1;
154  s2 = ctmp2;
155  if( ( s1.empty() ) && ( s2.empty() ) ) break;
156  }
157  num_verts++;
158  }
159 
160  file.seekg( top_of_block );
161 
162  std::vector< double* > arrays;
163  rval = readMeshIface->get_node_coords( 3, num_verts, 0, first_vertex, arrays );
164  if( MB_SUCCESS != rval ) return rval;
165 
166  Range verts;
167  verts.insert( first_vertex, first_vertex + num_verts - 1 );
168 
169  double* x = arrays[0];
170  double* y = arrays[1];
171  double* z = arrays[2];
172 
173  // For now, assume ids are sequential and begin with 1
174  Tag id_tag = MBI->globalId_tag();
175  const int beginning_node_id = 1;
176  int node_id = beginning_node_id;
177 
178  for( unsigned int i = 0; i < num_verts; i++ )
179  {
180  if( !file.getline( line1, 10000 ) ) return MB_FAILURE;
181  if( !file.getline( line2, 10000 ) ) return MB_FAILURE;
182 
183  // Get the id out of the 1st line. Check the assumption that node ids are
184  // sequential and begin with 1.
185  if( node_id != std::strtol( line1, &ctmp1, 10 ) )
186  MB_SET_ERR( MB_FAILURE, "node_id " << node_id << " line2:" << line2 << " ctmp1:" << ctmp1 );
187  else
188  ++node_id;
189 
190  // Get the doubles out of the 2nd line
191  x[i] = std::strtod( line2, &ctmp2 );
192  y[i] = std::strtod( ctmp2 + 1, &ctmp2 );
193  z[i] = std::strtod( ctmp2 + 1, NULL );
194  }
195 
196  if( !file.getline( line1, 10000 ) ) MB_SET_ERR( MB_FAILURE, " expect more lines" );
197  if( !file.getline( line2, 10000 ) ) MB_SET_ERR( MB_FAILURE, " expect more lines 2" );
198 
199  // Tag the nodes with ids
200  rval = readMeshIface->assign_ids( id_tag, verts, beginning_node_id );MB_CHK_SET_ERR( rval, "Failed to assign IDs" );
201  if( file_id_tag )
202  {
203  rval = readMeshIface->assign_ids( *file_id_tag, verts, beginning_node_id );MB_CHK_SET_ERR( rval, "Failed to assign file IDs" );
204  }
205 
206  return MB_SUCCESS;
207 }
208 
210 {
211  char line1[10000], line2[10000];
212  int il1, il2;
213  char *ctmp1, *ctmp2;
214  std::string s1, s2;
215  ErrorCode rval;
216  EntityHandle handle;
217 
218  Tag mat_tag, phys_tag, id_tag;
220  if( MB_SUCCESS != rval && MB_ALREADY_ALLOCATED != rval ) return rval;
222  if( MB_SUCCESS != rval && MB_ALREADY_ALLOCATED != rval ) return rval;
223  id_tag = MBI->globalId_tag();
224 
225  for( ;; )
226  {
227  if( !file.getline( line1, 10000 ) || !file.getline( line2, 10000 ) ) return MB_FAILURE;
228 
229  // Check if we are at the end of the block
230  il1 = std::strtol( line1, &ctmp1, 10 );
231  il2 = std::strtol( line2, &ctmp2, 10 );
232  if( ( il1 == -1 ) && ( il2 == -1 ) )
233  {
234  s1 = ctmp1;
235  s2 = ctmp2;
236  if( ( s1.empty() ) && ( s2.empty() ) ) return MB_SUCCESS;
237  }
238 
239  // The first line describes attributes of the element other than connectivity.
240  const int element_id = strtol( line1 + 1, &ctmp1, 10 );
241  const int ideas_type = strtol( line1 + 11, &ctmp1, 10 );
242  const int phys_table = strtol( line1 + 21, &ctmp1, 10 );
243  const int mat_table = strtol( line1 + 31, &ctmp1, 10 );
244 
245  // Determine the element type.
246  EntityType mb_type;
247  if( TRI0 == ideas_type || TRI1 == ideas_type )
248  mb_type = MBTRI;
249  else if( QUAD0 == ideas_type || QUAD1 == ideas_type )
250  mb_type = MBQUAD;
251  else if( TET == ideas_type )
252  mb_type = MBTET;
253  else if( HEX == ideas_type )
254  mb_type = MBHEX;
255  else if( WEDGE == ideas_type )
256  mb_type = MBPRISM;
257  else
258  {
259  std::cout << "IDEAS element type not yet added to MOAB reader." << std::endl;
260  return MB_NOT_IMPLEMENTED;
261  }
262 
263  // Get the connectivity out of the 2nd line
264  std::stringstream ss( line2 );
265  const int n_conn = CN::VerticesPerEntity( mb_type );
267  EntityHandle vert;
268  for( int i = 0; i < n_conn; ++i )
269  {
270  ss >> vert;
271  conn[i] = vstart + vert - 1;
272  }
273 
274  // Make the element. According to the Gmsh 2.2.3 source code, the IDEAS
275  // canonical numbering is the same as MBCN.
276  rval = MBI->create_element( mb_type, conn, n_conn, handle );MB_CHK_SET_ERR( rval, "can't create elements of type " << mb_type );
277 
278  // If the phys set does not already exist, create it.
279  Range phys_sets;
280  EntityHandle phys_set;
281  const void* const phys_set_id_val[] = { &phys_table };
282  rval = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &phys_tag, phys_set_id_val, 1, phys_sets );MB_CHK_SET_ERR( rval, "can't get phys sets" );
283  if( phys_sets.empty() )
284  {
285  rval = MBI->create_meshset( MESHSET_SET, phys_set );MB_CHK_SET_ERR( rval, "can't create phys set" );
286  rval = MBI->tag_set_data( phys_tag, &phys_set, 1, &phys_table );MB_CHK_SET_ERR( rval, "can't set tag to phys set" );
287  }
288  else if( 1 == phys_sets.size() )
289  {
290  phys_set = phys_sets.front();
291  }
292  else
293  {
295  }
296  rval = MBI->add_entities( phys_set, &handle, 1 );MB_CHK_SET_ERR( rval, "can't add entities to phys set" );
297 
298  // If the material set does not already exist, create it.
299  Range mat_sets;
300  EntityHandle mat_set;
301  const void* const mat_set_id_val[] = { &mat_table };
302  rval = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &mat_tag, mat_set_id_val, 1, mat_sets );
303  if( MB_SUCCESS != rval ) return rval;
304  if( mat_sets.empty() )
305  {
306  rval = MBI->create_meshset( MESHSET_SET, mat_set );
307  if( MB_SUCCESS != rval ) return rval;
308  rval = MBI->tag_set_data( mat_tag, &mat_set, 1, &mat_table );
309  if( MB_SUCCESS != rval ) return rval;
310  }
311  else if( 1 == mat_sets.size() )
312  {
313  mat_set = mat_sets.front();
314  }
315  else
316  {
318  }
319  rval = MBI->add_entities( mat_set, &handle, 1 );
320  if( MB_SUCCESS != rval ) return rval;
321 
322  // Tag the element with its id
323  rval = MBI->tag_set_data( id_tag, &handle, 1, &element_id );MB_CHK_SET_ERR( rval, "Failed to assign IDs" );
324  if( file_id_tag )
325  {
326  rval = MBI->tag_set_data( *file_id_tag, &handle, 1, &element_id );MB_CHK_SET_ERR( rval, "Failed to assign file IDs" );
327  }
328  }
329 
330  return MB_SUCCESS;
331 }
332 
333 } // namespace moab