EBMesher.hpp

Go to the documentation of this file.

#ifndef MESHKIT_EB_MESHER_HPP
#define MESHKIT_EB_MESHER_HPP

#include <vector>
#include <map>
#include <sys/resource.h>

#include "meshkit/Types.hpp"
#include "meshkit/Error.hpp"
#include "meshkit/MeshScheme.hpp"
#include "meshkit/ModelEnt.hpp"
#include "meshkit/iMesh.hpp"

#ifdef HAVE_MOAB
#include "moab/Interface.hpp"
#include "moab/GeomTopoTool.hpp"
#include "MBOrientedBoxTreeTool.hpp"
#endif

enum EdgeStatus {
  INSIDE,
  OUTSIDE,
  BOUNDARY
};

struct CutFraction {
  std::vector<double> fraction[3];

  CutFraction() {};

  CutFraction(int dir, const std::vector<double>& val) {
    add(dir, val);
  };
  
  void add(int dir, const std::vector<double>& val) {
    for (unsigned int i = 0; i < val.size(); i++) {
      fraction[dir].push_back(val[i]);
    }
  }
};

struct CutCellSurfEdgeKey {
  int i, j, k, l;

  CutCellSurfEdgeKey() {
    i = j = k = l = 0;
  };

  CutCellSurfEdgeKey(int ii, int jj, int kk, int ll) {
    i = ii;
    j = jj;
    k = kk;
    l = ll;
  };
};

struct IntersectDist {
  double distance;
  int index;

  IntersectDist() {};

  IntersectDist(double d, int i) {
    distance = d;
    index = i;
  };
};

struct VolFrac {
  double vol;
  bool closed;
  double ePnt[6];

  VolFrac() {};

  VolFrac(double f, int c) {
    vol = f;
    closed = c;
  };
};

struct LessThan
{
  bool operator() (const CutCellSurfEdgeKey key1, const CutCellSurfEdgeKey key2) const
  {
    if (key1.i < key2.i) return true;
    else if (key1.i > key2.i) return false;
    else {
      if (key1.j < key2.j) return true;
      else if (key1.j > key2.j) return false;
      else {
        if (key1.k < key2.k) return true;
        else if (key1.k > key2.k) return false;
        else {
          if (key1.l < key2.l) return true;
          else return false;
        }
      }
    }
  };
};

namespace MeshKit {

class MKCore;
    

class EBMesher : public MeshScheme
{
public:

  EBMesher(MKCore *mkcore, const MEntVector &me_vec,
           double size = -1., bool use_geom = true,
           int add_layer = 3);
  
  virtual ~EBMesher();
  
  static const char* name() 
    { return "EBMesher"; }

  static bool can_mesh(iBase_EntityType dim)
    { return iBase_REGION == dim; }
   
  static bool can_mesh(ModelEnt *me);
  
  static const moab::EntityType* output_types();
  
  virtual const moab::EntityType* mesh_types_arr() const
  { return output_types(); }
  
  
  virtual bool add_modelent(ModelEnt *model_ent);
  
  int set_division(double* min, double* max);
  
  virtual void setup_this();
  
  virtual void execute_this();
  
  void get_grid_and_edges_techX(double* boxMin, double* boxMax, int* nDiv,
                                std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& rmdCutCellSurfEdge,
                                std::vector<int>& rvnInsideCell, bool isCornerExterior = true);
  
  bool get_grid_and_edges(double* boxMin, double* boxMax, int* nDiv,
                          std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& rmdCutCellEdge,
                          std::vector<int>& rvnInsideCell, bool isCornerExterior = true);
  
  bool get_volume_fraction(int vol_frac_div);
  
  void export_mesh(const char* file_name, bool separate = false);

  void set_num_interval(int* n_interval);

  void increase_box(double box_increase);

  void use_whole_geom(int use);

  void use_mesh_geometry(bool use);

  void set_obb_tree_box_dimension();
  
protected:
  
private:

  EBMesher(const EBMesher &);
  
  MeshOp* get_scd_mesher();

  EBMesher &operator=(const EBMesher &);

