Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ReadRTT.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 Coroporation, 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 ReadRTT
18  * \brief ReadRTT based on ReadNASTRAN
19  *
20  * See:
21  *
22  * \author Andrew Davis
23  */
24 
25 #include "ReadRTT.hpp"
26 
27 #include <cassert>
28 #include <cmath>
29 #include <cstdlib>
30 #include <fstream>
31 #include <iostream>
32 #include <map>
33 #include <sstream>
34 #include <vector>
35 
36 #include "FileTokenizer.hpp"
37 #include "Internals.hpp" // for MB_START_ID
38 #include "MBTagConventions.hpp"
39 #include "moab/CN.hpp"
40 #include "moab/ErrorHandler.hpp"
41 #include "moab/FileOptions.hpp"
42 #include "moab/GeomTopoTool.hpp"
43 #include "moab/Interface.hpp"
44 #include "moab/Range.hpp"
45 #include "moab/ReadUtilIface.hpp"
46 
47 namespace moab
48 {
49 
51 {
52  return new ReadRTT( iface );
53 }
54 
55 // constructor
57  : MBI( impl ), geom_tag( 0 ), id_tag( 0 ), name_tag( 0 ), category_tag( 0 ), faceting_tol_tag( 0 )
58 {
59  assert( NULL != impl );
60  myGeomTool = new GeomTopoTool( impl );
62  assert( NULL != readMeshIface );
63 
64  // this section copied from ReadCGM initalisation
65  int negone = -1;
66  double zero = 0.;
67  ErrorCode rval;
69  &negone );
70  MB_CHK_ERR_CONT( rval );
71  id_tag = MBI->globalId_tag();
73  MB_CHK_ERR_CONT( rval );
76  MB_CHK_ERR_CONT( rval );
77  rval =
79  MB_CHK_ERR_CONT( rval );
80 }
81 
82 // destructor
84 {
85  if( readMeshIface )
86  {
88  readMeshIface = 0;
89  }
90 
91  delete myGeomTool;
92 }
93 
94 ErrorCode ReadRTT::read_tag_values( const char* /*file_name*/,
95  const char* /*tag_name*/,
96  const FileOptions& /*opts*/,
97  std::vector< int >& /*tag_values_out*/,
98  const SubsetList* /*subset_list*/ )
99 {
100  return MB_NOT_IMPLEMENTED;
101 }
102 
103 // load the file as called by the Interface function
104 ErrorCode ReadRTT::load_file( const char* filename,
105  const EntityHandle*,
106  const FileOptions&,
107  const ReaderIface::SubsetList* subset_list,
108  const Tag* /*file_id_tag*/ )
109 {
110  ErrorCode rval;
111 
112  // at this time there is no support for reading a subset of the file
113  if( subset_list )
114  {
115  std::cout << "Subset reading not supported for RTT meshes" << std::endl;
117  }
118 
119  // test to see if file exists
120  FILE* file = NULL;
121  file = fopen( filename, "r" );
122  if( file == NULL ) return MB_FILE_DOES_NOT_EXIST;
123  // otherwise close the file
124  fclose( file );
125 
126  // read the header
127  rval = ReadRTT::read_header( filename );
128  if( rval != MB_SUCCESS ) return rval;
129 
130  // read the side_flag data
131  rval = ReadRTT::read_side_flags( filename );
132  if( rval != MB_SUCCESS ) return rval;
133 
134  // read the cell data
135  rtt_flags cell_flags;
136  rval = ReadRTT::read_cell_flags( filename );
137  if( rval != MB_SUCCESS ) return rval;
138 
139  // read the node data
140  std::vector< node > node_data;
141  rval = ReadRTT::read_nodes( filename, node_data );
142  if( rval != MB_SUCCESS ) return rval;
143 
144  // read the facet data
145  std::vector< facet > facet_data;
146  rval = ReadRTT::read_facets( filename, facet_data );
147  if( rval != MB_SUCCESS ) return rval;
148 
149  // read the tetrahedra data
150  std::vector< tet > tet_data;
151  rval = ReadRTT::read_tets( filename, tet_data );
152  if( rval != MB_SUCCESS ) return rval;
153 
154  // make the map of surface number in the rttmesh to the surface meshset
155  std::map< int, EntityHandle > surface_map; // corrsespondance of surface number to entity handle
156  std::map< int, EntityHandle > volume_map; // corrsespondance of volume number to entity handle
157  rval = ReadRTT::generate_topology( side_data, cell_data, tet_data, surface_map, volume_map );
158  if( rval != MB_SUCCESS ) return rval;
159 
160  // generate the rest of the database, triangles to surface meshsets etc
161  rval = ReadRTT::build_moab( node_data, facet_data, tet_data, surface_map, volume_map );
162  if( rval != MB_SUCCESS ) return rval;
163 
164  return MB_SUCCESS;
165 }
166 
167 /*
168  * builds the topology of the problem
169  */
170 ErrorCode ReadRTT::generate_topology( std::vector< side > side_data,
171  std::vector< cell > cell_data,
172  std::vector< tet > tet_data,
173  std::map< int, EntityHandle >& surface_map,
174  std::map< int, EntityHandle >& volume_map )
175 {
176 
177  ErrorCode rval;
178  std::vector< EntityHandle > entmap[4];
179  int num_ents[4]; // number of entities in each dimension
180 
181  const char geom_categories[][CATEGORY_TAG_SIZE] = { "Vertex\0", "Curve\0", "Surface\0", "Volume\0", "Group\0" };
182 
183  std::vector< int > surface_numbers; // the surface numbers in the problem
184 
185  // corresponds to number of cad like surfaces and cad like volumes
186  num_ents[2] = side_data.size();
187  num_ents[3] = cell_data.size();
188 
189  // loop over surfaces & volumes
190  for( int dim = 2; dim <= 3; dim++ )
191  {
192  for( int i = 0; i != num_ents[dim]; i++ )
193  {
194  EntityHandle handle;
195  // create a meshset for each entity surface/volume
196  rval = MBI->create_meshset( dim == 1 ? MESHSET_ORDERED : MESHSET_SET, handle );
197  // if failure
198  if( rval != MB_SUCCESS ) return rval;
199 
200  // collect the entity handles into an
201  entmap[dim].push_back( handle );
202 
203  // set the dimension tag
204  rval = MBI->tag_set_data( geom_tag, &handle, 1, &dim );
205  // if fail
206  if( MB_SUCCESS != rval ) return rval;
207  // if we are a surface
208  if( dim == 2 )
209  {
210  // tag the id onto the surface meshset
211  rval = MBI->tag_set_data( id_tag, &handle, 1, &side_data[i].id );
212  // inesert entity into the map
213  surface_map[side_data[i].id] = handle;
214  }
215  else
216  {
217  // otherwise we set the volume tag data, loop is only 2 & 3 dim
218  rval = MBI->tag_set_data( id_tag, &handle, 1, &cell_data[i].id );
219  volume_map[cell_data[i].id] = handle;
220  }
221  // if fail
222  if( MB_SUCCESS != rval ) return rval;
223  // set the category tag
224  rval = MBI->tag_set_data( category_tag, &handle, 1, &geom_categories[dim] );
225  if( MB_SUCCESS != rval ) return rval;
226  }
227  }
228 
229  // generate parent child links
230  // best to loop over the surfaces and assign them to volumes, we can then
231  // assign facets to
232  // to each surface
233  generate_parent_child_links( num_ents, entmap, side_data, cell_data );
234 
235  // set the surface senses
236  set_surface_senses( num_ents, entmap, side_data, cell_data );
237 
238  // set the group data
239  rval = setup_group_data( entmap, tet_data, volume_map );
240 
241  return MB_SUCCESS;
242 }
243 
244 /*
245  * builds the moab representation of the mesh
246  */
247 ErrorCode ReadRTT::build_moab( std::vector< node > node_data,
248  std::vector< facet > facet_data,
249  std::vector< tet > tet_data,
250  const std::map< int, EntityHandle > surface_map,
251  const std::map< int, EntityHandle > volume_map )
252 {
253  ErrorCode rval;
254  EntityHandle file_set;
255 
256  rval = MBI->create_meshset( MESHSET_SET, file_set );
257  if( MB_SUCCESS != rval ) return rval;
258 
259  // adding vertex set to the file set
260  Range mb_coords;
261  for( const auto& n : node_data )
262  {
263  double coords[3] = { n.x, n.y, n.z };
264  EntityHandle v;
265  rval = MBI->create_vertex( coords, v );MB_CHK_ERR( rval );
266  mb_coords.insert( v );
267  }
268  rval = MBI->add_entities( file_set, mb_coords );MB_CHK_ERR( rval );
269 
270  // add facets to the file set
271  create_facets( facet_data, surface_map, mb_coords, file_set );
272 
273  // material number tag
274  Tag mat_num_tag;
275  rval = MBI->tag_get_handle( "MATERIAL_NUMBER", 1, MB_TYPE_INTEGER, mat_num_tag, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
276 
277  // adding material groups
278  std::string mat_flag = get_material_ref_flag();
279  std::string vol_flag = get_volume_ref_flag();
280 
281  // add tets to the file set
282  Range mb_tets;
283  for( const auto& t : tet_data )
284  {
285  EntityHandle tet_nodes[4] = { mb_coords[t.connectivity[0] - 1], mb_coords[t.connectivity[1] - 1],
286  mb_coords[t.connectivity[2] - 1], mb_coords[t.connectivity[3] - 1] };
287 
288  EntityHandle tet_h;
289  rval = MBI->create_element( MBTET, tet_nodes, 4, tet_h );MB_CHK_ERR( rval );
290 
291  int mat_no = t.flag_values[cell_flag_idx[mat_flag]];
292  rval = MBI->tag_set_data( mat_num_tag, &tet_h, 1, &mat_no );MB_CHK_ERR( rval );
293 
294  mb_tets.insert( tet_h );
295  }
296  rval = MBI->add_entities( file_set, mb_tets );MB_CHK_ERR( rval );
297 
298  return MB_SUCCESS;
299 }
300 
301 // Function to create a material group
302 ErrorCode ReadRTT::create_material_group( const std::string& material_name, int material_id, EntityHandle& handle )
303 {
304  ErrorCode rval = MBI->create_meshset( MESHSET_SET, handle );
305  if( rval != MB_SUCCESS ) return rval;
306 
307  // NAME
308  char name_val[NAME_TAG_SIZE] = { 0 };
309  std::strncpy( name_val, material_name.c_str(), NAME_TAG_SIZE - 1 );
310  rval = MBI->tag_set_data( name_tag, &handle, 1, name_val );MB_CHK_ERR( rval );
311 
312  // GLOBAL_ID
313  rval = MBI->tag_set_data( id_tag, &handle, 1, &material_id );MB_CHK_ERR( rval );
314 
315  // CATEGORY
316  char cat[CATEGORY_TAG_SIZE] = { 0 };
317  std::strncpy( cat, "Group", CATEGORY_TAG_SIZE - 1 );
318  rval = MBI->tag_set_data( category_tag, &handle, 1, cat );MB_CHK_ERR( rval );
319 
320  // GEOM_DIMENSION = 4
321  int dim4 = 4;
322  rval = MBI->tag_set_data( geom_tag, &handle, 1, &dim4 );MB_CHK_ERR( rval );
323 
324  return MB_SUCCESS;
325 }
326 
327 void ReadRTT::create_facets( const std::vector< facet >& facet_data,
328  const std::map< int, EntityHandle >& surface_map,
329  Range& mb_coords,
330  EntityHandle file_set )
331 {
332  ErrorCode rval;
333  Tag side_id_tag, surface_number_tag;
334  // Obtain or create tags for side IDs and surface numbers
335  rval = MBI->tag_get_handle( "SIDEID_TAG", 1, MB_TYPE_INTEGER, side_id_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
336  rval =
337  MBI->tag_get_handle( "SURFACE_NUMBER", 1, MB_TYPE_INTEGER, surface_number_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
338 
339  EntityHandle triangle;
340  Range mb_tris; // For storing triangles
341 
342  for( const auto& tmp : facet_data )
343  {
344  EntityHandle tri_nodes[3] = { mb_coords[tmp.connectivity[0] - 1], mb_coords[tmp.connectivity[1] - 1],
345  mb_coords[tmp.connectivity[2] - 1] };
346  rval = MBI->create_element( MBTRI, tri_nodes, 3, triangle );
347  // tag in side id on the triangle
348  rval = MBI->tag_set_data( side_id_tag, &triangle, 1, &tmp.side_id );
349  // tag the surface number on the triangle
350  rval = MBI->tag_set_data( surface_number_tag, &triangle, 1, &tmp.surface_number );
351  // insert vertices and triangles into the appropriate surface meshset
352  EntityHandle meshset_handle = surface_map.at( tmp.surface_number );
353  // also set surface tag
354  rval = MBI->tag_set_data( side_id_tag, &meshset_handle, 1, &tmp.side_id );
355  rval = MBI->tag_set_data( surface_number_tag, &meshset_handle, 1, &tmp.surface_number );
356  // add vertices to the mesh
357  rval = MBI->add_entities( meshset_handle, tri_nodes, 3 );
358  // add triangles to the meshset
359  rval = MBI->add_entities( meshset_handle, &triangle, 1 );
360  // insert triangles into mb_tris
361  mb_tris.insert( triangle );
362  }
363  rval = MBI->add_entities( file_set, mb_tris );
364 }
365 
367 {
368  moab::ErrorCode rval = MB_FAILURE;
369 
370  // Create CONTIGUITY tag and set its value
371  Tag contiguity_tag;
372  const char* contiguity_value = header_data.contiguity.c_str();
373  rval = MBI->tag_get_handle( "CONTIGUITY", strlen( contiguity_value ) + 1, MB_TYPE_OPAQUE, contiguity_tag,
375  if( rval != MB_SUCCESS ) return rval;
376  rval = MBI->tag_set_data( contiguity_tag, &file_set, 1, contiguity_value );
377 
378  return rval;
379 }
380 
381 /*
382  * read the header data from the filename pointed to
383  */
384 ErrorCode ReadRTT::read_header( const char* filename )
385 {
386  std::ifstream input_file( filename ); // filename for rtt file
387  // file ok?
388  if( !input_file.good() )
389  {
390  std::cout << "Problems reading file = " << filename << std::endl;
391  return MB_FAILURE;
392  }
393 
394  // if it works
395  std::string line;
396  moab::ErrorCode rval = MB_FAILURE;
397  if( input_file.is_open() )
398  {
399  while( std::getline( input_file, line ) )
400  {
401  if( line.compare( "header" ) == 0 )
402  {
403  rval = get_header_data( input_file );
404  }
405  else if( line.compare( "dims" ) == 0 )
406  {
407  rval = parse_dims( input_file );
408  }
409  else if( line.compare( "cell_defs" ) == 0 )
410  {
412  }
413  }
414  input_file.close();
415  }
416  return rval;
417 }
418 
419 // read all flags from a section
420 ErrorCode ReadRTT::read_all_flags( const char* filename,
421  std::vector< int > n_flags,
422  std::string flag_id,
423  rtt_flags& flags,
424  std::map< std::string, int >& flag_idx )
425 {
426  std::string start_flag = flag_id + "_flags";
427  std::string end_flag = "end_" + start_flag + "\0";
428  std::string line; // the current line being read
429  std::ifstream input_file( filename ); // filestream for rttfile
430  std::vector< std::string > flag_order; // order of the flags
431  // file ok?
432  if( !input_file.good() )
433  {
434  std::cout << "Problems reading file = " << filename << std::endl;
435  return MB_FAILURE;
436  }
437  // if it works
438  if( input_file.is_open() )
439  {
440  while( std::getline( input_file, line ) )
441  {
442  if( line.compare( start_flag ) == 0 )
443  {
444  while( std::getline( input_file, line ) )
445  {
446  // Read all the side block until we find the end
447  if( line.compare( end_flag ) == 0 ) break;
448  std::vector< std::string > token = ReadRTT::split_string( line, ' ' );
449  if( token.size() != 2 )
450  {
451  std::cout << "Error reading side flags" << std::endl;
452  return MB_FAILURE;
453  }
454  int flag_key = std::stoi( token[0] ) - 1;
455  std::string key = token[1];
456  flag_order.push_back( key );
457  for( int i = 0; i < n_flags[flag_key]; i++ )
458  {
459  std::getline( input_file, line );
460  flags[key].push_back( line );
461  }
462  }
463  }
464  }
465  input_file.close();
466  }
467  for( size_t i = 0; i < flag_order.size(); i++ )
468  {
469  std::string key = flag_order[i];
470  // fill the index
471  flag_idx[key] = i;
472  }
473  return MB_SUCCESS;
474 }
475 
476 /*
477  * reads the side data from the filename pointed to
478  */
479 ErrorCode ReadRTT::read_side_flags( const char* filename )
480 {
481  rtt_flags side_flags;
482 
483  ErrorCode rval = MB_FAILURE;
484  // read all the side data
485  rval = read_all_flags( filename, dim_data.nside_flags, "side", side_flags, side_flag_idx );
486 
487  //process sides
488  rval = ReadRTT::side_process_faces( side_flags, side_data );
489  if( rval != MB_SUCCESS ) return rval;
490 
491  return rval;
492 }
493 
494 /*
495  * process the FACES flag from the side_flags section
496  */
497 ErrorCode ReadRTT::side_process_faces( rtt_flags side_flags, std::vector< side >& side_data )
498 {
499  if( side_flags.find( "FACES" ) != side_flags.end() )
500  {
501  for( size_t i = 0; i < side_flags["FACES"].size(); i++ )
502  {
503  side data = ReadRTT::get_side_data( side_flags["FACES"][i] );
504  side_data.push_back( data );
505  }
506  }
507  if( side_data.size() == 0 ) return MB_FAILURE;
508  return MB_SUCCESS;
509 }
510 
512 {
513  std::string material_ref_flag = ""; // set defaul to REGIONS
514  if( cell_flag_datas.find( "MATERIAL" ) != cell_flag_datas.end() )
515  {
516  material_ref_flag = "MATERIAL";
517  }
518  return material_ref_flag;
519 }
520 
522 {
523  std::string part_ref_flag = "REGIONS"; // set defaul to REGIONS
524  if( cell_flag_datas.find( "MCNP_PSEUDO-CELLS" ) != cell_flag_datas.end() )
525  {
526  part_ref_flag = "MCNP_PSEUDO-CELLS";
527  }
528  else if( cell_flag_datas.find( "ABAQUS_PARTS" ) != cell_flag_datas.end() )
529  {
530  part_ref_flag = "ABAQUS_PARTS";
531  }
532  return part_ref_flag;
533 }
534 
535 /*
536  * reads the cell data from the filename pointed to
537  */
538 ErrorCode ReadRTT::read_cell_flags( const char* filename )
539 {
540  rtt_flags cell_flags;
541  ErrorCode rval = MB_FAILURE;
542  // read all the cell data
543  rval = read_all_flags( filename, dim_data.ncell_flags, "cell", cell_flags, cell_flag_idx );
544 
545  for( auto it = cell_flags.begin(); it != cell_flags.end(); ++it )
546  {
547  std::string key = it->first;
548  // fill the index
549  rval = ReadRTT::cell_process_flag( cell_flags, key );
550  if( rval != MB_SUCCESS ) return rval;
551  }
552 
553  std::string part_flag_name = get_volume_ref_flag();
554  cell_data = cell_flag_datas[part_flag_name];
555  cell_data_idx = cell_flag_indexes[part_flag_name];
556  return rval;
557 }
558 
559 /*
560  * process the standard flag from the cell_flags section
561  */
562 ErrorCode ReadRTT::cell_process_flag( rtt_flags cell_flags, std::string key )
563 {
564  std::vector< cell > cell_data;
565  // check if the key is in the cell_flags
566  if( cell_flags.find( key ) != cell_flags.end() )
567  {
568  for( size_t i = 0; i < cell_flags[key].size(); i++ )
569  {
570  cell data = ReadRTT::get_cell_data( cell_flags[key][i] );
571  cell_data.push_back( data );
572  }
573  if( cell_data.size() == 0 ) return MB_FAILURE;
574  }
575 
576  // fill the corresponding index
577  std::map< int, int > cell_data_idx;
578  for( size_t i = 0; i < cell_data.size(); ++i )
579  {
580  cell tmp = cell_data[i];
581  cell_data_idx[tmp.id] = i;
582  }
583  if( cell_data.size() > 0 ) cell_flag_datas[key] = cell_data;
585 
586  return MB_SUCCESS;
587 }
588 
589 /*
590  * Reads the node data fromt the filename pointed to
591  */
592 ErrorCode ReadRTT::read_nodes( const char* filename, std::vector< node >& node_data )
593 {
594  std::string line; // the current line being read
595  std::ifstream input_file( filename ); // filestream for rttfile
596  // file ok?
597  if( !input_file.good() )
598  {
599  std::cout << "Problems reading file = " << filename << std::endl;
600  return MB_FAILURE;
601  }
602 
603  // if it works
604  if( input_file.is_open() )
605  {
606  while( std::getline( input_file, line ) )
607  {
608  if( line.compare( "nodes\0" ) == 0 )
609  {
610  // read lines until find end nodes
611  while( std::getline( input_file, line ) )
612  {
613  if( line.compare( "end_nodes\0" ) == 0 ) break;
614  node data = ReadRTT::get_node_data( line );
615  node_data.push_back( data );
616  }
617  }
618  }
619  input_file.close();
620  }
621  if( node_data.size() == 0 ) return MB_FAILURE;
622  return MB_SUCCESS;
623 }
624 
625 /*
626  * Reads the facet data fromt the filename pointed to
627  */
628 ErrorCode ReadRTT::read_facets( const char* filename, std::vector< facet >& facet_data )
629 {
630  std::string line; // the current line being read
631  std::ifstream input_file( filename ); // filestream for rttfile
632  // file ok?
633  if( !input_file.good() )
634  {
635  std::cout << "Problems reading file = " << filename << std::endl;
636  return MB_FAILURE;
637  }
638 
639  // if it works
640  if( input_file.is_open() )
641  {
642  while( std::getline( input_file, line ) )
643  {
644  if( line.compare( "sides\0" ) == 0 )
645  {
646  // read lines until find end nodes
647  while( std::getline( input_file, line ) )
648  {
649  if( line.compare( "end_sides\0" ) == 0 ) break;
650  facet data = ReadRTT::get_facet_data( line );
651  facet_data.push_back( data );
652  }
653  }
654  }
655  input_file.close();
656  }
657  if( facet_data.size() == 0 ) return MB_FAILURE;
658  return MB_SUCCESS;
659 }
660 
661 /*
662  * Reads the facet data fromt the filename pointed to
663  */
664 ErrorCode ReadRTT::read_tets( const char* filename, std::vector< tet >& tet_data )
665 {
666  std::string line; // the current line being read
667  std::ifstream input_file( filename ); // filestream for rttfile
668  // file ok?
669  if( !input_file.good() )
670  {
671  std::cout << "Problems reading file = " << filename << std::endl;
672  return MB_FAILURE;
673  }
674  // if it works
675  if( input_file.is_open() )
676  {
677  while( std::getline( input_file, line ) )
678  {
679  if( line.compare( "cells\0" ) == 0 )
680  {
681  // read lines until find end nodes
682  while( std::getline( input_file, line ) )
683  {
684  if( line.compare( "end_cells\0" ) == 0 ) break;
685  tet data = ReadRTT::get_tet_data( line );
686  tet_data.push_back( data );
687  }
688  }
689  }
690  input_file.close();
691  }
692  if( tet_data.size() == 0 ) return MB_FAILURE;
693  return MB_SUCCESS;
694 }
695 
696 /*
697  * given the open file handle read until we find
698  */
700 {
701  std::string line;
702  while( std::getline( input_file, line ) )
703  {
704 
705  // tokenize the line
706  std::istringstream iss( line );
707  std::vector< std::string > split_string;
708  do
709  {
710  std::string sub_string;
711  iss >> sub_string;
712  split_string.push_back( sub_string );
713  } while( iss );
714 
715  // if we find version
716  if( line.find( "version" ) != std::string::npos )
717  {
718  if( split_string[1].find( "v" ) != std::string::npos &&
719  split_string[0].find( "version" ) != std::string::npos )
720  {
722  }
723  }
724  else if( line.find( "title" ) != std::string::npos )
725  {
727  }
728  else if( line.find( "date" ) != std::string::npos )
729  {
731  }
732  else if( line.find( "contiguity" ) != std::string::npos )
733  {
735  }
736  else if( line.find( "end_header" ) != std::string::npos )
737  {
738  return MB_SUCCESS;
739  }
740  }
741  // otherwise we never found the end_header keyword
742  return MB_FAILURE;
743 }
744 
746 {
747  if( !input_file.good() || !input_file.is_open() )
748  {
749  std::cout << "Problems reading file" << std::endl;
750  return MB_FAILURE;
751  }
752 
753  std::string line;
754  std::vector< std::string > tokens;
755  while( std::getline( input_file, line ) )
756  {
757  if( line == "" ) continue;
758  if( line.find( "end_dims" ) != std::string::npos ) break;
759 
760  tokens = ReadRTT::split_string( line, ' ' );
761  if( tokens[0] == "coor_units" )
762  {
763  dim_data.coor_units = tokens[1];
764  }
765  else if( tokens[0] == "prob_time_units" )
766  {
767  dim_data.prob_time_units = tokens[1];
768  }
769  else if( tokens[0] == "ncell_defs" )
770  {
771  dim_data.ncell_defs = std::atoi( tokens[1].c_str() );
772  }
773  else if( tokens[0] == "nnodes_max" )
774  {
775  dim_data.nnodes_max = std::atoi( tokens[1].c_str() );
776  }
777  else if( tokens[0] == "nsides_max" )
778  {
779  dim_data.nsides_max = std::atoi( tokens[1].c_str() );
780  }
781  else if( tokens[0] == "nnodes_sides_max" )
782  {
783  dim_data.nnodes_sides_max = std::atoi( tokens[1].c_str() );
784  }
785  else if( tokens[0] == "ndim" )
786  {
787  dim_data.ndim = std::atoi( tokens[1].c_str() );
788  }
789  else if( tokens[0] == "n_dim_topo" )
790  {
791  dim_data.n_dim_topo = std::atoi( tokens[1].c_str() );
792  }
793  else if( tokens[0] == "nnodes" )
794  {
795  dim_data.nnodes = std::atoi( tokens[1].c_str() );
796  }
797  else if( tokens[0] == "nnode_flag_types" )
798  {
799  dim_data.nnode_flag_types = std::atoi( tokens[1].c_str() );
800  }
801  else if( tokens[0] == "nnode_flags" )
802  {
803  for( size_t i = 1; i < tokens.size(); i++ )
804  {
805  dim_data.nnode_flags.push_back( std::atoi( tokens[i].c_str() ) );
806  }
807  }
808  else if( tokens[0] == "nnode_data" )
809  {
810  dim_data.nnode_data = std::atoi( tokens[1].c_str() );
811  }
812  else if( tokens[0] == "nsides" )
813  {
814  dim_data.nsides = std::atoi( tokens[1].c_str() );
815  }
816  else if( tokens[0] == "nside_types" )
817  {
818  dim_data.nside_types = std::atoi( tokens[1].c_str() );
819  }
820  else if( tokens[0] == "side_types" )
821  {
822  dim_data.side_types = std::atoi( tokens[1].c_str() );
823  }
824  else if( tokens[0] == "nside_flag_types" )
825  {
826  dim_data.nside_flag_types = std::atoi( tokens[1].c_str() );
827  }
828  else if( tokens[0] == "nside_flags" )
829  {
830  for( size_t i = 1; i < tokens.size(); i++ )
831  {
832  dim_data.nside_flags.push_back( std::atoi( tokens[i].c_str() ) );
833  }
834  }
835  else if( tokens[0] == "nside_data" )
836  {
837  dim_data.nside_data = std::atoi( tokens[1].c_str() );
838  }
839  else if( tokens[0] == "ncells" )
840  {
841  dim_data.ncells = std::atoi( tokens[1].c_str() );
842  }
843  else if( tokens[0] == "ncell_types" )
844  {
845  dim_data.ncell_types = std::atoi( tokens[1].c_str() );
846  }
847  else if( tokens[0] == "cell_types" )
848  {
849  dim_data.cell_types = std::atoi( tokens[1].c_str() );
850  }
851  else if( tokens[0] == "ncell_flag_types" )
852  {
853  dim_data.ncell_flag_types = std::atoi( tokens[1].c_str() );
854  }
855  else if( tokens[0] == "ncell_flags" )
856  {
857  for( size_t i = 1; i < tokens.size(); i++ )
858  {
859  dim_data.ncell_flags.push_back( std::atoi( tokens[i].c_str() ) );
860  }
861  }
862  else if( tokens[0] == "ncell_data" )
863  {
864  dim_data.ncell_data = std::atoi( tokens[1].c_str() );
865  }
866  }
867  // Check that the data is valid and has the expected number of entries
868  dim_data.validate();
869 
870  return MB_SUCCESS;
871 }
872 
873 /*
874  * parsing the cell deifinitions
875  */
877 {
878  if( !input_file.good() || !input_file.is_open() )
879  {
880  std::cout << "Problems reading file" << std::endl;
881  return MB_FAILURE;
882  }
883 
884  std::string line;
885  std::vector< std::string > tokens;
886  while( std::getline( input_file, line ) )
887  {
888  if( line == "" ) continue;
889  if( line.find( "end_cell_defs" ) != std::string::npos ) break;
890  cell_def new_cell_def;
891  // Tokenize the line
892  tokens = ReadRTT::split_string( line, ' ' );
893  new_cell_def.id = std::atoi( tokens[0].c_str() );
894  new_cell_def.name = tokens[1];
895  // Getting the number of nodes and sides
896  std::getline( input_file, line );
897  // Tokenize the line
898  tokens = ReadRTT::split_string( line, ' ' );
899  new_cell_def.nnodes = std::atoi( tokens[0].c_str() );
900  new_cell_def.nsides = std::atoi( tokens[1].c_str() );
901  // Side type index
902  std::getline( input_file, line );
903  // Tokenize the line
904  tokens = ReadRTT::split_string( line, ' ' );
905  for( int i = 0; i < new_cell_def.nsides; i++ )
906  {
907  int side_type = std::atoi( tokens[i].c_str() );
908  // Ensure side types exists
909  if( cell_def_data.find( side_type ) == cell_def_data.end() )
910  {
911  std::cout << "Error: side type " << side_type << " not found in cell definitions." << std::endl;
912  return MB_FAILURE;
913  }
914  new_cell_def.side_type.push_back( side_type );
915  }
916  // Read the nodes per side
917  for( int i = 0; i < new_cell_def.nsides; i++ )
918  {
919  std::vector< int > side_nodes;
920  std::getline( input_file, line );
921  // Tokenize the line
922  tokens = ReadRTT::split_string( line, ' ' );
923  for( int j = 0; j < new_cell_def.nnodes; j++ )
924  {
925  side_nodes.push_back( std::atoi( tokens[i].c_str() ) );
926  }
927  new_cell_def.sides_nodes.push_back( side_nodes );
928  }
929 
930  cell_def_data.insert( std::make_pair( new_cell_def.id, new_cell_def ) );
931  }
932 
933  return MB_SUCCESS;
934 }
935 
936 /*
937  * given the string sidedata, get the id number, senses and names of the sides
938  */
939 ReadRTT::side ReadRTT::get_side_data( std::string sidedata )
940 {
941  side new_side;
942  std::vector< std::string > tokens;
943  tokens = ReadRTT::split_string( sidedata, ' ' );
944 
945  // set the side id
946  if( tokens.size() != 2 )
947  {
948  MB_SET_ERR_RET_VAL( "Error, too many tokens found from side_data", new_side );
949  }
950  // create the new side
951  new_side.id = std::atoi( tokens[0].c_str() );
952 
953  std::vector< std::string > cell_names = ReadRTT::split_string( tokens[1], '/' );
954  // get the boundary
955  boundary new_bnd = ReadRTT::split_name( cell_names[0] );
956  // set the surface sense and name
957  new_side.senses[0] = new_bnd.sense;
958  new_side.names[0] = new_bnd.name;
959  //
960  if( cell_names.size() > 1 )
961  {
962  boundary bnd = ReadRTT::split_name( cell_names[1] );
963  new_side.senses[1] = bnd.sense;
964  new_side.names[1] = bnd.name;
965  }
966  else
967  {
968  new_side.senses[1] = 0;
969  new_side.names[1] = "\0";
970  }
971 
972  return new_side;
973 }
974 
975 /*
976  * given the string celldata, get the id number and name of each cell
977  */
978 ReadRTT::cell ReadRTT::get_cell_data( std::string celldata )
979 {
980  cell new_cell;
981  std::vector< std::string > tokens;
982  tokens = ReadRTT::split_string( celldata, ' ' );
983 
984  // set the side id
985  if( tokens.size() != 2 )
986  {
987  MB_SET_ERR_RET_VAL( "Error, too many tokens found from cell_data", new_cell );
988  }
989  // create the new side
990  new_cell.id = std::atoi( tokens[0].c_str() );
991  new_cell.name = tokens[1];
992 
993  return new_cell;
994 }
995 
996 /*
997  * given the string nodedata, get the id number and coordinates of the node
998  */
999 ReadRTT::node ReadRTT::get_node_data( std::string nodedata )
1000 {
1001  node new_node;
1002  std::vector< std::string > tokens;
1003  tokens = ReadRTT::split_string( nodedata, ' ' );
1004 
1005  // set the side id
1006  if( tokens.size() != 5 )
1007  {
1008  MB_SET_ERR_RET_VAL( "Error, too many tokens found from get_node_data", new_node );
1009  }
1010  new_node.id = std::atoi( tokens[0].c_str() );
1011  new_node.x = std::atof( tokens[1].c_str() );
1012  new_node.y = std::atof( tokens[2].c_str() );
1013  new_node.z = std::atof( tokens[3].c_str() );
1014  return new_node;
1015 }
1016 
1017 /*
1018  * given the string nodedata, get the id number, connectivity and sense data
1019  */
1020 ReadRTT::facet ReadRTT::get_facet_data( std::string facetdata )
1021 {
1022  facet new_facet;
1023  std::vector< std::string > tokens;
1024  tokens = ReadRTT::split_string( facetdata, ' ' );
1025 
1026  // ensure we have the correct number of tokens
1027  int base_token_size = 0;
1028  int idx_offset = 0;
1029  // branch on the rtt version number
1030  if( header_data.version == "v1.0.0" )
1031  {
1032  base_token_size = 4;
1033  }
1034  else if( header_data.version == "v1.0.1" )
1035  {
1036  base_token_size = 5;
1037  idx_offset = 1;
1038  }
1039  else
1040  {
1041  MB_SET_ERR_RET_VAL( "Error, version number not understood", new_facet );
1042  }
1043 
1044  if( (int)tokens.size() != base_token_size + dim_data.nside_flag_types )
1045  {
1046  std::cout << facetdata << std::endl;
1047  std::cout << header_data.version << " " << (int)tokens.size() << " " << base_token_size << " "
1048  << dim_data.nside_flag_types << std::endl;
1049  MB_SET_ERR_RET_VAL( "Error, too many tokens found from get_facet_data", new_facet );
1050  exit( 1 );
1051  }
1052 
1053  // set the side id
1054  new_facet.id = std::atoi( tokens[0].c_str() );
1055  new_facet.connectivity[0] = std::atoi( tokens[idx_offset + 1].c_str() );
1056  new_facet.connectivity[1] = std::atoi( tokens[idx_offset + 2].c_str() );
1057  new_facet.connectivity[2] = std::atoi( tokens[idx_offset + 3].c_str() );
1058  new_facet.side_id = std::atoi( tokens[idx_offset + 4].c_str() );
1059  new_facet.surface_number = std::atoi( tokens[idx_offset + 5].c_str() );
1060 
1061  return new_facet;
1062 }
1063 
1064 /*
1065  * given the string tetdata, get the id number, connectivity and mat num of the
1066  * tet
1067  */
1068 ReadRTT::tet ReadRTT::get_tet_data( std::string tetdata )
1069 {
1070  tet new_tet;
1071  std::vector< std::string > tokens;
1072  tokens = ReadRTT::split_string( tetdata, ' ' );
1073 
1074  // ensure we have the correct number of tokens
1075  int base_token_size = 0;
1076  int n_nodes = 0;
1077  // branch on the rtt version number
1078  if( header_data.version == "v1.0.0" )
1079  {
1080  base_token_size = 1;
1081  // for old format we assume 4 nodes (Tet)
1082  n_nodes = 4;
1083  }
1084  else if( header_data.version == "v1.0.1" )
1085  {
1086  base_token_size = 2;
1087 
1088  // get the cell type
1089  new_tet.type_id = std::atoi( tokens[1].c_str() );
1090  // check that the cell type exists
1091  if( cell_def_data.find( new_tet.type_id ) == cell_def_data.end() )
1092  {
1093  std::cout << "Error: cell type " << new_tet.type_id << " not found in cell definitions." << std::endl;
1094  MB_SET_ERR_RET_VAL( "Error, cell type not found in cell definitions", new_tet );
1095  }
1096  n_nodes = cell_def_data[new_tet.type_id].nnodes;
1097  int n_sides = cell_def_data[new_tet.type_id].nsides;
1098  if( n_nodes != 4 || n_sides != 4 )
1099  {
1100  std::cout << "Error: cell type " << new_tet.type_id << " does not have 4 nodes and 4 sides." << std::endl;
1101  MB_SET_ERR_RET_VAL( "Error, cell is not tetrahedricon", new_tet );
1102  }
1103  }
1104  else
1105  {
1106  MB_SET_ERR_RET_VAL( "Error, version number not understood", new_tet );
1107  }
1108  // set the side id
1109  new_tet.id = std::atoi( tokens[0].c_str() );
1110 
1111  if( (int)tokens.size() != base_token_size + n_nodes + dim_data.ncell_flag_types )
1112  {
1113  MB_SET_ERR_RET_VAL( "Error, unexpected number of tokens found from get_tet_data", new_tet );
1114  }
1115  for( int i = 0; i < n_nodes; i++ )
1116  {
1117  new_tet.connectivity[i] = std::atoi( tokens[i + base_token_size].c_str() );
1118  }
1119  for( int i = 0; i < dim_data.ncell_flag_types; i++ )
1120  {
1121  new_tet.flag_values.push_back( std::atoi( tokens[i + base_token_size + n_nodes].c_str() ) );
1122  }
1123 
1124  return new_tet;
1125 }
1126 
1127 /*
1128 * given the cell data, get the maximum name size
1129 */
1130 int ReadRTT::get_max_name_size( std::vector< cell > cell_data )
1131 {
1132  int max_size = 0;
1133  for( size_t i = 0; i < cell_data.size(); i++ )
1134  {
1135  if( (int)cell_data[i].name.length() > max_size ) max_size = cell_data[i].name.length();
1136  }
1137  return max_size;
1138 }
1139 
1140 /*
1141  * splits string into sense and name, to later facilitate the building
1142  * of sense data, strips off the tailing @ if it exists
1143  */
1144 ReadRTT::boundary ReadRTT::split_name( std::string atilla_cellname )
1145 {
1146  boundary new_boundary;
1147  // default initialisation
1148  new_boundary.sense = 0;
1149  new_boundary.name = "\0";
1150  // +ve sense
1151  if( atilla_cellname.find( "+" ) != std::string::npos )
1152  {
1153  new_boundary.sense = 1;
1154  // look for the @# we do not want it
1155  std::size_t found = atilla_cellname.find( "@" );
1156  if( found != std::string::npos )
1157  new_boundary.name = atilla_cellname.substr( 3, found );
1158  else
1159  new_boundary.name = atilla_cellname.substr( 3, atilla_cellname.length() );
1160  }
1161  else if( atilla_cellname.find( "-" ) != std::string::npos )
1162  {
1163  // negative sense
1164  new_boundary.sense = -1;
1165  new_boundary.name = atilla_cellname.substr( 3, atilla_cellname.length() );
1166  }
1167  return new_boundary;
1168 }
1169 
1170 /*
1171  * splits a string in a vector of strings split by spaces
1172  */
1173 std::vector< std::string > ReadRTT::split_string( std::string string_to_split, char split_char )
1174 {
1175  std::istringstream ss( string_to_split );
1176  std::vector< std::string > tokens;
1177  while( !ss.eof() )
1178  {
1179  std::string x; // here's a nice, empty string
1180  std::getline( ss, x, split_char ); // try to read the next field into it
1181  tokens.push_back( x );
1182  }
1183 
1184  // remove empty tokens
1185  std::vector< std::string >::iterator it;
1186  for( it = tokens.begin(); it != tokens.end(); )
1187  {
1188  std::string string = *it;
1189  if( string.compare( "\0" ) == 0 )
1190  it = tokens.erase( it );
1191  else
1192  ++it;
1193  }
1194  return tokens;
1195 }
1196 
1197 /*
1198  * Generate the parent-child links bwetween the cell and surface meshsets
1199  */
1201  std::vector< EntityHandle > entity_map[4],
1202  std::vector< side > side_data,
1203  std::vector< cell > cell_data )
1204 {
1205  ErrorCode rval; // return value
1206  // loop over the number of surfaces
1207  for( int i = 0; i < num_ents[2]; i++ )
1208  {
1209  // get the surface handle
1210  EntityHandle surf_handle = entity_map[2][i];
1211  // there are volumes that share this face
1212  for( unsigned int shared = 0; shared <= 1; shared++ )
1213  {
1214  std::string parent_name = side_data[i].names[shared];
1215  // find the @ sign
1216  unsigned pos = parent_name.find( "@" );
1217  parent_name = parent_name.substr( 0, pos );
1218 
1219  // loop over tets looking for matching name
1220  for( int j = 0; j < num_ents[3]; j++ )
1221  {
1222  // if match found
1223  if( cell_data[j].name.compare( parent_name ) == 0 )
1224  {
1225  EntityHandle cell_handle = entity_map[3][j];
1226  // parent
1227  rval = MBI->add_parent_child( cell_handle, surf_handle );
1228  if( rval != MB_SUCCESS )
1229  {
1230  std::cerr << "Failed to add parent child relationship" << std::endl;
1231  }
1232  }
1233  }
1234  }
1235  }
1236  return;
1237 }
1238 
1239 /*
1240  * sets the sense of the surfaces wrt to volumes using geom topo tool
1241  */
1242 void ReadRTT::set_surface_senses( int num_ents[4],
1243  std::vector< EntityHandle > entity_map[4],
1244  std::vector< side > side_data,
1245  std::vector< cell > cell_data )
1246 {
1247 
1248  ErrorCode rval; // return value
1249  // loop over the number of surfaces
1250  for( int i = 0; i < num_ents[2]; i++ )
1251  {
1252  EntityHandle surf_handle = entity_map[2][i];
1253  // there are 2 volumes that share this face
1254  for( unsigned int shared = 0; shared <= 1; shared++ )
1255  {
1256  std::string parent_name = side_data[i].names[shared];
1257  unsigned pos = parent_name.find( "@" );
1258  parent_name = parent_name.substr( 0, pos );
1259  // loop over tets looking for matching name
1260  for( int j = 0; j < num_ents[3]; j++ )
1261  {
1262  // if match found
1263  if( cell_data[j].name.compare( parent_name ) == 0 )
1264  {
1265  EntityHandle cell_handle = entity_map[3][j];
1266  // in rtt mesh +represents the inside and -represents outside
1267  // in moab reverse is outside and forward is inside
1268  if( side_data[i].senses[shared] == 1 )
1269  rval = myGeomTool->set_sense( surf_handle, cell_handle, SENSE_FORWARD );
1270  else if( side_data[i].senses[shared] == -1 )
1271  rval = myGeomTool->set_sense( surf_handle, cell_handle, SENSE_REVERSE );
1272  else
1273  rval = myGeomTool->set_sense( surf_handle, 0, SENSE_REVERSE );
1274 
1275  if( rval != MB_SUCCESS )
1276  {
1277  std::cerr << "Failed to set sense appropriately" << std::endl;
1278  }
1279  }
1280  }
1281  }
1282  }
1283  return;
1284 }
1285 
1286 /*
1287  * Add all entities that are to be part of the graveyard
1288  */
1289 ErrorCode ReadRTT::setup_group_data( std::vector< EntityHandle > entity_map[4],
1290  std::vector< tet > tet_data,
1291  std::map< int, EntityHandle >& volume_map )
1292 {
1293  ErrorCode rval; // error codes
1294  EntityHandle handle;
1295  handle = create_group( "graveyard_comp", 1 );
1296 
1297  // add any volume to group graveyard, it is ignored by dag
1298  EntityHandle vol_handle = entity_map[3][0];
1299  rval = MBI->add_entities( handle, &vol_handle, 1 );
1300 
1301  if( get_material_ref_flag() == "MATERIAL" )
1302  {
1303  std::string mat_flag = get_material_ref_flag();
1304  std::string vol_flag = get_volume_ref_flag();
1305  const std::vector< cell >& mat_cells = cell_flag_datas[mat_flag];
1306  const std::map< int, int >& mat_idx = cell_flag_indexes[mat_flag];
1307 
1308  std::map< int, int > volume2mat; // region → material
1309  std::map< int, EntityHandle > mat_groups; // material → group
1310 
1311  for( const auto& t : tet_data )
1312  {
1313  int mat_no = t.flag_values[cell_flag_idx[mat_flag]];
1314  int vol_no = t.flag_values[cell_flag_idx[vol_flag]];
1315 
1316  // record the material of this volume (consistency check)
1317  auto it = volume2mat.find( vol_no );
1318  if( it == volume2mat.end() )
1319  volume2mat[vol_no] = mat_no;
1320  else if( it->second != mat_no )
1321  {
1322  std::cerr << "Volume " << vol_no << " has conflicting material numbers: " << it->second << " and "
1323  << mat_no << std::endl;
1324  return MB_FAILURE;
1325  }
1326 
1327  // create material group the first time we meet this material
1328  if( mat_groups.find( mat_no ) == mat_groups.end() )
1329  {
1330  std::string name = mat_cells[mat_idx.at( mat_no )].name;
1331  if( name.rfind( "mat:", 0 ) != 0 ) name = "mat:" + name; // exactly one prefix
1332 
1333  EntityHandle mat_grp;
1334  rval = create_material_group( name, mat_no, mat_grp );MB_CHK_ERR( rval );
1335  mat_groups[mat_no] = mat_grp;
1336  }
1337  }
1338 
1339  // assigning volumes to material groups
1340  for( const auto& vp : volume2mat )
1341  {
1342  int vol_no = vp.first;
1343  int mat_no = vp.second;
1344  auto v_it = volume_map.find( vol_no );
1345  auto m_it = mat_groups.find( mat_no );
1346 
1347  if( v_it == volume_map.end() || m_it == mat_groups.end() )
1348  {
1349  std::cerr << "Missing handle while adding volume " << vol_no << " to material " << mat_no << std::endl;
1350  return MB_FAILURE;
1351  }
1352  EntityHandle vol_h = v_it->second;
1353  EntityHandle grp_h = m_it->second;
1354 
1355  rval = MBI->add_entities( grp_h, &vol_h, 1 );MB_CHK_ERR( rval );
1356  }
1357  }
1358 
1359  return rval;
1360 }
1361 
1362 /*
1363  * create a new group of with the name group name and id
1364  */
1365 EntityHandle ReadRTT::create_group( std::string group_name, int id )
1366 {
1367  ErrorCode rval;
1368  // category tags
1369  const char geom_categories[][CATEGORY_TAG_SIZE] = { "Vertex\0", "Curve\0", "Surface\0", "Volume\0", "Group\0" };
1370 
1371  EntityHandle handle;
1372  rval = MBI->create_meshset( MESHSET_SET, handle );
1373  if( MB_SUCCESS != rval ) return rval;
1374 
1375  rval = MBI->tag_set_data( name_tag, &handle, 1, group_name.c_str() );
1376  if( MB_SUCCESS != rval ) return MB_FAILURE;
1377 
1378  rval = MBI->tag_set_data( id_tag, &handle, 1, &id );
1379  if( MB_SUCCESS != rval ) return MB_FAILURE;
1380 
1381  rval = MBI->tag_set_data( category_tag, &handle, 1, &geom_categories[4] );
1382  if( MB_SUCCESS != rval ) return MB_FAILURE;
1383 
1384  return handle;
1385 }
1386 
1387 } // namespace moab