DagMC.hpp

Go to the documentation of this file.

#ifndef MOABMC_HPP
#define MOABMC_HPP

#include "moab/Interface.hpp"
#include "MBTagConventions.hpp"

#include <vector>
#include <map>
#include <string>
#include <assert.h>

#include "moab/OrientedBoxTreeTool.hpp"

class RefEntity;

struct DagmcVolData {
  int mat_id;
  double density, importance;
  std::string comp_name;
};
  

namespace moab {

class CartVect;

#define DAGMC_VERSION 1.0
#define DAGMC_VERSION_STRING "1.0"
#define DAGMC_INTERFACE_REVISION "$Rev$"


class DagMC 
{
public:
  static DagMC *instance(Interface *mb_impl = NULL);

  
  ~DagMC();

  static float version(std::string *version_string = NULL);
  static unsigned int interface_revision();


  /* SECTION I: Geometry Initialization */

  ErrorCode load_file(const char* cfile,
                      const double facet_tolerance = 0);

  /*\brief Use pre-loaded geometry set
   *
   * Works like load_file, but using data that has been externally 
   * loaded into DagMC's MOAB instance. 
   * Only one of the two functions should be called.
   * 
   * TODO: this function should accept a parameter, being the
   * entity set to use for DagMC's data.  Currently DagMC always
   * assumes that all the contents of its MOAB instance belong to it.
   */
  ErrorCode load_existing_contents();

  ErrorCode init_OBBTree();

private:
  ErrorCode finish_loading(); 

  bool have_obb_tree();

  ErrorCode get_impl_compl();
  
  ErrorCode build_obbs(Range &surfs, Range &vols);

  ErrorCode build_obb_impl_compl(Range &surfs);


  /* SECTION II: Fundamental Geometry Operations/Queries */
public:

  class RayHistory {

  public:
    void reset();

    void reset_to_last_intersection(); 

    void rollback_last_intersection();

    int size() const { return prev_facets.size(); }

  private:    
    std::vector<EntityHandle> prev_facets;

    friend class DagMC;
    
  };

  ErrorCode ray_fire(const EntityHandle volume, 
                     const double ray_start[3], const double ray_dir[3],
                     EntityHandle& next_surf, double& next_surf_dist,
                     RayHistory* history = NULL, double dist_limit = 0,
             int ray_orientation = 1, 
                     OrientedBoxTreeTool::TrvStats* stats = NULL );
  
  ErrorCode point_in_volume(const EntityHandle volume, 
                            const double xyz[3],
                            int& result,
                            const double* uvw = NULL,
                            const RayHistory* history = NULL );

  ErrorCode point_in_volume_slow( const EntityHandle volume, const double xyz[3], int& result ); 


  ErrorCode test_volume_boundary( const EntityHandle volume, const EntityHandle surface,
                                  const double xyz[3], const double uvw[3], int& result,
                                  const RayHistory* history = NULL );

  ErrorCode closest_to_location( EntityHandle volume, const double point[3], double& result);

  ErrorCode measure_volume( EntityHandle volume, double& result );

  ErrorCode measure_area( EntityHandle surface, double& result );

  ErrorCode surface_sense( EntityHandle volume, 
                             int num_surfaces,
                             const EntityHandle* surfaces,
                             int* senses_out );

  ErrorCode surface_sense( EntityHandle volume, EntityHandle surface, int& sense_out );

  ErrorCode get_angle(EntityHandle surf, const double xyz[3], double angle[3], 
                      const RayHistory* history = NULL );

  ErrorCode next_vol( EntityHandle surface, EntityHandle old_volume, 
                      EntityHandle& new_volume );

private:
  ErrorCode CAD_ray_intersect(const double *point, 
                                const double *dir, 
                                const double huge_val,
                                std::vector<double> &distances,
                                std::vector<EntityHandle> &surfaces, 
                                double &len);

  ErrorCode boundary_case( EntityHandle volume, int& result, 
                             double u, double v, double w,
                             EntityHandle facet,
                             EntityHandle surface);

  ErrorCode poly_solid_angle( EntityHandle face, const CartVect& point, double& area );

  /* SECTION III: Indexing & Cross-referencing */
public:
  /* Most calling apps refer to geometric entities with a combination of
   * base-1/0 ordinal index (or rank) and global ID (or name).
   * DagMC also has an internal EntityHandle reference to each geometric entity.
   * These method provide ways to translate from one to the other.
   */

  EntityHandle entity_by_index( int dimension, int index );
  int id_by_index( int dimension, int index );
  EntityHandle entity_by_id( int dimension, int id );
  int index_by_handle( EntityHandle handle );
  int get_entity_id(EntityHandle this_ent);

  int num_entities( int dimension );

private:
  ErrorCode build_indices(Range &surfs, Range &vols);
  

  /* SECTION IV: Handling DagMC settings */
public:

  double overlap_thickness() {return overlapThickness;}
  double numerical_precision() {return numericalPrecision;}
  double faceting_tolerance() {return facetingTolerance;}
  bool use_CAD() {return useCAD;}

  void set_overlap_thickness( double new_overlap_thickness );

  void set_numerical_precision( double new_precision );

  void set_use_CAD( bool use_cad ); 

  /* SECTION V: Metadata handling */
  ErrorCode detect_available_props( std::vector<std::string>& keywords_out );

  ErrorCode parse_properties( const std::vector<std::string>& keywords, 
                              const std::map<std::string,std::string>& synonyms = no_synonyms );

