Mesh Oriented datABase  (version 5.5.0)
An array-based unstructured mesh library
ReadCGM.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 #ifdef WIN32
16 #pragma warning( disable : 4786 )
17 #endif
18 
19 #include "CGMConfig.h"
20 #include "GeometryQueryTool.hpp"
21 #include "ModelQueryEngine.hpp"
22 #include "RefEntityName.hpp"
23 #include "GMem.hpp"
24 
25 #include "RefGroup.hpp"
26 #include "RefVolume.hpp"
27 #include "RefFace.hpp"
28 #include "RefEdge.hpp"
29 #include "RefVertex.hpp"
30 
31 #include "SenseEntity.hpp"
32 #include "Surface.hpp"
33 #include "Curve.hpp"
34 #include "Body.hpp"
35 #include "InitCGMA.hpp"
36 
37 #include "moab/Core.hpp"
38 #include "moab/Interface.hpp"
39 #include "moab/Range.hpp"
40 #include "MBTagConventions.hpp"
41 #include "moab/FileOptions.hpp"
42 
43 #include "moab/GeomTopoTool.hpp"
44 
45 #include "CubitCompat.hpp"
46 
47 #include <cstdio>
48 #include <algorithm>
49 #include <cassert>
50 #include <iostream>
51 
52 #include "ReadCGM.hpp"
53 #include "moab/ReadUtilIface.hpp"
54 
55 namespace moab
56 {
57 
58 #define GF_CUBIT_FILE_TYPE "CUBIT"
59 #define GF_STEP_FILE_TYPE "STEP"
60 #define GF_IGES_FILE_TYPE "IGES"
61 #define GF_OCC_BREP_FILE_TYPE "OCC"
62 #define GF_FACET_FILE_TYPE "FACET"
63 
64 ReaderIface* ReadCGM::factory( Interface* iface )
65 {
66  return new ReadCGM( iface );
67 }
68 
69 ReadCGM::ReadCGM( Interface* impl )
70  : geom_tag( 0 ), id_tag( 0 ), name_tag( 0 ), category_tag( 0 ), faceting_tol_tag( 0 ), geometry_resabs_tag( 0 )
71 {
72  assert( NULL != impl );
73  mdbImpl = impl;
74  myGeomTool = new GeomTopoTool( impl );
75  impl->query_interface( readUtilIface );
76  assert( NULL != readUtilIface );
77 
78  // initialise counters
79  failed_curve_count = 0;
80  failed_surface_count = 0;
81 
82  ErrorCode rval;
83 
84  // get some tag handles
85  int negone = -1, zero = 0 /*, negonearr[] = {-1, -1, -1, -1}*/;
86  rval = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag, MB_TAG_SPARSE | MB_TAG_CREAT,
87  &negone );
88  assert( !rval );
89  id_tag = mdbImpl->globalId_tag();
90 
91  rval =
92  mdbImpl->tag_get_handle( NAME_TAG_NAME, NAME_TAG_SIZE, MB_TYPE_OPAQUE, name_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
93  assert( !rval );
94 
95  rval = mdbImpl->tag_get_handle( CATEGORY_TAG_NAME, CATEGORY_TAG_SIZE, MB_TYPE_OPAQUE, category_tag,
96  MB_TAG_SPARSE | MB_TAG_CREAT );
97  assert( !rval );
98  rval = mdbImpl->tag_get_handle( "FACETING_TOL", 1, MB_TYPE_DOUBLE, faceting_tol_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
99  assert( !rval );
100  rval = mdbImpl->tag_get_handle( "GEOMETRY_RESABS", 1, MB_TYPE_DOUBLE, geometry_resabs_tag,
101  MB_TAG_SPARSE | MB_TAG_CREAT );
102  assert( !rval );
103 #ifdef NDEBUG
104  if( !rval )
105  {
106  }; // Line to avoid compiler warning about variable set but not used
107 #endif
108 }
109 
110 ReadCGM::~ReadCGM()
111 {
112  mdbImpl->release_interface( readUtilIface );
113  delete myGeomTool;
114 }
115 
116 ErrorCode ReadCGM::read_tag_values( const char* /* file_name */,
117  const char* /* tag_name */,
118  const FileOptions& /* opts */,
119  std::vector< int >& /* tag_values_out */,
120  const SubsetList* /* subset_list */ )
121 {
122  return MB_NOT_IMPLEMENTED;
123 }
124 
125 // Sets options passed into ReadCGM::load_file
126 ErrorCode ReadCGM::set_options( const FileOptions& opts,
127  int& norm_tol,
128  double& faceting_tol,
129  double& len_tol,
130  bool& act_att,
131  bool& verbose_warnings,
132  bool& fatal_on_curves )
133 {
134  ErrorCode rval;
135 
136  // Default Values
137  int DEFAULT_NORM = 5;
138  double DEFAULT_FACET_TOL = 0.001;
139  double DEFAULT_LEN_TOL = 0.0;
140  act_att = true;
141 
142  // Check for the options.
143  if( MB_SUCCESS != opts.get_int_option( "FACET_NORMAL_TOLERANCE", norm_tol ) ) norm_tol = DEFAULT_NORM;
144 
145  if( MB_SUCCESS != opts.get_real_option( "FACET_DISTANCE_TOLERANCE", faceting_tol ) )
146  faceting_tol = DEFAULT_FACET_TOL;
147 
148  if( MB_SUCCESS != opts.get_real_option( "MAX_FACET_EDGE_LENGTH", len_tol ) ) len_tol = DEFAULT_LEN_TOL;
149 
150  if( MB_SUCCESS == opts.get_null_option( "VERBOSE_CGM_WARNINGS" ) ) verbose_warnings = true;
151 
152  if( MB_SUCCESS == opts.get_null_option( "FATAL_ON_CURVES" ) ) fatal_on_curves = true;
153 
154  const char* name = "CGM_ATTRIBS";
155  const char* value = "no";
156  rval = opts.match_option( name, value );
157  if( MB_SUCCESS == rval ) act_att = false;
158 
159  return MB_SUCCESS;
160 }
161 
162 ErrorCode ReadCGM::create_entity_sets( std::map< RefEntity*, EntityHandle > ( &entmap )[5] )
163 {
164  ErrorCode rval;
165  const char geom_categories[][CATEGORY_TAG_SIZE] = { "Vertex\0", "Curve\0", "Surface\0", "Volume\0", "Group\0" };
166  const char* const names[] = { "Vertex", "Curve", "Surface", "Volume" };
167  DLIList< RefEntity* > entlist;
168 
169  for( int dim = 0; dim < 4; dim++ )
170  {
171  entlist.clean_out();
172  GeometryQueryTool::instance()->ref_entity_list( names[dim], entlist, true );
173  entlist.reset();
174 
175  for( int i = entlist.size(); i--; )
176  {
177  RefEntity* ent = entlist.get_and_step();
178  EntityHandle handle;
179  // Create the new meshset
180  rval = mdbImpl->create_meshset( dim == 1 ? MESHSET_ORDERED : MESHSET_SET, handle );
181  if( MB_SUCCESS != rval ) return rval;
182 
183  // Map the geom reference entity to the corresponding moab meshset
184  entmap[dim][ent] = handle;
185 
186  // Create tags for the new meshset
187  rval = mdbImpl->tag_set_data( geom_tag, &handle, 1, &dim );
188  if( MB_SUCCESS != rval ) return rval;
189 
190  int id = ent->id();
191  rval = mdbImpl->tag_set_data( id_tag, &handle, 1, &id );
192  if( MB_SUCCESS != rval ) return rval;
193 
194  rval = mdbImpl->tag_set_data( category_tag, &handle, 1, &geom_categories[dim] );
195  if( MB_SUCCESS != rval ) return rval;
196  }
197  }
198 
199  return MB_SUCCESS;
200 }
201 
202 ErrorCode ReadCGM::create_topology( std::map< RefEntity*, EntityHandle > ( &entitymap )[5] )
203 {
204  ErrorCode rval;
205  DLIList< RefEntity* > entitylist;
206  std::map< RefEntity*, EntityHandle >::iterator ci;
207 
208  for( int dim = 1; dim < 4; ++dim )
209  {
210  for( ci = entitymap[dim].begin(); ci != entitymap[dim].end(); ++ci )
211  {
212  entitylist.clean_out();
213  ci->first->get_child_ref_entities( entitylist );
214 
215  entitylist.reset();
216  for( int i = entitylist.size(); i--; )
217  {
218  RefEntity* ent = entitylist.get_and_step();
219  EntityHandle h = entitymap[dim - 1][ent];
220  rval = mdbImpl->add_parent_child( ci->second, h );
221  if( MB_SUCCESS != rval ) return rval;
222  }
223  }
224  }
225 
226  return MB_SUCCESS;
227 }
228 
229 ErrorCode ReadCGM::store_surface_senses( std::map< RefEntity*, EntityHandle >& surface_map,
230  std::map< RefEntity*, EntityHandle >& volume_map )
231 {
232  ErrorCode rval;
233  std::map< RefEntity*, EntityHandle >::iterator ci;
234 
235  for( ci = surface_map.begin(); ci != surface_map.end(); ++ci )
236  {
237  RefFace* face = (RefFace*)( ci->first );
238  BasicTopologyEntity *forward = 0, *reverse = 0;
239  for( SenseEntity* cf = face->get_first_sense_entity_ptr(); cf; cf = cf->next_on_bte() )
240  {
241  BasicTopologyEntity* vol = cf->get_parent_basic_topology_entity_ptr();
242  // Allocate vol to the proper topology entity (forward or reverse)
243  if( cf->get_sense() == CUBIT_UNKNOWN || cf->get_sense() != face->get_surface_ptr()->bridge_sense() )
244  {
245  // Check that each surface has a sense for only one volume
246  if( reverse )
247  {
248  std::cout << "Surface " << face->id() << " has reverse sense "
249  << "with multiple volume " << reverse->id() << " and "
250  << "volume " << vol->id() << std::endl;
251  return MB_FAILURE;
252  }
253  reverse = vol;
254  }
255  if( cf->get_sense() == CUBIT_UNKNOWN || cf->get_sense() == face->get_surface_ptr()->bridge_sense() )
256  {
257  // Check that each surface has a sense for only one volume
258  if( forward )
259  {
260  std::cout << "Surface " << face->id() << " has forward sense "
261  << "with multiple volume " << forward->id() << " and "
262  << "volume " << vol->id() << std::endl;
263  return MB_FAILURE;
264  }
265  forward = vol;
266  }
267  }
268 
269  if( forward )
270  {
271  rval = myGeomTool->set_sense( ci->second, volume_map[forward], SENSE_FORWARD );
272  if( MB_SUCCESS != rval ) return rval;
273  }
274  if( reverse )
275  {
276  rval = myGeomTool->set_sense( ci->second, volume_map[reverse], SENSE_REVERSE );
277  if( MB_SUCCESS != rval ) return rval;
278  }
279  }
280 
281  return MB_SUCCESS;
282 }
283 
284 ErrorCode ReadCGM::store_curve_senses( std::map< RefEntity*, EntityHandle >& curve_map,
285  std::map< RefEntity*, EntityHandle >& surface_map )
286 {
287  ErrorCode rval;
288  std::vector< EntityHandle > ents;
289  std::vector< int > senses;
290  std::map< RefEntity*, EntityHandle >::iterator ci;
291  for( ci = curve_map.begin(); ci != curve_map.end(); ++ci )
292  {
293  RefEdge* edge = (RefEdge*)( ci->first );
294  ents.clear();
295  senses.clear();
296  for( SenseEntity* ce = edge->get_first_sense_entity_ptr(); ce; ce = ce->next_on_bte() )
297  {
298  BasicTopologyEntity* fac = ce->get_parent_basic_topology_entity_ptr();
299  EntityHandle face = surface_map[fac];
300  if( ce->get_sense() == CUBIT_UNKNOWN || ce->get_sense() != edge->get_curve_ptr()->bridge_sense() )
301  {
302  ents.push_back( face );
303  senses.push_back( SENSE_REVERSE );
304  }
305  if( ce->get_sense() == CUBIT_UNKNOWN || ce->get_sense() == edge->get_curve_ptr()->bridge_sense() )
306  {
307  ents.push_back( face );
308  senses.push_back( SENSE_FORWARD );
309  }
310  }
311 
312  rval = myGeomTool->set_senses( ci->second, ents, senses );
313  if( MB_SUCCESS != rval ) return rval;
314  }
315  return MB_SUCCESS;
316 }
317 
318 ErrorCode ReadCGM::store_groups( std::map< RefEntity*, EntityHandle > ( &entitymap )[5] )
319 {
320  ErrorCode rval;
321 
322  // Create entity sets for all ref groups
323  rval = create_group_entsets( entitymap[4] );
324  if( rval != MB_SUCCESS ) return rval;
325 
326  // Store group names and entities in the mesh
327  rval = store_group_content( entitymap );
328  if( rval != MB_SUCCESS ) return rval;
329 
330  return MB_SUCCESS;
331 }
332 
333 ErrorCode ReadCGM::create_group_entsets( std::map< RefEntity*, EntityHandle >& group_map )
334 {
335  ErrorCode rval;
336  const char geom_categories[][CATEGORY_TAG_SIZE] = { "Vertex\0", "Curve\0", "Surface\0", "Volume\0", "Group\0" };
337  DLIList< RefEntity* > entitylist;
338  // Create entity sets for all ref groups
339  std::vector< Tag > extra_name_tags;
340 #if CGM_MAJOR_VERSION > 13
341  DLIList< CubitString > name_list;
342 #else
343  DLIList< CubitString* > name_list;
344 #endif
345  entitylist.clean_out();
346  // Get all entity groups from the CGM model
347  GeometryQueryTool::instance()->ref_entity_list( "group", entitylist );
348  entitylist.reset();
349  // Loop over all groups
350  for( int i = entitylist.size(); i--; )
351  {
352  // Take the next group
353  RefEntity* grp = entitylist.get_and_step();
354  name_list.clean_out();
355 // Get the names of all entities in this group from the solid model
356 #if CGM_MAJOR_VERSION > 13
357  RefEntityName::instance()->get_refentity_name( grp, name_list );
358 #else
359  // True argument is optional, but for large multi-names situation, it should save
360  // some cpu time
361  RefEntityName::instance()->get_refentity_name( grp, name_list, true );
362 #endif
363  if( name_list.size() == 0 ) continue;
364  // Set pointer to first name of the group and set the first name to name1
365  name_list.reset();
366 #if CGM_MAJOR_VERSION > 13
367  CubitString name1 = name_list.get();
368 #else
369  CubitString name1 = *name_list.get();
370 #endif
371  // Create entity handle for the group
372  EntityHandle h;
373  rval = mdbImpl->create_meshset( MESHSET_SET, h );
374  if( MB_SUCCESS != rval ) return rval;
375  // Set tag data for the group
376  char namebuf[NAME_TAG_SIZE];
377  memset( namebuf, '\0', NAME_TAG_SIZE );
378  strncpy( namebuf, name1.c_str(), NAME_TAG_SIZE - 1 );
379  if( name1.length() >= (unsigned)NAME_TAG_SIZE )
380  std::cout << "WARNING: group name '" << name1.c_str() << "' truncated to '" << namebuf << "'" << std::endl;
381  rval = mdbImpl->tag_set_data( name_tag, &h, 1, namebuf );
382  if( MB_SUCCESS != rval ) return MB_FAILURE;
383 
384  int id = grp->id();
385  rval = mdbImpl->tag_set_data( id_tag, &h, 1, &id );
386  if( MB_SUCCESS != rval ) return MB_FAILURE;
387 
388  rval = mdbImpl->tag_set_data( category_tag, &h, 1, &geom_categories[4] );
389  if( MB_SUCCESS != rval ) return MB_FAILURE;
390  // Check for extra group names
391  if( name_list.size() > 1 )
392  {
393  for( int j = extra_name_tags.size(); j < name_list.size(); ++j )
394  {
395  sprintf( namebuf, "EXTRA_%s%d", NAME_TAG_NAME, j );
396  Tag t;
397  rval =
398  mdbImpl->tag_get_handle( namebuf, NAME_TAG_SIZE, MB_TYPE_OPAQUE, t, MB_TAG_SPARSE | MB_TAG_CREAT );
399  assert( !rval );
400  extra_name_tags.push_back( t );
401  }
402  // Add extra group names to the group handle
403  for( int j = 0; j < name_list.size(); ++j )
404  {
405 #if CGM_MAJOR_VERSION > 13
406  name1 = name_list.get_and_step();
407 #else
408  name1 = *name_list.get_and_step();
409 #endif
410  memset( namebuf, '\0', NAME_TAG_SIZE );
411  strncpy( namebuf, name1.c_str(), NAME_TAG_SIZE - 1 );
412  if( name1.length() >= (unsigned)NAME_TAG_SIZE )
413  std::cout << "WARNING: group name '" << name1.c_str() << "' truncated to '" << namebuf << "'"
414  << std::endl;
415  rval = mdbImpl->tag_set_data( extra_name_tags[j], &h, 1, namebuf );
416  if( MB_SUCCESS != rval ) return MB_FAILURE;
417  }
418  }
419  // Add the group handle
420  group_map[grp] = h;
421  }
422 
423  return MB_SUCCESS;
424 }
425 
426 ErrorCode ReadCGM::store_group_content( std::map< RefEntity*, EntityHandle > ( &entitymap )[5] )
427 {
428  ErrorCode rval;
429  DLIList< RefEntity* > entlist;
430  std::map< RefEntity*, EntityHandle >::iterator ci;
431  // Store contents for each group
432  entlist.reset();
433  for( ci = entitymap[4].begin(); ci != entitymap[4].end(); ++ci )
434  {
435  RefGroup* grp = (RefGroup*)( ci->first );
436  entlist.clean_out();
437  grp->get_child_ref_entities( entlist );
438 
439  Range entities;
440  while( entlist.size() )
441  {
442  RefEntity* ent = entlist.pop();
443  int dim = ent->dimension();
444 
445  if( dim < 0 )
446  {
447  Body* body;
448  if( entitymap[4].find( ent ) != entitymap[4].end() )
449  {
450  // Child is another group; examine its contents
451  entities.insert( entitymap[4][ent] );
452  }
453  else if( ( body = dynamic_cast< Body* >( ent ) ) != NULL )
454  {
455  // Child is a CGM Body, which presumably comprises some volumes--
456  // extract volumes as if they belonged to group.
457  DLIList< RefVolume* > vols;
458  body->ref_volumes( vols );
459  for( int vi = vols.size(); vi--; )
460  {
461  RefVolume* vol = vols.get_and_step();
462  if( entitymap[3].find( vol ) != entitymap[3].end() )
463  {
464  entities.insert( entitymap[3][vol] );
465  }
466  else
467  {
468  std::cerr << "Warning: CGM Body has orphan RefVolume" << std::endl;
469  }
470  }
471  }
472  else
473  {
474  // Otherwise, warn user.
475  std::cerr << "Warning: A dim<0 entity is being ignored by ReadCGM." << std::endl;
476  }
477  }
478  else if( dim < 4 )
479  {
480  if( entitymap[dim].find( ent ) != entitymap[dim].end() ) entities.insert( entitymap[dim][ent] );
481  }
482  }
483 
484  if( !entities.empty() )
485  {
486  rval = mdbImpl->add_entities( ci->second, entities );
487  if( MB_SUCCESS != rval ) return MB_FAILURE;
488  }
489  }
490 
491  return MB_SUCCESS;
492 }
493 
494 void ReadCGM::set_cgm_attributes( bool const act_attributes, bool const verbose )
495 {
496  if( act_attributes )
497  {
498  CGMApp::instance()->attrib_manager()->set_all_auto_read_flags( act_attributes );
499  CGMApp::instance()->attrib_manager()->set_all_auto_actuate_flags( act_attributes );
500  }
501 
502  if( !verbose )
503  {
504  CGMApp::instance()->attrib_manager()->silent_flag( true );
505  }
506 }
507 
508 ErrorCode ReadCGM::create_vertices( std::map< RefEntity*, EntityHandle >& vertex_map )
509 {
510  ErrorCode rval;
511  std::map< RefEntity*, EntityHandle >::iterator ci;
512  for( ci = vertex_map.begin(); ci != vertex_map.end(); ++ci )
513  {
514  CubitVector pos = dynamic_cast< RefVertex* >( ci->first )->coordinates();
515  double coords[3] = { pos.x(), pos.y(), pos.z() };
516  EntityHandle vh;
517  rval = mdbImpl->create_vertex( coords, vh );
518  if( MB_SUCCESS != rval ) return MB_FAILURE;
519 
520  // Add the vertex to its tagged meshset
521  rval = mdbImpl->add_entities( ci->second, &vh, 1 );
522  if( MB_SUCCESS != rval ) return MB_FAILURE;
523 
524  // Replace the meshset handle with the vertex handle
525  // This makes adding the vertex to higher dim sets easier
526  ci->second = vh;
527  }
528 
529  return MB_SUCCESS;
530 }
531 
532 ErrorCode ReadCGM::create_curve_facets( std::map< RefEntity*, EntityHandle >& curve_map,
533  std::map< RefEntity*, EntityHandle >& vertex_map,
534 #if CGM_MAJOR_VERSION > 12
535  int norm_tol,
536 #else
537  int /* norm_tol */,
538 #endif
539  double faceting_tol,
540  bool verbose_warn,
541  bool fatal_on_curves )
542 {
543  ErrorCode rval;
544  CubitStatus s;
545  // Maximum allowable curve-endpoint proximity warnings
546  // If this integer becomes negative, then abs(curve_warnings) is the
547  // number of warnings that were suppressed.
548  int curve_warnings = 0;
549 
550  // Map iterator
551  std::map< RefEntity*, EntityHandle >::iterator ci;
552 
553  // Create geometry for all curves
554  GMem data;
555  for( ci = curve_map.begin(); ci != curve_map.end(); ++ci )
556  {
557  // Get the start and end points of the curve in the form of a reference edge
558  RefEdge* edge = dynamic_cast< RefEdge* >( ci->first );
559  // Get the edge's curve information
560  Curve* curve = edge->get_curve_ptr();
561  // Clean out previous curve information
562  data.clean_out();
563  // Facet curve according to parameters and CGM version
564 #if CGM_MAJOR_VERSION > 12
565  s = edge->get_graphics( data, norm_tol, faceting_tol );
566 #else
567  s = edge->get_graphics( data, faceting_tol );
568 #endif
569 
570  if( s != CUBIT_SUCCESS )
571  {
572  // if we fatal on curves
573  if( fatal_on_curves )
574  {
575  std::cout << "Failed to facet the curve " << edge->id() << std::endl;
576  return MB_FAILURE;
577  }
578  // otherwise record them
579  else
580  {
581  failed_curve_count++;
582  failed_curves.push_back( edge->id() );
583  }
584  continue;
585  }
586 
587  std::vector< CubitVector > points;
588  for( int i = 0; i < data.pointListCount; ++i )
589  // Add Cubit vertext points to a list
590  points.push_back( CubitVector( data.point_list()[i].x, data.point_list()[i].y, data.point_list()[i].z ) );
591 
592  // Need to reverse data?
593  if( curve->bridge_sense() == CUBIT_REVERSED ) std::reverse( points.begin(), points.end() );
594 
595  // Check for closed curve
596  RefVertex *start_vtx, *end_vtx;
597  start_vtx = edge->start_vertex();
598  end_vtx = edge->end_vertex();
599 
600  // Special case for point curve
601  if( points.size() < 2 )
602  {
603  if( start_vtx != end_vtx || curve->measure() > GEOMETRY_RESABS )
604  {
605  std::cerr << "Warning: No facetting for curve " << edge->id() << std::endl;
606  continue;
607  }
608  EntityHandle h = vertex_map[start_vtx];
609  rval = mdbImpl->add_entities( ci->second, &h, 1 );
610  if( MB_SUCCESS != rval ) return MB_FAILURE;
611  continue;
612  }
613  // Check to see if the first and last interior vertices are considered to be
614  // coincident by CUBIT
615  const bool closed = ( points.front() - points.back() ).length() < GEOMETRY_RESABS;
616  if( closed != ( start_vtx == end_vtx ) )
617  {
618  std::cerr << "Warning: topology and geometry inconsistant for possibly closed curve " << edge->id()
619  << std::endl;
620  }
621 
622  // Check proximity of vertices to end coordinates
623  if( ( start_vtx->coordinates() - points.front() ).length() > GEOMETRY_RESABS ||
624  ( end_vtx->coordinates() - points.back() ).length() > GEOMETRY_RESABS )
625  {
626 
627  curve_warnings--;
628  if( curve_warnings >= 0 || verbose_warn )
629  {
630  std::cerr << "Warning: vertices not at ends of curve " << edge->id() << std::endl;
631  if( curve_warnings == 0 && !verbose_warn )
632  {
633  std::cerr << " further instances of this warning will be suppressed..." << std::endl;
634  }
635  }
636  }
637  // Create interior points
638  std::vector< EntityHandle > verts, edges;
639  verts.push_back( vertex_map[start_vtx] );
640  for( size_t i = 1; i < points.size() - 1; ++i )
641  {
642  double coords[] = { points[i].x(), points[i].y(), points[i].z() };
643  EntityHandle h;
644  // Create vertex entity
645  rval = mdbImpl->create_vertex( coords, h );
646  if( MB_SUCCESS != rval ) return MB_FAILURE;
647  verts.push_back( h );
648  }
649  verts.push_back( vertex_map[end_vtx] );
650 
651  // Create edges
652  for( size_t i = 0; i < verts.size() - 1; ++i )
653  {
654  EntityHandle h;
655  rval = mdbImpl->create_element( MBEDGE, &verts[i], 2, h );
656  if( MB_SUCCESS != rval ) return MB_FAILURE;
657  edges.push_back( h );
658  }
659 
660  // If closed, remove duplicate
661  if( verts.front() == verts.back() ) verts.pop_back();
662  // Add entities to the curve meshset from entitymap
663  rval = mdbImpl->add_entities( ci->second, &verts[0], verts.size() );
664  if( MB_SUCCESS != rval ) return MB_FAILURE;
665  rval = mdbImpl->add_entities( ci->second, &edges[0], edges.size() );
666  if( MB_SUCCESS != rval ) return MB_FAILURE;
667  }
668 
669  if( !verbose_warn && curve_warnings < 0 )
670  {
671  std::cerr << "Suppressed " << -curve_warnings << " 'vertices not at ends of curve' warnings." << std::endl;
672  std::cerr << "To see all warnings, use reader param VERBOSE_CGM_WARNINGS." << std::endl;
673  }
674 
675  return MB_SUCCESS;
676 }
677 
678 ErrorCode ReadCGM::create_surface_facets( std::map< RefEntity*, EntityHandle >& surface_map,
679  std::map< RefEntity*, EntityHandle >& vertex_map,
680  int norm_tol,
681  double facet_tol,
682  double length_tol )
683 {
684  ErrorCode rval;
685  std::map< RefEntity*, EntityHandle >::iterator ci;
686  CubitStatus s;
687 #if( ( CGM_MAJOR_VERSION == 14 && CGM_MINOR_VERSION > 2 ) || CGM_MAJOR_VERSION >= 15 )
688  DLIList< TopologyEntity* > me_list;
689 #else
690  DLIList< ModelEntity* > me_list;
691 #endif
692 
693  GMem data;
694  // Create geometry for all surfaces
695  for( ci = surface_map.begin(); ci != surface_map.end(); ++ci )
696  {
697  RefFace* face = dynamic_cast< RefFace* >( ci->first );
698 
699  data.clean_out();
700  s = face->get_graphics( data, norm_tol, facet_tol, length_tol );
701 
702  if( CUBIT_SUCCESS != s ) return MB_FAILURE;
703 
704  // Declare array of all vertex handles
705  std::vector< EntityHandle > verts( data.pointListCount, 0 );
706 
707  // Get list of geometric vertices in surface
708  me_list.clean_out();
709  ModelQueryEngine::instance()->query_model( *face, DagType::ref_vertex_type(), me_list );
710 
711  // For each geometric vertex, find a single coincident point in facets
712  // Otherwise, print a warning
713  for( int i = me_list.size(); i--; )
714  {
715  // Assign geometric vertex
716  RefVertex* vtx = dynamic_cast< RefVertex* >( me_list.get_and_step() );
717  CubitVector pos = vtx->coordinates();
718 
719  for( int j = 0; j < data.pointListCount; ++j )
720  {
721  // Assign facet vertex
722  CubitVector vpos( data.point_list()[j].x, data.point_list()[j].y, data.point_list()[j].z );
723  // Check to see if they are considered coincident
724  if( ( pos - vpos ).length_squared() < GEOMETRY_RESABS * GEOMETRY_RESABS )
725  {
726  // If this facet vertex has already been found coincident, print warning
727  if( verts[j] ) std::cerr << "Warning: Coincident vertices in surface " << face->id() << std::endl;
728  // If a coincidence is found, keep track of it in the verts vector
729  verts[j] = vertex_map[vtx];
730  break;
731  }
732  }
733  }
734 
735  // Now create vertices for the remaining points in the facetting
736  for( int i = 0; i < data.pointListCount; ++i )
737  {
738  if( verts[i] ) // If a geometric vertex
739  continue;
740  double coords[] = { data.point_list()[i].x, data.point_list()[i].y, data.point_list()[i].z };
741  // Return vertex handle to verts to fill in all remaining facet
742  // vertices
743  rval = mdbImpl->create_vertex( coords, verts[i] );
744  if( MB_SUCCESS != rval ) return rval;
745  }
746 
747  // record the failures for information
748  if( data.fListCount == 0 )
749  {
750  failed_surface_count++;
751  failed_surfaces.push_back( face->id() );
752  }
753 
754  // Now create facets
755  Range facets;
756  std::vector< EntityHandle > corners;
757  for( int i = 0; i < data.fListCount; i += data.facet_list()[i] + 1 )
758  {
759  // Get number of facet verts
760  int* facet = data.facet_list() + i;
761  corners.resize( *facet );
762  for( int j = 1; j <= *facet; ++j )
763  {
764  if( facet[j] >= (int)verts.size() )
765  {
766  std::cerr << "ERROR: Invalid facet data for surface " << face->id() << std::endl;
767  return MB_FAILURE;
768  }
769  corners[j - 1] = verts[facet[j]];
770  }
771  EntityType type;
772  if( *facet == 3 )
773  type = MBTRI;
774  else
775  {
776  std::cerr << "Warning: non-triangle facet in surface " << face->id() << std::endl;
777  std::cerr << " entity has " << *facet << " edges" << std::endl;
778  if( *facet == 4 )
779  type = MBQUAD;
780  else
781  type = MBPOLYGON;
782  }
783 
784  // if (surf->bridge_sense() == CUBIT_REVERSED)
785  // std::reverse(corners.begin(), corners.end());
786 
787  EntityHandle h;
788  rval = mdbImpl->create_element( type, &corners[0], corners.size(), h );
789  if( MB_SUCCESS != rval ) return MB_FAILURE;
790 
791  facets.insert( h );
792  }
793 
794  // Add vertices and facets to surface set
795  rval = mdbImpl->add_entities( ci->second, &verts[0], verts.size() );
796  if( MB_SUCCESS != rval ) return MB_FAILURE;
797  rval = mdbImpl->add_entities( ci->second, facets );
798  if( MB_SUCCESS != rval ) return MB_FAILURE;
799  }
800 
801  return MB_SUCCESS;
802 }
803 
804 // Copy geometry into mesh database
805 ErrorCode ReadCGM::load_file( const char* cgm_file_name,
806  const EntityHandle* file_set,
807  const FileOptions& opts,
808  const ReaderIface::SubsetList* subset_list,
809  const Tag* /*file_id_tag*/ )
810 {
811  // Blocks_to_load and num_blocks are ignored.
812  ErrorCode rval;
813 
814  if( subset_list )
815  {
816  MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for CGM data" );
817  }
818 
819  int norm_tol;
820  double faceting_tol;
821  double len_tol;
822  bool act_att = true;
823  bool verbose_warnings = false;
824  bool fatal_on_curves = false;
825 
826  rval = set_options( opts, norm_tol, faceting_tol, len_tol, act_att, verbose_warnings, fatal_on_curves );
827  if( MB_SUCCESS != rval ) return rval;
828 
829  // Always tag with the faceting_tol and geometry absolute resolution
830  // If file_set is defined, use that, otherwise (file_set == NULL) tag the interface
831  EntityHandle set = file_set ? *file_set : 0;
832  rval = mdbImpl->tag_set_data( faceting_tol_tag, &set, 1, &faceting_tol );
833  if( MB_SUCCESS != rval ) return rval;
834 
835  rval = mdbImpl->tag_set_data( geometry_resabs_tag, &set, 1, &GEOMETRY_RESABS );
836  if( MB_SUCCESS != rval ) return rval;
837 
838  // Initialize CGM
839  InitCGMA::initialize_cgma();
840 
841  // Determine CGM settings and amount of output
842  set_cgm_attributes( act_att, verbose_warnings );
843 
844  CubitStatus s;
845 
846  // Get CGM file type
847  const char* file_type = 0;
848  file_type = get_geom_file_type( cgm_file_name );
849  if( !file_type || !strcmp( file_type, "CUBIT" ) ) return MB_FAILURE;
850 
851  s = CubitCompat_import_solid_model( cgm_file_name, file_type );
852  if( CUBIT_SUCCESS != s )
853  {
854  MB_SET_ERR( MB_FAILURE, cgm_file_name << ": Failed to read file of type \"" << file_type << "\"" );
855  }
856 
857  // Create entity sets for all geometric entities
858  std::map< RefEntity*, EntityHandle > entmap[5]; // One for each dim, and one for groups
859 
860  rval = create_entity_sets( entmap );
861  if( rval != MB_SUCCESS ) return rval;
862 
863  // Create topology for all geometric entities
864  rval = create_topology( entmap );
865  if( rval != MB_SUCCESS ) return rval;
866 
867  // Store CoFace senses
868  rval = store_surface_senses( entmap[2], entmap[3] );
869  if( rval != MB_SUCCESS ) return rval;
870 
871  // Store CoEdge senses
872  rval = store_curve_senses( entmap[1], entmap[2] );
873  if( rval != MB_SUCCESS ) return rval;
874 
875  // Get group information and store it in the mesh
876  rval = store_groups( entmap );
877  if( rval != MB_SUCCESS ) return rval;
878 
879  // Done with volumes and groups
880  entmap[3].clear();
881  entmap[4].clear();
882 
883  // Create geometry for all vertices and replace
884  rval = create_vertices( entmap[0] );
885  if( rval != MB_SUCCESS ) return rval;
886 
887  // Create facets for all curves
888  rval = create_curve_facets( entmap[1], entmap[0], norm_tol, faceting_tol, verbose_warnings, fatal_on_curves );
889  if( rval != MB_SUCCESS ) return rval;
890 
891  // Create facets for surfaces
892  rval = create_surface_facets( entmap[2], entmap[0], norm_tol, faceting_tol, len_tol );
893  if( rval != MB_SUCCESS ) return rval;
894 
895  // print the fail information
896  dump_fail_counts();
897 
898  return MB_SUCCESS;
899 }
900 
901 // return the number of curves that failed to facet
902 int ReadCGM::get_failed_curve_count()
903 {
904  return failed_curve_count;
905 }
906 
907 // return the number of surfaces that failed to facet
908 int ReadCGM::get_failed_surface_count()
909 {
910  return failed_surface_count;
911 }
912 
913 void ReadCGM::dump_fail_counts()
914 {
915  std::cout << "***** Faceting Summary Information *****" << std::endl;
916  std::cout << "----- Curve Fail Information -----" << std::endl;
917  std::cout << "There were " << failed_curve_count << " curves that could not be faceted." << std::endl;
918 
919  if( failed_curve_count > 0 )
920  {
921  std::cout << "The curves were ";
922  for( int i = 0; i < failed_curve_count; i++ )
923  {
924  std::cout << failed_curves[i] << " ";
925  if( ( i % 10 == 0 ) & ( i > 0 ) ) std::cout << std::endl;
926  }
927  }
928  std::cout << std::endl;
929  std::cout << "----- Facet Fail Information -----" << std::endl;
930  std::cout << "There were " << failed_surface_count << " surfaces that could not be faceted." << std::endl;
931  if( failed_surface_count > 0 )
932  {
933  std::cout << "The surfaces were ";
934  for( int i = 0; i < failed_surface_count; i++ )
935  {
936  std::cout << failed_surfaces[i] << " ";
937  if( ( i % 10 == 0 ) & ( i > 0 ) ) std::cout << std::endl;
938  }
939  }
940  std::cout << std::endl;
941  std::cout << "***** End of Faceting Summary Information *****" << std::endl;
942  return;
943 }
944 
945 const char* ReadCGM::get_geom_file_type( const char* name )
946 {
947  FILE* file;
948  const char* result = 0;
949 
950  file = fopen( name, "r" );
951  if( file )
952  {
953  result = get_geom_fptr_type( file );
954  fclose( file );
955  }
956 
957  return result;
958 }
959 
960 const char* ReadCGM::get_geom_fptr_type( FILE* file )
961 {
962  static const char* CUBIT_NAME = GF_CUBIT_FILE_TYPE;
963  static const char* STEP_NAME = GF_STEP_FILE_TYPE;
964  static const char* IGES_NAME = GF_IGES_FILE_TYPE;
965  static const char* BREP_NAME = GF_OCC_BREP_FILE_TYPE;
966  static const char* FACET_NAME = GF_FACET_FILE_TYPE;
967 
968  if( is_cubit_file( file ) )
969  return CUBIT_NAME;
970  else if( is_step_file( file ) )
971  return STEP_NAME;
972  else if( is_iges_file( file ) )
973  return IGES_NAME;
974  else if( is_occ_brep_file( file ) )
975  return BREP_NAME;
976  else if( is_facet_file( file ) )
977  return FACET_NAME;
978  else
979  return NULL;
980 }
981 
982 int ReadCGM::is_cubit_file( FILE* file )
983 {
984  unsigned char buffer[4];
985  return !fseek( file, 0, SEEK_SET ) && fread( buffer, 4, 1, file ) && !memcmp( buffer, "CUBE", 4 );
986 }
987 
988 int ReadCGM::is_step_file( FILE* file )
989 {
990  unsigned char buffer[9];
991  return !fseek( file, 0, SEEK_SET ) && fread( buffer, 9, 1, file ) && !memcmp( buffer, "ISO-10303", 9 );
992 }
993 
994 int ReadCGM::is_iges_file( FILE* file )
995 {
996  unsigned char buffer[10];
997  return !fseek( file, 72, SEEK_SET ) && fread( buffer, 10, 1, file ) && !memcmp( buffer, "S 1", 8 );
998 }
999 
1000 int ReadCGM::is_occ_brep_file( FILE* file )
1001 {
1002  unsigned char buffer[6];
1003  return !fseek( file, 0, SEEK_SET ) && fread( buffer, 6, 1, file ) && !memcmp( buffer, "DBRep_", 6 );
1004 }
1005 int ReadCGM::is_facet_file( FILE* file )
1006 {
1007  unsigned char buffer[10];
1008  return !fseek( file, 0, SEEK_SET ) && fread( buffer, 10, 1, file ) && !memcmp( buffer, "MESH_BASED", 10 );
1009 }
1010 
1011 void ReadCGM::tokenize( const std::string& str, std::vector< std::string >& tokens, const char* delimiters )
1012 {
1013  std::string::size_type last = str.find_first_not_of( delimiters, 0 );
1014  std::string::size_type pos = str.find_first_of( delimiters, last );
1015  while( std::string::npos != pos && std::string::npos != last )
1016  {
1017  tokens.push_back( str.substr( last, pos - last ) );
1018  last = str.find_first_not_of( delimiters, pos );
1019  pos = str.find_first_of( delimiters, last );
1020  if( std::string::npos == pos ) pos = str.size();
1021  }
1022 }
1023 
1024 } // namespace moab