1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
#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"

#include "moab/Interface.hpp"
#include "moab/GeomTopoTool.hpp"
#include "moab/OrientedBoxTreeTool.hpp"

//! edge status
enum EdgeStatus {
  INSIDE,
  OUTSIDE,
  BOUNDARY
};

//! stores boundary element edge cut fraction
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]);
    }
  }
};

//! boundary cut-cell surface edge key
//! made for TechX
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;
  };
};

//! stores intersection distances
struct IntersectDist {
  double distance;
  int index;

  IntersectDist() {};

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

//! stores volume fractions
struct VolFrac {
  double vol;
  bool closed;
  double ePnt[6];

  VolFrac() {};

  VolFrac(double f, int c) {<--- Member variable 'VolFrac::ePnt' is not initialized in the constructor.
    vol = f;
    closed = c;
  };
};

struct LessThan
{
  bool operator() (const CutCellSurfEdgeKey key1, const CutCellSurfEdgeKey key2) const<--- Function parameter 'key1' should be passed by reference.<--- Function parameter 'key2' should be passed by reference.
  {
    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 EBMesher.hpp "meshkit/EBMesher.hpp"
 * \brief A meshing geometry as Cartesian structured mesh
 * It makes constructs tree structure with  triangles and
 * makes hexes bounding geometry
 * With ray-tracing, find intersections and determine element inside/outside/boundary.
 * Intersection fraction is stored to boundary elements.
 * Element inside/outside/boundary status are stored as tag.
 */
class EBMesher : public MeshScheme
{
public:

    //! Bare constructor
  EBMesher(MKCore *mkcore, const MEntVector &me_vec,
           double size = -1., bool use_geom = true,
           int add_layer = 3);
  
    //! Destructor
  virtual ~EBMesher();
  
  /**\brief Get class name */
  static const char* name() 
    { return "EBMesher"; }

  /**\brief Function returning whether this scheme can mesh entities of t
   *        the specified dimension.
   *\param dim entity dimension
   */
  static bool can_mesh(iBase_EntityType dim)
    { return iBase_REGION == dim; }
   
  /** \brief Function returning whether this scheme can mesh the specified entity
   * 
   * Used by MeshOpFactory to find scheme for an entity.
   * \param me ModelEnt being queried
   * \return If true, this scheme can mesh the specified ModelEnt
   */
  static bool can_mesh(ModelEnt *me);
  
  /**\brief Get list of mesh entity types that can be generated.
   *\return array terminated with \c moab::MBMAXTYPE
   */
  static const moab::EntityType* output_types();
  
  /** \brief Return the mesh entity types operated on by this scheme
   * \return array terminated with \c moab::MBMAXTYPE
   */
  virtual const moab::EntityType* mesh_types_arr() const
  { return output_types(); }
  
  
  virtual bool add_modelent(ModelEnt *model_ent);
  
  /** \brief set # of divisions (x/y/z directions) of Cartesian box to use SCDMesh output
   * \param min Cartesian box min coordinates
   * \param max Cartesian box max coordinates
   * \return int if is working correctly
   */
  int set_division(double* min, double* max);
  
  /**\brief Setup is a no-op, but must be provided since it's pure virtual
   */  
  virtual void setup_this();
  
  /**\ The only setup/execute function we need, since meshing vertices is trivial
   */
  virtual void execute_this();
  
  /** \brief query function for techX
   * \param boxMin Cartesian box min coordinates returned
   * \param boxMax Cartesian box max coordinates returned
   * \param nDiv nDiv(x/yz directions) returned
   * \param rmdCutCellSurfEdge map of cut-cell surface index and edge cut information returned
   * \param rvnInsideCell Inside elements returned
   * \param isCornerExterior if box corner is exterior returned
   */
  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);
  
  /** \brief query function to get multiple cut-cell edges
   * \param boxMin Cartesian box min coordinates returned
   * \param boxMax Cartesian box max coordinates returned
   * \param nDiv nDiv(x/yz directions) returned
   * \param rmdCutCellEdge map of cut-cell surface index and edge cut information returned
   * \param rvnInsideCell Inside elements returned
   * \param isCornerExterior if box corner is exterior returned
   * \return if this function is working correctly
   */
  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);
  
  /** \brief get volume fraction for each material
   * \param vol_frac_div resolution to get volume fraction
   * \return if this function is working correctly
   */
  bool get_volume_fraction(int vol_frac_div);
  
  /** \brief export mesh to file
   * \param file_name export file name
   * \param separate export file separately(inside/outside/boundary)
   */
  void export_mesh(const char* file_name, bool separate = false);

  /** \brief set number of intervals
   * \param n_interval number of interval array for x/y/z directions
   */
  void set_num_interval(int* n_interval);

  /** \brief set number of intervals
   * \param box_increase number of interval array for x/y/z directions
   */
  void increase_box(double box_increase);

  /** \brief set if mesh for whole geometry at once or individually
   * \param use if mesh for whole geometry at once
   */
  void use_whole_geom(int use);

  /** \brief set if mesh based geometry is used
   * \param use if mesh based geometry is used
   */
  void use_mesh_geometry(bool use);

  /** \brief construct obb tree with faceted geometry and set the box max/min coordinates
   */
  void set_obb_tree_box_dimension();
  
protected:
  
private:

  //! No copy constructor, since there's only meant to be one of these
  EBMesher(const EBMesher &);
  
  MeshOp* get_scd_mesher();

  //! No operator=, since there's only meant to be one of these
  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];
  
  /** \brief get hex edge status (inside/outside/boundary)
   * \param dZ edge end coordinate
   * \param iSkip how many index skipped for next intersection checking
   * \return EdgeStatus edge status
   */
  EdgeStatus get_edge_status(const double dZ, int& iSkip);

  /** \brief set hex status for neighboring elements
   * \param dir current fired ray direction
   * \param i index for i direction
   * \param j index for j direction
   * \param k index for k direction
   * \return bool if is working correctly
   */
  bool set_neighbor_hex_status(int dir, int i, int j, int k);

  /** \ brief set hex status by edge status
   * \param index hex index in m_vnHexStatus vector
   * \param value edge status
   * \param dir ray direction
   * \return bool if is working correctly
   */
  bool set_hex_status(int index, EdgeStatus value, int dir);

  /** \brief set edge status
   * \param dir ray direction
   * \return bool if is working correctly
   */
  bool set_edge_status(int dir);

  /** \brief set all produced mesh information as tag
   * \ hex status, edge-cut information.....
   */
  void set_tag_info();

  /** \brief wirte mesh
   * \param file_name
   * \param type element type (inside:0, outside:1, boundary:2)
   * \param handles element handles
   * \param n_elem # of elements
   * \return int if is working correctly
   */
  void write_mesh(const char* file_name, int type,
                  iBase_EntityHandle* handles, int& n_elem);

  /** \brief get edge fraction information
   * \param idHex index in m_mdCutFraction
   * \param dir ray direction
   * \return double edge fraction
   */
  double get_edge_fraction(int idHex, int dir);

  /** \brief get if the edge is fully inside(returns 1) or outside(returns 0)
   * \param i index for i direction
   * \param j index for j direction
   * \param k index for k direction
   * \param dir ray direction
   * \return double edge fraction
   */
  double get_uncut_edge_fraction(int i, int j, int k, int dir);

  /** \brief check if the intersected surface geometry is shared or overlapped
   * \param index intersection index in m_vIntersection vector
   * \return bool if is working correctly
   */
  bool is_shared_overlapped_surf(int index);

  /** \brief move intersection pairs to check element status
   * \param n_dir ray direction
   * \param n_inter index of intersection points
   * \param start_pnt ray starting point
   * \return bool if is working correctly
   */
  bool move_intersections(int n_dir, int n_inter, double start_pnt[3]);

  /** \brief get inside status elements
   * \param rvnInsideCell cell indices (i,j,k triple)
   * \return bool if is working correctly
   */
  bool get_inside_hex(std::vector<int>& rvnInsideCell);

  /** \brief check if ray is passing shared vertices or edges
   * \ And, check if the passing surface is overlapped one
   * \param bMoveOnce if ray is already moved before
   * \return bool if is working correctly
   */
  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 variable, used in registration
  static int init;

  /** \brief get MOAB instance
   * \return MOAB instance
   */
  moab::Interface* moab_instance() {return mk_core()->moab_instance();}

  /** \brief get MOAB's Tag
   * \param name Tag name
   * \param size Tag size
   * \param store Tag type
   * \param type data type stored
   * \param def_value default value
   * \param create_if_missing create Tag if it is missed (flag)
   * \return Tag handle
   */
  iBase_TagHandle get_tag(const char* name, int size, unsigned flags, moab::DataType type,
                          const void* def_value = NULL, bool create_if_missing = true);
  
  /** \brief get MOAB's various length Tag
   * \param name Tag name
   * \param store Tag type
   * \param type data type stored
   * \return Tag handle
   */
  iBase_TagHandle get_various_length_tag(const char* name,
                                         unsigned store, moab::DataType type);

  /** \brief construct hexes in MOAB structured mesh format
   * \return int if is working correctly
   */
  int make_scd_hexes();

  /** \brief construct hexes in MOAB unstructured mesh format
   * \return int if is working correctly
   */
  int make_uscd_hexes();

  /** \brief construct OBB tree
   * \input geometry is faceted which are constructed to OBB tree
   * \return int if is working correctly
   */
  int construct_obb_tree();

  /** \brief construct and set OBB tree
   * \input Model entity pointer for tree construction
   */
  void set_tree_root(ModelEnt* me);

  /** \brief set # of divisions
   * \calcalate the best ones from the size of faceted triangles
   */
  void set_division();

  /** \brief find intersections by firing rays
   */
  void find_intersections();

  /** \brief fires rays to 3 directions
   * \param dir give the ray direction
   * \it calls fire_ray
   */
  void fire_rays(int dir);

  /** \brief fires ray
   * \param nIntersect # of intersections returned
   * \param startPnt ray starting point
   * \param endPnt ray ending point
   * \param tol tolerance to find intersection
   * \param dir ray direction
   * \param rayLength ray length
   * \return bool if is working correctly
   */
  bool fire_ray(int& nIntersect, double startPnt[3],
                double endPnt[3], double tol, int dir,
                double rayLength);

  /** \brief moves ray which passes any singluar point
   * \param nIntersect # of intersections returned
   * \param startPnt ray starting point
   * \param endPnt ray ending point
   * \param tol tolerance to find intersection
   * \param dir ray direction
   * \param bMoveOnce if ray is already moved before
   * \return bool if is working correctly
   */
  bool move_ray(int& nIntersect, double* startPnt, double* endPnt,
                double tol, int dir, bool bMoveOnce);

  /** \brief if the facet has the same direction to the ray
   * \param i index in m_vIntersection (intersection vector)
   * \param dir ray direction
   * \return bool if is working correctly
   */
  bool is_same_direct_to_ray(int i, int dir);

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

  // ! Tree root
  moab::EntityHandle m_hTreeRoot;

  // ! OBB tree tool instance
  moab::OrientedBoxTreeTool* m_hObbTree;

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

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

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

  std::map<moab::EntityHandle, moab::EntityHandle>  m_mRootSets;

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

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