  iBase_TagHandle m_elemStatusTag, m_edgeCutFracLengthTag,
    m_edgeCutFracTag, m_volFracLengthTag, m_volFracHandleTag, m_volFracTag, m_bbTag;
  iMesh_Instance m_mesh;
  iBase_EntitySetHandle m_hRootSet;
  std::vector<IntersectDist> m_vIntersection;
  int m_nTri, m_nHex, m_iInter, m_nFraction, m_iStartHex, m_nMove, m_nAddLayer,
    m_nIteration, m_iOverlap, m_iElem, m_nNode[3], m_nDiv[3],
    m_iFirstP, m_iSecondP;
  double m_dFirstP, m_dSecondP, m_curPnt, m_prevPnt, m_boxIncrease,
    m_dIntervalSize[3], m_origin_coords[3], m_dInputSize,
    m_min[3], m_max[3];
  EdgeStatus m_nStatus;
  bool m_bMovedOnce, m_bUseGeom, m_bUseWholeGeom, m_bUseMeshGeom;
  std::vector<iBase_EntityHandle> m_vhVertex;
  std::vector<int> m_vnHexStatus;
  std::map<int, CutFraction> m_mdCutFraction;
  std::vector<EdgeStatus> m_vnEdgeStatus[3];
  
  EdgeStatus get_edge_status(const double dZ, int& iSkip);

  bool set_neighbor_hex_status(int dir, int i, int j, int k);

  bool set_hex_status(int index, EdgeStatus value, int dir);

  bool set_edge_status(int dir);

  void set_tag_info();

  void write_mesh(const char* file_name, int type,
                  iBase_EntityHandle* handles, int& n_elem);

  double get_edge_fraction(int idHex, int dir);

  double get_uncut_edge_fraction(int i, int j, int k, int dir);

  bool is_shared_overlapped_surf(int index);

  bool move_intersections(int n_dir, int n_inter, double start_pnt[3]);

  bool get_inside_hex(std::vector<int>& rvnInsideCell);

  bool is_ray_move_and_set_overlap_surf(bool& bMoveOnce);

  void remove_intersection_duplicates();

  // test function 1 for debugging
  bool export_fraction_edges(std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& rmdCutCellSurfEdge);

  // test functions 2 for debugging
  bool export_fraction_points(std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& mdCutCellEdge);

  // test functions 3 for debugging
  bool make_edge(double ePnt[6], std::vector<iBase_EntityHandle>& edge_handles);

  static int init;

#ifdef HAVE_MOAB

  MBInterface* moab_instance() {return mk_core()->moab_instance();}

  iBase_TagHandle get_tag(const char* name, int size, unsigned flags, MBDataType type,
                          const void* def_value = NULL, bool create_if_missing = true);
  
  iBase_TagHandle get_various_length_tag(const char* name,
                                         unsigned store, MBDataType type);

  int make_scd_hexes();

  int make_uscd_hexes();

  int construct_obb_tree();

  void set_tree_root(ModelEnt* me);

  void set_division();

  void find_intersections();

  void fire_rays(int dir);

  bool fire_ray(int& nIntersect, double startPnt[3],
                double endPnt[3], double tol, int dir,
                double rayLength);

  bool move_ray(int& nIntersect, double* startPnt, double* endPnt,
                double tol, int dir, bool bMoveOnce);

  bool is_same_direct_to_ray(int i, int dir);

  // ! GeomTopoTool instance
  moab::GeomTopoTool* m_GeomTopoTool;

  // ! Tree root
  MBEntityHandle m_hTreeRoot;

  // ! OBB tree tool instance
  MBOrientedBoxTreeTool* m_hObbTree;

  // ! intersected surface geometry list
  std::vector<MBEntityHandle> m_vhInterSurf;

  // ! intersected facet list
  std::vector<MBEntityHandle> m_vhInterFacet;

  // ! overlapped surface list
  std::map<MBEntityHandle, int> m_mhOverlappedSurf;

  std::map<MBEntityHandle, MBEntityHandle>  m_mRootSets;

  double m_minCoord[3], m_maxCoord[3];
#endif
};

inline void EBMesher::increase_box(double box_increase)
{
  m_boxIncrease = box_increase;
}

inline void EBMesher::use_whole_geom(int use)
{
  m_bUseWholeGeom = use;
}

inline void EBMesher::use_mesh_geometry(bool use)
{
  m_bUseMeshGeom = use;
}
}

#endif