  ErrorCode prop_value( EntityHandle eh, const std::string& prop, std::string& value );

  ErrorCode prop_values( EntityHandle eh, const std::string& prop, 
                         std::vector< std::string >& value );

  bool has_prop( EntityHandle eh, const std::string& prop );

  ErrorCode get_all_prop_values( const std::string& prop, std::vector<std::string>& return_list );

  ErrorCode entities_by_property( const std::string& prop, std::vector<EntityHandle>& return_list,
                                  int dimension = 0, const std::string* value = NULL );

  bool is_implicit_complement(EntityHandle volume);

  Tag name_tag() {return nameTag;}

    // Get the tag used to associate OBB trees with geometry in load_file(..).
  Tag obb_tag() {return obbTag;}
  Tag geom_tag() {return geomTag;}
  Tag id_tag() {return idTag;}
  Tag sense_tag() { return senseTag; }

private:
  void tokenize( const std::string& str,
                 std::vector<std::string>& tokens,
                 const char* delimiters ) const;

  // a common type within the property and group name functions
  typedef std::map<std::string, std::string> prop_map;

  ErrorCode get_group_name( EntityHandle group_set, std::string& name );
  ErrorCode parse_group_name( EntityHandle group_set, prop_map& result );
  ErrorCode append_packed_string( Tag, EntityHandle, std::string& );
  ErrorCode unpack_packed_string( Tag tag, EntityHandle eh, 
                                  std::vector< std::string >& values );
  
  std::vector<EntityHandle>& surf_handles() {return entHandles[2];}
  std::vector<EntityHandle>& vol_handles() {return entHandles[3];}
  std::vector<EntityHandle>& group_handles() {return entHandles[4];}

  Tag get_tag( const char* name, int size, TagType store, DataType type,
                 const void* def_value = NULL, bool create_if_missing = true);
  
  /* SECTION VI: Other */
public:
  OrientedBoxTreeTool *obb_tree() {return &obbTree;}

  ErrorCode write_mesh(const char* ffile,
                       const int flen);

    // get the corners of the OBB for a given volume
  ErrorCode getobb(EntityHandle volume, double minPt[3], double maxPt[3]);

    // get the center point and three vectors for the OBB of a given volume
  ErrorCode getobb(EntityHandle volume, double center[3], 
                     double axis1[3], double axis2[3], double axis3[3]);

    // get the root of the obbtree for a given entity
  ErrorCode get_root(EntityHandle vol_or_surf, EntityHandle &root);

    // Get the instance of MOAB used by functions in this file.
  Interface* moab_instance() {return mbImpl;}

  
private:

  DagMC(Interface *mb_impl);
  
  static void create_instance(Interface *mb_impl = NULL);

  /* PRIVATE MEMBER DATA */

  static DagMC *instance_;
  Interface *mbImpl;

  OrientedBoxTreeTool obbTree;
  EntityHandle impl_compl_handle;
  Tag obbTag, geomTag, idTag, nameTag, senseTag, facetingTolTag;

  std::vector<EntityHandle> entHandles[5];
    // store some lists indexed by handle
    // this is the lowest-valued handle among entity sets representing
    // surfaces & volumes
  EntityHandle setOffset;
    // list of obbTree root sets for surfaces and volumes, 
    // indexed by [surf_or_vol_handle - setOffset]
  std::vector<EntityHandle> rootSets;
    // entity index (contiguous 1-N indices) indexed like rootSets are
  std::vector<int> entIndices;

    // corresponding geometric entities indexed like rootSets are
  std::vector<RefEntity *> geomEntities;
  
  // metadata
  // an empty synonym map to provide as a default argument to parse_properties()
  static const std::map<std::string,std::string> no_synonyms;
  // a map from the canonical property names to the tags representing them
  std::map<std::string, Tag> property_tagmap;

  char implComplName[NAME_TAG_SIZE];

  double overlapThickness;
  double numericalPrecision;
  double facetingTolerance, defaultFacetingTolerance;
  bool useCAD;         
  bool have_cgm_geom;  

  // temporary storage so functions don't have to reallocate vectors
  // for ray_fire:
  std::vector<double> distList;
  std::vector<EntityHandle> prevFacetList, surfList, facetList;
  // for point_in_volume:
  std::vector<double> disList;
  std::vector<int>    dirList;
  std::vector<EntityHandle> surList, facList;

  // for (optional) counting
  long long int n_pt_in_vol_calls, n_ray_fire_calls;

};

inline DagMC *DagMC::instance(Interface *mb_impl) 
{
  if (NULL == instance_) create_instance(mb_impl);
  
  return instance_;
}

inline EntityHandle DagMC::entity_by_index( int dimension, int index )
{
  assert(2 <= dimension && 3 >= dimension && (unsigned) index < entHandles[dimension].size());
  return entHandles[dimension][index];
}

inline int DagMC::index_by_handle( EntityHandle handle )
{
  assert(handle-setOffset < entIndices.size());
  return entIndices[handle-setOffset];
}

inline int DagMC::num_entities( int dimension )
{
  assert(0 <= dimension && 3 >= dimension);
  
  return entHandles[dimension].size() - 1;
}

    // get the root of the obbtree for a given entity
inline ErrorCode DagMC::get_root(EntityHandle vol_or_surf, EntityHandle &root) 
{
  unsigned int index = vol_or_surf - setOffset;
  root = (index < rootSets.size() ? rootSets[index] : 0);
  return (root ? MB_SUCCESS : MB_INDEX_OUT_OF_RANGE);
}

} // namespace moab

#endif