FuncViz lib
Generalized bitree/quadtree/octree library
Loading...
Searching...
No Matches
MR_cell_cplx.hpp
Go to the documentation of this file.
1// -*- Mode:C++; Coding:us-ascii-unix; fill-column:158 -*-
2/*******************************************************************************************************************************************************.H.S.**/
3/**
4 @file MR_cell_cplx.hpp
5 @author Mitch Richling http://www.mitchr.me/
6 @date 2024-07-13
7 @brief Implementation for the MR_cell_cplx class.@EOL
8 @keywords VTK polydata PLY file MR_rect_tree polygon triangulation cell complex tessellation
9 @std C++23
10 @see MR_rect_tree.hpp, MR_rt_to_cc.hpp
11 @copyright
12 @parblock
13 Copyright (c) 2024, Mitchell Jay Richling <http://www.mitchr.me/> All rights reserved.
14
15 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
16
17 1. Redistributions of source code must retain the above copyright notice, this list of conditions, and the following disclaimer.
18
19 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions, and the following disclaimer in the documentation
20 and/or other materials provided with the distribution.
21
22 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software
23 without specific prior written permission.
24
25 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
26 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
28 OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
30 DAMAGE.
31 @endparblock
32*/
33/*******************************************************************************************************************************************************.H.E.**/
34
35////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
36#ifndef MJR_INCLUDE_MR_cell_cplx
37
38////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
39#include <algorithm> /* STL algorithm C++11 */
40#include <array> /* array template C++11 */
41#include <cmath> /* std:: C math.h C++11 */
42#include <fstream> /* C++ fstream C++98 */
43#include <functional> /* STL funcs C++98 */
44#include <iomanip> /* C++ stream formatting C++11 */
45#include <iostream> /* C++ iostream C++11 */
46#include <map> /* STL map C++11 */
47#include <set> /* STL set C++98 */
48#include <sstream> /* C++ string stream C++ */
49#include <string> /* C++ strings C++11 */
50#include <unordered_map> /* STL hash map C++11 */
51#include <variant> /* C++ variant type C++17 */
52#include <vector> /* STL vector C++11 */
53
54////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
55#include "MR_math.hpp" /* My Simple Math Utilities */
56
57////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
58// Put everything in the mjr namespace
59namespace mjr {
60 /** @brief Template class used to store and transform cell complexes (mesh/triangluation/etc...) and data sets generated from MR_rect_tree sample data.
61
62 The primary use case is to store and manipulate geometric (i.e. mesh) data derived from MR_rect_tree objects. On the mesh manipulation front the focus
63 is on computations that require knowledge related to the function being approximated.
64
65 This code is a quick-n-dirty hack job. I simply wanted to quickly get to a point where I could visually test MR_rect_tree objects. Logically this
66 class should be split into several classes (point, point list, cell, cell list, 3D vectors, 3D analytic geometry, etc...). It's a bit of a mess.
67
68 Design Constraints:
69 - Geometry
70 - Only four types of cells are supported: points, segments, triangles, quads, pyramids, & hexahedrons.
71 - No tetrahedrons! This may seem bazaar; however, they are simply not a natural product of MR_rect_tree tessellation.
72 - Unlike MR_rect_tree, this class can't operate on data of dimension higher than 3!
73 - Data
74 - Only point data is supported (i.e. no cell data)
75 - Data elements are scalar and/or vector. No tensors, fields, etc...
76 - uft_t types are used for atomic data components (scalars are uft_t, and components of vectors are uft_t).
77 - Computations/Manipulation
78 - Supported is focused on manipulation that benefits from knowledge of the approximated function
79 - The idea is that one can use one of the many mesh packages for standard types of mesh computations
80 - I/O
81 - Write only.
82 - VTK files
83 - Only Unstructured_Grid files are supported
84 - Both legacy & XML files are supported, but they are ASCII only.
85 - XML files are serial and self contained
86 - Note that NaN's in ASCII legacy files are not properly handled by many VTK applications
87 - Error Checking
88 - It's pretty limited
89 - Memory allocation -- you run out, the thing crashes
90 - Cells that are a part of other cells may be added (i.e. a segment that is part of an existing triangle)
91 - Geometry checks are quite basic, and I don't even bother with things like convexity. ;)
92
93 Cells supported:
94
95 - Point
96 - Zero-dimensional cell.
97 - Defined by a single point.
98
99 - Segment
100 - One-dimensional cell.
101 - Defined an ordered list of two points.
102 - The direction along the line is from the first point to the second point.
103
104 - Triangle
105 - Two-dimensional cell.
106 - Defined by an ordered list of three points.
107 - The points are ordered counterclockwise, defining a surface normal using the right-hand rule.
108
109 - Quadrilateral
110 - Two-dimensional cell.
111 - Defined by an ordered list of four points.
112 - The points are ordered counterclockwise, defining a surface normal using the right-hand rule.
113 - The points must be coplainer
114 - The quadrilateral is convex and its edges must not intersect.
115
116 - Hexahedron
117 - Three-dimensional cell consisting of six quadrilateral faces, twelve edges, and eight vertices.
118 - The hexahedron is defined by an ordered list of eight points.
119 - The faces and edges must not intersect any other faces and edges, and the hexahedron must be convex.
120
121 - Pyramid
122 - Three-dimensional cell consisting of one quadrilateral face, four triangular faces, eight edges, and five vertices.
123 - The pyramid is defined by an ordered list of five points.
124 - The four points defining the quadrilateral base plane must be convex
125 - The fifth point, the apex, must not be co-planar with the base points.
126
127 \verbatim
128 HEXAHEDRON 4
129 PYRAMID .^.
130 3---------2 . / \ .
131 TRIANGLE QUAD /| /| . / \ .
132 / | / | . / \ .
133 POINT SEGMENT 2 3---------2 7---------6 | 1...../.......\.....0 Back
134 / \ | | | | | | | / \ |
135 0 0-------1 / \ | | | 0------|--1 Back | / \ |
136 / \ | | | / | / | / \ |
137 0-------1 | | |/ |/ | / \ |
138 0---------1 4---------5 Front |/ \|
139 2 ----------------- 3 Front
140 \endverbatim
141
142 A number of quality checks may be performed on points and cells before they are added to the object. These checks can slow down execution of add_cell()
143 by an order of magnitude, and almost double the RAM required for this class. These checks are entirely optional.
144
145 - chk_point_unique:
146 - Many 3D file formats store the master list of potential points, and then use an integer index into this list when defining geometric objects.
147 - Most visualization software is pretty tolerant of having duplicate points on this list when just doing rendering.
148 - Duplicate points can break many software packages pretty badly when it comes to computation.
149 - Depending on how points are added to this object, this check can avoid quite a bit of wasted RAM and produce *much* smaller output files.
150 - I normally leave this check turned on.
151 - chk_cell_unique:
152 - Most packages can render objects with duplicate cells, but you might see graphical artifacts or glitching.
153 - When using a mesh with duplicate cells, computations can sometimes fail or take a *very* long time
154 - I normally leave this check turned on.
155 - chk_cell_vertexes:
156 - This check makes sure that cells have the correct number of vertexes and that those vertexes are valid (on the master point list)
157 - This check also makes sure cells have no duplicate vertexes.
158 - In general this check is pretty cheap, and catches a great many common mesh problems.
159 - This is the only check required for cell_kind_t::POINT and cell_kind_t::SEG cells
160 - chk_cell_dimension:
161 - This check makes sure that PYRAMIDS & HEXAHEDRONS are not co-planar and that TRIANGLES & QUADS are not co-linear.
162 - Combined with chk_cell_vertexes this check is all that is required for cell_kind_t::TRIANGLE cells
163 - chk_cell_edges:
164 - This check makes sure that no cell edges intersect in disallowed ways
165 - This check when combined with the above checks, usually produce cells good enough for most software packages.
166 - cell_kind_t::QUAD cells might be concave or non-plainer.
167 - cell_kind_t::HEXAHEDRON cells might be concave or degenerate (bad edge-face intersections)
168 - cell_kind_t::PYRAMID cells might be concave.
169
170 Levels Of Cell Goodness
171
172 \verbatim
173 +-------------+-------------------+--------------------+--------------------+-------------------------------+-------------------------------+
174 | | chk_point_unique | chk_point_unique | chk_point_unique | chk_point_unique | chk_point_unique |
175 | | chk_cell_unique | chk_cell_unique | chk_cell_unique | chk_cell_unique | chk_cell_unique |
176 | | chk_cell_vertexes | chk_cell_vertexes | chk_cell_vertexes | chk_cell_vertexes | chk_cell_vertexes |
177 | | | chk_cell_dimension | chk_cell_dimension | chk_cell_dimension | chk_cell_dimension |
178 | | | | chk_cell_edges | chk_cell_edges | chk_cell_edges |
179 | | | | | check_cell_face_intersections | check_cell_face_intersections |
180 | | | | | | check_cell_faces_plainer |
181 | | | | | | check_cell_convex |
182 +-------------+-------------------+--------------------+--------------------+-------------------------------+-------------------------------+
183 | vertex | Perfect | | | | |
184 | segment | Perfect | | | | |
185 | triangle | | Perfect | | | |
186 | quad | | | Good Enough | | Perfect |
187 | hexahedron | | | Mostly OK | Good Enough | Perfect |
188 | pyramid | | | Good Enough | | Perfect |
189 +-------------+-------------------+--------------------+--------------------+-------------------------------+-------------------------------+
190 \endverbatim
191
192
193 @tparam chk_point_unique Only allow unique points on the master point list
194 @tparam chk_cell_unique Only allow unique cells on the master cell list. Won't prevent cells that are sub-cells of existing cells.
195 @tparam chk_cell_vertexes Do cell vertex checks (See: cell_stat_t)
196 @tparam chk_cell_dimension Do cell dimension checks (See: cell_stat_t)
197 @tparam chk_cell_edges Do cell edge checks (See: cell_stat_t)
198 @tparam flt_type Type for floating point values (
199 @tparam eps Epsilon used to detect zero */
200 template <bool chk_point_unique,
201 bool chk_cell_unique,
202 bool chk_cell_vertexes,
203 bool chk_cell_dimension,
204 bool chk_cell_edges,
205 typename flt_type,
206 flt_type eps
207 // double eps
208 >
209 requires ((std::is_floating_point<flt_type>::value) &&
210 (eps > std::numeric_limits<flt_type>::epsilon()))
211 class MR_cell_cplx {
212
213 public:
214
215 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
216 /** @name Global Types & Constants. */
217 //@{
218 //--------------------------------------------------------------------------------------------------------------------------------------------------------
219 /** Universal floating point type -- used for point components, scalar data values, vector components, all computation, etc.... */
220 typedef flt_type uft_t;
221 //--------------------------------------------------------------------------------------------------------------------------------------------------------
222 /** epsilon */
223 constexpr static uft_t epsilon = eps;
224 constexpr static uft_t uft_epsilon = std::numeric_limits<uft_t>::epsilon();
225 //@}
226
227 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
228 /** @name Node Data.
229
230 Cell complex objects have a central repository of nodes. Each node has a simple data set consisting of a vector of uft_t values. This data may be
231 used for various purposes; however, three of the values in the vector are always required to specify the node's location in R^3 (i.e. the node's
232 "point"). Other uses include normal vectors and RGB colors. In addition subsets of these values may be named, and extracted by name. Lists of
233 indexes are used to define cells. */
234 //@{
235 //--------------------------------------------------------------------------------------------------------------------------------------------------------
236 /** Type to hold a data for a single node. */
237 typedef std::vector<uft_t> node_data_t;
238 //--------------------------------------------------------------------------------------------------------------------------------------------------------
239 /** Type used to hold an index into a node_data_t or a constant MR_cell_cplx::uft_t. */
240 typedef std::variant<int, uft_t> node_data_idx_t;
241 //--------------------------------------------------------------------------------------------------------------------------------------------------------
242 /** Type to hold a list of data indexes -- one element for scalars, 3 for vectors, etc.... */
243 typedef std::vector<node_data_idx_t> node_data_idx_lst_t;
244 //--------------------------------------------------------------------------------------------------------------------------------------------------------
245 /** Type to hold a list values extracted from node data via a MR_cell_cplx::node_data_idx_lst_t (also see MR_cell_cplx::fvec3_t). */
246 typedef node_data_t vdat_t;
247 //--------------------------------------------------------------------------------------------------------------------------------------------------------
248 /** How to construct a point from a node dataset */
249 node_data_idx_lst_t node_data_to_pnt;
250 //--------------------------------------------------------------------------------------------------------------------------------------------------------
251 /** Integer type used to indentify/index nodes. */
252 typedef int node_idx_t;
253 //--------------------------------------------------------------------------------------------------------------------------------------------------------
254 /** list of node indexes. */
255 typedef std::vector<node_idx_t> node_idx_list_t;
256 //--------------------------------------------------------------------------------------------------------------------------------------------------------
257 /** Hold a tuple of real values defining a point/vector in R^3 (also see MR_cell_cplx::vdat_t). */
258 typedef std::array<uft_t, 3> fvec3_t;
259 //--------------------------------------------------------------------------------------------------------------------------------------------------------
260 /** Hold a list of MR_cell_cplx::fvec3_t objects. */
261 typedef std::vector<fvec3_t> fvec3_list_t;
262 //@}
263
264 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
265 /** @name Named Node Datasets.
266
267 node data may be grouped into named datasets, and written to geometry files. */
268 //@{
269 //--------------------------------------------------------------------------------------------------------------------------------------------------------
270 /** Type to hold a dataset name */
271 typedef std::string named_data_name_t;
272 //--------------------------------------------------------------------------------------------------------------------------------------------------------
273 /** Type to hold a list of dataset names */
274 typedef std::vector<named_data_name_t> named_data_name_list_t;
275 //--------------------------------------------------------------------------------------------------------------------------------------------------------
276 /** Type to map names to named data sets (index lists). */
277 typedef std::map<named_data_name_t, node_data_idx_lst_t> data_name_to_node_data_idx_lst_t;
278 //@}
279
280 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
281 /** @name Cells. */
282 //@{
283 //--------------------------------------------------------------------------------------------------------------------------------------------------------
284 /** Cell Status. */
285 enum class cell_stat_t { GOOD, //!< Looks like a good cell N/A
286 TOO_FEW_PNT, //!< List of points was empty check_cell_vertexes
287 TOO_MANY_PNT, //!< List of points was too long (>5) check_cell_vertexes
288 NEG_PNT_IDX, //!< Negative point index (i.e. not on the central node data) check_cell_vertexes
289 BIG_PNT_IDX, //!< Point index was too big (i.e. not on the central node data) check_cell_vertexes
290 DUP_PNT, //!< At least two points have the same index check_cell_vertexes
291 DIM_LOW, //!< Dimension low (degenerate) check_cell_dimension
292 BAD_EDGEI, //!< Bad edge-edge intersection check_cell_edge_intersections
293 BAD_FACEI, //!< Bad face-edge intersection check_cell_face_intersections
294 FACE_BENT, //!< A face (QUAD, HEXAHEDRON, or PYRAMID) was not plainer check_cell_faces_plainer
295 CONCAVE, //!< (QUAD, HEXAHEDRON, or PYRAMID) was concave TBD
296 };
297 //--------------------------------------------------------------------------------------------------------------------------------------------------------
298 /** Cell Types. */
299 enum class cell_kind_t { POINT,
300 SEGMENT,
301 TRIANGLE,
302 QUAD,
303 // TETRAHEDRON
304 HEXAHEDRON,
305 PYRAMID,
306 };
307 //--------------------------------------------------------------------------------------------------------------------------------------------------------
308 /** Type to hold the vertexes of a poly cell */
309 typedef node_idx_list_t cell_verts_t;
310 //--------------------------------------------------------------------------------------------------------------------------------------------------------
311 /** Type to hold a poly cell -- a list of point indexes */
312 // struct cell_t {
313 // cell_kind_t type;
314 // cell_verts_t verts;
315 // };
316 //@}
317
318 private:
319
320 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
321 /** @name Node Data. */
322 //@{
323 //--------------------------------------------------------------------------------------------------------------------------------------------------------
324 /** The index of the last point added via the add_node() method. */
325 node_idx_t last_point_idx = -1;
326 //--------------------------------------------------------------------------------------------------------------------------------------------------------
327 /** True if the last point given to the add_node() method was new -- i.e. not on the central node data.
328 Only updated if chk_point_unique is true. See: last_point_added_was_new() */
329 bool last_point_new = true;
330 /** @cond */ // Leaving this out of doxygen because I don't like how it shows up...
331 //--------------------------------------------------------------------------------------------------------------------------------------------------------
332 /** Less operator for pnt_lst_uniq_t.
333 @return If a & b are close in space, then return false. Otherwise uses lexicographic ordering. */
334 struct pnt_less {
335 bool operator()(const fvec3_t& a, const fvec3_t& b) const { return (((std::abs(a[0]-b[0]) > eps) ||
336 (std::abs(a[1]-b[1]) > eps) ||
337 (std::abs(a[2]-b[2]) > eps)) &&
338 ((a[0] < b[0]) ||
339 ((a[0] == b[0]) &&
340 ((a[1] < b[1]) ||
341 ((a[1] == b[1]) && (a[2] < b[2])))))); }
342 };
343 /** @endcond */
344 //--------------------------------------------------------------------------------------------------------------------------------------------------------
345 /** Used to uniquify points and assign point index values */
346 typedef std::map<fvec3_t, node_idx_t, pnt_less> pnt_to_node_idx_map_t;
347 //--------------------------------------------------------------------------------------------------------------------------------------------------------
348 /** Maps points to point index -- used to detect points with nearly identical geometric points in R^3 */
349 pnt_to_node_idx_map_t pnt_to_node_idx_map;
350 //--------------------------------------------------------------------------------------------------------------------------------------------------------
351 /** Type to hold a all node data. */
352 typedef std::vector<node_data_t> node_idx_to_node_data_t;
353 //--------------------------------------------------------------------------------------------------------------------------------------------------------
354 /** Hold all node data sets */
355 node_idx_to_node_data_t node_idx_to_node_data;
356 //@}
357
358 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
359 /** @name Named Node Datasets. */
360 //@{
361 //--------------------------------------------------------------------------------------------------------------------------------------------------------
362 /** Hold named descriptions of data (scalars & vectors) */
363 data_name_to_node_data_idx_lst_t data_name_to_data_idx_lst;
364 //@}
365
366 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
367 /** @name Cells. */
368 //@{
369 //--------------------------------------------------------------------------------------------------------------------------------------------------------
370 /** Type for a list of poly cells */
371 typedef std::vector<cell_verts_t> cell_lst_t;
372 //--------------------------------------------------------------------------------------------------------------------------------------------------------
373 /** House all poly cells */
374 cell_lst_t cell_lst;
375 //--------------------------------------------------------------------------------------------------------------------------------------------------------
376 /** Used to uniquify cells */
377 typedef std::set<cell_verts_t> uniq_cell_lst_t;
378 //--------------------------------------------------------------------------------------------------------------------------------------------------------
379 /** Unique cell list. */
380 uniq_cell_lst_t uniq_cell_lst;
381 //--------------------------------------------------------------------------------------------------------------------------------------------------------
382 /** True if the last cell given to the add_cell() method was new -- i.e. not on the master cell list.
383 Only updated if chk_cell_unique is true. Value is invalid if last_cell_stat is NOT cell_stat_t::GOOD.
384 See: last_cell_added_was_new() */
385 bool last_cell_new = true;
386 //--------------------------------------------------------------------------------------------------------------------------------------------------------
387 /** Status of the last cell added via add_cell() method.
388 Only updated if (chk_cell_vertexes || chk_cell_dimension | chk_cell_edges) is true. See: status_of_last_cell_added() */
389 cell_stat_t last_cell_stat = cell_stat_t::GOOD;
390 //@}
391
392 public:
393
394 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
395 /** @name Node Data. */
396 //@{
397 //--------------------------------------------------------------------------------------------------------------------------------------------------------
398 /** Return the number of points in the central node data list.
399 Note the return type is node_idx_t (a signed integer type) and not a size_t. */
400 inline node_idx_t node_count() const {
401 return static_cast<node_idx_t>(node_idx_to_node_data.size()); // Yes. We mean node_idx_to_node_data.
402 }
403 //--------------------------------------------------------------------------------------------------------------------------------------------------------
404 /** Return point given a node index */
405 inline fvec3_t get_pnt(node_idx_t pnt_idx) const {
406 return fvec3_from_node_data(node_data_to_pnt, node_idx_to_node_data[pnt_idx]);
407 }
408 //--------------------------------------------------------------------------------------------------------------------------------------------------------
409 /** Retruns the index of the last point given to the add_node() method. */
410 inline node_idx_t idx_of_last_point_added() const {
411 return last_point_idx;
412 }
413 //--------------------------------------------------------------------------------------------------------------------------------------------------------
414 /** Retruns true if the last point given to the add_node() method was a new point. */
415 inline node_idx_t last_point_added_was_new() const {
416 return last_point_new;
417 }
418 //--------------------------------------------------------------------------------------------------------------------------------------------------------
419 /** Add a point to the central node data store.
420
421 The geometric point coordinates are extracted, and tested before the node data is added leading to three cases:
422
423 Case 1: If the geometric point contains a NaN, then the point is NOT added to the central node data store, last_point_idx=-1, and
424 last_point_new=false.
425
426 Case 2: If the geometric point is already on the list and chk_point_unique is true, then the point is NOT added to the central node data store,
427 last_point_idx is stet to the existing point's index, and last_point_new=false. Only the geometric point components are used to determine if a
428 point exists in the store.
429
430 Case 3: If the geometric point is not on the list, then the point is added to the central node data store, last_point_idx is set to the new point's
431 index, and last_point_new=true Note that last_point_idx is always the return value.
432
433 @param node_data The point to add! */
434 inline node_idx_t add_node(node_data_t node_data) {
435 if (node_data_to_pnt.empty())
436 node_data_to_pnt = {0, 1, 2};
437 fvec3_t new_pnt = fvec3_from_node_data(node_data_to_pnt, node_data);
438 if (pnt_has_nan(new_pnt)) {
439 last_point_idx = -1;
440 last_point_new = false;
441 } else {
442 if constexpr (chk_point_unique) {
443 if (pnt_to_node_idx_map.contains(new_pnt)) {
444 /* Point is already in list */
445 last_point_idx = pnt_to_node_idx_map[new_pnt];
446 last_point_new = false;
447 } else {
448 /* Point is not already in list */
449 last_point_idx = node_count();
450 pnt_to_node_idx_map[new_pnt] = last_point_idx;
451 //node_idx_to_pnt.push_back(new_pnt);
452 node_idx_to_node_data.push_back(node_data);
453 last_point_new = true;
454 }
455 } else {
456 /* Add point without regard to uniqueness */
457 last_point_idx = node_count();
458 //node_idx_to_pnt.push_back(new_pnt);
459 node_idx_to_node_data.push_back(node_data);
460 last_point_new = true;
461 }
462 }
463 return last_point_idx;
464 }
465 //--------------------------------------------------------------------------------------------------------------------------------------------------------
466 /** Construct a fvec3_t from a node_data_t and a node_data_idx_lst_t.
467 This is a speciazlied form of vector_from_node_data() that produces fvec3_t objects. Primarly used to extract the location coordinates from node data.
468 Note data_idx_lst must contain at least three values, but only the first three are used. If the node_data_idx_lst_t has fewer than three elements, then
469 this function will throw.
470 @param data_idx_lst Index list used to identify the data to pull
471 @param node_data node data vector from which to pull data */
472 inline fvec3_t fvec3_from_node_data(node_data_idx_lst_t data_idx_lst, node_data_t node_data) const {
473 fvec3_t v;
474 for(int i=0; i<3; ++i)
475 if (data_idx_lst[i].index() == 0)
476 v[i] = node_data[std::get<int>(data_idx_lst.at(i))];
477 else
478 v[i] = std::get<uft_t>(data_idx_lst.at(i));
479 return v;
480 }
481 //--------------------------------------------------------------------------------------------------------------------------------------------------------
482 /** Construct a vdat_t from a node_data_t and a node_data_idx_lst_t.
483 @param data_idx_lst Index list used to identify the data to pull
484 @param node_data node data vector from which to pull data */
485 inline vdat_t vector_from_node_data(node_data_idx_lst_t data_idx_lst, node_data_t node_data) const {
486 vdat_t v;
487 for(decltype(data_idx_lst.size()) i=0; i<data_idx_lst.size(); ++i)
488 if (data_idx_lst[i].index() == 0)
489 v.push_back(node_data[std::get<int>(data_idx_lst.at(i))]);
490 else
491 v.push_back(std::get<uft_t>(data_idx_lst.at(i)));
492 return v;
493 }
494 //--------------------------------------------------------------------------------------------------------------------------------------------------------
495 /** Extract a uft_t from a node_data_t and a node_data_idx_t.
496 @param data_idx Index list used to identify the data to pull
497 @param node_data node data vector from which to pull data */
498 inline uft_t scalar_from_node_data(node_data_idx_t data_idx, node_data_t node_data) const {
499 if (data_idx.index() == 0)
500 return (node_data[std::get<int>(data_idx)]);
501 else
502 return (std::get<int>(data_idx));
503 }
504 //--------------------------------------------------------------------------------------------------------------------------------------------------------
505 /** Convert a fvec3_t to a string representation
506 @param pnt_idx The point to convert.*/
507 inline std::string node_to_string(node_idx_t pnt_idx) const {
508 std::ostringstream convert;
509 if (pnt_idx >= 0) {
510 convert << "[ ";
511 for(auto c: node_idx_to_node_data[pnt_idx])
512 convert << std::setprecision(5) << static_cast<uft_t>(c) << " ";
513 convert << "]";
514 } else {
515 convert << "[ DNE ]";
516 }
517 return(convert.str());
518 }
519 //@}
520
521 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
522 /** @name Named Node Datasets. */
523 //@{
524 //--------------------------------------------------------------------------------------------------------------------------------------------------------
525 /** Return the number of named scalar datasets */
526 inline int named_scalar_datasets_count() const {
527 return static_cast<int>(std::count_if(data_name_to_data_idx_lst.cbegin(), data_name_to_data_idx_lst.cend(), [](auto x) { return x.second.size()==1; }));
528 }
529 //--------------------------------------------------------------------------------------------------------------------------------------------------------
530 /** Return the number of named vector datasets */
531 inline int named_vector_datasets_count() const {
532 return static_cast<int>(std::count_if(data_name_to_data_idx_lst.cbegin(), data_name_to_data_idx_lst.cend(), [](auto x) { return x.second.size()!=1; }));
533 }
534 //--------------------------------------------------------------------------------------------------------------------------------------------------------
535 /** Return the number of named (both vector and scalar) datasets */
536 inline int named_datasets_count() const {
537 return static_cast<int>(data_name_to_data_idx_lst.size());
538 }
539 //--------------------------------------------------------------------------------------------------------------------------------------------------------
540 /** Create named data sets by setting raw member.
541 @param names value to set. */
545 //--------------------------------------------------------------------------------------------------------------------------------------------------------
546 /** Create named data sets.
547 Simple way to create named scalar datasets for the first n elements in each point's data.
548 @param scalar_name_strings Scalar data set names to apply, in order, to node data elements. */
549 inline void create_named_datasets(named_data_name_list_t scalar_name_strings) {
550 data_name_to_data_idx_lst.clear();
551 for(int i=0; i<static_cast<int>(scalar_name_strings.size()); ++i)
552 data_name_to_data_idx_lst[scalar_name_strings[i]] = {i};
553 }
554 //--------------------------------------------------------------------------------------------------------------------------------------------------------
555 /** @overload
556 @param scalar_name_strings Scalar data set names to apply, in order, to node data elements.
557 @param names data_name_to_data_idx entries to create -- frequently used to add a named vector dataset. */
558 inline void create_named_datasets(named_data_name_list_t scalar_name_strings, data_name_to_node_data_idx_lst_t names) {
559 data_name_to_data_idx_lst.clear();
560 for(int i=0; i<static_cast<int>(scalar_name_strings.size()); ++i)
561 data_name_to_data_idx_lst[scalar_name_strings[i]] = {i};
562 for(auto kv : names)
563 data_name_to_data_idx_lst[kv.first] = kv.second;
564 }
565 //--------------------------------------------------------------------------------------------------------------------------------------------------------
566 /** Create mapping beteen node data and geometric pionts.
567 Geometric points are constructed from node data via a node_data_idx_lst_t -- just like named data sets. Each time a geometric point is required,
568 it is dynamically pulled from the point's data via this mapping.
569 @param point_data_index A node_data_idx_lst_t defining how to construct a geometric point. */
570 inline void create_dataset_to_point_mapping(node_data_idx_lst_t point_data_index) {
571 node_data_to_pnt = point_data_index;
572 }
573 //@}
574
575 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
576 /** @name 3D Vector Computations.
577
578 These functions operate on fvec3_t objects (arrays of 3 real numbers), but treat them as 3-vectors. */
579 //@{
580 //--------------------------------------------------------------------------------------------------------------------------------------------------------
581 /** Compute the magnitude */
582 inline uft_t vec3_two_norm(const fvec3_t& pnt) const {
583 return std::sqrt(vec3_self_dot_product(pnt));
584 }
585 //--------------------------------------------------------------------------------------------------------------------------------------------------------
586 /** Compute the self dot product -- i.e. magnitude squared */
587 inline uft_t vec3_self_dot_product(const fvec3_t& pnt) const {
588 uft_t tmp = 0.0;
589 for(int i=0; i<3; ++i)
590 tmp += pnt[i]*pnt[i];
591 return tmp;
592 }
593 //--------------------------------------------------------------------------------------------------------------------------------------------------------
594 /** Compute the cross prodcut */
595 inline uft_t vec3_dot_product(const fvec3_t& pnt1, const fvec3_t& pnt2) const {
596 uft_t tmp = 0.0;
597 for(int i=0; i<3; ++i)
598 tmp += pnt1[i]*pnt2[i];
599 return tmp;
600 }
601 //--------------------------------------------------------------------------------------------------------------------------------------------------------
602 /** Compute the cross prodcut */
603 inline fvec3_t vec3_cross_product(const fvec3_t& pnt1, const fvec3_t& pnt2) const {
604 fvec3_t tmp;
605 tmp[0] = pnt1[1]*pnt2[2]-pnt1[2]*pnt2[1];
606 tmp[1] = pnt1[2]*pnt2[0]-pnt1[0]*pnt2[2];
607 tmp[2] = pnt1[0]*pnt2[1]-pnt1[1]*pnt2[0];
608 return tmp;
609 }
610 //--------------------------------------------------------------------------------------------------------------------------------------------------------
611 /** Compute the vector diffrence */
612 inline fvec3_t vec3_diff(const fvec3_t& pnt1, const fvec3_t& pnt2) const {
613 fvec3_t tmp;
614 for(int i=0; i<3; ++i)
615 tmp[i] = pnt1[i] - pnt2[i];
616 return tmp;
617 }
618 //--------------------------------------------------------------------------------------------------------------------------------------------------------
619 /** Compute the vector diffrence */
620 inline uft_t vec3_scalar_triple_product(const fvec3_t& pnt1, const fvec3_t& pnt2, const fvec3_t& pnt3) const {
621 return vec3_dot_product(pnt1, vec3_cross_product(pnt2, pnt3));
622 }
623 //--------------------------------------------------------------------------------------------------------------------------------------------------------
624 /** Unitize the given point in place. Return true if the result is valid, and false otherwise */
625 inline bool vec3_unitize(fvec3_t& pnt) const {
626 uft_t length = vec3_two_norm(pnt);
627 if (std::abs(length) > eps) {
628 for(int i=0; i<3; ++i)
629 pnt[i] = pnt[i]/length;
630 return true;
631 } else {
632 return false;
633 }
634 }
635 //--------------------------------------------------------------------------------------------------------------------------------------------------------
636 /** Compute the linear combination */
637 inline fvec3_t vec3_linear_combination(uft_t scl1, const fvec3_t& pnt1, uft_t scl2, const fvec3_t& pnt2) const {
638 fvec3_t tmp;
639 for(int i=0; i<3; ++i)
640 tmp[i] = scl1 * pnt1[i] + scl2 * pnt2[i];
641 return tmp;
642 }
643 //--------------------------------------------------------------------------------------------------------------------------------------------------------
644 /** Determinant of 3x3 matrix with given vectors as columns/rows. */
645 inline uft_t vec3_det3(const fvec3_t& pnt1, const fvec3_t& pnt2, const fvec3_t& pnt3) const {
646 // MJR TODO NOTE <2024-07-22T15:48:31-0500> vec3_det3: UNTESTED! UNTESTED! UNTESTED! UNTESTED!
647 static_assert(false, "vec3_det3: Under active development. Untested!");
648 return (pnt1[0] * pnt2[1] * pnt3[2] -
649 pnt1[0] * pnt2[2] * pnt3[1] -
650 pnt1[1] * pnt2[0] * pnt3[2] +
651 pnt1[1] * pnt2[2] * pnt3[0] +
652 pnt1[2] * pnt2[0] * pnt3[1] -
653 pnt1[2] * pnt2[1] * pnt3[0]);
654 }
655 //@}
656
657 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
658 /** @name 3D Geometry Computations.
659
660 Methods with geomi_ prefix take point index types while methods with a geomr_ prefix take double tuples. */
661 //@{
662 //--------------------------------------------------------------------------------------------------------------------------------------------------------
663 /** Segment intersection types */
664 enum class seg_isect_t { C0_EMPTY, //!< Cardinality=0. Intersection is the empty set.
665 C1_VERTEX1, //!< Cardinality=1. Intersection is a single, shared vertex.
666 C1_INTERIOR, //!< Cardinality=1. Intersection is a single point in interior of at least one segment.
667 CI_VERTEX2, //!< Cardinality=infinite. Intersection is a segment -- equal to input segments.
668 CI_VERTEX1, //!< Cardinality=infinite. Intersection is a segment -- exactially one vertex in the intersection.
669 CI_VERTEX0, //!< Cardinality=infinite. Intersection is a segment -- no vertexes included
670 BAD_SEGMENT, //!< At least one of the segments was degenerate
671 };
672 //--------------------------------------------------------------------------------------------------------------------------------------------------------
673 /** Convert an seg_isect_t to a string */
674 inline std::string seg_isect_to_string(seg_isect_t seg_isect) const {
675 switch(seg_isect) {
676 case seg_isect_t::C0_EMPTY: return std::string("C0_EMPTY"); break;
677 case seg_isect_t::C1_VERTEX1: return std::string("C1_VERTEX1"); break;
678 case seg_isect_t::C1_INTERIOR: return std::string("C1_INTERIOR"); break;
679 case seg_isect_t::CI_VERTEX2: return std::string("CI_VERTEX2"); break;
680 case seg_isect_t::CI_VERTEX1: return std::string("CI_VERTEX1"); break;
681 case seg_isect_t::CI_VERTEX0: return std::string("CI_VERTEX0"); break;
682 case seg_isect_t::BAD_SEGMENT: return std::string("BAD_SEGMENT"); break;
683 }
684 return std::string(""); // Never get here, but some compilers can't figure that out. ;)
685 }
686 //--------------------------------------------------------------------------------------------------------------------------------------------------------
687 /** Determine the nature of the intersection between two line segments */
688 seg_isect_t geomi_seg_isect_type(node_idx_t ilin1pnt1, node_idx_t ilin1pnt2, node_idx_t ilin2pnt1, node_idx_t ilin2pnt2) const {
689 // MJR TODO NOTE <2024-08-02T09:41:48-0500> geomi_seg_isect_type: Repeated point look-up slows things down
690 // MJR TODO NOTE <2024-08-02T09:41:48-0500> geomi_seg_isect_type: Add unit tests for each code branch
691 // Check for degenerate segments
692 if (ilin1pnt1 == ilin1pnt2)
693 return seg_isect_t::BAD_SEGMENT;
694 if (ilin2pnt1 == ilin2pnt2)
695 return seg_isect_t::BAD_SEGMENT;
696 // Count unique points & break into cases
697 std::set<node_idx_t> points_sorted;
698 points_sorted.insert(ilin1pnt1);
699 points_sorted.insert(ilin1pnt2);
700 points_sorted.insert(ilin2pnt1);
701 points_sorted.insert(ilin2pnt2);
702 if (points_sorted.size() == 4) { // ...................................................... REMAINING CASES: C0_EMPTY, C1_INTERIOR, CI_VERTEX0
703 if (geomi_pts_colinear(ilin1pnt1, ilin1pnt2, ilin2pnt1, ilin2pnt2)) { // ............... REMAINING CASES: C0_EMPTY, CI_VERTEX0
704 if ( (geomi_pnt_line_distance(ilin1pnt1, ilin1pnt2, ilin2pnt1, true) < eps) ||
705 (geomi_pnt_line_distance(ilin1pnt1, ilin1pnt2, ilin2pnt2, true) < eps) ||
706 (geomi_pnt_line_distance(ilin2pnt1, ilin2pnt2, ilin1pnt1, true) < eps) ||
707 (geomi_pnt_line_distance(ilin2pnt1, ilin2pnt2, ilin1pnt2, true) < eps) ) { // .. REMAINING CASES: CI_VERTEX0
708 return seg_isect_t::CI_VERTEX0;
709 } else { // ......................................................................... REMAINING CASES: C0_EMPTY
710 return seg_isect_t::C0_EMPTY;
711 }
712 } else { // ........................................................................... REMAINING CASES: C0_EMPTY, C1_INTERIOR
713 if (geomi_seg_isect1(ilin1pnt1, ilin1pnt2, ilin2pnt1, ilin2pnt2)) { // .............. REMAINING CASES: C1_INTERIOR
714 return seg_isect_t::C1_INTERIOR;
715 } else { // ......................................................................... REMAINING CASES: C0_EMPTY
716 return seg_isect_t::C0_EMPTY;
717 }
718 }
719 return seg_isect_t::C0_EMPTY;
720 } else if (points_sorted.size() == 3) { // .............................................. REMAINING CASES: C1_VERTEX1, CI_VERTEX1
721 node_idx_t ipnt1, ipnt2, ipntc;
722 if (ilin1pnt1 == ilin2pnt1) {
723 ipntc = ilin1pnt1; ipnt1 = ilin1pnt2; ipnt2 = ilin2pnt2;
724 } else if (ilin1pnt1 == ilin2pnt2) {
725 ipntc = ilin1pnt1; ipnt1 = ilin1pnt2; ipnt2 = ilin2pnt1;
726 } else if (ilin1pnt2 == ilin2pnt1) {
727 ipntc = ilin1pnt2; ipnt1 = ilin1pnt1; ipnt2 = ilin2pnt2;
728 } else if (ilin1pnt2 == ilin2pnt2) {
729 ipntc = ilin1pnt2; ipnt1 = ilin1pnt1; ipnt2 = ilin2pnt1;
730 } else { // Never get here. Silences compiler warnings.
731 ipntc = 0; ipnt1 = 0; ipnt2 = 0;
732 }
733 if (geomi_pts_colinear(ipnt1, ipnt2, ipntc)) { // ..................................... REMAINING CASES: C1_VERTEX1, CI_VERTEX1
734 if ( (geomi_pnt_line_distance(ipnt1, ipntc, ipnt2, true) < eps) ||
735 (geomi_pnt_line_distance(ipnt2, ipntc, ipnt1, true) < eps) ) { // .............. REMAINING CASES: CI_VERTEX1
736 return seg_isect_t::CI_VERTEX1;
737 } else { // ......................................................................... REMAINING CASES: CI_VERTEX1
738 return seg_isect_t::C1_VERTEX1;
739 }
740 } else { // ........................................................................... REMAINING CASES: C1_VERTEX1
741 return seg_isect_t::C1_VERTEX1;
742 }
743 } else if (points_sorted.size() == 2) { // .............................................. REMAINING CASES: CI_VERTEX2
744 return seg_isect_t::CI_VERTEX2;
745 } else { // (points_sorted.size() == 1) which is an error. // ........................... REMAINING CASES: CI_VERTEX2
746 return seg_isect_t::CI_VERTEX2;
747 }
748 return seg_isect_t::C0_EMPTY; // Should never get here...
749 }
750 //--------------------------------------------------------------------------------------------------------------------------------------------------------
751 /** Check if two line segments intersect in a single point */
752 inline bool geomi_seg_isect1(node_idx_t ilin1pnt1, node_idx_t ilin1pnt2, node_idx_t ilin2pnt1, node_idx_t ilin2pnt2) const {
753 return geomr_seg_isect1(get_pnt(ilin1pnt1), get_pnt(ilin1pnt2), get_pnt(ilin2pnt1), get_pnt(ilin2pnt2));
754 }
755 //--------------------------------------------------------------------------------------------------------------------------------------------------------
756 /** Check if two line segments intersect in a single point */
757 inline bool geomr_seg_isect1(const fvec3_t& lin1pnt1, const fvec3_t& lin1pnt2, const fvec3_t& lin2pnt1, const fvec3_t& lin2pnt2) const {
758 uft_t denom =
759 lin1pnt1[0] * lin2pnt1[1] - lin1pnt1[0] * lin2pnt2[1] - lin1pnt1[1] * lin2pnt1[0] + lin1pnt1[1] * lin2pnt2[0] -
760 lin1pnt2[0] * lin2pnt1[1] + lin1pnt2[0] * lin2pnt2[1] + lin1pnt2[1] * lin2pnt1[0] - lin1pnt2[1] * lin2pnt2[0];
761 if (std::abs(denom) < eps) // Lines are parallel
762 return false;
763 uft_t numera =
764 lin1pnt1[0]*lin2pnt1[1] - lin1pnt1[0]*lin2pnt2[1] -
765 lin1pnt1[1]*lin2pnt1[0] + lin1pnt1[1]*lin2pnt2[0] +
766 lin2pnt1[0]*lin2pnt2[1] - lin2pnt1[1]*lin2pnt2[0];
767 uft_t numerb =
768 -(lin1pnt1[0]*lin1pnt2[1] - lin1pnt1[0]*lin2pnt1[1] -
769 lin1pnt1[1]*lin1pnt2[0] + lin1pnt1[1]*lin2pnt1[0] +
770 lin1pnt2[0]*lin2pnt1[1] - lin1pnt2[1]*lin2pnt1[0]);
771 uft_t ua = numera/denom;
772 uft_t ub = numerb/denom;
773 if ( (ua < 0) || (ub < 0) || (ua > 1) || (ub > 1) )
774 return false;
775 uft_t eq3 = numera/denom * (lin1pnt2[2] - lin1pnt1[2]) - numerb/denom * (lin2pnt2[2] - lin2pnt1[2]) + lin2pnt1[2] + lin1pnt1[2];
776 if (std::abs(eq3) < eps) // Equation in third coordinate is satisfied
777 return true;
778 else
779 return false;
780 }
781 //--------------------------------------------------------------------------------------------------------------------------------------------------------
782 /** Distance between a point and a line.
783 See: geomr_pnt_line_distance(). */
784 inline uft_t geomi_pnt_line_distance(node_idx_t ilinpnt1, node_idx_t ilinpnt2, node_idx_t ipnt, bool seg_distance) const {
785 return geomr_pnt_line_distance(get_pnt(ilinpnt1), get_pnt(ilinpnt2), get_pnt(ipnt), seg_distance);
786 }
787 //--------------------------------------------------------------------------------------------------------------------------------------------------------
788 /** Distance between a point and a line.
789 The result depends on the value of seg_distance:
790 - seg_distance==false: Distance between the point pnt and the line containing linpnt1 & linpnt2.
791 - seg_distance==true: Distance between the point pnt and the line segment with endpoints linpnt1 & linpnt2. */
792 uft_t geomr_pnt_line_distance(const fvec3_t& linpnt1, const fvec3_t& linpnt2, const fvec3_t& pnt, bool seg_distance) const {
793 uft_t segd = geomr_pnt_pnt_distance(linpnt1, linpnt2);
794 fvec3_t d, p;
795 uft_t t = 0;
796 for(int i=0; i<3; i++) {
797 d[i] = (linpnt2[i] - linpnt1[i]) / segd;
798 t += (pnt[i] - linpnt1[i])*d[i];
799 }
800 for(int i=0; i<3; i++) {
801 p[i] = linpnt1[i] + t * d[i];
802 }
803 // p is the point on the line nearest pnt
804 if (seg_distance) {
805 uft_t dp1 = geomr_pnt_pnt_distance(linpnt1, p);
806 uft_t dp2 = geomr_pnt_pnt_distance(linpnt2, p);
807 // MJR TODO NOTE <2024-08-02T09:42:19-0500> geomr_pnt_line_distance: check logic -- dp1>segd || dp2>segd?
808 if (std::abs((dp1+dp2)-segd) > eps)
809 return std::min(geomr_pnt_pnt_distance(linpnt1, pnt), geomr_pnt_pnt_distance(linpnt2, pnt));
810 }
811 return geomr_pnt_pnt_distance(p, pnt);
812 }
813 //--------------------------------------------------------------------------------------------------------------------------------------------------------
814 /** Compute the euculedian (2 norm) distance between two points */
815 inline uft_t geomr_pnt_pnt_distance(const fvec3_t& pnt1, const fvec3_t& pnt2) const {
816 return vec3_two_norm(vec3_diff(pnt1, pnt2));
817 }
818 //--------------------------------------------------------------------------------------------------------------------------------------------------------
819 /** Compute the normal of a triangle (or a plane defined by 3 points) */
820 inline fvec3_t geomr_tri_normal(const fvec3_t& tripnt1, const fvec3_t& tripnt2, const fvec3_t& tripnt3, bool unit) const {
821 fvec3_t basisv1 = vec3_diff(tripnt1, tripnt2); // basis vectors for pln containing triagnel
822 fvec3_t basisv2 = vec3_diff(tripnt3, tripnt2); // basis vectors for pln containing triagnel
823 fvec3_t normal = vec3_cross_product(basisv1, basisv2); // normal vector for tri. n=pld1xpld2
824 if (unit)
825 vec3_unitize(normal);
826 return normal;
827 }
828 //--------------------------------------------------------------------------------------------------------------------------------------------------------
829 /** Compute the distance between a point and a plane */
830 inline uft_t geomr_pnt_pln_distance(const fvec3_t& plnpnt1, const fvec3_t& plnpnt2, const fvec3_t& plnpnt3, const fvec3_t& pnt) const {
831 fvec3_t n = geomr_tri_normal(plnpnt1, plnpnt2, plnpnt3, true);
832 return std::abs(vec3_dot_product(n, pnt) - vec3_dot_product(n, plnpnt2));
833 }
834 //--------------------------------------------------------------------------------------------------------------------------------------------------------
835 /** Compute the distance between a point and a triangle */
836 inline uft_t geomr_pnt_tri_distance(const fvec3_t& tripnt1, const fvec3_t& tripnt2, const fvec3_t& tripnt3, const fvec3_t& pnt) const {
837 // MJR TODO NOTE <2024-07-22T15:48:31-0500> geomr_pnt_tri_distance: UNTESTED! UNTESTED! UNTESTED! UNTESTED! UNTESTED!
838 fvec3_t basisv1 = vec3_diff(tripnt1, tripnt2); // basis vectors for pln containing triagnel
839 fvec3_t basisv2 = vec3_diff(tripnt3, tripnt2); // basis vectors for pln containing triagnel
840 fvec3_t normal = vec3_cross_product(basisv1, basisv2); // normal vector for tri. ax+by+cz+d=0, a=normal[0], b=normal[1], c=normal[2]
841 vec3_unitize(normal);
842 uft_t d = -vec3_dot_product(normal, tripnt2); // ax+by+cz+d=0
843 uft_t lambda = vec3_dot_product(normal, pnt) + d;
844 fvec3_t q = vec3_diff(vec3_linear_combination(1.0, pnt, lambda, normal), tripnt2); // q is the point in the plane closest to pnt
845 uft_t denom = basisv1[0] * basisv2[1] - basisv2[0] * basisv1[1]; // If zero, then triangle is broken!
846 uft_t u = (q[0] * basisv2[1] - basisv2[0] * q[1]) / denom;
847 uft_t v = (basisv1[0] * q[1] - q[0] * basisv1[1]) / denom;
848 uft_t dd = std::abs(u*basisv1[2] + v*basisv2[2] - q[2]);
849 if ( (u>=0) && (v>=0) && ((u+v)<=1) && (dd<eps) ) { // q is in the triangle => Use the plane distance
850 return std::abs(lambda);
851 } else { // q is not in the triangle => Distance will be minimum distance to an edge.
852 uft_t d1 = geomr_pnt_line_distance(tripnt1, tripnt2, pnt, true);
853 uft_t d2 = geomr_pnt_line_distance(tripnt2, tripnt3, pnt, true);
854 uft_t d3 = geomr_pnt_line_distance(tripnt3, tripnt1, pnt, true);
855 return std::min(std::min(d1, d2), d3);
856 }
857 }
858 //--------------------------------------------------------------------------------------------------------------------------------------------------------
859 /** Check if points (just one point for this function) are zero.
860 Checks if the infinity norm is less than epsilon.*/
861 inline bool geomr_pnt_zero(const fvec3_t& p1) const {
862 return ((std::abs(p1[0]) < eps) &&
863 (std::abs(p1[1]) < eps) &&
864 (std::abs(p1[2]) < eps));
865 }
866 //--------------------------------------------------------------------------------------------------------------------------------------------------------
867 /** Check if points are colinear */
868 inline bool geomi_pts_colinear(node_idx_t pi1, node_idx_t pi2, node_idx_t pi3, node_idx_t pi4) const {
869 return ( geomr_pts_colinear(get_pnt(pi1), get_pnt(pi2), get_pnt(pi3)) &&
870 geomr_pts_colinear(get_pnt(pi1), get_pnt(pi2), get_pnt(pi4)) );
871 }
872 //--------------------------------------------------------------------------------------------------------------------------------------------------------
873 /** Check if points are colinear. */
874 inline bool geomi_pts_colinear(node_idx_t pi1, node_idx_t pi2, node_idx_t pi3) const {
875 return geomr_pts_colinear(get_pnt(pi1), get_pnt(pi2), get_pnt(pi3));
876 }
877 //--------------------------------------------------------------------------------------------------------------------------------------------------------
878 /** Check if points are colinear */
879 inline bool geomr_pts_colinear(fvec3_t p1, fvec3_t p2, fvec3_t p3) const {
880 return geomr_pnt_zero(vec3_cross_product(vec3_diff(p1, p2), vec3_diff(p1, p3)));
881 }
882 //--------------------------------------------------------------------------------------------------------------------------------------------------------
883 /** Check if points are coplanar */
884 inline bool geomi_pts_coplanar(const node_idx_list_t& pnt_list) const {
885 if (pnt_list.size() > 3) {
886 if ( !(geomi_pts_coplanar(pnt_list[0], pnt_list[1], pnt_list[2], pnt_list[3])))
887 return false;
888 for(decltype(pnt_list.size()) i=4; i<pnt_list.size(); i++)
889 if ( !(geomi_pts_coplanar(pnt_list[0], pnt_list[1], pnt_list[2], pnt_list[i])))
890 return false;
891 }
892 return true;
893 }
894 //--------------------------------------------------------------------------------------------------------------------------------------------------------
895 /** Check if points are coplanar */
896 inline bool geomi_pts_coplanar(node_idx_t pi1, node_idx_t pi2, node_idx_t pi3, node_idx_t pi4) const {
897 return geomr_pts_coplanar(get_pnt(pi1), get_pnt(pi2), get_pnt(pi3), get_pnt(pi4));
898 }
899 //--------------------------------------------------------------------------------------------------------------------------------------------------------
900 /** Check if points are coplanar */
901 inline bool geomr_pts_coplanar(fvec3_t p1, fvec3_t p2, fvec3_t p3, fvec3_t p4) const {
902 return (std::abs(vec3_scalar_triple_product(vec3_diff(p3, p1), vec3_diff(p2, p1), vec3_diff(p4, p1))) < eps);
903 }
904 //--------------------------------------------------------------------------------------------------------------------------------------------------------
905 /** Compute intersection of a segment and a triangle */
906 inline bool geomr_seg_tri_intersection(fvec3_t tp1, fvec3_t tp2, fvec3_t tp3, fvec3_t sp1, fvec3_t sp2) const {
907 // MJR TODO NOTE <2024-07-11T16:08:16-0500> geomr_seg_tri_intersection: implement
908 // MJR TODO NOTE <2024-07-11T16:08:27-0500> geomr_seg_tri_intersection: Should this be a bool or an enum?
909 static_assert(false, "geomr_seg_tri_intersection: Not yet implemented!");
910 return (tp1[0]+tp2[0]+tp3[0]+sp1[0]+sp2[0]>1.0);
911 }
912 //@}
913
914 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
915 /** @name Class utilities. */
916 //@{
917 void clear() {
918 last_point_idx = -1;
919 last_point_new = true;
920 pnt_to_node_idx_map.clear();
921 node_idx_to_node_data.clear();
922 data_name_to_data_idx_lst.clear();
923 node_data_to_pnt.clear();
924 cell_lst.clear();
925 uniq_cell_lst.clear();
926 last_cell_new = true;
927 last_cell_stat = cell_stat_t::GOOD;
928 }
929 //@}
930
931 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
932 /** @name Cells. */
933 //@{
934 //--------------------------------------------------------------------------------------------------------------------------------------------------------
935 /** Cell segment/face/etc structure. */
936 typedef cell_lst_t cell_structure_t;
937 //--------------------------------------------------------------------------------------------------------------------------------------------------------
938 /** Number of cells. */
939 inline int num_cells() const {
940 return static_cast<int>(cell_lst.size());
941 }
942 //--------------------------------------------------------------------------------------------------------------------------------------------------------
943 /** Convert cell_kind_t value to a list of vertexes. */
944 inline cell_structure_t& cell_kind_to_structure(cell_kind_t cell_kind, int dimension) const {
945 // MJR TODO NOTE <2024-07-11T16:06:42-0500> cell_kind_to_structure: make sure polygons are all oriented correctly
946 int logical_dim = cell_kind_to_dimension(cell_kind);
947 if ( (dimension < 0) || (dimension > logical_dim) )
948 dimension = logical_dim;
949 if (logical_dim > 3) {
950 std::cerr << "MR_cell_cplx API USAGE ERROR: Maximum supported dimension is 3" << std::endl;
951 exit(0);
952 }
953 int idx = 0;
954 switch(cell_kind) {
955 case cell_kind_t::POINT: idx = 0; break;
956 case cell_kind_t::SEGMENT: idx = 1; break;
957 case cell_kind_t::TRIANGLE: idx = 2; break;
958 case cell_kind_t::QUAD: idx = 3; break;
959 case cell_kind_t::PYRAMID: idx = 4; break;
960 case cell_kind_t::HEXAHEDRON: idx = 5; break;
961 }
962 static std::vector<std::vector<cell_structure_t>> cst = {{{{0}}, // vertex
963 {{0},{1}}, // segment
964 {{0},{1},{2}}, // triangle
965 {{0},{1},{2},{3}}, // Quad
966 {{0},{1},{2},{3}, // Pyramid: Base vertexes
967 {4}}, // Pyramid: Tip vertex
968 {{0},{1},{2},{3}, // Hexahedron: Back vertexes
969 {4},{5},{6},{7}}}, // Hexahedron: Front vertexes
970 {{}, // vertex
971 {{0,1}}, // segment
972 {{0,1},{1,2},{2,0}}, // triangle
973 {{0,1},{1,2},{2,3},{3,0}}, // Quad
974 {{0,1},{1,2},{2,3},{3,0}, // Pyramid: Base segments
975 {0,4},{1,4},{2,4},{3,4}}, // Pyramid: Side sets
976 {{0,1},{1,2},{2,3},{3,0}, // Hexahedron: Back segments
977 {4,5},{5,6},{6,7},{7,4}, // Hexahedron: Front segments
978 {0,4},{1,5},{2,6},{3,7}}}, // Hexahedron: Back to front segments
979 {{}, // vertex
980 {}, // segment
981 {{0,1,2}}, // triangle
982 {{0,1,2,3}}, // Quad
983 {{0,1,2,3}, // Pyramid: Base face
984 {0,1,4}, // Pyramid: Back face
985 {1,2,4}, // Pyramid: Left face
986 {2,3,4}, // Pyramid: Front face
987 {3,0,4}}, // Pyramid: Right faces
988 {{0,1,2,3}, // Hexahedron: Back face
989 {4,5,6,7}, // Hexahedron: Front face
990 {0,3,7,4}, // Hexahedron: Left face
991 {2,3,7,6}, // Hexahedron: Top face
992 {1,2,6,5}, // Hexahedron: Right face
993 {0,1,4,5}}}, // Hexahedron: Bottom face
994 {{}, // vertex
995 {}, // segment
996 {}, // triangle
997 {}, // Quad
998 {{0,1,2,3,4}}, // Pyramid
999 {{0,1,2,3,4,5,6,7}}}}; // Hexahedron
1000 return (cst[dimension][idx]);
1001 }
1002 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1003 /** Convert cell_kind_t value to the logical dimension of the cell.
1004
1005 Note that the vertexes of a cell_kind_t::QUAD cell might not be coplainer; however, the cell is logically a 2D cell. Similarly, the vertexes of a
1006 cell_kind_t::PYRAMID might be coplainer; however, the cell is logically a 3D cell. */
1007 inline int cell_kind_to_dimension(cell_kind_t cell_kind) const {
1008 switch(cell_kind) {
1009 case cell_kind_t::POINT: return (0); break;
1010 case cell_kind_t::SEGMENT: return (1); break;
1011 case cell_kind_t::TRIANGLE: return (2); break;
1012 case cell_kind_t::QUAD: return (2); break;
1013 case cell_kind_t::PYRAMID: return (3); break;
1014 case cell_kind_t::HEXAHEDRON: return (3); break;
1015 }
1016 return -1;
1017 }
1018 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1019 /** Convert cell_kind_t value to the number of points required for cell type. */
1020 inline int cell_kind_to_req_pt_cnt(cell_kind_t cell_kind) const {
1021 switch(cell_kind) {
1022 case cell_kind_t::POINT: return (1); break;
1023 case cell_kind_t::SEGMENT: return (2); break;
1024 case cell_kind_t::TRIANGLE: return (3); break;
1025 case cell_kind_t::QUAD: return (4); break;
1026 case cell_kind_t::PYRAMID: return (5); break;
1027 case cell_kind_t::HEXAHEDRON: return (8); break;
1028 }
1029 return -1;
1030 }
1031 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1032 /** Convert cell_kind_t value to the VTK cell type integer. */
1033 inline int cell_kind_to_vtk_type(cell_kind_t cell_kind) const {
1034 switch(cell_kind) {
1035 case cell_kind_t::POINT: return ( 1); break;
1036 case cell_kind_t::SEGMENT: return ( 3); break;
1037 case cell_kind_t::TRIANGLE: return ( 5); break;
1038 case cell_kind_t::QUAD: return ( 9); break;
1039 case cell_kind_t::HEXAHEDRON: return (12); break;
1040 case cell_kind_t::PYRAMID: return (14); break;
1041 }
1042 return -1;
1043 }
1044 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1045 /** Convert cell_kind_t value to the VTK cell type integer. */
1046 inline std::string cell_kind_to_string(cell_kind_t cell_kind) const {
1047 switch(cell_kind) {
1048 case cell_kind_t::POINT: return ("POINT" ); break;
1049 case cell_kind_t::SEGMENT: return ("SEGMENT" ); break;
1050 case cell_kind_t::TRIANGLE: return ("TRIANGLE" ); break;
1051 case cell_kind_t::QUAD: return ("QUAD" ); break;
1052 case cell_kind_t::HEXAHEDRON: return ("HEXAHEDRON"); break;
1053 case cell_kind_t::PYRAMID: return ("PYRAMID" ); break;
1054 }
1055 return ""; // Never get here.
1056 }
1057 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1058 /** Convert number of points in a cell to cell_kind_t. */
1059 inline cell_kind_t req_pt_cnt_to_cell_kind(std::vector<int>::size_type node_count) const {
1060 switch(node_count) {
1061 case 1: return (cell_kind_t::POINT); break;
1062 case 2: return (cell_kind_t::SEGMENT); break;
1063 case 3: return (cell_kind_t::TRIANGLE); break;
1064 case 4: return (cell_kind_t::QUAD); break;
1065 case 5: return (cell_kind_t::PYRAMID); break;
1066 case 8: return (cell_kind_t::HEXAHEDRON); break;
1067 }
1068 return (cell_kind_t::POINT); // Never get here
1069 }
1070 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1071 /** Convert cell_stat_t to a bool (true if GOOD). */
1072 inline bool cell_stat_is_good(cell_stat_t cell_stat) const {
1073 return (cell_stat == cell_stat_t::GOOD);
1074 }
1075 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1076 /** Convert cell_stat_t to a bool (true if BAD). */
1077 inline bool cell_stat_is_bad(cell_stat_t cell_stat) const {
1078 return (cell_stat != cell_stat_t::GOOD);
1079 }
1080 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1081 /** Convert cell_stat_t enum value to a string. */
1082 inline std::string cell_stat_to_string(cell_stat_t cell_stat) const {
1083 switch(cell_stat) {
1084 case cell_stat_t::GOOD: return (std::string("GOOD")); break;
1085 case cell_stat_t::TOO_FEW_PNT: return (std::string("TOO_FEW_PNT")); break;
1086 case cell_stat_t::TOO_MANY_PNT: return (std::string("TOO_MANY_PNT")); break;
1087 case cell_stat_t::NEG_PNT_IDX: return (std::string("NEG_PNT_IDX")); break;
1088 case cell_stat_t::BIG_PNT_IDX: return (std::string("BIG_PNT_IDX")); break;
1089 case cell_stat_t::DUP_PNT: return (std::string("DUP_PNT")); break;
1090 case cell_stat_t::DIM_LOW: return (std::string("DIM_LOW")); break;
1091 case cell_stat_t::BAD_EDGEI: return (std::string("BAD_EDGEI")); break;
1092 case cell_stat_t::BAD_FACEI: return (std::string("BAD_FACEI")); break;
1093 case cell_stat_t::FACE_BENT: return (std::string("FACE_BENT")); break;
1094 case cell_stat_t::CONCAVE: return (std::string("CONCAVE")); break;
1095 }
1096 return std::string("ERROR");
1097 }
1098 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1099 /** Perform cell vertex checks
1100
1101 @param cell_kind The type of the potential cell.
1102 @param cell_verts The vertexes of the potential cell. */
1103 inline cell_stat_t check_cell_vertexes(cell_kind_t cell_kind, cell_verts_t cell_verts) const {
1104 // Check number of points
1105 std::vector<int>::size_type num_verts = cell_verts.size();
1106 std::vector<int>::size_type req_num_pts = cell_kind_to_req_pt_cnt(cell_kind);
1107 if (num_verts < req_num_pts)
1108 return cell_stat_t::TOO_FEW_PNT;
1109 if (num_verts > req_num_pts)
1110 return cell_stat_t::TOO_MANY_PNT;
1111 // Check for negative point index
1112 if (std::any_of(cell_verts.cbegin(), cell_verts.cend(), [](node_idx_t i) { return (i < 0); }))
1113 return cell_stat_t::NEG_PNT_IDX;
1114 // Check for too big point index
1115 if (std::any_of(cell_verts.cbegin(), cell_verts.cend(), [this](node_idx_t i) { return (i >= node_count()); }))
1116 return cell_stat_t::BIG_PNT_IDX;
1117 // Check for duplicate points
1118 if (num_verts > 1) {
1119 std::set<node_idx_t> cell_verts_sorted;
1120 for(node_idx_t pnt_idx: cell_verts) {
1121 if (cell_verts_sorted.contains(pnt_idx))
1122 return cell_stat_t::DUP_PNT;
1123 cell_verts_sorted.insert(pnt_idx);
1124 }
1125 }
1126 // Return GOOD
1127 return cell_stat_t::GOOD;
1128 }
1129 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1130 /** Perform cell dimension (2D cells must not be colinear & 3D cells must not be coplainer)
1131
1132 @warning This function assumes that check_cell_vertexes() would have returned GOOD for the cell being checked!
1133
1134 @warning Will not detect degenerate cell_kind_t::SEGMENT as that is handled by check_cell_vertexes.
1135
1136 @param cell_kind The type of the potential cell.
1137 @param cell_verts The vertexes of the potential cell. */
1138 inline cell_stat_t check_cell_dimension(cell_kind_t cell_kind, cell_verts_t cell_verts) const {
1139 if (cell_kind == cell_kind_t::TRIANGLE) {
1140 if (geomi_pts_colinear(cell_verts[0], cell_verts[1], cell_verts[2]))
1141 return cell_stat_t::DIM_LOW;
1142 } else if (cell_kind == cell_kind_t::QUAD) {
1143 if (geomi_pts_colinear(cell_verts[0], cell_verts[1], cell_verts[2], cell_verts[3]))
1144 return cell_stat_t::DIM_LOW;
1145 } else if (cell_kind == cell_kind_t::HEXAHEDRON) {
1146 if ( geomi_pts_coplanar(cell_verts))
1147 return cell_stat_t::DIM_LOW;
1148 } else if (cell_kind == cell_kind_t::PYRAMID) {
1149 if ( geomi_pts_coplanar(cell_verts))
1150 return cell_stat_t::DIM_LOW;
1151 }
1152 return cell_stat_t::GOOD;
1153 }
1154 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1155 /** Checks that cell edges have expected intersections.
1156
1157 Checks that every pair of cell edges has the correct intersection type (i.e. intersect in a single vertex or an empty intersection). If a bad
1158 intersection is detected, then a cell_stat_t of cell_stat_t::BAD_EDGEI will be returned. Otherwise cell_stat_t::GOOD will be returned.
1159
1160 @warning This function assumes that check_cell_vertexes() would have returned GOOD for the cell being checked!
1161
1162 Note that this function will detect some conditions caught by other checks. For example, this function will return cell_stat_t::BAD_EDGEI for
1163 cell_kind_t::SEG cells for which check_cell_dimension() returns cell_stat_t::DUP_PNT. Note the cell_stat_t values are different depending upon which
1164 check function is called.
1165
1166 @param cell_kind The type of the potential cell.
1167 @param cell_verts The vertexes of the potential cell. */
1168 inline cell_stat_t check_cell_edge_intersections(cell_kind_t cell_kind, cell_verts_t cell_verts) const {
1169 cell_structure_t& segs = cell_kind_to_structure(cell_kind, 1);
1170 if ( !(segs.empty())) {
1171 for(decltype(segs.size()) i=0; i<segs.size()-1; i++) {
1172 for(decltype(segs.size()) j=i+1; j<segs.size(); j++) {
1173 std::set<node_idx_t> points_sorted;
1174 points_sorted.insert(segs[i][0]);
1175 points_sorted.insert(segs[i][1]);
1176 points_sorted.insert(segs[j][0]);
1177 points_sorted.insert(segs[j][1]);
1178 auto it = geomi_seg_isect_type(cell_verts[segs[i][0]], cell_verts[segs[i][1]], cell_verts[segs[j][0]], cell_verts[segs[j][1]]);
1179 if(points_sorted.size() == 4) {
1180 if (it != seg_isect_t::C0_EMPTY)
1181 return cell_stat_t::BAD_EDGEI;
1182 } else if(points_sorted.size() == 3) {
1183 if (it != seg_isect_t::C1_VERTEX1)
1184 return cell_stat_t::BAD_EDGEI;
1185 } else {
1186 return cell_stat_t::BAD_EDGEI;
1187 }
1188 }
1189 }
1190 }
1191 return cell_stat_t::GOOD;
1192 }
1193 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1194 /** Checks that cell faces have expected intersections.
1195
1196 @param cell_kind The type of the potential cell.
1197 @param cell_verts The vertexes of the potential cell. */
1198 inline cell_stat_t check_cell_face_intersections(cell_kind_t cell_kind, cell_verts_t cell_verts) const {
1199 // MJR TODO NOTE <2024-08-02T09:42:38-0500> check_cell_face_intersections: Implement
1200 if (cell_kind == cell_kind_t::HEXAHEDRON) {
1201 if ( geomi_pts_coplanar(cell_verts))
1202 return cell_stat_t::DIM_LOW;
1203 } else if (cell_kind == cell_kind_t::PYRAMID) {
1204 if ( geomi_pts_coplanar(cell_verts))
1205 return cell_stat_t::DIM_LOW;
1206 }
1207 return cell_stat_t::GOOD;
1208 }
1209 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1210 /** Check if the vertexes of a cell are coplainer.
1211
1212 Returns cell_stat_t::FACE_BENT if the cell is *not* plainer, and cell_stat_t::GOOD if it is.
1213
1214 Note that most visualization packages are tolerant of non-plainer faces, and can render them just fine. For example, most applications will
1215 automatically split a non-plainer cell_kind_t::QUAD into two triangles. That said such faces can be an issue when using a tessellation for
1216 computation. Such issues are most severe for things like FEM, but they can also show up for common visualization computations like the extraction
1217 of level curves/surfaces.
1218
1219 @param cell_kind The type of the potential cell.
1220 @param cell_verts The vertexes of the potential cell. */
1221 inline cell_stat_t check_cell_faces_plainer(cell_kind_t cell_kind, cell_verts_t cell_verts) const {
1222 const cell_structure_t& face_structures = cell_kind_to_structure(cell_kind, 2);
1223 for(auto face_structure: face_structures) {
1224 cell_verts_t face;
1225 for(auto idx: face_structure)
1226 face.push_back(cell_verts[idx]);
1227 if ( !(geomi_pts_coplanar(face)))
1228 return cell_stat_t::FACE_BENT;
1229 }
1230 return cell_stat_t::GOOD;
1231 }
1232 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1233 /** Add parts of a cell of the specified dimension.
1234
1235 Example 1: a valid cell_kind_t::PYRAMID cell added wth dimension == 2 results in three cell_kind_t::TRIANGLE cells and one cell_kind_t::QUAD cell.
1236
1237 Example 2: a valid cell_kind_t::PYRAMID cell added wth dimension == 3 results in one cell_kind_t::PYRAMID cell added
1238
1239 Example 3: a valid cell_kind_t::PYRAMID cell added wth dimension == 1 results in 8 cell_kind_t::SEG cells added
1240
1241 @param new_cell_kind The type of the new cell
1242 @param new_cell_verts The vertexes of the new cell
1243 @param dimension The dimension of the parts to add.
1244 @return Number of cells added */
1245 inline int add_cell(cell_kind_t new_cell_kind, cell_verts_t new_cell_verts, int dimension) {
1246 int num_added = 0;
1247 if ( (dimension < 0) || (dimension >= cell_kind_to_dimension(new_cell_kind)) ) {
1248 if (add_cell(new_cell_kind, new_cell_verts))
1249 num_added++;
1250 } else { // We need to break the cell up into lower dimensional bits, and add the bits.
1251 const cell_structure_t& cell_parts = cell_kind_to_structure(new_cell_kind, dimension);
1252 for(auto cell_part: cell_parts) {
1253 cell_verts_t newer_cell_verts;
1254 for(auto idx: cell_part)
1255 newer_cell_verts.push_back(new_cell_verts[idx]);
1256 if (add_cell(req_pt_cnt_to_cell_kind(newer_cell_verts.size()), newer_cell_verts))
1257 num_added++;
1258 }
1259 }
1260 return num_added;
1261 }
1262 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1263 /** Add a cell
1264 @param new_cell_kind The type of the new cell
1265 @param new_cell_verts The vertexes of the new cell
1266 @return A boolean indicateing success
1267 @retval true The cell was added or had been added previously
1268 @retval false The cell could not be added (because of a failed geometric check) */
1269 inline bool add_cell(cell_kind_t new_cell_kind, cell_verts_t new_cell_verts) {
1270 // Check vertexes if required
1271 if constexpr (chk_cell_vertexes) {
1272 last_cell_stat = check_cell_vertexes(new_cell_kind, new_cell_verts);
1273 if (cell_stat_is_bad(last_cell_stat))
1274 return false;
1275 }
1276 // Check dimension if required
1277 if constexpr (chk_cell_dimension) {
1278 last_cell_stat = check_cell_dimension(new_cell_kind, new_cell_verts);
1279 if (cell_stat_is_bad(last_cell_stat))
1280 return false;
1281 }
1282 // Check edges
1283 if constexpr (chk_cell_edges) {
1284 last_cell_stat = check_cell_edge_intersections(new_cell_kind, new_cell_verts);
1285 if (cell_stat_is_bad(last_cell_stat))
1286 return false;
1287 }
1288 // Geom was good or we didn't need to check.
1289 if constexpr (chk_cell_unique) {
1290 cell_verts_t new_cell_sorted_verts = new_cell_verts;
1291 std::sort(new_cell_sorted_verts.begin(), new_cell_sorted_verts.end());
1292 if (uniq_cell_lst.contains(new_cell_sorted_verts)) {
1293 last_cell_new = false;
1294 } else {
1295 last_cell_new = true;
1296 cell_lst.push_back(new_cell_verts);
1297 uniq_cell_lst.insert(new_cell_sorted_verts);
1298 }
1299 } else {
1300 cell_lst.push_back(new_cell_verts);
1301 }
1302 return true;
1303 }
1304 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1305 /** Retruns the status of the last cell given to the add_cell() method.
1306 If (chk_cell_vertexes || chk_cell_dimension | chk_cell_edges) is true, this value is updated each time add_cell() is called.
1307 Otherwise its value is always cell_stat_t::GOOD. */
1308 inline cell_stat_t status_of_last_cell_added() const {
1309 return last_cell_stat;
1310 }
1311 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1312 /** Retruns true if the last cell given to the add_cell() was new and the return value of status_of_last_cell_added() is cell_stat_t::GOOD.
1313 If chk_cell_unique is true, this value is updated each time add_cell() is called. Otherwise its value is always true. */
1314 inline bool last_cell_added_was_new() const {
1315 return last_cell_new;
1316 }
1317 //@}
1318
1319 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1320 /** @name I/O. */
1321 //@{
1322 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1323 /** Type returned by I/O functions */
1324 typedef int io_result;
1325 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1326 /** Dump to an XML VTK unstructured grid file
1327
1328 Note: A point vector data set named "NORMALS" will be used for normal vectors.
1329
1330 @param file_name The name of the output file
1331 @param description This is included as a file comment at the start of the file.
1332 @return 0 if everything worked, and non-zero otherwise */
1333 io_result write_xml_vtk(std::string file_name, std::string description) {
1334 /* Check that we have data */
1335 if (node_count() <= 0) {
1336 std::cout << "ERROR(write_xml_vtk): No points!" << std::endl;
1337 return 1;
1338 }
1339 if (num_cells() <= 0) {
1340 std::cout << "ERROR(write_xml_vtk): No cells!" << std::endl;
1341 return 2;
1342 }
1343 /* Looks like we have data. Let's open our file */
1344 std::ofstream out_stream;
1345 out_stream.open(file_name, std::ios::out | std::ios::binary | std::ios::trunc);
1346 if (out_stream.is_open()) {
1347 out_stream.imbue(std::locale::classic());
1348 } else {
1349 std::cout << "ERROR(write_xml_vtk): Could not open file!" << std::endl;
1350 return 3;
1351 }
1352 out_stream << "<VTKFile type='UnstructuredGrid' version='0.1' byte_order='LittleEndian'>" << std::endl;
1353 out_stream << "<!-- " << description << " -->" << std::endl;
1354 out_stream << " <UnstructuredGrid>" << std::endl;
1355 out_stream << " <Piece NumberOfPoints='" << node_count() << "' NumberOfCells='" << num_cells() << "'>" << std::endl;
1356 if ( !(data_name_to_data_idx_lst.empty())) {
1357 std::string scalars_attr_value, vectors_attr_value, normals_attr_value;
1358 for (auto& kv : data_name_to_data_idx_lst)
1359 if (kv.second.size()==1)
1360 scalars_attr_value += (scalars_attr_value.empty() ? "" : " ") + kv.first;
1361 else
1362 if (kv.first == "NORMALS")
1363 normals_attr_value = "NORMALS";
1364 else
1365 vectors_attr_value += (vectors_attr_value.empty() ? "" : " ") + kv.first;
1366 out_stream << " <PointData";
1367 if ( !(scalars_attr_value.empty()))
1368 out_stream << " Scalars='" << scalars_attr_value << "'";
1369 if ( !(normals_attr_value.empty()))
1370 out_stream << " Normals='" << normals_attr_value << "'";
1371 if ( !(vectors_attr_value.empty()))
1372 out_stream << " Vectors='" << vectors_attr_value << "'";
1373 out_stream << ">" << std::endl;
1374 for (auto& kv : data_name_to_data_idx_lst) {
1375 out_stream << " <DataArray Name='" << kv.first << "' type='Float64' format='ascii' NumberOfComponents='" << kv.second.size() << "'>" << std::endl;
1376 out_stream << " ";
1377 for (const auto& dv : node_idx_to_node_data) {
1378 for (auto& idx : kv.second)
1379 out_stream << std::setprecision(10) << scalar_from_node_data(idx, dv) << " ";
1380 }
1381 out_stream << std::endl << " </DataArray>" << std::endl;
1382 }
1383 out_stream << " </PointData>" << std::endl;
1384 }
1385 out_stream << " <Points>" << std::endl;
1386 out_stream << " <DataArray Name='Points' type='Float64' format='ascii' NumberOfComponents='3'>" << std::endl;
1387 for(node_idx_t pnt_idx=0; pnt_idx<static_cast<node_idx_t>(node_idx_to_node_data.size()); pnt_idx++) {
1388 fvec3_t pnt = get_pnt(pnt_idx);
1389 out_stream << " " << std::setprecision(10) << pnt[0] << " " << pnt[1] << " " << pnt[2] << std::endl;
1390 }
1391 out_stream << " </DataArray>" << std::endl;
1392 out_stream << " </Points>" << std::endl;
1393 out_stream << " <Cells>" << std::endl;
1394 out_stream << " <DataArray type='Int32' Name='connectivity' format='ascii'>" << std::endl;
1395 for(auto& poly: cell_lst) {
1396 out_stream << " " ;
1397 for(auto& vert: poly)
1398 out_stream << vert << " ";
1399 out_stream << std::endl;
1400 }
1401 out_stream << " </DataArray>" << std::endl;
1402 out_stream << " <DataArray type='Int32' Name='offsets' format='ascii'>" << std::endl;
1403 out_stream << " ";
1404 std::vector<int>::size_type j = 0;
1405 for(auto& poly: cell_lst) {
1406 j += poly.size();
1407 out_stream << j << " ";
1408 }
1409 out_stream << std::endl;
1410 out_stream << " </DataArray>" << std::endl;
1411 out_stream << " <DataArray type='Int8' Name='types' format='ascii'>" << std::endl;
1412 out_stream << " ";
1413 for(auto& poly: cell_lst)
1414 out_stream << cell_kind_to_vtk_type(req_pt_cnt_to_cell_kind(poly.size())) << " ";
1415 out_stream << std::endl;
1416 out_stream << " </DataArray>" << std::endl;
1417 out_stream << " </Cells>" << std::endl;
1418 out_stream << " </Piece>" << std::endl;
1419 out_stream << " </UnstructuredGrid>" << std::endl;
1420 out_stream << "</VTKFile>" << std::endl;
1421
1422 /* Final newline */
1423 out_stream << std::endl;
1424 out_stream.close();
1425 return 0;
1426 }
1427 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1428 /** Dump to a legacy VTK UNSTRUCTUREDGRID file.
1429
1430 Note: A point vector data set named "NORMALS" will be used for normal vectors.
1431 Note: A point vector data set named "COLORS" will be used as COLOR_SCALARS. Components should be in [0.0, 1.0].
1432
1433 @param file_name The name of the output file
1434 @param description This is the file description.
1435 @return 0 if everything worked, and non-zero otherwise */
1436 io_result write_legacy_vtk(std::string file_name, std::string description) {
1437 /* Check that we have data */
1438 if (node_count() <= 0) {
1439 std::cout << "ERROR(write_legacy_vtk): No points!" << std::endl;
1440 return 1;
1441 }
1442 if (num_cells() <= 0) {
1443 std::cout << "ERROR(write_legacy_vtk): No cells!" << std::endl;
1444 return 2;
1445 }
1446 /* Looks like we have data. Let's open our file */
1447 std::ofstream out_stream;
1448 out_stream.open(file_name, std::ios::out | std::ios::binary | std::ios::trunc);
1449 if (out_stream.is_open()) {
1450 out_stream.imbue(std::locale::classic());
1451 } else {
1452 std::cout << "ERROR(write_legacy_vtk): Could not open file!" << std::endl;
1453 return 3;
1454 }
1455 /* Dump the header */
1456 out_stream << "# vtk DataFile Version 3.0" << std::endl;
1457 out_stream << description << std::endl;
1458 out_stream << "ASCII" << std::endl;
1459 out_stream << "DATASET UNSTRUCTURED_GRID" << std::endl;
1460 /* Dump the points */
1461 out_stream << "POINTS " << node_count() << " double" << std::endl;
1462 //for (const auto& pnt : node_idx_to_pnt)
1463 for(node_idx_t pnt_idx=0; pnt_idx<static_cast<node_idx_t>(node_idx_to_node_data.size()); pnt_idx++) {
1464 fvec3_t pnt = get_pnt(pnt_idx);
1465 out_stream << std::setprecision(10) << pnt[0] << " " << pnt[1] << " " << pnt[2] << std::endl;
1466 }
1467 /* Dump the cell data */
1468 std::vector<int>::size_type total_cells_ints = 0;
1469 for(auto& cell: cell_lst)
1470 total_cells_ints += (1+cell.size());
1471 out_stream << "CELLS " << num_cells() << " " << total_cells_ints << std::endl;
1472 for(auto& poly: cell_lst) {
1473 out_stream << poly.size() << " ";
1474 for(auto& vert: poly) {
1475 out_stream << vert << " ";
1476 }
1477 out_stream << std::endl;
1478 }
1479 out_stream << "CELL_KINDS " << num_cells() << std::endl;
1480 for(auto& poly: cell_lst)
1481 out_stream << cell_kind_to_vtk_type(req_pt_cnt_to_cell_kind(poly.size())) << std::endl;
1482 /* Dump point scalar data */
1483 if (data_name_to_data_idx_lst.size() > 0) {
1484 out_stream << "POINT_DATA " << node_count() << std::endl;
1485 if (named_scalar_datasets_count() > 0) {
1486 for (auto& kv : data_name_to_data_idx_lst) {
1487 if (kv.second.size() == 1) {
1488 out_stream << "SCALARS " << kv.first << " double 1" << std::endl;
1489 out_stream << "LOOKUP_TABLE default" << std::endl;
1490 for (const auto& dv : node_idx_to_node_data) {
1491 uft_t v = scalar_from_node_data(kv.second[0], dv);
1492 out_stream << std::setprecision(10) << v << std::endl;
1493 }
1494 }
1495 }
1496 }
1497 if (named_vector_datasets_count() > 0) {
1498 for (auto& kv : data_name_to_data_idx_lst) {
1499 if (kv.second.size() == 3) { // Legacy format only supports three element vectors
1500 if ("NORMALS" == kv.first) {
1501 out_stream << "NORMALS " << kv.first << " double" << std::endl;
1502 } else if ("COLORS" == kv.first) {
1503 out_stream << "COLOR_SCALARS " << kv.first << " 3" << std::endl;
1504 } else {
1505 out_stream << "VECTORS " << kv.first << " double" << std::endl;
1506 }
1507 for (const auto& dv : node_idx_to_node_data) {
1508 vdat_t v = vector_from_node_data(kv.second, dv);
1509 out_stream << std::setprecision(10) << v[0] << " " << v[1] << " " << v[2] << std::endl;
1510 }
1511 }
1512 }
1513 }
1514 }
1515 /* Final newline */
1516 out_stream << std::endl;
1517 out_stream.close();
1518 return 0;
1519 }
1520 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1521 /** Dump object to STDOUT.
1522
1523 @param max_num_print Maximum number of points/cells to print. Use 0 to print all points/cells. */
1524 void dump_cplx(int max_num_print) const {
1525 std::cout << "Meta Data" << std::endl;
1526 std::cout << " Points .................. " << node_count() << std::endl;
1527 std::cout << " Data Scalars Per Point .. " << node_idx_to_node_data.size() << std::endl;
1528 std::cout << " Named Data Sets ......... " << named_datasets_count() << std::endl;
1529 std::cout << " Scalar Data Sets ...... " << named_scalar_datasets_count() << std::endl;
1530 std::cout << " Vector Data Sets ...... " << named_vector_datasets_count() << std::endl;
1531 std::cout << " Cells ................... " << num_cells() << std::endl;
1532 if (node_count() > 0) {
1533 int num_nodes_printed = 0;
1534 std::cout << "NODES BEGIN (" << node_count() << ")" << std::endl;
1535 for(node_idx_t pnt_idx = 0; pnt_idx<node_count(); ++pnt_idx) {
1536 std::cout << " " << pnt_idx << ": " << node_to_string(pnt_idx) << std::endl;
1537 num_nodes_printed++;
1538 if ((max_num_print > 0) && (num_nodes_printed >= max_num_print)) {
1539 std::cout << " Maximum number of nodes reached. Halting tree dump." << std::endl;
1540 break;
1541 }
1542 }
1543 std::cout << "NODES END" << std::endl;
1544 }
1545 if (num_cells() > 0) {
1546 int num_cells_printed = 0;
1547 std::cout << "CELLS BEGIN (" << num_cells() << ")" << std::endl;
1548 for(int i=0; i<num_cells(); i++) {
1549 std::cout << " ";
1550 for(auto& vert: cell_lst[i])
1551 std::cout << vert << " ";
1552 std::cout << " " << cell_kind_to_string(req_pt_cnt_to_cell_kind(cell_lst[i].size())) << std::endl;
1553 num_cells_printed++;
1554 if ((max_num_print > 0) && (num_cells_printed >= max_num_print)) {
1555 std::cout << " Maximum number of cells reached. Halting tree dump." << std::endl;
1556 break;
1557 }
1558 }
1559 std::cout << "CELLS END" << std::endl;
1560 }
1561 }
1562 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1563 /** Dump to a PLY file
1564
1565 Note: If one of the point vector data sets is named "NORMALS", then it will be used for normal vectors.
1566 Note: If one of the point vector data sets is named "COLORS", then it will be used for color data. Components should be in [0.0, 1.0].
1567
1568 @param file_name The name of the output file
1569 @param description For legacy files, this is the file description. For XML files this is included as a file comment at the start of the file.
1570 @return 0 if everything worked, and non-zero otherwise */
1571 io_result write_ply(std::string file_name, std::string description) {
1572 /* Check that we have data */
1573 if (node_count() <= 0) {
1574 std::cout << "ERROR(write_ply): No points!" << std::endl;
1575 return 1;
1576 }
1577 if (num_cells() <= 0) {
1578 std::cout << "ERROR(write_ply): No cells!" << std::endl;
1579 return 2;
1580 }
1581 for(auto& poly: cell_lst) {
1582 cell_kind_t cell_kind = req_pt_cnt_to_cell_kind(poly.size());
1583 if ( !((cell_kind == cell_kind_t::TRIANGLE) || (cell_kind == cell_kind_t::QUAD))) {
1584 std::cout << "ERROR(write_ply): Cells must all be 2D (triangles or quads)!" << std::endl;
1585 return 2;
1586 }
1587 }
1588 /* Looks like we have data. Let's open our file */
1589 std::ofstream out_stream;
1590 out_stream.open(file_name, std::ios::out | std::ios::binary | std::ios::trunc);
1591 if (out_stream.is_open()) {
1592 out_stream.imbue(std::locale::classic());
1593 } else {
1594 std::cout << "ERROR(write_ply): Could not open file!" << std::endl;
1595 return 3;
1596 }
1597 bool have_colors_data = data_name_to_data_idx_lst.contains("COLORS");
1598 bool have_normal_data = data_name_to_data_idx_lst.contains("NORMALS");
1599 out_stream << "ply" << std::endl;
1600 out_stream << "format ascii 1.0" << std::endl;
1601 out_stream << "comment software: Mitch Richling's MR_rect_tree package" << std::endl;
1602 out_stream << "comment note: " << description << std::endl;
1603 out_stream << "element vertex " << node_count() << std::endl;
1604 out_stream << "property float x" << std::endl;
1605 out_stream << "property float y" << std::endl;
1606 out_stream << "property float z" << std::endl;
1607 if (have_colors_data) {
1608 out_stream << "property uchar red" << std::endl;
1609 out_stream << "property uchar green" << std::endl;
1610 out_stream << "property uchar blue" << std::endl;
1611 }
1612 if (have_normal_data) {
1613 out_stream << "property float nx" << std::endl;
1614 out_stream << "property float ny" << std::endl;
1615 out_stream << "property float nz" << std::endl;
1616 }
1617 out_stream << "element face " << num_cells() << std::endl; // May need to be adjusted if cells are not triangles..
1618 out_stream << "property list uchar int vertex_index" << std::endl;
1619 out_stream << "end_header" << std::endl;
1620 // Dump Vertex Data
1621 for (int i=0; i<node_count(); i++) {
1622 fvec3_t pnt = get_pnt(i);
1623 out_stream << std::setprecision(10) << pnt[0] << " " << pnt[1] << " " << pnt[2];
1624 if (have_colors_data) {
1625 fvec3_t clr = fvec3_from_node_data(data_name_to_data_idx_lst["COLORS"], node_idx_to_node_data[i]);
1626 out_stream << " " << static_cast<int>(255*clr[0]) << " " << static_cast<int>(255*clr[1]) << " " << static_cast<int>(255*clr[2]);
1627 }
1628 if (have_normal_data) {
1629 fvec3_t nml = fvec3_from_node_data(data_name_to_data_idx_lst["NORMALS"], node_idx_to_node_data[i]);
1630 vec3_unitize(nml);
1631 out_stream << " " << std::setprecision(10) << nml[0] << " " << nml[1] << " " << nml[2];
1632 }
1633 out_stream << std::endl;
1634 }
1635 // Dump Cells
1636 for(auto& poly: cell_lst) {
1637 out_stream << poly.size() << " ";
1638 for(auto& vert: poly) {
1639 out_stream << vert << " ";
1640 }
1641 out_stream << std::endl;
1642 }
1643 /* Final newline */
1644 out_stream << std::endl;
1645 out_stream.close();
1646 return 0;
1647 }
1648 //@}
1649
1650 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1651 /** @name Point Tests & Predicates */
1652 //@{
1653 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1654 /** Return true if any component of the point has a NaN.
1655 @param test_pnt The point to test */
1656 inline bool pnt_has_nan(const fvec3_t& test_pnt) {
1657 return (std::isnan(test_pnt[0]) || std::isnan(test_pnt[1]) || std::isnan(test_pnt[2]));
1658 }
1659 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1660 /** Return 0 if point is close to the level, -1 if it is below the level, and 1 if it is above the level.
1661 @param test_pnt The point to test
1662 @param axis_index Which axis to compare to the level
1663 @param level Level to test aginst
1664 @param close_epsilon Epsilon used to check for "closeness". */
1665 inline int pnt_vs_level(const fvec3_t& test_pnt, int axis_index, uft_t level, uft_t close_epsilon=epsilon) {
1666 if (std::abs(test_pnt[axis_index]-level) < close_epsilon)
1667 return 0;
1668 else if(test_pnt[axis_index] < level)
1669 return -1;
1670 else
1671 return 1;
1672 }
1673 //@}
1674
1675 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1676 /** @name Complex Computation. */
1677 //@{
1678 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1679 /** Function that takes and returns a node_data_t (Node Data Transform) */
1680 typedef std::function<node_data_t(const node_data_t&)> p2data_func_t; // pd2pd
1681 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1682 /** Function that takes node data and returns a float (Node Data SDF) */
1683 typedef std::function<uft_t(const node_data_t&)> p2real_func_t; // pd2r
1684 /** Function that takes node data and returns a bool (Node Data Predicate) */
1685 typedef std::function<bool(const node_data_t&)> p2bool_func_t; // pd2b
1686 /** Function that takes a cell and returns a bool (Cell Predicate) */
1687 typedef std::function<bool(const cell_verts_t&)> c2bool_func_t; // c2b
1688 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1689 /** Delete cells matching a predicate function.
1690
1691 @param func Predicate function -- if true we get rid of the cell */
1692 int cull_cells(c2bool_func_t func) {
1693 int idx_last_good = -1;
1694 int start_size = num_cells();
1695 for(int i=0; i<start_size; i++)
1696 if ( !(func(cell_lst[i]))) {
1697 idx_last_good++;
1698 if (idx_last_good != i)
1699 cell_lst[idx_last_good] = cell_lst[i];
1700 }
1701 idx_last_good++;
1702 cell_lst.resize(idx_last_good);
1703 return (start_size-num_cells());
1704 }
1705 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1706 /** Add new cells with points from existing cells with the given coordinates negated.
1707
1708 When zero_epsilon is positive, some existing points may be adjusted such that components on the flipped axis near zero will become zero. When
1709 chk_point_unique is true, both the original point coordinates and the new point coordinates will be mapped to the same point index -- that is to say if
1710 you add a new point to the complex with the old coordinates or the new coordinates, you will get the same index. Note that a point collision may
1711 occur if the adjusted point is equal-ish to another, existing point -- a message is printed when this occurs.
1712
1713 @param flip_list A vector of booleans (0 or 1). 1 indicates the coordinate in a node_data_t vector should be negated.
1714 @param zero_epsilon If non-negative, will collapse flipped coordinates near zero to be precisely zero.
1715 @param reverse_vertex_order Reverse the order of vertexes in each cell -- useful if vertex order impacts rendering. */
1716 void mirror(std::vector<int> flip_list, uft_t zero_epsilon=epsilon*1000, bool reverse_vertex_order=true) {
1717 int num_start_cells = num_cells();
1718 for(int cell_idx=0; cell_idx<num_start_cells; ++cell_idx) {
1719 cell_verts_t new_cell;
1720 for(auto pidx: cell_lst[cell_idx]) {
1721 node_data_t od = node_idx_to_node_data[pidx];
1722 for(int flip_list_idx=0; flip_list_idx<static_cast<int>(flip_list.size()); ++flip_list_idx)
1723 if ((flip_list[flip_list_idx]) && (std::abs(od[flip_list_idx]) < zero_epsilon))
1724 od[flip_list_idx] = 0;
1725 node_idx_to_node_data[pidx] = od;
1726 fvec3_t new_old_pnt = fvec3_from_node_data(node_data_to_pnt, od);
1727 if constexpr (chk_point_unique) {
1728 if (pnt_to_node_idx_map.contains(new_old_pnt))
1729 if (pnt_to_node_idx_map[new_old_pnt] != pidx)
1730 std::cout << "ERROR(mirror): Collapse caused collision!" << std::endl;
1731 pnt_to_node_idx_map[new_old_pnt] = pidx;
1732 }
1733 node_data_t nd = od;
1734 for(int flip_list_idx=0; flip_list_idx<static_cast<int>(flip_list.size()); ++flip_list_idx)
1735 if (flip_list[flip_list_idx])
1736 nd[flip_list_idx] = -nd[flip_list_idx];
1737 node_idx_t p = add_node(nd);
1738 new_cell.push_back(p);
1739 }
1740 if (reverse_vertex_order)
1741 std::reverse(new_cell.begin(), new_cell.end());
1742 add_cell(req_pt_cnt_to_cell_kind(new_cell.size()), new_cell);
1743 }
1744 }
1745 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1746 /** Clear the solution cache stored in edge_solver_sdf(). */
1747 void clear_cache_edge_solver_sdf() {
1748 edge_solver_sdf(nullptr, 0, 0, nullptr, 0.0);
1749 }
1750 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1751 /** Solve an SDF function for zero on a line between two node data sets, and add the solution to the central node data.
1752
1753 The SDF must have different signs on the given points. If this is not the case, then the endpoint with SDF closest to zero will be returned.
1754
1755 Uses bisection. Halting conditions:
1756 - The SDF value at the most recent guess is within solve_epsilon of zero
1757 - The absolute difference in SDF values at the end points is less than solve_epsilon
1758
1759 This function stores a cache of previous solution results. This cache may be cleared by clear_cache_edge_solver_sdf().
1760
1761 @param pnt_idx1 First edge vertex
1762 @param pnt_idx2 Second edge vertex
1763 @param sdf_func Data function (if nullptr, then clear solver cache and return immediately)
1764 @param solve_epsilon Used to detect SDF value near zero
1765 @param dat_func Produce the node data vector for the newly solved point. */
1766 node_idx_t edge_solver_sdf(p2data_func_t dat_func, node_idx_t pnt_idx1, node_idx_t pnt_idx2, p2real_func_t sdf_func, uft_t solve_epsilon=epsilon/10) {
1767 // Solver cache.
1768 static std::unordered_map<node_idx_t, std::unordered_map<node_idx_t, node_idx_t>> edge_solver_cache;
1769 if (sdf_func == nullptr) {
1770 edge_solver_cache.clear();
1771 return pnt_idx1;
1772 }
1773 node_idx_t cache_key1 = std::min(pnt_idx1, pnt_idx2);
1774 node_idx_t cache_key2 = std::max(pnt_idx1, pnt_idx2);
1775 // Check to see if we solved this one before
1776 if (edge_solver_cache.contains(cache_key1))
1777 if (edge_solver_cache[cache_key1].contains(cache_key2))
1778 return edge_solver_cache[cache_key1][cache_key2];
1779 // Gotta do the work and solve...
1780 // Init neg&pos with the points that have eng&pos sdf values
1781 node_data_t neg_node_data;
1782 uft_t neg_pnt_sdfv;
1783 node_data_t pos_node_data = node_idx_to_node_data[pnt_idx1];
1784 uft_t pos_pnt_sdfv = sdf_func(pos_node_data);
1785 if (pos_pnt_sdfv > 0) {
1786 neg_node_data = node_idx_to_node_data[pnt_idx2];
1787 neg_pnt_sdfv = sdf_func(neg_node_data);
1788 } else {
1789 neg_node_data = pos_node_data;
1790 neg_pnt_sdfv = pos_pnt_sdfv;
1791 pos_node_data = node_idx_to_node_data[pnt_idx2];
1792 pos_pnt_sdfv = sdf_func(pos_node_data);
1793 }
1794 // Init sol_node_data to end point with sdf value nearest zero
1795 node_data_t sol_node_data;
1796 uft_t sol_pnt_sdfv;
1797 if (std::abs(pos_pnt_sdfv) < std::abs(neg_pnt_sdfv)) {
1798 sol_node_data = pos_node_data;
1799 sol_pnt_sdfv = pos_pnt_sdfv;
1800 } else {
1801 sol_node_data = neg_node_data;
1802 sol_pnt_sdfv = neg_pnt_sdfv;
1803 }
1804 if (neg_pnt_sdfv < 0) { // Just to make sure pos&neg are pos&neg...
1805 while ((std::abs(sol_pnt_sdfv) > solve_epsilon) && (pos_pnt_sdfv-neg_pnt_sdfv) > solve_epsilon) {
1806 for(decltype(pos_node_data.size()) i=0; i<pos_node_data.size(); i++)
1807 sol_node_data[i] = (pos_node_data[i] + neg_node_data[i])/static_cast<uft_t>(2.0);
1808 sol_pnt_sdfv = sdf_func(sol_node_data);
1809 if (sol_pnt_sdfv > 0) {
1810 pos_node_data = sol_node_data;
1811 pos_pnt_sdfv = sol_pnt_sdfv;
1812 } else {
1813 neg_node_data = sol_node_data;
1814 neg_pnt_sdfv = sol_pnt_sdfv;
1815 }
1816 }
1817 }
1818 // Add out point, update solver cache, and return index
1819 node_idx_t sol_pnt_idx = add_node(dat_func(sol_node_data));
1820 edge_solver_cache[cache_key1][cache_key2] = sol_pnt_idx;
1821 return sol_pnt_idx;
1822 }
1823 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1824 /** Fold triangles that cross over an SDF function.
1825 @param data_func Data function
1826 @param sdf_func SDF function
1827 @param solve_epsilon Used to detect SDF value near zero */
1828 void triangle_folder(p2data_func_t data_func, p2real_func_t sdf_func, uft_t solve_epsilon=epsilon/10) {
1829 clear_cache_edge_solver_sdf();
1830 int num_start_cells = num_cells();
1831 for(int cell_idx=0; cell_idx<num_start_cells; ++cell_idx) {
1832 if (static_cast<int>(cell_lst[cell_idx].size()) == cell_kind_to_req_pt_cnt(cell_kind_t::TRIANGLE)) {
1833 auto& cur_cell = cell_lst[cell_idx];
1834 // Find and count zeros, positive, and negative vertexes
1835 int zero_cnt= 0, plus_cnt= 0, negv_cnt= 0;
1836 int zero_loc=-1, plus_loc=-1, negv_loc=-1;
1837 for(int i=0; i<3; i++) {
1838 uft_t sdf_val = sdf_func(node_idx_to_node_data[cur_cell[i]]);
1839 if (std::abs(sdf_val) <= solve_epsilon) {
1840 zero_cnt++;
1841 zero_loc = i;
1842 } else {
1843 if (sdf_val < static_cast<uft_t>(0.0)) {
1844 plus_cnt++;
1845 plus_loc = i;
1846 } else {
1847 negv_cnt++;
1848 negv_loc = i;
1849 }
1850 }
1851 }
1852 const std::vector<std::array<int, 3>> pmat { {0, 1, 2}, {1, 2, 0}, {2, 0, 1}};
1853 std::array<int, 3> p;
1854 if ((zero_cnt == 0) && (plus_cnt > 0) && (negv_cnt >0)) { // three triangles
1855 if (plus_cnt == 1)
1856 p = pmat[plus_loc];
1857 else
1858 p = pmat[negv_loc];
1859 auto orgv0 = cur_cell[p[0]];
1860 auto orgv1 = cur_cell[p[1]];
1861 auto orgv2 = cur_cell[p[2]];
1862 auto newv1 = edge_solver_sdf(data_func, orgv0, orgv1, sdf_func, solve_epsilon);
1863 auto newv2 = edge_solver_sdf(data_func, orgv0, orgv2, sdf_func, solve_epsilon);
1864 if ((newv1 != orgv0) && (newv1 != orgv1) && (newv2 != orgv0) && (newv2 != orgv2)) {
1865 cur_cell[p[1]] = newv1; // Modify current triangle in place
1866 cur_cell[p[2]] = newv2; // Modify current triangle in place
1867 add_cell(cell_kind_t::TRIANGLE, {newv1, orgv1, newv2}); // Add new triangle
1868 add_cell(cell_kind_t::TRIANGLE, {orgv1, orgv2, newv2}); // Add new triangle
1869 }
1870 } else if ((zero_cnt == 1) && (plus_cnt == 1) && (negv_cnt == 1)) { // two triangles
1871 p = pmat[zero_loc];
1872 auto orgv0 = cur_cell[p[0]];
1873 auto orgv1 = cur_cell[p[1]];
1874 auto orgv2 = cur_cell[p[2]];
1875 auto newv0 = edge_solver_sdf(data_func, orgv1, orgv2, sdf_func, solve_epsilon);
1876 if ((newv0 != orgv1) && (newv0 != orgv2)) {
1877 cur_cell[2] = newv0; // Modify current triangle in place
1878 add_cell(cell_kind_t::TRIANGLE, {orgv0, newv0, orgv2}); // Add new triangle
1879 }
1880 }
1881 }
1882 }
1883 }
1884 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1885 /** Fold segments that cross over an SDF function.
1886 @param data_func Data function
1887 @param sdf_func SDF function
1888 @param solve_epsilon Used to detect SDF value near zero */
1889 void segment_folder(p2data_func_t data_func, p2real_func_t sdf_func, uft_t solve_epsilon=epsilon/10) {
1890 clear_cache_edge_solver_sdf();
1891 int num_start_cells = num_cells();
1892 for(int cell_idx=0; cell_idx<num_start_cells; ++cell_idx) {
1893 if (static_cast<int>(cell_lst[cell_idx].size()) == cell_kind_to_req_pt_cnt(cell_kind_t::SEGMENT)) {
1894 auto& cur_cell = cell_lst[cell_idx];
1895 int plus_cnt=0, negv_cnt=0;
1896 for(int i=0; i<2; i++) {
1897 uft_t sdf_val = sdf_func(node_idx_to_node_data[cur_cell[i]]);
1898 if (sdf_val < static_cast<uft_t>(0.0))
1899 plus_cnt++;
1900 else if (sdf_val > static_cast<uft_t>(0.0))
1901 negv_cnt++;
1902 }
1903 if ((plus_cnt == 1) && (negv_cnt == 1)) {
1904 auto orgv0 = cur_cell[0];
1905 auto orgv1 = cur_cell[1];
1906 auto newv = edge_solver_sdf(data_func, orgv0, orgv1, sdf_func, solve_epsilon);
1907 if ((newv != orgv0) && (newv != orgv1)) {
1908 cur_cell[1] = newv; // Modify current segment in place
1909 add_cell(cell_kind_t::SEGMENT, {newv, orgv1}); // Add new segment
1910 }
1911 }
1912 }
1913 }
1914 }
1915 //@}
1916
1917 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1918 /** @name Cell Predicates */
1919 //@{
1920 //--------------------------------------------------------------------------------------------------------------------------------------------------------
1921 /** Return true if ALL vertexes are above the level by more than epsilon */
1922 bool cell_above_level(const cell_verts_t cell_verts, int level_index, uft_t level, uft_t level_epsilon=epsilon) {
1923 return std::all_of(cell_verts.cbegin(), cell_verts.cend(), [this, level_index, level, level_epsilon](int v) { return (node_idx_to_node_data[v][level_index] > level+level_epsilon); });
1924 }
1925 /** Return true if ALL vertexes are below the level by more than epsilon */
1926 bool cell_below_level(const cell_verts_t cell_verts, int level_index, uft_t level, uft_t level_epsilon=epsilon) {
1927 return std::all_of(cell_verts.cbegin(), cell_verts.cend(), [this, level_index, level, level_epsilon](int v) { return (node_idx_to_node_data[v][level_index] < level+level_epsilon); });
1928 }
1929 /** Return true if cell is near SDF boundry */
1930 bool cell_near_sdf_boundry(const cell_verts_t cell_verts, p2real_func_t sdf_function, uft_t sdf_epsilon=epsilon) {
1931 int pos_cnt=0, neg_cnt=0;
1932 for(auto v: cell_verts) {
1933 uft_t sv = sdf_function(node_idx_to_node_data[v]);
1934 if (mjr::math::fnear_zero(sv, sdf_epsilon))
1935 return true;
1936 if (sv < 0) {
1937 neg_cnt++;
1938 } else {
1939 pos_cnt++;
1940 }
1941 if ((neg_cnt > 0) && (pos_cnt > 0))
1942 return true;
1943 }
1944 return false;
1945 }
1946 //@}
1947
1948 };
1949
1950 typedef MR_cell_cplx<false, false, false, false, false, double, 1.0e-3> MRccF3;
1951 typedef MR_cell_cplx<true, true, true, true, true, double, 1.0e-3> MRccT3;
1952
1953 typedef MR_cell_cplx<false, false, false, false, false, double, 1.0e-5> MRccF5;
1954 typedef MR_cell_cplx<true, true, true, true, true, double, 1.0e-5> MRccT5;
1955
1956 typedef MR_cell_cplx<false, false, false, false, false, double, 1.0e-9> MRccF9;
1957 typedef MR_cell_cplx<true, true, true, true, true, double, 1.0e-9> MRccT9;
1958}
1959
1960#define MJR_INCLUDE_MR_cell_cplx
1961#endif
mjr::MR_cell_cplx
Template class used to store and transform cell complexes (mesh/triangluation/etc....
Definition MR_cell_cplx.hpp:211
mjr::MR_cell_cplx::fvec3_list_t
std::vector< fvec3_t > fvec3_list_t
Hold a list of MR_cell_cplx::fvec3_t objects.
Definition MR_cell_cplx.hpp:261
mjr::MR_cell_cplx::p2data_func_t
std::function< node_data_t(const node_data_t &)> p2data_func_t
Function that takes and returns a node_data_t (Node Data Transform)
Definition MR_cell_cplx.hpp:1680
mjr::MR_cell_cplx::vec3_scalar_triple_product
uft_t vec3_scalar_triple_product(const fvec3_t &pnt1, const fvec3_t &pnt2, const fvec3_t &pnt3) const
Compute the vector diffrence.
Definition MR_cell_cplx.hpp:620
mjr::MR_cell_cplx::uniq_cell_lst_t
std::set< cell_verts_t > uniq_cell_lst_t
Used to uniquify cells.
Definition MR_cell_cplx.hpp:377
mjr::MR_cell_cplx::check_cell_vertexes
cell_stat_t check_cell_vertexes(cell_kind_t cell_kind, cell_verts_t cell_verts) const
Perform cell vertex checks.
Definition MR_cell_cplx.hpp:1103
mjr::MR_cell_cplx::create_named_datasets
void create_named_datasets(named_data_name_list_t scalar_name_strings, data_name_to_node_data_idx_lst_t names)
This is an overloaded member function, provided for convenience. It differs from the above function o...
Definition MR_cell_cplx.hpp:558
mjr::MR_cell_cplx::idx_of_last_point_added
node_idx_t idx_of_last_point_added() const
Retruns the index of the last point given to the add_node() method.
Definition MR_cell_cplx.hpp:410
mjr::MR_cell_cplx::vec3_linear_combination
fvec3_t vec3_linear_combination(uft_t scl1, const fvec3_t &pnt1, uft_t scl2, const fvec3_t &pnt2) const
Compute the linear combination.
Definition MR_cell_cplx.hpp:637
mjr::MR_cell_cplx::cell_verts_t
node_idx_list_t cell_verts_t
Type to hold the vertexes of a poly cell.
Definition MR_cell_cplx.hpp:309
mjr::MR_cell_cplx::check_cell_dimension
cell_stat_t check_cell_dimension(cell_kind_t cell_kind, cell_verts_t cell_verts) const
Perform cell dimension (2D cells must not be colinear & 3D cells must not be coplainer)
Definition MR_cell_cplx.hpp:1138
mjr::MR_cell_cplx::clear
void clear()
Definition MR_cell_cplx.hpp:917
mjr::MR_cell_cplx::seg_isect_to_string
std::string seg_isect_to_string(seg_isect_t seg_isect) const
Convert an seg_isect_t to a string.
Definition MR_cell_cplx.hpp:674
mjr::MR_cell_cplx::vec3_cross_product
fvec3_t vec3_cross_product(const fvec3_t &pnt1, const fvec3_t &pnt2) const
Compute the cross prodcut.
Definition MR_cell_cplx.hpp:603
mjr::MR_cell_cplx::node_to_string
std::string node_to_string(node_idx_t pnt_idx) const
Convert a fvec3_t to a string representation.
Definition MR_cell_cplx.hpp:507
mjr::MR_cell_cplx::cell_stat_is_good
bool cell_stat_is_good(cell_stat_t cell_stat) const
Convert cell_stat_t to a bool (true if GOOD).
Definition MR_cell_cplx.hpp:1072
mjr::MR_cell_cplx::cell_kind_to_structure
cell_structure_t & cell_kind_to_structure(cell_kind_t cell_kind, int dimension) const
Convert cell_kind_t value to a list of vertexes.
Definition MR_cell_cplx.hpp:944
mjr::MR_cell_cplx::geomr_pts_coplanar
bool geomr_pts_coplanar(fvec3_t p1, fvec3_t p2, fvec3_t p3, fvec3_t p4) const
Check if points are coplanar.
Definition MR_cell_cplx.hpp:901
mjr::MR_cell_cplx::c2bool_func_t
std::function< bool(const cell_verts_t &)> c2bool_func_t
Function that takes a cell and returns a bool (Cell Predicate)
Definition MR_cell_cplx.hpp:1687
mjr::MR_cell_cplx::cell_near_sdf_boundry
bool cell_near_sdf_boundry(const cell_verts_t cell_verts, p2real_func_t sdf_function, uft_t sdf_epsilon=epsilon)
Return true if cell is near SDF boundry.
Definition MR_cell_cplx.hpp:1930
mjr::MR_cell_cplx::fvec3_t
std::array< uft_t, 3 > fvec3_t
Hold a tuple of real values defining a point/vector in R^3 (also see MR_cell_cplx::vdat_t).
Definition MR_cell_cplx.hpp:258
mjr::MR_cell_cplx::status_of_last_cell_added
cell_stat_t status_of_last_cell_added() const
Retruns the status of the last cell given to the add_cell() method.
Definition MR_cell_cplx.hpp:1308
mjr::MR_cell_cplx::cell_stat_to_string
std::string cell_stat_to_string(cell_stat_t cell_stat) const
Convert cell_stat_t enum value to a string.
Definition MR_cell_cplx.hpp:1082
mjr::MR_cell_cplx::set_data_name_to_data_idx_lst
void set_data_name_to_data_idx_lst(data_name_to_node_data_idx_lst_t names)
Create named data sets by setting raw member.
Definition MR_cell_cplx.hpp:542
mjr::MR_cell_cplx::last_cell_added_was_new
bool last_cell_added_was_new() const
Retruns true if the last cell given to the add_cell() was new and the return value of status_of_last_...
Definition MR_cell_cplx.hpp:1314
mjr::MR_cell_cplx::geomr_pnt_pln_distance
uft_t geomr_pnt_pln_distance(const fvec3_t &plnpnt1, const fvec3_t &plnpnt2, const fvec3_t &plnpnt3, const fvec3_t &pnt) const
Compute the distance between a point and a plane.
Definition MR_cell_cplx.hpp:830
mjr::MR_cell_cplx::node_idx_list_t
std::vector< node_idx_t > node_idx_list_t
list of node indexes.
Definition MR_cell_cplx.hpp:255
mjr::MR_cell_cplx::node_data_to_pnt
node_data_idx_lst_t node_data_to_pnt
How to construct a point from a node dataset.
Definition MR_cell_cplx.hpp:249
mjr::MR_cell_cplx::io_result
int io_result
Type returned by I/O functions.
Definition MR_cell_cplx.hpp:1324
mjr::MR_cell_cplx::cell_below_level
bool cell_below_level(const cell_verts_t cell_verts, int level_index, uft_t level, uft_t level_epsilon=epsilon)
Return true if ALL vertexes are below the level by more than epsilon.
Definition MR_cell_cplx.hpp:1926
mjr::MR_cell_cplx::node_data_t
std::vector< uft_t > node_data_t
Type to hold a data for a single node.
Definition MR_cell_cplx.hpp:237
mjr::MR_cell_cplx::dump_cplx
void dump_cplx(int max_num_print) const
Dump object to STDOUT.
Definition MR_cell_cplx.hpp:1524
mjr::MR_cell_cplx::cell_lst_t
std::vector< cell_verts_t > cell_lst_t
Type for a list of poly cells.
Definition MR_cell_cplx.hpp:371
mjr::MR_cell_cplx::vdat_t
node_data_t vdat_t
Type to hold a list values extracted from node data via a MR_cell_cplx::node_data_idx_lst_t (also see...
Definition MR_cell_cplx.hpp:246
mjr::MR_cell_cplx::num_cells
int num_cells() const
Number of cells.
Definition MR_cell_cplx.hpp:939
mjr::MR_cell_cplx::geomi_pts_colinear
bool geomi_pts_colinear(node_idx_t pi1, node_idx_t pi2, node_idx_t pi3, node_idx_t pi4) const
Check if points are colinear.
Definition MR_cell_cplx.hpp:868
mjr::MR_cell_cplx::geomr_pts_colinear
bool geomr_pts_colinear(fvec3_t p1, fvec3_t p2, fvec3_t p3) const
Check if points are colinear.
Definition MR_cell_cplx.hpp:879
mjr::MR_cell_cplx::scalar_from_node_data
uft_t scalar_from_node_data(node_data_idx_t data_idx, node_data_t node_data) const
Extract a uft_t from a node_data_t and a node_data_idx_t.
Definition MR_cell_cplx.hpp:498
mjr::MR_cell_cplx::last_cell_stat
cell_stat_t last_cell_stat
Status of the last cell added via add_cell() method.
Definition MR_cell_cplx.hpp:389
mjr::MR_cell_cplx::vec3_dot_product
uft_t vec3_dot_product(const fvec3_t &pnt1, const fvec3_t &pnt2) const
Compute the cross prodcut.
Definition MR_cell_cplx.hpp:595
mjr::MR_cell_cplx::seg_isect_t
seg_isect_t
Segment intersection types.
Definition MR_cell_cplx.hpp:664
mjr::MR_cell_cplx::seg_isect_t::CI_VERTEX0
@ CI_VERTEX0
Cardinality=infinite. Intersection is a segment – no vertexes included.
mjr::MR_cell_cplx::seg_isect_t::BAD_SEGMENT
@ BAD_SEGMENT
At least one of the segments was degenerate.
mjr::MR_cell_cplx::seg_isect_t::C1_INTERIOR
@ C1_INTERIOR
Cardinality=1. Intersection is a single point in interior of at least one segment.
mjr::MR_cell_cplx::seg_isect_t::C0_EMPTY
@ C0_EMPTY
Cardinality=0. Intersection is the empty set.
mjr::MR_cell_cplx::seg_isect_t::CI_VERTEX1
@ CI_VERTEX1
Cardinality=infinite. Intersection is a segment – exactially one vertex in the intersection.
mjr::MR_cell_cplx::seg_isect_t::CI_VERTEX2
@ CI_VERTEX2
Cardinality=infinite. Intersection is a segment – equal to input segments.
mjr::MR_cell_cplx::seg_isect_t::C1_VERTEX1
@ C1_VERTEX1
Cardinality=1. Intersection is a single, shared vertex.
mjr::MR_cell_cplx::geomi_pts_colinear
bool geomi_pts_colinear(node_idx_t pi1, node_idx_t pi2, node_idx_t pi3) const
Check if points are colinear.
Definition MR_cell_cplx.hpp:874
mjr::MR_cell_cplx::geomi_pts_coplanar
bool geomi_pts_coplanar(const node_idx_list_t &pnt_list) const
Check if points are coplanar.
Definition MR_cell_cplx.hpp:884
mjr::MR_cell_cplx::named_vector_datasets_count
int named_vector_datasets_count() const
Return the number of named vector datasets.
Definition MR_cell_cplx.hpp:531
mjr::MR_cell_cplx::named_datasets_count
int named_datasets_count() const
Return the number of named (both vector and scalar) datasets.
Definition MR_cell_cplx.hpp:536
mjr::MR_cell_cplx::cell_stat_is_bad
bool cell_stat_is_bad(cell_stat_t cell_stat) const
Convert cell_stat_t to a bool (true if BAD).
Definition MR_cell_cplx.hpp:1077
mjr::MR_cell_cplx::geomr_tri_normal
fvec3_t geomr_tri_normal(const fvec3_t &tripnt1, const fvec3_t &tripnt2, const fvec3_t &tripnt3, bool unit) const
Compute the normal of a triangle (or a plane defined by 3 points)
Definition MR_cell_cplx.hpp:820
mjr::MR_cell_cplx::vec3_det3
uft_t vec3_det3(const fvec3_t &pnt1, const fvec3_t &pnt2, const fvec3_t &pnt3) const
Determinant of 3x3 matrix with given vectors as columns/rows.
Definition MR_cell_cplx.hpp:645
mjr::MR_cell_cplx::geomr_pnt_pnt_distance
uft_t geomr_pnt_pnt_distance(const fvec3_t &pnt1, const fvec3_t &pnt2) const
Compute the euculedian (2 norm) distance between two points.
Definition MR_cell_cplx.hpp:815
mjr::MR_cell_cplx::geomi_seg_isect_type
seg_isect_t geomi_seg_isect_type(node_idx_t ilin1pnt1, node_idx_t ilin1pnt2, node_idx_t ilin2pnt1, node_idx_t ilin2pnt2) const
Determine the nature of the intersection between two line segments.
Definition MR_cell_cplx.hpp:688
mjr::MR_cell_cplx::last_point_added_was_new
node_idx_t last_point_added_was_new() const
Retruns true if the last point given to the add_node() method was a new point.
Definition MR_cell_cplx.hpp:415
mjr::MR_cell_cplx::segment_folder
void segment_folder(p2data_func_t data_func, p2real_func_t sdf_func, uft_t solve_epsilon=epsilon/10)
Fold segments that cross over an SDF function.
Definition MR_cell_cplx.hpp:1889
mjr::MR_cell_cplx::req_pt_cnt_to_cell_kind
cell_kind_t req_pt_cnt_to_cell_kind(std::vector< int >::size_type node_count) const
Convert number of points in a cell to cell_kind_t.
Definition MR_cell_cplx.hpp:1059
mjr::MR_cell_cplx::geomi_seg_isect1
bool geomi_seg_isect1(node_idx_t ilin1pnt1, node_idx_t ilin1pnt2, node_idx_t ilin2pnt1, node_idx_t ilin2pnt2) const
Check if two line segments intersect in a single point.
Definition MR_cell_cplx.hpp:752
mjr::MR_cell_cplx::last_point_idx
node_idx_t last_point_idx
The index of the last point added via the add_node() method.
Definition MR_cell_cplx.hpp:325
mjr::MR_cell_cplx::named_data_name_list_t
std::vector< named_data_name_t > named_data_name_list_t
Type to hold a list of dataset names.
Definition MR_cell_cplx.hpp:274
mjr::MR_cell_cplx::node_idx_to_node_data_t
std::vector< node_data_t > node_idx_to_node_data_t
Type to hold a all node data.
Definition MR_cell_cplx.hpp:352
mjr::MR_cell_cplx::check_cell_faces_plainer
cell_stat_t check_cell_faces_plainer(cell_kind_t cell_kind, cell_verts_t cell_verts) const
Check if the vertexes of a cell are coplainer.
Definition MR_cell_cplx.hpp:1221
mjr::MR_cell_cplx::uft_t
flt_type uft_t
Universal floating point type – used for point components, scalar data values, vector components,...
Definition MR_cell_cplx.hpp:220
mjr::MR_cell_cplx::vec3_self_dot_product
uft_t vec3_self_dot_product(const fvec3_t &pnt) const
Compute the self dot product – i.e.
Definition MR_cell_cplx.hpp:587
mjr::MR_cell_cplx::cell_structure_t
cell_lst_t cell_structure_t
Cell segment/face/etc structure.
Definition MR_cell_cplx.hpp:936
mjr::MR_cell_cplx::create_dataset_to_point_mapping
void create_dataset_to_point_mapping(node_data_idx_lst_t point_data_index)
Create mapping beteen node data and geometric pionts.
Definition MR_cell_cplx.hpp:570
mjr::MR_cell_cplx::vec3_two_norm
uft_t vec3_two_norm(const fvec3_t &pnt) const
Compute the magnitude.
Definition MR_cell_cplx.hpp:582
mjr::MR_cell_cplx::node_count
node_idx_t node_count() const
Return the number of points in the central node data list.
Definition MR_cell_cplx.hpp:400
mjr::MR_cell_cplx::get_pnt
fvec3_t get_pnt(node_idx_t pnt_idx) const
Return point given a node index.
Definition MR_cell_cplx.hpp:405
mjr::MR_cell_cplx::add_cell
bool add_cell(cell_kind_t new_cell_kind, cell_verts_t new_cell_verts)
Add a cell.
Definition MR_cell_cplx.hpp:1269
mjr::MR_cell_cplx::epsilon
static constexpr uft_t epsilon
epsilon
Definition MR_cell_cplx.hpp:223
mjr::MR_cell_cplx::geomi_pts_coplanar
bool geomi_pts_coplanar(node_idx_t pi1, node_idx_t pi2, node_idx_t pi3, node_idx_t pi4) const
Check if points are coplanar.
Definition MR_cell_cplx.hpp:896
mjr::MR_cell_cplx::add_node
node_idx_t add_node(node_data_t node_data)
Add a point to the central node data store.
Definition MR_cell_cplx.hpp:434
mjr::MR_cell_cplx::data_name_to_data_idx_lst
data_name_to_node_data_idx_lst_t data_name_to_data_idx_lst
Hold named descriptions of data (scalars & vectors)
Definition MR_cell_cplx.hpp:363
mjr::MR_cell_cplx::write_xml_vtk
io_result write_xml_vtk(std::string file_name, std::string description)
Dump to an XML VTK unstructured grid file.
Definition MR_cell_cplx.hpp:1333
mjr::MR_cell_cplx::clear_cache_edge_solver_sdf
void clear_cache_edge_solver_sdf()
Clear the solution cache stored in edge_solver_sdf().
Definition MR_cell_cplx.hpp:1747
mjr::MR_cell_cplx::cell_kind_to_string
std::string cell_kind_to_string(cell_kind_t cell_kind) const
Convert cell_kind_t value to the VTK cell type integer.
Definition MR_cell_cplx.hpp:1046
mjr::MR_cell_cplx::last_cell_new
bool last_cell_new
True if the last cell given to the add_cell() method was new – i.e.
Definition MR_cell_cplx.hpp:385
mjr::MR_cell_cplx::create_named_datasets
void create_named_datasets(named_data_name_list_t scalar_name_strings)
Create named data sets.
Definition MR_cell_cplx.hpp:549
mjr::MR_cell_cplx::cell_kind_t
cell_kind_t
Cell Types.
Definition MR_cell_cplx.hpp:299
mjr::MR_cell_cplx::cell_kind_t::SEGMENT
@ SEGMENT
mjr::MR_cell_cplx::cell_kind_t::HEXAHEDRON
@ HEXAHEDRON
mjr::MR_cell_cplx::cell_kind_t::TRIANGLE
@ TRIANGLE
mjr::MR_cell_cplx::cell_kind_t::POINT
@ POINT
mjr::MR_cell_cplx::cell_kind_t::QUAD
@ QUAD
mjr::MR_cell_cplx::cell_kind_t::PYRAMID
@ PYRAMID
mjr::MR_cell_cplx::node_idx_t
int node_idx_t
Integer type used to indentify/index nodes.
Definition MR_cell_cplx.hpp:252
mjr::MR_cell_cplx::pnt_to_node_idx_map_t
std::map< fvec3_t, node_idx_t, pnt_less > pnt_to_node_idx_map_t
Used to uniquify points and assign point index values.
Definition MR_cell_cplx.hpp:346
mjr::MR_cell_cplx::geomr_seg_isect1
bool geomr_seg_isect1(const fvec3_t &lin1pnt1, const fvec3_t &lin1pnt2, const fvec3_t &lin2pnt1, const fvec3_t &lin2pnt2) const
Check if two line segments intersect in a single point.
Definition MR_cell_cplx.hpp:757
mjr::MR_cell_cplx::geomr_seg_tri_intersection
bool geomr_seg_tri_intersection(fvec3_t tp1, fvec3_t tp2, fvec3_t tp3, fvec3_t sp1, fvec3_t sp2) const
Compute intersection of a segment and a triangle.
Definition MR_cell_cplx.hpp:906
mjr::MR_cell_cplx::geomr_pnt_line_distance
uft_t geomr_pnt_line_distance(const fvec3_t &linpnt1, const fvec3_t &linpnt2, const fvec3_t &pnt, bool seg_distance) const
Distance between a point and a line.
Definition MR_cell_cplx.hpp:792
mjr::MR_cell_cplx::cell_stat_t
cell_stat_t
Cell Status.
Definition MR_cell_cplx.hpp:285
mjr::MR_cell_cplx::cell_stat_t::TOO_FEW_PNT
@ TOO_FEW_PNT
List of points was empty check_cell_vertexes.
mjr::MR_cell_cplx::cell_stat_t::GOOD
@ GOOD
Looks like a good cell N/A.
mjr::MR_cell_cplx::cell_stat_t::DUP_PNT
@ DUP_PNT
At least two points have the same index check_cell_vertexes.
mjr::MR_cell_cplx::cell_stat_t::BAD_FACEI
@ BAD_FACEI
Bad face-edge intersection check_cell_face_intersections.
mjr::MR_cell_cplx::cell_stat_t::NEG_PNT_IDX
@ NEG_PNT_IDX
Negative point index (i.e. not on the central node data) check_cell_vertexes.
mjr::MR_cell_cplx::cell_stat_t::TOO_MANY_PNT
@ TOO_MANY_PNT
List of points was too long (>5) check_cell_vertexes.
mjr::MR_cell_cplx::cell_stat_t::BIG_PNT_IDX
@ BIG_PNT_IDX
Point index was too big (i.e. not on the central node data) check_cell_vertexes.
mjr::MR_cell_cplx::cell_stat_t::FACE_BENT
@ FACE_BENT
A face (QUAD, HEXAHEDRON, or PYRAMID) was not plainer check_cell_faces_plainer.
mjr::MR_cell_cplx::cell_stat_t::BAD_EDGEI
@ BAD_EDGEI
Bad edge-edge intersection check_cell_edge_intersections.
mjr::MR_cell_cplx::cell_stat_t::DIM_LOW
@ DIM_LOW
Dimension low (degenerate) check_cell_dimension.
mjr::MR_cell_cplx::cell_stat_t::CONCAVE
@ CONCAVE
(QUAD, HEXAHEDRON, or PYRAMID) was concave TBD
mjr::MR_cell_cplx::last_point_new
bool last_point_new
True if the last point given to the add_node() method was new – i.e.
Definition MR_cell_cplx.hpp:329
mjr::MR_cell_cplx::cull_cells
int cull_cells(c2bool_func_t func)
Delete cells matching a predicate function.
Definition MR_cell_cplx.hpp:1692
mjr::MR_cell_cplx::named_data_name_t
std::string named_data_name_t
Type to hold a dataset name.
Definition MR_cell_cplx.hpp:271
mjr::MR_cell_cplx::cell_kind_to_dimension
int cell_kind_to_dimension(cell_kind_t cell_kind) const
Convert cell_kind_t value to the logical dimension of the cell.
Definition MR_cell_cplx.hpp:1007
mjr::MR_cell_cplx::vec3_unitize
bool vec3_unitize(fvec3_t &pnt) const
Unitize the given point in place.
Definition MR_cell_cplx.hpp:625
mjr::MR_cell_cplx::mirror
void mirror(std::vector< int > flip_list, uft_t zero_epsilon=epsilon *1000, bool reverse_vertex_order=true)
Add new cells with points from existing cells with the given coordinates negated.
Definition MR_cell_cplx.hpp:1716
mjr::MR_cell_cplx::add_cell
int add_cell(cell_kind_t new_cell_kind, cell_verts_t new_cell_verts, int dimension)
Add parts of a cell of the specified dimension.
Definition MR_cell_cplx.hpp:1245
mjr::MR_cell_cplx::edge_solver_sdf
node_idx_t edge_solver_sdf(p2data_func_t dat_func, node_idx_t pnt_idx1, node_idx_t pnt_idx2, p2real_func_t sdf_func, uft_t solve_epsilon=epsilon/10)
Solve an SDF function for zero on a line between two node data sets, and add the solution to the cent...
Definition MR_cell_cplx.hpp:1766
mjr::MR_cell_cplx::write_ply
io_result write_ply(std::string file_name, std::string description)
Dump to a PLY file.
Definition MR_cell_cplx.hpp:1571
mjr::MR_cell_cplx::pnt_has_nan
bool pnt_has_nan(const fvec3_t &test_pnt)
Return true if any component of the point has a NaN.
Definition MR_cell_cplx.hpp:1656
mjr::MR_cell_cplx::uniq_cell_lst
uniq_cell_lst_t uniq_cell_lst
Unique cell list.
Definition MR_cell_cplx.hpp:380
mjr::MR_cell_cplx::data_name_to_node_data_idx_lst_t
std::map< named_data_name_t, node_data_idx_lst_t > data_name_to_node_data_idx_lst_t
Type to map names to named data sets (index lists).
Definition MR_cell_cplx.hpp:277
mjr::MR_cell_cplx::cell_kind_to_vtk_type
int cell_kind_to_vtk_type(cell_kind_t cell_kind) const
Convert cell_kind_t value to the VTK cell type integer.
Definition MR_cell_cplx.hpp:1033
mjr::MR_cell_cplx::uft_epsilon
static constexpr uft_t uft_epsilon
Definition MR_cell_cplx.hpp:224
mjr::MR_cell_cplx::check_cell_face_intersections
cell_stat_t check_cell_face_intersections(cell_kind_t cell_kind, cell_verts_t cell_verts) const
Checks that cell faces have expected intersections.
Definition MR_cell_cplx.hpp:1198
mjr::MR_cell_cplx::node_data_idx_t
std::variant< int, uft_t > node_data_idx_t
Type used to hold an index into a node_data_t or a constant MR_cell_cplx::uft_t.
Definition MR_cell_cplx.hpp:240
mjr::MR_cell_cplx::geomi_pnt_line_distance
uft_t geomi_pnt_line_distance(node_idx_t ilinpnt1, node_idx_t ilinpnt2, node_idx_t ipnt, bool seg_distance) const
Distance between a point and a line.
Definition MR_cell_cplx.hpp:784
mjr::MR_cell_cplx::pnt_vs_level
int pnt_vs_level(const fvec3_t &test_pnt, int axis_index, uft_t level, uft_t close_epsilon=epsilon)
Return 0 if point is close to the level, -1 if it is below the level, and 1 if it is above the level.
Definition MR_cell_cplx.hpp:1665
mjr::MR_cell_cplx::named_scalar_datasets_count
int named_scalar_datasets_count() const
Return the number of named scalar datasets.
Definition MR_cell_cplx.hpp:526
mjr::MR_cell_cplx::cell_lst
cell_lst_t cell_lst
House all poly cells.
Definition MR_cell_cplx.hpp:374
mjr::MR_cell_cplx::check_cell_edge_intersections
cell_stat_t check_cell_edge_intersections(cell_kind_t cell_kind, cell_verts_t cell_verts) const
Checks that cell edges have expected intersections.
Definition MR_cell_cplx.hpp:1168
mjr::MR_cell_cplx::geomr_pnt_tri_distance
uft_t geomr_pnt_tri_distance(const fvec3_t &tripnt1, const fvec3_t &tripnt2, const fvec3_t &tripnt3, const fvec3_t &pnt) const
Compute the distance between a point and a triangle.
Definition MR_cell_cplx.hpp:836
mjr::MR_cell_cplx::vector_from_node_data
vdat_t vector_from_node_data(node_data_idx_lst_t data_idx_lst, node_data_t node_data) const
Construct a vdat_t from a node_data_t and a node_data_idx_lst_t.
Definition MR_cell_cplx.hpp:485
mjr::MR_cell_cplx::p2bool_func_t
std::function< bool(const node_data_t &)> p2bool_func_t
Function that takes node data and returns a bool (Node Data Predicate)
Definition MR_cell_cplx.hpp:1685
mjr::MR_cell_cplx::node_idx_to_node_data
node_idx_to_node_data_t node_idx_to_node_data
Hold all node data sets.
Definition MR_cell_cplx.hpp:355
mjr::MR_cell_cplx::cell_above_level
bool cell_above_level(const cell_verts_t cell_verts, int level_index, uft_t level, uft_t level_epsilon=epsilon)
Return true if ALL vertexes are above the level by more than epsilon.
Definition MR_cell_cplx.hpp:1922
mjr::MR_cell_cplx::triangle_folder
void triangle_folder(p2data_func_t data_func, p2real_func_t sdf_func, uft_t solve_epsilon=epsilon/10)
Fold triangles that cross over an SDF function.
Definition MR_cell_cplx.hpp:1828
mjr::MR_cell_cplx::geomr_pnt_zero
bool geomr_pnt_zero(const fvec3_t &p1) const
Check if points (just one point for this function) are zero.
Definition MR_cell_cplx.hpp:861
mjr::MR_cell_cplx::cell_kind_to_req_pt_cnt
int cell_kind_to_req_pt_cnt(cell_kind_t cell_kind) const
Convert cell_kind_t value to the number of points required for cell type.
Definition MR_cell_cplx.hpp:1020
mjr::MR_cell_cplx::node_data_idx_lst_t
std::vector< node_data_idx_t > node_data_idx_lst_t
Type to hold a list of data indexes – one element for scalars, 3 for vectors, etc....
Definition MR_cell_cplx.hpp:243
mjr::MR_cell_cplx::p2real_func_t
std::function< uft_t(const node_data_t &)> p2real_func_t
Function that takes node data and returns a float (Node Data SDF)
Definition MR_cell_cplx.hpp:1683
mjr::MR_cell_cplx::write_legacy_vtk
io_result write_legacy_vtk(std::string file_name, std::string description)
Dump to a legacy VTK UNSTRUCTUREDGRID file.
Definition MR_cell_cplx.hpp:1436
mjr::MR_cell_cplx::pnt_to_node_idx_map
pnt_to_node_idx_map_t pnt_to_node_idx_map
Maps points to point index – used to detect points with nearly identical geometric points in R^3.
Definition MR_cell_cplx.hpp:349
mjr::MR_cell_cplx::fvec3_from_node_data
fvec3_t fvec3_from_node_data(node_data_idx_lst_t data_idx_lst, node_data_t node_data) const
Construct a fvec3_t from a node_data_t and a node_data_idx_lst_t.
Definition MR_cell_cplx.hpp:472
mjr::MR_cell_cplx::vec3_diff
fvec3_t vec3_diff(const fvec3_t &pnt1, const fvec3_t &pnt2) const
Compute the vector diffrence.
Definition MR_cell_cplx.hpp:612
mjr::MRccF9
MR_cell_cplx< false, false, false, false, false, double, 1.0e-9 > MRccF9
Definition MR_cell_cplx.hpp:1956
mjr::MRccF5
MR_cell_cplx< false, false, false, false, false, double, 1.0e-5 > MRccF5
Definition MR_cell_cplx.hpp:1953
mjr::MRccT3
MR_cell_cplx< true, true, true, true, true, double, 1.0e-3 > MRccT3
Definition MR_cell_cplx.hpp:1951
mjr::MRccT5
MR_cell_cplx< true, true, true, true, true, double, 1.0e-5 > MRccT5
Definition MR_cell_cplx.hpp:1954
mjr::MRccT9
MR_cell_cplx< true, true, true, true, true, double, 1.0e-9 > MRccT9
Definition MR_cell_cplx.hpp:1957
mjr::MRccF3
MR_cell_cplx< false, false, false, false, false, double, 1.0e-3 > MRccF3
Definition MR_cell_cplx.hpp:1950
Generated by doxygen 1.11.0