MeshKit  1.0
EBMesher.hpp
Go to the documentation of this file.
00001 #ifndef MESHKIT_EB_MESHER_HPP
00002 #define MESHKIT_EB_MESHER_HPP
00003 
00004 #include <vector>
00005 #include <map>
00006 
00007 #include "meshkit/Types.hpp"
00008 #include "meshkit/Error.hpp"
00009 #include "meshkit/MeshScheme.hpp"
00010 #include "meshkit/ModelEnt.hpp"
00011 #include "meshkit/iMesh.hpp"
00012 
00013 #include "moab/Interface.hpp"
00014 #include "moab/GeomTopoTool.hpp"
00015 #include "moab/OrientedBoxTreeTool.hpp"
00016 
00018 enum EdgeStatus {
00019   INSIDE,
00020   OUTSIDE,
00021   BOUNDARY
00022 };
00023 
00025 struct CutFraction {
00026   std::vector<double> fraction[3];
00027 
00028   CutFraction() {};
00029 
00030   CutFraction(int dir, const std::vector<double>& val) {
00031     add(dir, val);
00032   };
00033   
00034   void add(int dir, const std::vector<double>& val) {
00035     for (unsigned int i = 0; i < val.size(); i++) {
00036       fraction[dir].push_back(val[i]);
00037     }
00038   }
00039 };
00040 
00043 struct CutCellSurfEdgeKey {
00044   int i, j, k, l;
00045 
00046   CutCellSurfEdgeKey() {
00047     i = j = k = l = 0;
00048   };
00049 
00050   CutCellSurfEdgeKey(int ii, int jj, int kk, int ll) {
00051     i = ii;
00052     j = jj;
00053     k = kk;
00054     l = ll;
00055   };
00056 };
00057 
00059 struct IntersectDist {
00060   double distance;
00061   int index;
00062 
00063   IntersectDist() {};
00064 
00065   IntersectDist(double d, int i) {
00066     distance = d;
00067     index = i;
00068   };
00069 };
00070 
00072 struct VolFrac {
00073   double vol;
00074   bool closed;
00075   double ePnt[6];
00076 
00077   VolFrac() {};
00078 
00079   VolFrac(double f, int c) {
00080     vol = f;
00081     closed = c;
00082   };
00083 };
00084 
00085 struct LessThan
00086 {
00087   bool operator() (const CutCellSurfEdgeKey key1, const CutCellSurfEdgeKey key2) const
00088   {
00089     if (key1.i < key2.i) return true;
00090     else if (key1.i > key2.i) return false;
00091     else {
00092       if (key1.j < key2.j) return true;
00093       else if (key1.j > key2.j) return false;
00094       else {
00095         if (key1.k < key2.k) return true;
00096         else if (key1.k > key2.k) return false;
00097         else {
00098           if (key1.l < key2.l) return true;
00099           else return false;
00100         }
00101       }
00102     }
00103   };
00104 };
00105 
00106 namespace MeshKit {
00107 
00108 class MKCore;
00109     
00110 
00119 class EBMesher : public MeshScheme
00120 {
00121 public:
00122 
00124   EBMesher(MKCore *mkcore, const MEntVector &me_vec,
00125            double size = -1., bool use_geom = true,
00126            int add_layer = 3);
00127   
00129   virtual ~EBMesher();
00130   
00132   static const char* name() 
00133     { return "EBMesher"; }
00134 
00139   static bool can_mesh(iBase_EntityType dim)
00140     { return iBase_REGION == dim; }
00141    
00148   static bool can_mesh(ModelEnt *me);
00149   
00153   static const moab::EntityType* output_types();
00154   
00158   virtual const moab::EntityType* mesh_types_arr() const
00159   { return output_types(); }
00160   
00161   
00162   virtual bool add_modelent(ModelEnt *model_ent);
00163   
00169   int set_division(double* min, double* max);
00170   
00173   virtual void setup_this();
00174   
00177   virtual void execute_this();
00178   
00187   void get_grid_and_edges_techX(double* boxMin, double* boxMax, int* nDiv,
00188                                 std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& rmdCutCellSurfEdge,
00189                                 std::vector<int>& rvnInsideCell, bool isCornerExterior = true);
00190   
00200   bool get_grid_and_edges(double* boxMin, double* boxMax, int* nDiv,
00201                           std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& rmdCutCellEdge,
00202                           std::vector<int>& rvnInsideCell, bool isCornerExterior = true);
00203   
00208   bool get_volume_fraction(int vol_frac_div);
00209   
00214   void export_mesh(const char* file_name, bool separate = false);
00215 
00219   void set_num_interval(int* n_interval);
00220 
00224   void increase_box(double box_increase);
00225 
00229   void use_whole_geom(int use);
00230 
00234   void use_mesh_geometry(bool use);
00235 
00238   void set_obb_tree_box_dimension();
00239   
00240 protected:
00241   
00242 private:
00243 
00245   EBMesher(const EBMesher &);
00246   
00247   MeshOp* get_scd_mesher();
00248 
00250   EBMesher &operator=(const EBMesher &);
00251 
00252   iBase_TagHandle m_elemStatusTag, m_edgeCutFracLengthTag,
00253     m_edgeCutFracTag, m_volFracLengthTag, m_volFracHandleTag, m_volFracTag, m_bbTag;
00254   iMesh_Instance m_mesh;
00255   iBase_EntitySetHandle m_hRootSet;
00256   std::vector<IntersectDist> m_vIntersection;
00257   int m_nTri, m_nHex, m_iInter, m_nFraction, m_iStartHex, m_nMove, m_nAddLayer,
00258     m_nIteration, m_iOverlap, m_iElem, m_nNode[3], m_nDiv[3],
00259     m_iFirstP, m_iSecondP;
00260   double m_dFirstP, m_dSecondP, m_curPnt, m_prevPnt, m_boxIncrease,
00261     m_dIntervalSize[3], m_origin_coords[3], m_dInputSize,
00262     m_min[3], m_max[3];
00263   EdgeStatus m_nStatus;
00264   bool m_bMovedOnce, m_bUseGeom, m_bUseWholeGeom, m_bUseMeshGeom;
00265   std::vector<iBase_EntityHandle> m_vhVertex;
00266   std::vector<int> m_vnHexStatus;
00267   std::map<int, CutFraction> m_mdCutFraction;
00268   std::vector<EdgeStatus> m_vnEdgeStatus[3];
00269   
00275   EdgeStatus get_edge_status(const double dZ, int& iSkip);
00276 
00284   bool set_neighbor_hex_status(int dir, int i, int j, int k);
00285 
00292   bool set_hex_status(int index, EdgeStatus value, int dir);
00293 
00298   bool set_edge_status(int dir);
00299 
00303   void set_tag_info();
00304 
00312   void write_mesh(const char* file_name, int type,
00313                   iBase_EntityHandle* handles, int& n_elem);
00314 
00320   double get_edge_fraction(int idHex, int dir);
00321 
00329   double get_uncut_edge_fraction(int i, int j, int k, int dir);
00330 
00335   bool is_shared_overlapped_surf(int index);
00336 
00343   bool move_intersections(int n_dir, int n_inter, double start_pnt[3]);
00344 
00349   bool get_inside_hex(std::vector<int>& rvnInsideCell);
00350 
00356   bool is_ray_move_and_set_overlap_surf(bool& bMoveOnce);
00357 
00358   void remove_intersection_duplicates();
00359 
00360   // test function 1 for debugging
00361   bool export_fraction_edges(std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& rmdCutCellSurfEdge);
00362 
00363   // test functions 2 for debugging
00364   bool export_fraction_points(std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& mdCutCellEdge);
00365 
00366   // test functions 3 for debugging
00367   bool make_edge(double ePnt[6], std::vector<iBase_EntityHandle>& edge_handles);
00368 
00370   static int init;
00371 
00375   moab::Interface* moab_instance() {return mk_core()->moab_instance();}
00376 
00386   iBase_TagHandle get_tag(const char* name, int size, unsigned flags, moab::DataType type,
00387                           const void* def_value = NULL, bool create_if_missing = true);
00388   
00395   iBase_TagHandle get_various_length_tag(const char* name,
00396                                          unsigned store, moab::DataType type);
00397 
00401   int make_scd_hexes();
00402 
00406   int make_uscd_hexes();
00407 
00412   int construct_obb_tree();
00413 
00417   void set_tree_root(ModelEnt* me);
00418 
00422   void set_division();
00423 
00426   void find_intersections();
00427 
00432   void fire_rays(int dir);
00433 
00443   bool fire_ray(int& nIntersect, double startPnt[3],
00444                 double endPnt[3], double tol, int dir,
00445                 double rayLength);
00446 
00456   bool move_ray(int& nIntersect, double* startPnt, double* endPnt,
00457                 double tol, int dir, bool bMoveOnce);
00458 
00464   bool is_same_direct_to_ray(int i, int dir);
00465 
00466   // ! GeomTopoTool instance
00467   moab::GeomTopoTool* m_GeomTopoTool;
00468 
00469   // ! Tree root
00470   moab::EntityHandle m_hTreeRoot;
00471 
00472   // ! OBB tree tool instance
00473   moab::OrientedBoxTreeTool* m_hObbTree;
00474 
00475   // ! intersected surface geometry list
00476   std::vector<moab::EntityHandle> m_vhInterSurf;
00477 
00478   // ! intersected facet list
00479   std::vector<moab::EntityHandle> m_vhInterFacet;
00480 
00481   // ! overlapped surface list
00482   std::map<moab::EntityHandle, int> m_mhOverlappedSurf;
00483 
00484   std::map<moab::EntityHandle, moab::EntityHandle>  m_mRootSets;
00485 
00486   double m_minCoord[3], m_maxCoord[3];
00487 };
00488 
00489 inline void EBMesher::increase_box(double box_increase)
00490 {
00491   m_boxIncrease = box_increase;
00492 }
00493 
00494 inline void EBMesher::use_whole_geom(int use)
00495 {
00496   m_bUseWholeGeom = use;
00497 }
00498 
00499 inline void EBMesher::use_mesh_geometry(bool use)
00500 {
00501   m_bUseMeshGeom = use;
00502 }
00503 }
00504 
00505 #endif
00506 
00507   
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines