Branch data Line data Source code
1 : : #ifndef MESHKIT_EB_MESHER_HPP
2 : : #define MESHKIT_EB_MESHER_HPP
3 : :
4 : : #include <vector>
5 : : #include <map>
6 : :
7 : : #include "meshkit/Types.hpp"
8 : : #include "meshkit/Error.hpp"
9 : : #include "meshkit/MeshScheme.hpp"
10 : : #include "meshkit/ModelEnt.hpp"
11 : : #include "meshkit/iMesh.hpp"
12 : :
13 : : #include "moab/Interface.hpp"
14 : : #include "moab/GeomTopoTool.hpp"
15 : : #include "moab/OrientedBoxTreeTool.hpp"
16 : :
17 : : //! edge status
18 : : enum EdgeStatus {
19 : : INSIDE,
20 : : OUTSIDE,
21 : : BOUNDARY
22 : : };
23 : :
24 : : //! stores boundary element edge cut fraction
25 [ + - ][ + + ]: 9112 : struct CutFraction {
[ # # # # ]
[ + + ][ + - ]
[ + + ]
26 : : std::vector<double> fraction[3];
27 : :
28 [ + - ][ + + ]: 1072 : CutFraction() {};
[ # # # # ]
29 : :
30 [ + - ][ + + : 1340 : CutFraction(int dir, const std::vector<double>& val) {
# # # # #
# ]
31 [ + - ]: 268 : add(dir, val);
32 [ # # ]: 268 : };
33 : :
34 : 366 : void add(int dir, const std::vector<double>& val) {
35 [ + + ]: 732 : for (unsigned int i = 0; i < val.size(); i++) {
36 : 366 : fraction[dir].push_back(val[i]);
37 : : }
38 : 366 : }
39 : : };
40 : :
41 : : //! boundary cut-cell surface edge key
42 : : //! made for TechX
43 : : struct CutCellSurfEdgeKey {
44 : : int i, j, k, l;
45 : :
46 : : CutCellSurfEdgeKey() {
47 : : i = j = k = l = 0;
48 : : };
49 : :
50 : 1489 : CutCellSurfEdgeKey(int ii, int jj, int kk, int ll) {
51 : 1489 : i = ii;
52 : 1489 : j = jj;
53 : 1489 : k = kk;
54 : 1489 : l = ll;
55 : 1489 : };
56 : : };
57 : :
58 : : //! stores intersection distances
59 : : struct IntersectDist {
60 : : double distance;
61 : : int index;
62 : :
63 : 523 : IntersectDist() {};
64 : :
65 : 523 : IntersectDist(double d, int i) {
66 : 523 : distance = d;
67 : 523 : index = i;
68 : 523 : };
69 : : };
70 : :
71 : : //! stores volume fractions
72 : : struct VolFrac {
73 : : double vol;
74 : : bool closed;
75 : : double ePnt[6];
76 : :
77 : 0 : VolFrac() {};
78 : :
79 : 0 : VolFrac(double f, int c) {
80 : 0 : vol = f;
81 : 0 : closed = c;
82 : 0 : };
83 : : };
84 : :
85 : : struct LessThan
86 : : {
87 : 19495 : bool operator() (const CutCellSurfEdgeKey key1, const CutCellSurfEdgeKey key2) const
88 : : {
89 [ + + ]: 19495 : if (key1.i < key2.i) return true;
90 [ + + ]: 15747 : else if (key1.i > key2.i) return false;
91 : : else {
92 [ + + ]: 12829 : if (key1.j < key2.j) return true;
93 [ + + ]: 8780 : else if (key1.j > key2.j) return false;
94 : : else {
95 [ + + ]: 6394 : if (key1.k < key2.k) return true;
96 [ + + ]: 3190 : else if (key1.k > key2.k) return false;
97 : : else {
98 [ + + ]: 2931 : if (key1.l < key2.l) return true;
99 : 726 : else return false;
100 : : }
101 : : }
102 : : }
103 : : };
104 : : };
105 : :
106 : : namespace MeshKit {
107 : :
108 : : class MKCore;
109 : :
110 : :
111 : : /** \class EBMesher EBMesher.hpp "meshkit/EBMesher.hpp"
112 : : * \brief A meshing geometry as Cartesian structured mesh
113 : : * It makes constructs tree structure with triangles and
114 : : * makes hexes bounding geometry
115 : : * With ray-tracing, find intersections and determine element inside/outside/boundary.
116 : : * Intersection fraction is stored to boundary elements.
117 : : * Element inside/outside/boundary status are stored as tag.
118 : : */
119 : : class EBMesher : public MeshScheme
120 : : {
121 : : public:
122 : :
123 : : //! Bare constructor
124 : : EBMesher(MKCore *mkcore, const MEntVector &me_vec,
125 : : double size = -1., bool use_geom = true,
126 : : int add_layer = 3);
127 : :
128 : : //! Destructor
129 : : virtual ~EBMesher();
130 : :
131 : : /**\brief Get class name */
132 : 701 : static const char* name()
133 : 701 : { return "EBMesher"; }
134 : :
135 : : /**\brief Function returning whether this scheme can mesh entities of t
136 : : * the specified dimension.
137 : : *\param dim entity dimension
138 : : */
139 : 160 : static bool can_mesh(iBase_EntityType dim)
140 : 160 : { return iBase_REGION == dim; }
141 : :
142 : : /** \brief Function returning whether this scheme can mesh the specified entity
143 : : *
144 : : * Used by MeshOpFactory to find scheme for an entity.
145 : : * \param me ModelEnt being queried
146 : : * \return If true, this scheme can mesh the specified ModelEnt
147 : : */
148 : : static bool can_mesh(ModelEnt *me);
149 : :
150 : : /**\brief Get list of mesh entity types that can be generated.
151 : : *\return array terminated with \c moab::MBMAXTYPE
152 : : */
153 : : static const moab::EntityType* output_types();
154 : :
155 : : /** \brief Return the mesh entity types operated on by this scheme
156 : : * \return array terminated with \c moab::MBMAXTYPE
157 : : */
158 : 0 : virtual const moab::EntityType* mesh_types_arr() const
159 : 0 : { return output_types(); }
160 : :
161 : :
162 : : virtual bool add_modelent(ModelEnt *model_ent);
163 : :
164 : : /** \brief set # of divisions (x/y/z directions) of Cartesian box to use SCDMesh output
165 : : * \param min Cartesian box min coordinates
166 : : * \param max Cartesian box max coordinates
167 : : * \return int if is working correctly
168 : : */
169 : : int set_division(double* min, double* max);
170 : :
171 : : /**\brief Setup is a no-op, but must be provided since it's pure virtual
172 : : */
173 : : virtual void setup_this();
174 : :
175 : : /**\ The only setup/execute function we need, since meshing vertices is trivial
176 : : */
177 : : virtual void execute_this();
178 : :
179 : : /** \brief query function for techX
180 : : * \param boxMin Cartesian box min coordinates returned
181 : : * \param boxMax Cartesian box max coordinates returned
182 : : * \param nDiv nDiv(x/yz directions) returned
183 : : * \param rmdCutCellSurfEdge map of cut-cell surface index and edge cut information returned
184 : : * \param rvnInsideCell Inside elements returned
185 : : * \param isCornerExterior if box corner is exterior returned
186 : : */
187 : : void get_grid_and_edges_techX(double* boxMin, double* boxMax, int* nDiv,
188 : : std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& rmdCutCellSurfEdge,
189 : : std::vector<int>& rvnInsideCell, bool isCornerExterior = true);
190 : :
191 : : /** \brief query function to get multiple cut-cell edges
192 : : * \param boxMin Cartesian box min coordinates returned
193 : : * \param boxMax Cartesian box max coordinates returned
194 : : * \param nDiv nDiv(x/yz directions) returned
195 : : * \param rmdCutCellEdge map of cut-cell surface index and edge cut information returned
196 : : * \param rvnInsideCell Inside elements returned
197 : : * \param isCornerExterior if box corner is exterior returned
198 : : * \return if this function is working correctly
199 : : */
200 : : bool get_grid_and_edges(double* boxMin, double* boxMax, int* nDiv,
201 : : std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& rmdCutCellEdge,
202 : : std::vector<int>& rvnInsideCell, bool isCornerExterior = true);
203 : :
204 : : /** \brief get volume fraction for each material
205 : : * \param vol_frac_div resolution to get volume fraction
206 : : * \return if this function is working correctly
207 : : */
208 : : bool get_volume_fraction(int vol_frac_div);
209 : :
210 : : /** \brief export mesh to file
211 : : * \param file_name export file name
212 : : * \param separate export file separately(inside/outside/boundary)
213 : : */
214 : : void export_mesh(const char* file_name, bool separate = false);
215 : :
216 : : /** \brief set number of intervals
217 : : * \param n_interval number of interval array for x/y/z directions
218 : : */
219 : : void set_num_interval(int* n_interval);
220 : :
221 : : /** \brief set number of intervals
222 : : * \param box_increase number of interval array for x/y/z directions
223 : : */
224 : : void increase_box(double box_increase);
225 : :
226 : : /** \brief set if mesh for whole geometry at once or individually
227 : : * \param use if mesh for whole geometry at once
228 : : */
229 : : void use_whole_geom(int use);
230 : :
231 : : /** \brief set if mesh based geometry is used
232 : : * \param use if mesh based geometry is used
233 : : */
234 : : void use_mesh_geometry(bool use);
235 : :
236 : : /** \brief construct obb tree with faceted geometry and set the box max/min coordinates
237 : : */
238 : : void set_obb_tree_box_dimension();
239 : :
240 : : protected:
241 : :
242 : : private:
243 : :
244 : : //! No copy constructor, since there's only meant to be one of these
245 : : EBMesher(const EBMesher &);
246 : :
247 : : MeshOp* get_scd_mesher();
248 : :
249 : : //! No operator=, since there's only meant to be one of these
250 : : EBMesher &operator=(const EBMesher &);
251 : :
252 : : iBase_TagHandle m_elemStatusTag, m_edgeCutFracLengthTag,
253 : : m_edgeCutFracTag, m_volFracLengthTag, m_volFracHandleTag, m_volFracTag, m_bbTag;
254 : : iMesh_Instance m_mesh;
255 : : iBase_EntitySetHandle m_hRootSet;
256 : : std::vector<IntersectDist> m_vIntersection;
257 : : int m_nTri, m_nHex, m_iInter, m_nFraction, m_iStartHex, m_nMove, m_nAddLayer,
258 : : m_nIteration, m_iOverlap, m_iElem, m_nNode[3], m_nDiv[3],
259 : : m_iFirstP, m_iSecondP;
260 : : double m_dFirstP, m_dSecondP, m_curPnt, m_prevPnt, m_boxIncrease,
261 : : m_dIntervalSize[3], m_origin_coords[3], m_dInputSize,
262 : : m_min[3], m_max[3];
263 : : EdgeStatus m_nStatus;
264 : : bool m_bMovedOnce, m_bUseGeom, m_bUseWholeGeom, m_bUseMeshGeom;
265 : : std::vector<iBase_EntityHandle> m_vhVertex;
266 : : std::vector<int> m_vnHexStatus;
267 : : std::map<int, CutFraction> m_mdCutFraction;
268 : : std::vector<EdgeStatus> m_vnEdgeStatus[3];
269 : :
270 : : /** \brief get hex edge status (inside/outside/boundary)
271 : : * \param dZ edge end coordinate
272 : : * \param iSkip how many index skipped for next intersection checking
273 : : * \return EdgeStatus edge status
274 : : */
275 : : EdgeStatus get_edge_status(const double dZ, int& iSkip);
276 : :
277 : : /** \brief set hex status for neighboring elements
278 : : * \param dir current fired ray direction
279 : : * \param i index for i direction
280 : : * \param j index for j direction
281 : : * \param k index for k direction
282 : : * \return bool if is working correctly
283 : : */
284 : : bool set_neighbor_hex_status(int dir, int i, int j, int k);
285 : :
286 : : /** \ brief set hex status by edge status
287 : : * \param index hex index in m_vnHexStatus vector
288 : : * \param value edge status
289 : : * \param dir ray direction
290 : : * \return bool if is working correctly
291 : : */
292 : : bool set_hex_status(int index, EdgeStatus value, int dir);
293 : :
294 : : /** \brief set edge status
295 : : * \param dir ray direction
296 : : * \return bool if is working correctly
297 : : */
298 : : bool set_edge_status(int dir);
299 : :
300 : : /** \brief set all produced mesh information as tag
301 : : * \ hex status, edge-cut information.....
302 : : */
303 : : void set_tag_info();
304 : :
305 : : /** \brief wirte mesh
306 : : * \param file_name
307 : : * \param type element type (inside:0, outside:1, boundary:2)
308 : : * \param handles element handles
309 : : * \param n_elem # of elements
310 : : * \return int if is working correctly
311 : : */
312 : : void write_mesh(const char* file_name, int type,
313 : : iBase_EntityHandle* handles, int& n_elem);
314 : :
315 : : /** \brief get edge fraction information
316 : : * \param idHex index in m_mdCutFraction
317 : : * \param dir ray direction
318 : : * \return double edge fraction
319 : : */
320 : : double get_edge_fraction(int idHex, int dir);
321 : :
322 : : /** \brief get if the edge is fully inside(returns 1) or outside(returns 0)
323 : : * \param i index for i direction
324 : : * \param j index for j direction
325 : : * \param k index for k direction
326 : : * \param dir ray direction
327 : : * \return double edge fraction
328 : : */
329 : : double get_uncut_edge_fraction(int i, int j, int k, int dir);
330 : :
331 : : /** \brief check if the intersected surface geometry is shared or overlapped
332 : : * \param index intersection index in m_vIntersection vector
333 : : * \return bool if is working correctly
334 : : */
335 : : bool is_shared_overlapped_surf(int index);
336 : :
337 : : /** \brief move intersection pairs to check element status
338 : : * \param n_dir ray direction
339 : : * \param n_inter index of intersection points
340 : : * \param start_pnt ray starting point
341 : : * \return bool if is working correctly
342 : : */
343 : : bool move_intersections(int n_dir, int n_inter, double start_pnt[3]);
344 : :
345 : : /** \brief get inside status elements
346 : : * \param rvnInsideCell cell indices (i,j,k triple)
347 : : * \return bool if is working correctly
348 : : */
349 : : bool get_inside_hex(std::vector<int>& rvnInsideCell);
350 : :
351 : : /** \brief check if ray is passing shared vertices or edges
352 : : * \ And, check if the passing surface is overlapped one
353 : : * \param bMoveOnce if ray is already moved before
354 : : * \return bool if is working correctly
355 : : */
356 : : bool is_ray_move_and_set_overlap_surf(bool& bMoveOnce);
357 : :
358 : : void remove_intersection_duplicates();
359 : :
360 : : // test function 1 for debugging
361 : : bool export_fraction_edges(std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& rmdCutCellSurfEdge);
362 : :
363 : : // test functions 2 for debugging
364 : : bool export_fraction_points(std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& mdCutCellEdge);
365 : :
366 : : // test functions 3 for debugging
367 : : bool make_edge(double ePnt[6], std::vector<iBase_EntityHandle>& edge_handles);
368 : :
369 : : //! Static variable, used in registration
370 : : static int init;
371 : :
372 : : /** \brief get MOAB instance
373 : : * \return MOAB instance
374 : : */
375 : 1242 : moab::Interface* moab_instance() {return mk_core()->moab_instance();}
376 : :
377 : : /** \brief get MOAB's Tag
378 : : * \param name Tag name
379 : : * \param size Tag size
380 : : * \param store Tag type
381 : : * \param type data type stored
382 : : * \param def_value default value
383 : : * \param create_if_missing create Tag if it is missed (flag)
384 : : * \return Tag handle
385 : : */
386 : : iBase_TagHandle get_tag(const char* name, int size, unsigned flags, moab::DataType type,
387 : : const void* def_value = NULL, bool create_if_missing = true);
388 : :
389 : : /** \brief get MOAB's various length Tag
390 : : * \param name Tag name
391 : : * \param store Tag type
392 : : * \param type data type stored
393 : : * \return Tag handle
394 : : */
395 : : iBase_TagHandle get_various_length_tag(const char* name,
396 : : unsigned store, moab::DataType type);
397 : :
398 : : /** \brief construct hexes in MOAB structured mesh format
399 : : * \return int if is working correctly
400 : : */
401 : : int make_scd_hexes();
402 : :
403 : : /** \brief construct hexes in MOAB unstructured mesh format
404 : : * \return int if is working correctly
405 : : */
406 : : int make_uscd_hexes();
407 : :
408 : : /** \brief construct OBB tree
409 : : * \input geometry is faceted which are constructed to OBB tree
410 : : * \return int if is working correctly
411 : : */
412 : : int construct_obb_tree();
413 : :
414 : : /** \brief construct and set OBB tree
415 : : * \input Model entity pointer for tree construction
416 : : */
417 : : void set_tree_root(ModelEnt* me);
418 : :
419 : : /** \brief set # of divisions
420 : : * \calcalate the best ones from the size of faceted triangles
421 : : */
422 : : void set_division();
423 : :
424 : : /** \brief find intersections by firing rays
425 : : */
426 : : void find_intersections();
427 : :
428 : : /** \brief fires rays to 3 directions
429 : : * \param dir give the ray direction
430 : : * \it calls fire_ray
431 : : */
432 : : void fire_rays(int dir);
433 : :
434 : : /** \brief fires ray
435 : : * \param nIntersect # of intersections returned
436 : : * \param startPnt ray starting point
437 : : * \param endPnt ray ending point
438 : : * \param tol tolerance to find intersection
439 : : * \param dir ray direction
440 : : * \param rayLength ray length
441 : : * \return bool if is working correctly
442 : : */
443 : : bool fire_ray(int& nIntersect, double startPnt[3],
444 : : double endPnt[3], double tol, int dir,
445 : : double rayLength);
446 : :
447 : : /** \brief moves ray which passes any singluar point
448 : : * \param nIntersect # of intersections returned
449 : : * \param startPnt ray starting point
450 : : * \param endPnt ray ending point
451 : : * \param tol tolerance to find intersection
452 : : * \param dir ray direction
453 : : * \param bMoveOnce if ray is already moved before
454 : : * \return bool if is working correctly
455 : : */
456 : : bool move_ray(int& nIntersect, double* startPnt, double* endPnt,
457 : : double tol, int dir, bool bMoveOnce);
458 : :
459 : : /** \brief if the facet has the same direction to the ray
460 : : * \param i index in m_vIntersection (intersection vector)
461 : : * \param dir ray direction
462 : : * \return bool if is working correctly
463 : : */
464 : : bool is_same_direct_to_ray(int i, int dir);
465 : :
466 : : // ! GeomTopoTool instance
467 : : moab::GeomTopoTool* m_GeomTopoTool;
468 : :
469 : : // ! Tree root
470 : : moab::EntityHandle m_hTreeRoot;
471 : :
472 : : // ! OBB tree tool instance
473 : : moab::OrientedBoxTreeTool* m_hObbTree;
474 : :
475 : : // ! intersected surface geometry list
476 : : std::vector<moab::EntityHandle> m_vhInterSurf;
477 : :
478 : : // ! intersected facet list
479 : : std::vector<moab::EntityHandle> m_vhInterFacet;
480 : :
481 : : // ! overlapped surface list
482 : : std::map<moab::EntityHandle, int> m_mhOverlappedSurf;
483 : :
484 : : std::map<moab::EntityHandle, moab::EntityHandle> m_mRootSets;
485 : :
486 : : double m_minCoord[3], m_maxCoord[3];
487 : : };
488 : :
489 : 1 : inline void EBMesher::increase_box(double box_increase)
490 : : {
491 : 1 : m_boxIncrease = box_increase;
492 : 1 : }
493 : :
494 : 1 : inline void EBMesher::use_whole_geom(int use)
495 : : {
496 : 1 : m_bUseWholeGeom = use;
497 : 1 : }
498 : :
499 : 1 : inline void EBMesher::use_mesh_geometry(bool use)
500 : : {
501 : 1 : m_bUseMeshGeom = use;
502 : 1 : }
503 : : }
504 : :
505 : : #endif
506 : :
507 : :
|
