From f2aecfac644aa9ce82efc837782a916533cbc041 Mon Sep 17 00:00:00 2001 From: dengwirda Date: Mon, 4 Mar 2019 13:11:06 -0500 Subject: [PATCH] 0.9.9 updates: "Skinny" off-centres, arcs on ellipsoid, etc --- src/geo_load.hpp | 167 +- src/hfn_load.hpp | 23 +- src/ini_load.hpp | 27 +- src/jig_load.hpp | 18 +- src/jig_read.hpp | 12 +- src/jigsaw.cpp | 6 +- src/libcpp/aabb_tree/aabb_pred_k.hpp | 757 +++++++-- src/libcpp/geom_base/intersect_k.hpp | 131 +- src/libcpp/geom_base/vect_base_k.hpp | 2 + .../geom_type/geom_mesh_ellipsoid_3.hpp | 1437 ++++++++++++++--- .../geom_type/geom_mesh_euclidean_2.hpp | 140 +- .../geom_type/geom_mesh_euclidean_3.hpp | 320 ++-- src/libcpp/iter_mesh/iter_divs_2.inc | 36 +- src/libcpp/iter_mesh/iter_dual_2.inc | 2 +- src/libcpp/iter_mesh/iter_mesh_2.hpp | 181 ++- src/libcpp/iter_mesh/iter_node_2.inc | 8 +- src/libcpp/iter_mesh/iter_params.hpp | 55 +- src/libcpp/iter_mesh/iter_timers.hpp | 107 ++ src/libcpp/iter_mesh/iter_zips_2.inc | 36 +- src/libcpp/itermesh.hpp | 5 +- src/libcpp/libbasic.hpp | 28 +- .../mesh_func/hfun_mesh_euclidean_2.hpp | 9 +- .../mesh_func/hfun_mesh_euclidean_3.hpp | 8 +- src/libcpp/rdel_mesh/rdel_base_3.hpp | 102 +- src/libcpp/rdel_mesh/rdel_cost_delfront_2.inc | 2 +- src/libcpp/rdel_mesh/rdel_cost_delfront_3.inc | 4 +- src/libcpp/rdel_mesh/rdel_mesh_2.hpp | 33 +- src/libcpp/rdel_mesh/rdel_mesh_3.hpp | 56 +- src/libcpp/rdel_mesh/rdel_offh_delfront_2.inc | 27 +- src/libcpp/rdel_mesh/rdel_offh_delfront_3.inc | 52 +- src/libcpp/rdel_mesh/rdel_params.hpp | 121 +- src/libcpp/rdel_mesh/rdel_pred_delaunay_2.hpp | 8 +- src/libcpp/rdel_mesh/rdel_pred_delaunay_3.hpp | 8 +- src/libcpp/rdel_mesh/rdel_pred_delfront_2.hpp | 76 +- src/libcpp/rdel_mesh/rdel_pred_delfront_3.hpp | 133 +- src/libcpp/rdel_mesh/rdel_refine_ball_3.inc | 155 +- src/libcpp/rdel_mesh/rdel_refine_base_2.inc | 12 +- src/libcpp/rdel_mesh/rdel_refine_base_3.inc | 16 +- src/libcpp/rdel_mesh/rdel_refine_face_3.inc | 2 +- src/libcpp/rdel_mesh/rdel_refine_topo_2.inc | 172 +- src/libcpp/rdel_mesh/rdel_refine_topo_3.inc | 384 ++++- src/libcpp/rdel_mesh/rdel_timers.hpp | 110 ++ src/libcpp/rdel_mesh/rdel_update_face_2.inc | 6 +- src/libcpp/rdel_mesh/rdel_update_face_3.inc | 8 +- src/libcpp/rdelmesh.hpp | 5 +- src/libcpp/tessellate/del_tri_euclidean_2.hpp | 28 +- src/libcpp/tessellate/del_tri_euclidean_3.hpp | 24 +- src/libcpp/tessellate/delaunay_tri_k.hpp | 2 +- src/libcpp/tessellate/delaunay_walk_mesh.inc | 89 +- 49 files changed, 3956 insertions(+), 1194 deletions(-) create mode 100644 src/libcpp/iter_mesh/iter_timers.hpp create mode 100644 src/libcpp/rdel_mesh/rdel_timers.hpp diff --git a/src/geo_load.hpp b/src/geo_load.hpp index 8b96fea..b453310 100644 --- a/src/geo_load.hpp +++ b/src/geo_load.hpp @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 19 January, 2019 + * Last updated: 28 February, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -154,7 +154,17 @@ if (this->_kind == jmsh_kind::ellipsoid_mesh) { - //!! do things here... + typename + geom_data::ellipsoid_mesh_3d + ::node_type _ndat ; + _ndat.pval(0) = _pval[0]; + _ndat.pval(1) = _pval[1]; + _ndat.pval(2) = + 0.0 ; + _ndat.itag () = _itag ; + + this->_geom-> + _ellipsoid_mesh_3d. + _mesh.push_node(_ndat, false) ; } } /*---------------------------------- parse EDGE2 data */ @@ -201,7 +211,16 @@ if (this->_kind == jmsh_kind::ellipsoid_mesh) { - //!! do things here... + typename + geom_data::ellipsoid_mesh_3d + ::edge_type _edat ; + _edat.node(0) = _node[0]; + _edat.node(1) = _node[1]; + _edat.itag () = _itag ; + + this->_geom-> + _ellipsoid_mesh_3d. + _mesh.push_edge(_edat, false) ; } } /*---------------------------------- parse TRIA3 data */ @@ -308,7 +327,10 @@ _file, geom_reader(&_geom)); } else - { + { + _jlog.push( + "**parse error: file not found!\n" ) ; + _errv = __file_not_located ; } _file.close (); @@ -531,6 +553,43 @@ _geom._ellipsoid_mesh_3d. _radC = _gmsh._radii._data[0] ; } + + for (auto _ipos = (size_t) +0 ; + _ipos != _gmsh._vert3._size ; + ++_ipos ) + { + typename + geom_data::ellipsoid_mesh_3d + ::node_type _ndat ; + _ndat.pval(0) = _gmsh. + _vert3._data[_ipos]._ppos[0]; + _ndat.pval(1) = _gmsh. + _vert3._data[_ipos]._ppos[1]; + _ndat.pval(2) = + 0.0 ; + _ndat.itag () = _gmsh. + _vert3._data[_ipos]._itag ; + + _geom._ellipsoid_mesh_3d. + _mesh.push_node(_ndat , false) ; + } + + for (auto _ipos = (size_t) +0 ; + _ipos != _gmsh._edge2._size ; + ++_ipos ) + { + typename + geom_data::ellipsoid_mesh_3d + ::edge_type _edat ; + _edat.node(0) = _gmsh. + _edge2._data[_ipos]._node[0]; + _edat.node(1) = _gmsh. + _edge2._data[_ipos]._node[1]; + _edat.itag () = _gmsh. + _edge2._data[_ipos]._itag ; + + _geom._ellipsoid_mesh_3d. + _mesh.push_edge(_edat , false) ; + } } } @@ -557,26 +616,8 @@ geom_data &_geom ) { - iptr_type _errv = __no_error ; - - std::string _path ; - std::string _name ; - std::string _fext ; - file_part ( - _jcfg._geom_file, - _path, _name, _fext ) ; - - if (_fext.find("msh") == +0 ) - { return geom_from_jmsh ( _jcfg, _jlog, _geom ) ; - } - else - { - _errv =__file_not_located ; - } - - return ( _errv ) ; } /* @@ -869,6 +910,59 @@ _errv = __invalid_argument ; } + + iptr_type _imin = + std::numeric_limits::max() ; + iptr_type _imax = + std::numeric_limits::min() ; + + iptr_type _nnPT = +0 ; + iptr_type _nnE2 = +0 ; + + for (auto _iter = _geom. + _ellipsoid_mesh_3d._mesh._set1.head() ; + _iter != _geom. + _ellipsoid_mesh_3d._mesh._set1.tend() ; + ++_iter ) + { + if (_iter->mark() < 0) continue; + + _nnPT += +1 ; + } + + for (auto _iter = _geom. + _ellipsoid_mesh_3d._mesh._set2.head() ; + _iter != _geom. + _ellipsoid_mesh_3d._mesh._set2.tend() ; + ++_iter ) + { + if (_iter->mark() < 0) continue; + + _imin = std::min( + _imin, _iter->node(0)) ; + _imin = std::min( + _imin, _iter->node(1)) ; + _imax = std::max( + _imax, _iter->node(0)) ; + _imax = std::max( + _imax, _iter->node(1)) ; + + _nnE2 += +1 ; + + if (_imin < +0 || + _imax >= _nnPT) + { + _errv = __invalid_argument ; + } + } + + if (_errv != __no_error) + { + _jlog. push ( + "**input error: GEOM. EDGE2 indexing is incorrect.\n") ; + + return _errv ; + } } return ( _errv ) ; @@ -1020,7 +1114,34 @@ _geom._ellipsoid_mesh_3d._radB) ; __dumpREAL("|3-RAD.|", - _geom._ellipsoid_mesh_3d._radC) ; + _geom._ellipsoid_mesh_3d._radC) ; + + _jlog.push("\n") ; + + iptr_type _nnPT = +0 ; + iptr_type _nnE2 = +0 ; + + for (auto _iter = _geom. + _ellipsoid_mesh_3d._mesh._set1.head() ; + _iter != _geom. + _ellipsoid_mesh_3d._mesh._set1.tend() ; + ++_iter ) + { + if (_iter->mark()>=+0) _nnPT += +1 ; + } + + __dumpINTS("|COORD.|", _nnPT) + + for (auto _iter = _geom. + _ellipsoid_mesh_3d._mesh._set2.head() ; + _iter != _geom. + _ellipsoid_mesh_3d._mesh._set2.tend() ; + ++_iter ) + { + if (_iter->mark()>=+0) _nnE2 += +1 ; + } + + __dumpINTS("|EDGE-2|", _nnE2) } _jlog.push("\n") ; diff --git a/src/hfn_load.hpp b/src/hfn_load.hpp index b719ec0..845e228 100644 --- a/src/hfn_load.hpp +++ b/src/hfn_load.hpp @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 19 January, 2019 + * Last updated: 02 February, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -372,6 +372,9 @@ } else { + _jlog.push( + "**parse error: file not found!\n" ) ; + _errv = __file_not_located ; } _file.close (); @@ -670,26 +673,8 @@ hfun_data &_hfun ) { - iptr_type _errv = __no_error ; - - std::string _path ; - std::string _name ; - std::string _fext ; - file_part ( - _jcfg._hfun_file , - _path, _name, _fext ) ; - - if (_fext.find("msh") == +0 ) - { return hfun_from_jmsh ( _jcfg, _jlog, _hfun ) ; - } - else - { - _errv =__file_not_located ; - } - - return ( _errv ) ; } /* diff --git a/src/ini_load.hpp b/src/ini_load.hpp index 541d3d2..2a45494 100644 --- a/src/ini_load.hpp +++ b/src/ini_load.hpp @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 19 January, 2019 + * Last updated: 02 February, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -287,7 +287,7 @@ jmsh_reader _jmsh ; std::ifstream _file ; _file. open( - _jcfg._init_file, std::ifstream::in); + _jcfg._init_file, std::ifstream::in) ; if (_file.is_open() ) { @@ -296,6 +296,9 @@ } else { + _jlog.push( + "**parse error: file not found!\n" ) ; + _errv = __file_not_located ; } _file.close (); @@ -307,7 +310,7 @@ ++_iter ) { _jlog.push( - "**parse error: " + * _iter + "\n") ; + "**parse error: " + * _iter + "\n" ) ; } } catch (...) @@ -547,26 +550,8 @@ mesh_data &_init ) { - iptr_type _errv = __no_error ; - - std::string _path ; - std::string _name ; - std::string _fext ; - file_part ( - _jcfg._init_file, - _path, _name, _fext ) ; - - if (_fext.find("msh") == +0 ) - { return init_from_jmsh ( _jcfg, _jlog, _init ) ; - } - else - { - _errv =__file_not_located ; - } - - return ( _errv ) ; } /* diff --git a/src/jig_load.hpp b/src/jig_load.hpp index 4644bbc..2321f5f 100644 --- a/src/jig_load.hpp +++ b/src/jig_load.hpp @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 02 January, 2019 + * Last updated: 20 February, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -231,6 +231,13 @@ this->_jjig-> _rdel_opts.iter() = _iter; } + __normal_call void_type push_mesh_rule ( + std::int32_t _rule + ) + { + this->_jjig-> + _rdel_opts.rule() = _rule; + } __normal_call void_type push_mesh_siz1 ( double _siz1 ) @@ -391,7 +398,7 @@ jcfg_reader _read; std::ifstream _file; _file. open( - _jcfg._jcfg_file, std::ifstream::in); + _jcfg._jcfg_file, std::ifstream::in) ; if (_file.is_open()) { @@ -399,7 +406,10 @@ _file, jcfg_loader(&_jcfg)); } else - { + { + _jlog.push( + "**parse error: file not found!\n" ) ; + _errv = __file_not_located ; } _file.close (); @@ -411,7 +421,7 @@ ++_iter ) { _jlog.push( - "**parse error: " + * _iter + "\n") ; + "**parse error: " + * _iter + "\n" ) ; } } catch (...) diff --git a/src/jig_read.hpp b/src/jig_read.hpp index 1733cdc..69b81b6 100644 --- a/src/jig_read.hpp +++ b/src/jig_read.hpp @@ -31,9 +31,9 @@ * -------------------------------------------------------- * - * Last updated: 31 July, 2018 + * Last updated: 20 February, 2019 * - * Copyright 2013-2018 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -114,6 +114,9 @@ __normal_call void_type push_mesh_iter ( std::int32_t /*_iter*/ ) { } + __normal_call void_type push_mesh_rule ( + std::int32_t /*_rule*/ + ) { } __normal_call void_type push_mesh_siz1 ( double /*_siz1*/ ) { } @@ -535,6 +538,11 @@ __putINTS(push_mesh_iter, _stok) ; } else + if (_stok[0] == "MESH_RULE") + { + __putINTS(push_mesh_rule, _stok) ; + } + else /*---------------------------- read OPTM keywords */ if (_stok[0] == "OPTM_ITER") { diff --git a/src/jigsaw.cpp b/src/jigsaw.cpp index a2a78d8..dd59eef 100644 --- a/src/jigsaw.cpp +++ b/src/jigsaw.cpp @@ -33,8 +33,8 @@ * JIGSAW: an unstructured mesh generation library. -------------------------------------------------------- * - * JIGSAW release 0.9.8.x - * Last updated: 22 January, 2019 + * JIGSAW release 0.9.9.x + * Last updated: 20 February, 2019 * * Copyright 2013 -- 2019 * Darren Engwirda @@ -200,7 +200,7 @@ # endif -# define __JGSWVSTR "JIGSAW VERSION 0.9.8" +# define __JGSWVSTR "JIGSAW VERSION 0.9.9" /*---------------------------------- for i/o on files */ diff --git a/src/libcpp/aabb_tree/aabb_pred_k.hpp b/src/libcpp/aabb_tree/aabb_pred_k.hpp index fadf42d..b652a8d 100644 --- a/src/libcpp/aabb_tree/aabb_pred_k.hpp +++ b/src/libcpp/aabb_tree/aabb_pred_k.hpp @@ -31,9 +31,9 @@ * ------------------------------------------------------------ * - * Last updated: 31 December, 2018 + * Last updated: 02 March, 2019 * - * Copyright 2013-2018 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -47,230 +47,721 @@ # define __AABB_PRED_K__ namespace geom_tree { + - template < typename R , - typename I , - size_t K + typename I > - class aabb_pred_node_k + class aabb_pred_node_2 { /*---------------------- node-tree intersection predicate */ public : - typedef R real_type ; - typedef I iptr_type ; + typedef R real_type ; + typedef I iptr_type ; - iptr_type static constexpr _dims = K; + iptr_type static constexpr _dims = 2; private : - real_type _ppos [ K] ; + real_type _ppos [ 2] ; public : /*---------------------------------- construct from _src. */ - __normal_call aabb_pred_node_k ( - __write_ptr (real_type) _psrc + __normal_call aabb_pred_node_2 ( + __const_ptr (real_type) _psrc ) { - /*----------------------------------- make local copy */ - for (auto _idim = _dims; _idim-- != +0; ) + this->_ppos[0] = _psrc[0]; + this->_ppos[1] = _psrc[1]; + } + +/*---------------------------------- TRUE if intersection */ + __inline_call bool_type operator() ( + __const_ptr (real_type) _bmin , + __const_ptr (real_type) _bmax + ) const + { + /*------------------------------ test if bbox overlap */ + if(!geometry::node_rect_2d ( + this->_ppos , + _bmin, _bmax)) + { + /*------------------------------ bbox can't intersect */ + return false ; + } + else { - this->_ppos[_idim] = _psrc[_idim] ; + return true ; } } -/*----------------------- TRUE for node-aabb intersection */ - __normal_call bool_type operator() ( - __const_ptr (real_type) _bmin, + } ; + + template < + typename R , + typename I + > + class aabb_pred_node_3 + { +/*---------------------- node-tree intersection predicate */ + public : + + typedef R real_type ; + typedef I iptr_type ; + + iptr_type static constexpr _dims = 3; + + private : + + real_type _ppos [ 3] ; + + public : +/*---------------------------------- construct from _src. */ + __normal_call aabb_pred_node_3 ( + __const_ptr (real_type) _psrc + ) + { + this->_ppos[0] = _psrc[0]; + this->_ppos[1] = _psrc[1]; + this->_ppos[2] = _psrc[2]; + } + +/*---------------------------------- TRUE if intersection */ + __inline_call bool_type operator() ( + __const_ptr (real_type) _bmin , __const_ptr (real_type) _bmax ) const - { - /*----------------------- test bounds along each axis */ - for (auto _idim = _dims; _idim-- != +0; ) + { + /*------------------------------ test if bbox overlap */ + if(!geometry::node_rect_3d ( + this->_ppos , + _bmin, _bmax)) { - if (_bmin[_idim] > this->_ppos[_idim] || - _bmax[_idim] < this->_ppos[_idim] ) - return ( false ) ; + /*------------------------------ bbox can't intersect */ + return false ; + } + else + { + return true ; } - - /*----------------------- intersecting if we got here */ - return ( true ) ; } } ; - template < typename R , - typename I , - size_t K + typename I > - class aabb_pred_rect_k + class aabb_pred_rect_2 { /*---------------------- rect-tree intersection predicate */ public : - typedef R real_type ; - typedef I iptr_type ; + typedef R real_type ; + typedef I iptr_type ; - iptr_type static constexpr _dims = K; + iptr_type static constexpr _dims = 2; private : - real_type _rmin [ K] ; - real_type _rmax [ K] ; + real_type _rmin [ 2] ; + real_type _rmax [ 2] ; public : /*---------------------------------- construct from _src. */ - __normal_call aabb_pred_rect_k ( - __write_ptr (real_type) _lsrc , - __write_ptr (real_type) _rsrc + __normal_call aabb_pred_rect_2 ( + __const_ptr (real_type) _asrc , + __const_ptr (real_type) _bsrc ) { - /*----------------------------------- make local copy */ - for (auto _idim = _dims; _idim-- != +0; ) + this->_rmin[0] = _asrc[0]; + this->_rmin[1] = _asrc[1]; + + this->_rmax[0] = _bsrc[0]; + this->_rmax[1] = _bsrc[1]; + } + +/*---------------------------------- TRUE if intersection */ + __inline_call bool_type operator() ( + __const_ptr (real_type) _bmin , + __const_ptr (real_type) _bmax + ) const + { + /*------------------------------ test if bbox overlap */ + if(!geometry::rect_rect_2d ( + this->_rmin , + this->_rmax , + _bmin, _bmax)) + { + /*------------------------------ bbox can't intersect */ + return false ; + } + else { - this->_rmin[_idim] = _lsrc[_idim] ; - this->_rmax[_idim] = _rsrc[_idim] ; + return true ; } } -/*----------------------- TRUE for rect-aabb intersection */ - __normal_call bool_type operator() ( - __const_ptr (real_type) _bmin, + } ; + + template < + typename R , + typename I + > + class aabb_pred_rect_3 + { +/*---------------------- rect-tree intersection predicate */ + public : + + typedef R real_type ; + typedef I iptr_type ; + + iptr_type static constexpr _dims = 3; + + private : + + real_type _rmin [ 3] ; + real_type _rmax [ 3] ; + + public : +/*---------------------------------- construct from _src. */ + __normal_call aabb_pred_rect_3 ( + __const_ptr (real_type) _asrc , + __const_ptr (real_type) _bsrc + ) + { + this->_rmin[0] = _asrc[0]; + this->_rmin[1] = _asrc[1]; + this->_rmin[2] = _asrc[2]; + + this->_rmax[0] = _bsrc[0]; + this->_rmax[1] = _bsrc[1]; + this->_rmax[2] = _bsrc[2]; + } + +/*---------------------------------- TRUE if intersection */ + __inline_call bool_type operator() ( + __const_ptr (real_type) _bmin , __const_ptr (real_type) _bmax ) const - { - /*----------------------- test bounds along each axis */ - for (auto _idim = _dims; _idim-- != +0; ) + { + /*------------------------------ test if bbox overlap */ + if(!geometry::rect_rect_3d ( + this->_rmin , + this->_rmax , + _bmin, _bmax)) { - if (_bmin[_idim] > this->_rmax[_idim] || - _bmax[_idim] < this->_rmin[_idim] ) - return ( false ) ; + /*------------------------------ bbox can't intersect */ + return false ; + } + else + { + return true ; } - - /*----------------------- intersecting if we got here */ - return ( true ) ; } } ; - template < typename R , - typename I , - size_t K + typename I > - class aabb_pred_line_k + class aabb_pred_line_2 { /*---------------------- line-tree intersection predicate */ public : - typedef R real_type ; - typedef I iptr_type ; + typedef R real_type ; + typedef I iptr_type ; - iptr_type static constexpr _dims = K; + iptr_type static constexpr _dims = 2; private : - real_type _ipos [ K] ; - real_type _jpos [ K] ; - - real_type _xden [ K] ; + real_type _ipos [ 2] ; + real_type _jpos [ 2] ; - real_type _rmin [ K] ; - real_type _rmax [ K] ; + real_type _xden [ 2] ; + + real_type _rmin [ 2] ; + real_type _rmax [ 2] ; public : /*---------------------------------- construct from _src. */ - __normal_call aabb_pred_line_k ( - __write_ptr (real_type) _isrc , - __write_ptr (real_type) _jsrc + __normal_call aabb_pred_line_2 ( + __const_ptr (real_type) _isrc , + __const_ptr (real_type) _jsrc ) { - /*------------------------------ init bbox at -+ inf. */ - for (auto _idim = _dims; _idim-- != +0; ) - { - _rmin[_idim] = - +std::numeric_limits - ::infinity() ; - - _rmax[_idim] = - -std::numeric_limits - ::infinity() ; - } + this->_ipos[0] = _isrc[0]; + this->_ipos[1] = _isrc[1]; + + this->_jpos[0] = _jsrc[0]; + this->_jpos[1] = _jsrc[1]; + + this->_rmin[0] = std::min( + _isrc[0] , _jsrc[0]) ; + this->_rmin[1] = std::min( + _isrc[1] , _jsrc[1]) ; + + this->_rmax[0] = std::max( + _isrc[0] , _jsrc[0]) ; + this->_rmax[1] = std::max( + _isrc[1] , _jsrc[1]) ; + + this->_xden[0] = _jsrc[0]- + _isrc[0]; + this->_xden[1] = _jsrc[1]- + _isrc[1]; + } - /*------------------------------ copy and calc. bbox. */ - for (auto _idim = _dims; _idim-- != +0; ) - { - this->_ipos[_idim] = _isrc[_idim] ; - this->_jpos[_idim] = _jsrc[_idim] ; +/*---------------------------------- TRUE if intersection */ + __inline_call bool_type operator() ( + __const_ptr (real_type) _bmin , + __const_ptr (real_type) _bmax + ) const + { + /*------------------------------ test if bbox overlap */ + if(!geometry::rect_rect_2d ( + this->_rmin , + this->_rmax , + _bmin, _bmax)) + /*------------------------------ bbox can't intersect */ + return false ; + + /*------------------------------ test if line overlap */ + real_type _aval, _bval ; + _aval = (_bmin[0]-_ipos[0]) + /_xden[0]; + _bval = (_bmax[0]-_ipos[0]) + /_xden[0]; + + real_type _tmin, _tmax ; + _tmin = + std::min( _aval, _bval ) ; + _tmax = + std::max( _aval, _bval ) ; + + if (_tmax<_tmin) return false ; + + _aval = (_bmin[1]-_ipos[1]) + /_xden[1]; + _bval = (_bmax[1]-_ipos[1]) + /_xden[1]; + + _tmin = std::max(_tmin, + std::min( _aval, _bval )) ; + _tmax = std::min(_tmax, + std::max( _aval, _bval )) ; + + if (_tmax<_tmin) return false ; + + return true ; + } + + } ; + + template < + typename R , + typename I + > + class aabb_pred_line_3 + { +/*---------------------- line-tree intersection predicate */ + public : + + typedef R real_type ; + typedef I iptr_type ; + + iptr_type static constexpr _dims = 3; + + private : + + real_type _ipos [ 3] ; + real_type _jpos [ 3] ; + + real_type _xden [ 3] ; + + real_type _rmin [ 3] ; + real_type _rmax [ 3] ; - this->_xden[_idim] = _jsrc[_idim] - - _isrc[_idim] ; + public : +/*---------------------------------- construct from _src. */ + __normal_call aabb_pred_line_3 ( + __const_ptr (real_type) _isrc , + __const_ptr (real_type) _jsrc + ) + { + this->_ipos[0] = _isrc[0]; + this->_ipos[1] = _isrc[1]; + this->_ipos[2] = _isrc[2]; + + this->_jpos[0] = _jsrc[0]; + this->_jpos[1] = _jsrc[1]; + this->_jpos[2] = _jsrc[2]; + + this->_rmin[0] = std::min( + _isrc[0] , _jsrc[0]) ; + this->_rmin[1] = std::min( + _isrc[1] , _jsrc[1]) ; + this->_rmin[2] = std::min( + _isrc[2] , _jsrc[2]) ; + + this->_rmax[0] = std::max( + _isrc[0] , _jsrc[0]) ; + this->_rmax[1] = std::max( + _isrc[1] , _jsrc[1]) ; + this->_rmax[2] = std::max( + _isrc[2] , _jsrc[2]) ; + + this->_xden[0] = _jsrc[0]- + _isrc[0]; + this->_xden[1] = _jsrc[1]- + _isrc[1]; + this->_xden[2] = _jsrc[2]- + _isrc[2]; + } - real_type _xmin = std::min ( - _isrc[_idim], - _jsrc[_idim]) ; +/*---------------------------------- TRUE if intersection */ + __inline_call bool_type operator() ( + __const_ptr (real_type) _bmin , + __const_ptr (real_type) _bmax + ) const + { + /*------------------------------ test if bbox overlap */ + if(!geometry::rect_rect_3d ( + this->_rmin , + this->_rmax , + _bmin, _bmax)) + /*------------------------------ bbox can't intersect */ + return false ; + + /*------------------------------ test if line overlap */ + real_type _aval, _bval ; + _aval = (_bmin[0]-_ipos[0]) + /_xden[0]; + _bval = (_bmax[0]-_ipos[0]) + /_xden[0]; - real_type _xmax = std::max ( - _isrc[_idim], - _jsrc[_idim]) ; + real_type _tmin, _tmax ; + _tmin = + std::min( _aval, _bval ) ; + _tmax = + std::max( _aval, _bval ) ; + + if (_tmax<_tmin) return false ; + + _aval = (_bmin[1]-_ipos[1]) + /_xden[1]; + _bval = (_bmax[1]-_ipos[1]) + /_xden[1]; + + _tmin = std::max(_tmin, + std::min( _aval, _bval )) ; + _tmax = std::min(_tmax, + std::max( _aval, _bval )) ; + + if (_tmax<_tmin) return false ; + + _aval = (_bmin[2]-_ipos[2]) + /_xden[2]; + _bval = (_bmax[2]-_ipos[2]) + /_xden[2]; + + _tmin = std::max(_tmin, + std::min( _aval, _bval )) ; + _tmax = std::min(_tmax, + std::max( _aval, _bval )) ; + + if (_tmax<_tmin) return false ; + + return true ; + } + + } ; + + template < + typename R , + typename I + > + class aabb_pred_flat_3 + { +/*---------------------- flat-tree intersection predicate */ + public : + + typedef R real_type ; + typedef I iptr_type ; + + iptr_type static constexpr _dims = 3; - if (this->_rmin[_idim] > _xmin) - this->_rmin[_idim] = _xmin; + private : + + real_type _ppos [ 3] ; + real_type _vnrm [ 3] ; + real_type _vabs [ 3] ; + real_type _dval ; - if (this->_rmax[_idim] < _xmax) - this->_rmax[_idim] = _xmax; - } + real_type _rmin [ 3] ; + real_type _rmax [ 3] ; + + public : +/*---------------------------------- construct from _src. */ + __normal_call aabb_pred_flat_3 ( + __const_ptr (real_type) _psrc , + __const_ptr (real_type) _vsrc , + __const_ptr (real_type) _asrc , + __const_ptr (real_type) _bsrc + ) + { + this->_dval = + geometry::dot_3d(_vsrc, _psrc) ; + + this->_ppos[0] = _psrc[0]; + this->_ppos[1] = _psrc[1]; + this->_ppos[2] = _psrc[2]; + + this->_vnrm[0] = _vsrc[0]; + this->_vnrm[1] = _vsrc[1]; + this->_vnrm[2] = _vsrc[2]; + + this->_vabs[0] = + std::abs (_vsrc[0]) ; + this->_vabs[1] = + std::abs (_vsrc[1]) ; + this->_vabs[2] = + std::abs (_vsrc[2]) ; + + this->_rmin[0] = _asrc[0]; + this->_rmin[1] = _asrc[1]; + this->_rmin[2] = _asrc[2]; + + this->_rmax[0] = _bsrc[0]; + this->_rmax[1] = _bsrc[1]; + this->_rmax[2] = _bsrc[2]; + } + +/*---------------------------------- TRUE if intersection */ + __normal_call bool_type operator() ( + __const_ptr (real_type) _bmin , + __const_ptr (real_type) _bmax + ) const + { + /*------------------------------ test if bbox overlap */ + if(!geometry::rect_rect_3d ( + this->_rmin , + this->_rmax , + _bmin, _bmax)) + /*------------------------------ bbox can't intersect */ + return false ; + + /*------------------------------ check flat intersect */ + real_type _bmid[3] = { + (real_type) +.5 * _bmin [0] + + (real_type) +.5 * _bmax [0] , + (real_type) +.5 * _bmin [1] + + (real_type) +.5 * _bmax [1] , + (real_type) +.5 * _bmin [2] + + (real_type) +.5 * _bmax [2] , + } ; + real_type _blen[3] = { + (real_type) +1. * _bmax [0] - + (real_type) +1. * _bmid [0] , + (real_type) +1. * _bmax [1] - + (real_type) +1. * _bmid [1] , + (real_type) +1. * _bmax [2] - + (real_type) +1. * _bmid [2] , + } ; + + real_type _rval = + _blen[0] * this->_vabs [0] + + _blen[1] * this->_vabs [1] + + _blen[2] * this->_vabs [2] ; + + real_type _sval = + geometry::dot_3d(_vnrm, _bmid) ; + + _sval -= this->_dval; + + return std::abs(_sval)<=_rval; + } + + } ; + + template < + typename R , + typename I + > + class aabb_pred_ball_2 + { +/*---------------------- ball-tree intersection predicate */ + public : + + typedef R real_type ; + typedef I iptr_type ; + + iptr_type static constexpr _dims = 2; + + private : + + real_type _ppos [ 2] ; + real_type _rsqr ; + + public : +/*---------------------------------- construct from _src. */ + __normal_call aabb_pred_ball_2 ( + __const_ptr (real_type) _psrc , + real_type _rsrc + ) + { + this->_ppos[0] = _psrc[0]; + this->_ppos[1] = _psrc[1]; + + this->_rsqr = _rsrc * + _rsrc ; } -/*----------------------- TRUE for line-aabb intersection */ +/*---------------------------------- TRUE if intersection */ __normal_call bool_type operator() ( - __const_ptr (real_type) _bmin, + __const_ptr (real_type) _bmin , __const_ptr (real_type) _bmax ) const { - /*--------------------------------- test bbox overlap */ - for (auto _idim = _dims; _idim-- != +0; ) + real_type _dsqr = (real_type)+0. ; + real_type _diff = (real_type)+0. ; + + /*------------------------------ check ball intersect */ + if (this->_ppos[0] < _bmin[0]) + { + _diff = _bmin[0] - + this->_ppos[0] ; + + _dsqr += _diff * _diff ; + } + else + if (this->_ppos[0] > _bmax[0]) + { + _diff = _bmax[0] - + this->_ppos[0] ; + + _dsqr += _diff * _diff ; + } + + if (this->_ppos[1] < _bmin[1]) + { + _diff = _bmin[1] - + this->_ppos[1] ; + + _dsqr += _diff * _diff ; + } + else + if (this->_ppos[1] > _bmax[1]) { - if (_bmin[_idim] > this->_rmax[_idim] || - _bmax[_idim] < this->_rmin[_idim] ) - /*----------------------------------- can't intersect */ - return ( false ) ; + _diff = _bmax[1] - + this->_ppos[1] ; + + _dsqr += _diff * _diff ; } + + return _dsqr<= this->_rsqr ; + } + + } ; + + template < + typename R , + typename I + > + class aabb_pred_ball_3 + { +/*---------------------- ball-tree intersection predicate */ + public : + + typedef R real_type ; + typedef I iptr_type ; + + iptr_type static constexpr _dims = 3; - /*--------------------------------- line-face overlap */ - real_type _tmin = - -std::numeric_limits - ::infinity() ; - - real_type _tmax = - +std::numeric_limits - ::infinity() ; + private : + + real_type _ppos [ 3] ; + real_type _rsqr ; + + public : +/*---------------------------------- construct from _src. */ + __normal_call aabb_pred_ball_3 ( + __const_ptr (real_type) _psrc , + real_type _rsrc + ) + { + this->_ppos[0] = _psrc[0]; + this->_ppos[1] = _psrc[1]; + this->_ppos[2] = _psrc[2]; + + this->_rsqr = _rsrc * + _rsrc ; + } - for (auto _idim = _dims; _idim-- != +0; ) +/*---------------------------------- TRUE if intersection */ + __normal_call bool_type operator() ( + __const_ptr (real_type) _bmin , + __const_ptr (real_type) _bmax + ) const + { + real_type _dsqr = (real_type)+0. ; + real_type _diff = (real_type)+0. ; + + /*------------------------------ check ball intersect */ + if (this->_ppos[0] < _bmin[0]) + { + _diff = _bmin[0] - + this->_ppos[0] ; + + _dsqr += _diff * _diff ; + } + else + if (this->_ppos[0] > _bmax[0]) + { + _diff = _bmax[0] - + this->_ppos[0] ; + + _dsqr += _diff * _diff ; + } + + if (this->_ppos[1] < _bmin[1]) + { + _diff = _bmin[1] - + this->_ppos[1] ; + + _dsqr += _diff * _diff ; + } + else + if (this->_ppos[1] > _bmax[1]) + { + _diff = _bmax[1] - + this->_ppos[1] ; + + _dsqr += _diff * _diff ; + } + + if (this->_ppos[2] < _bmin[2]) + { + _diff = _bmin[2] - + this->_ppos[2] ; + + _dsqr += _diff * _diff ; + } + else + if (this->_ppos[2] > _bmax[2]) { - real_type _ival=(_bmin[_idim] - - this->_ipos[_idim])/ - this->_xden[_idim] ; - - real_type _jval=(_bmax[_idim] - - this->_ipos[_idim])/ - this->_xden[_idim] ; - - _tmin = std::max(_tmin, - std::min(_ival, _jval)) ; - - _tmax = std::min(_tmax, - std::max(_ival, _jval)) ; - - if (_tmax < _tmin) return false ; + _diff = _bmax[2] - + this->_ppos[2] ; + + _dsqr += _diff * _diff ; } - return ( true ) ; + return _dsqr<= this->_rsqr ; } } ; diff --git a/src/libcpp/geom_base/intersect_k.hpp b/src/libcpp/geom_base/intersect_k.hpp index 0619630..21f4fa3 100644 --- a/src/libcpp/geom_base/intersect_k.hpp +++ b/src/libcpp/geom_base/intersect_k.hpp @@ -44,7 +44,7 @@ * -------------------------------------------------------- * - * Last updated: 07 January, 2019 + * Last updated: 02 March, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -347,6 +347,80 @@ return ( _nn ) ; // return num roots } + /* + -------------------------------------------------------- + * RECT-RECT-KD: rect-rect intersections + -------------------------------------------------------- + */ + + template < + typename data_type + > + __inline_call bool node_rect_2d ( + __const_ptr (data_type) _pp, + __const_ptr (data_type) _b0, + __const_ptr (data_type) _b1 + ) + { return _pp[0] >= _b0[0] && + _pp[0] <= _b1[0] && + _pp[1] >= _b0[1] && + _pp[1] <= _b1[1] ; + } + + template < + typename data_type + > + __inline_call bool node_rect_3d ( + __const_ptr (data_type) _pp, + __const_ptr (data_type) _b0, + __const_ptr (data_type) _b1 + ) + { return _pp[0] >= _b0[0] && + _pp[0] <= _b1[0] && + _pp[1] >= _b0[1] && + _pp[1] <= _b1[1] && + _pp[2] >= _b0[2] && + _pp[2] <= _b1[2] ; + } + + template < + typename data_type + > + __inline_call bool rect_rect_2d ( + __const_ptr (data_type) _a0, + __const_ptr (data_type) _a1, + __const_ptr (data_type) _b0, + __const_ptr (data_type) _b1 + ) + { return _a0[0] <= _b1[0] && + _b0[0] <= _a1[0] && + _a0[1] <= _b1[1] && + _b0[1] <= _a1[1] ; + } + + template < + typename data_type + > + __inline_call bool rect_rect_3d ( + __const_ptr (data_type) _a0, + __const_ptr (data_type) _a1, + __const_ptr (data_type) _b0, + __const_ptr (data_type) _b1 + ) + { return _a0[0] <= _b1[0] && + _b0[0] <= _a1[0] && + _a0[1] <= _b1[1] && + _b0[1] <= _a1[1] && + _a0[2] <= _b1[2] && + _b0[2] <= _a1[2] ; + } + + /* + -------------------------------------------------------- + * LINE-FLAT-KD: line-flat intersections + -------------------------------------------------------- + */ + template < typename data_type > @@ -466,6 +540,12 @@ } } + /* + -------------------------------------------------------- + * TRIA-FLAT-KD: tria-flat intersections + -------------------------------------------------------- + */ + template < typename data_type > @@ -500,6 +580,12 @@ return ( _ni ) ; } + /* + -------------------------------------------------------- + * PROJ-LINE-KD: node-line projections + -------------------------------------------------------- + */ + template < typename real_type > @@ -564,6 +650,49 @@ return ( true ) ; } + + template < + typename real_type + > + __normal_call bool proj_flat_3d ( + __const_ptr (real_type) _pi, + __const_ptr (real_type) _pp, // node on flat + __const_ptr (real_type) _nv, // nrm. of flat + __write_ptr (real_type) _qq + ) + { + real_type _ip[3]; + _ip[0] = _pp[0] - _pi[0] ; + _ip[1] = _pp[1] - _pi[1] ; + _ip[2] = _pp[2] - _pi[2] ; + + _qq[0] = _pi[0] ; + _qq[1] = _pi[1] ; + _qq[2] = _pi[2] ; + + real_type _ep = +100. * + std::numeric_limits::epsilon(); + + real_type _tt ; + real_type _d1 = + geometry::dot_3d(_ip, _nv) ; + real_type _d2 = + geometry::dot_3d(_nv, _nv) ; + + if (std::abs(_d2) <= _ep * std::abs(_d1) ) + return ( false ) ; + + _tt = _d1/_d2 ; + + _qq[0] = + _pi[0] + _tt * _nv[0] ; + _qq[1] = + _pi[1] + _tt * _nv[1] ; + _qq[2] = + _pi[2] + _tt * _nv[2] ; + + return ( true ) ; + } /* -------------------------------------------------------- diff --git a/src/libcpp/geom_base/vect_base_k.hpp b/src/libcpp/geom_base/vect_base_k.hpp index 42b9338..8804d2a 100644 --- a/src/libcpp/geom_base/vect_base_k.hpp +++ b/src/libcpp/geom_base/vect_base_k.hpp @@ -194,8 +194,10 @@ { _cp[0] = _v1[1] * _v2[2] - _v1[2] * _v2[1] ; + _cp[1] = _v1[2] * _v2[0] - _v1[0] * _v2[2] ; + _cp[2] = _v1[0] * _v2[1] - _v1[1] * _v2[0] ; } diff --git a/src/libcpp/geom_type/geom_mesh_ellipsoid_3.hpp b/src/libcpp/geom_type/geom_mesh_ellipsoid_3.hpp index 1ce719d..2673082 100644 --- a/src/libcpp/geom_type/geom_mesh_ellipsoid_3.hpp +++ b/src/libcpp/geom_type/geom_mesh_ellipsoid_3.hpp @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 16 January, 2019 + * Last updated: 02 March, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -50,7 +50,8 @@ template < typename R , - typename I + typename I , + typename A = allocators::basic_alloc > class geom_mesh_ellipsoid_3d : public geom_base_3d @@ -61,10 +62,12 @@ typedef R real_type ; typedef I iptr_type ; + typedef A allocator ; typedef geom_mesh_ellipsoid_3d < real_type , - iptr_type > geom_type ; + iptr_type , + allocator > geom_type ; typedef geom_base_3d < real_type , @@ -78,6 +81,113 @@ base_type::disc_type disc_type ; typedef typename base_type::ball_type ball_type ; + + + class node_type: public tria_complex_node_3 + { + /*------------------------------------ loc. node type */ + public : + iptr_type _itag ; + + char_type _fdim ; + char_type _feat ; + char_type _topo ; + + public : + /*------------------------------------ "write" access */ + __inline_call iptr_type& itag ( + ) + { return this->_itag ; + } + __inline_call char_type& feat ( + ) + { return this->_feat ; + } + __inline_call char_type& fdim ( + ) + { return this->_fdim ; + } + __inline_call char_type& topo ( + ) + { return this->_topo ; + } + /*------------------------------------ "const" access */ + __inline_call iptr_type const& itag ( + ) const + { return this->_itag ; + } + __inline_call char_type const& feat ( + ) const + { return this->_feat ; + } + __inline_call char_type const& fdim ( + ) const + { return this->_fdim ; + } + __inline_call char_type const& topo ( + ) const + { return this->_topo ; + } + + } ; + + class edge_type: public tria_complex_edge_2 + { + /*------------------------------------ loc. edge type */ + public : + iptr_type _itag ; + + char_type _feat ; + char_type _topo ; + + public : + /*------------------------------------ "write" access */ + __inline_call iptr_type& itag ( + ) + { return this->_itag ; + } + __inline_call char_type& topo ( + ) + { return this->_topo ; + } + __inline_call char_type& feat ( + ) + { return this->_feat ; + } + /*------------------------------------ "const" access */ + __inline_call iptr_type const& itag ( + ) const + { return this->_itag ; + } + __inline_call char_type const& topo ( + ) const + { return this->_topo ; + } + __inline_call char_type const& feat ( + ) const + { return this->_feat ; + } + + } ; + + typedef mesh::tria_complex_1< + node_type, + edge_type, + allocator > mesh_type ; + + typedef geom_tree::aabb_node_base_k + tree_node ; + + typedef geom_tree::aabb_item_rect_k < + real_type, + iptr_type, +3 > tree_item ; + + typedef geom_tree::aabb_tree< + tree_item, +3 , + tree_node, + allocator > tree_type ; + + iptr_type static constexpr _nbox = +4 ; public : @@ -85,12 +195,18 @@ fixed_array _bmin ; containers:: fixed_array _bmax ; + + mesh_type _mesh ; + + tree_type _ebox ; /*--------------- (x/a)**2 + (y/b)**2 + (z/c)**2 = 1. */ real_type _radA ; real_type _radB ; real_type _radC ; + + real_type _rEPS ; public : @@ -109,16 +225,92 @@ { __unreferenced(_opts) ; + /*--------------------------- form AABB for full obj. */ this->_bmin[0] = -this->_radA ; this->_bmin[1] = -this->_radB ; this->_bmin[2] = -this->_radC ; this->_bmax[0] = +this->_radA ; this->_bmax[1] = +this->_radB ; - this->_bmax[2] = +this->_radC ; + this->_bmax[2] = +this->_radC ; + + /*--------------------------- convert to R^3 coord.'s */ + for (auto _iter = + this->_mesh._set1.head() ; + _iter != + this->_mesh._set1.tend() ; + ++_iter ) + { + if (_iter->mark() < 0) + continue ; + + real_type _alon, _alat ; + _alon = _iter->pval(0) ; + _alat = _iter->pval(1) ; + + _iter->pval(0) = + this-> _radA* + std::cos(_alon) * + std::cos(_alat) ; + _iter->pval(1) = + this-> _radB* + std::sin(_alon) * + std::cos(_alat) ; + _iter->pval(2) = + this-> _radC* + std::sin(_alat) ; + } + + /*--------------------------- form rel.-tol. for prj. */ + this->_rEPS = std::sqrt ( + std::numeric_limits + ::epsilon()) ; + + real_type _rBAR ; + _rBAR = (real_type) +0. ; + _rBAR += this->_radA ; + _rBAR += this->_radB ; + _rBAR += this->_radC ; + _rBAR /= (real_type) +3. ; + + this->_rEPS *= _rBAR ; + + /*--------------------------- init. AABB for arc-seg. */ + containers:: + block_array _bbox; + + iptr_type _inum = +0; + + for (auto _iter = + this->_mesh._set2.head() ; + _iter != + this->_mesh._set2.tend() ; + ++_iter, ++_inum ) + { + if (_iter->mark() < 0) + continue ; + + iptr_type _enod[ +2] ; + _enod[0] = _iter->node(0) ; + _enod[1] = _iter->node(1) ; + + tree_item _tdat ; + _tdat.ipos() = _inum ; - //!! blah something about coastlines... + make_aabb ( + &this->_mesh. + _set1[_enod[0]].pval(0) , + &this->_mesh. + _set1[_enod[1]].pval(0) , + &_tdat .pmin(0), + &_tdat .pmax(0)) ; + + _bbox.push_tail(_tdat) ; + } + this->_ebox.load ( + _bbox.head(), + _bbox.tend(),this->_nbox) ; } /* @@ -132,12 +324,12 @@ ) { if (_fdim == +2) - return ( true ) ; + return ( true ) ; else if (_fdim == +1) - return ( false ) ; + return ! this->_ebox.empty () ; else - return ( false ) ; + return ( false ) ; } /* @@ -192,119 +384,294 @@ real_type _pi = (real_type)std::atan(1.0) * 4. ; + real_type _la = + (real_type)std::atan(0.5) * 1. ; + + real_type _lo = 2.*_pi / 10. ; + if (_mesh._tria. _nset.count() <= +8 ) { + /*--------------------------- init. reg.-icosahedron */ real_type _ppos[3] ; iptr_type _inod; _ppos[0] = this->_radA * + std::cos(_pi*(real_type)+0.0) * + std::cos(_pi*(real_type)+0.5) ; + _ppos[1] = this->_radB * std::sin(_pi*(real_type)+0.0) * - std::cos(_pi*(real_type)+0.0) ; + std::cos(_pi*(real_type)+0.5) ; + _ppos[2] = this->_radC * + std::sin(_pi*(real_type)+0.5) ; + _mesh. + _tria.push_node(_ppos, _inod) ; + _mesh. + _tria.node(_inod)->fdim() = +0; + _mesh. + _tria.node(_inod)->feat() = +0; + _mesh. + _tria.node(_inod)->topo() = +2; + + _ppos[0] = this->_radA * + std::cos(_pi*(real_type)+0.0) * + std::cos(_pi*(real_type)-0.5) ; _ppos[1] = this->_radB * std::sin(_pi*(real_type)+0.0) * - std::sin(_pi*(real_type)+0.0) ; + std::cos(_pi*(real_type)-0.5) ; _ppos[2] = this->_radC * - std::cos(_pi*(real_type)+0.0) ; + std::sin(_pi*(real_type)-0.5) ; _mesh. _tria.push_node(_ppos, _inod) ; _mesh. - _tria.node(_inod)->fdim() = +2; + _tria.node(_inod)->fdim() = +0; _mesh. _tria.node(_inod)->feat() = +0; _mesh. _tria.node(_inod)->topo() = +2; _ppos[0] = this->_radA * - std::sin(_pi*(real_type)+1.0) * - std::cos(_pi*(real_type)+0.0) ; + std::cos(_lo*(real_type)+0.0) * + std::cos(_la*(real_type)+1.0) ; _ppos[1] = this->_radB * - std::sin(_pi*(real_type)+1.0) * - std::sin(_pi*(real_type)+0.0) ; + std::sin(_lo*(real_type)+0.0) * + std::cos(_la*(real_type)+1.0) ; _ppos[2] = this->_radC * - std::cos(_pi*(real_type)+1.0) ; + std::sin(_la*(real_type)+1.0) ; _mesh. _tria.push_node(_ppos, _inod) ; _mesh. - _tria.node(_inod)->fdim() = +2; + _tria.node(_inod)->fdim() = +0; _mesh. _tria.node(_inod)->feat() = +0; _mesh. _tria.node(_inod)->topo() = +2; _ppos[0] = this->_radA * - std::sin(_pi*(real_type)+0.5) * - std::cos(_pi*(real_type)+0.0) ; + std::cos(_lo*(real_type)+1.0) * + std::cos(_la*(real_type)-1.0) ; _ppos[1] = this->_radB * - std::sin(_pi*(real_type)+0.5) * - std::sin(_pi*(real_type)+0.0) ; + std::sin(_lo*(real_type)+1.0) * + std::cos(_la*(real_type)-1.0) ; _ppos[2] = this->_radC * - std::cos(_pi*(real_type)+0.5) ; + std::sin(_la*(real_type)-1.0) ; _mesh. _tria.push_node(_ppos, _inod) ; _mesh. - _tria.node(_inod)->fdim() = +2; + _tria.node(_inod)->fdim() = +0; _mesh. _tria.node(_inod)->feat() = +0; _mesh. - _tria.node(_inod)->topo() = +2; + _tria.node(_inod)->topo() = +2; _ppos[0] = this->_radA * - std::sin(_pi*(real_type)+0.5) * - std::cos(_pi*(real_type)+0.5) ; + std::cos(_lo*(real_type)+2.0) * + std::cos(_la*(real_type)+1.0) ; _ppos[1] = this->_radB * - std::sin(_pi*(real_type)+0.5) * - std::sin(_pi*(real_type)+0.5) ; + std::sin(_lo*(real_type)+2.0) * + std::cos(_la*(real_type)+1.0) ; _ppos[2] = this->_radC * - std::cos(_pi*(real_type)+0.5) ; + std::sin(_la*(real_type)+1.0) ; + _mesh. + _tria.push_node(_ppos, _inod) ; + _mesh. + _tria.node(_inod)->fdim() = +0; + _mesh. + _tria.node(_inod)->feat() = +0; + _mesh. + _tria.node(_inod)->topo() = +2; + + _ppos[0] = this->_radA * + std::cos(_lo*(real_type)+3.0) * + std::cos(_la*(real_type)-1.0) ; + _ppos[1] = this->_radB * + std::sin(_lo*(real_type)+3.0) * + std::cos(_la*(real_type)-1.0) ; + _ppos[2] = this->_radC * + std::sin(_la*(real_type)-1.0) ; + _mesh. + _tria.push_node(_ppos, _inod) ; + _mesh. + _tria.node(_inod)->fdim() = +0; + _mesh. + _tria.node(_inod)->feat() = +0; + _mesh. + _tria.node(_inod)->topo() = +2; + + _ppos[0] = this->_radA * + std::cos(_lo*(real_type)+4.0) * + std::cos(_la*(real_type)+1.0) ; + _ppos[1] = this->_radB * + std::sin(_lo*(real_type)+4.0) * + std::cos(_la*(real_type)+1.0) ; + _ppos[2] = this->_radC * + std::sin(_la*(real_type)+1.0) ; + _mesh. + _tria.push_node(_ppos, _inod) ; + _mesh. + _tria.node(_inod)->fdim() = +0; + _mesh. + _tria.node(_inod)->feat() = +0; + _mesh. + _tria.node(_inod)->topo() = +2; + + _ppos[0] = this->_radA * + std::cos(_lo*(real_type)+5.0) * + std::cos(_la*(real_type)-1.0) ; + _ppos[1] = this->_radB * + std::sin(_lo*(real_type)+5.0) * + std::cos(_la*(real_type)-1.0) ; + _ppos[2] = this->_radC * + std::sin(_la*(real_type)-1.0) ; _mesh. _tria.push_node(_ppos, _inod) ; _mesh. - _tria.node(_inod)->fdim() = +2; + _tria.node(_inod)->fdim() = +0; _mesh. _tria.node(_inod)->feat() = +0; _mesh. - _tria.node(_inod)->topo() = +2; + _tria.node(_inod)->topo() = +2; _ppos[0] = this->_radA * - std::sin(_pi*(real_type)+0.5) * - std::cos(_pi*(real_type)+1.0) ; + std::cos(_lo*(real_type)+6.0) * + std::cos(_la*(real_type)+1.0) ; _ppos[1] = this->_radB * - std::sin(_pi*(real_type)+0.5) * - std::sin(_pi*(real_type)+1.0) ; + std::sin(_lo*(real_type)+6.0) * + std::cos(_la*(real_type)+1.0) ; _ppos[2] = this->_radC * - std::cos(_pi*(real_type)+0.5) ; + std::sin(_la*(real_type)+1.0) ; _mesh. _tria.push_node(_ppos, _inod) ; _mesh. - _tria.node(_inod)->fdim() = +2; + _tria.node(_inod)->fdim() = +0; _mesh. _tria.node(_inod)->feat() = +0; _mesh. - _tria.node(_inod)->topo() = +2; - + _tria.node(_inod)->topo() = +2; + _ppos[0] = this->_radA * - std::sin(_pi*(real_type)+0.5) * - std::cos(_pi*(real_type)+1.5) ; + std::cos(_lo*(real_type)+7.0) * + std::cos(_la*(real_type)-1.0) ; _ppos[1] = this->_radB * - std::sin(_pi*(real_type)+0.5) * - std::sin(_pi*(real_type)+1.5) ; + std::sin(_lo*(real_type)+7.0) * + std::cos(_la*(real_type)-1.0) ; _ppos[2] = this->_radC * - std::cos(_pi*(real_type)+0.5) ; + std::sin(_la*(real_type)-1.0) ; + _mesh. + _tria.push_node(_ppos, _inod) ; + _mesh. + _tria.node(_inod)->fdim() = +0; + _mesh. + _tria.node(_inod)->feat() = +0; + _mesh. + _tria.node(_inod)->topo() = +2; + + _ppos[0] = this->_radA * + std::cos(_lo*(real_type)+8.0) * + std::cos(_la*(real_type)+1.0) ; + _ppos[1] = this->_radB * + std::sin(_lo*(real_type)+8.0) * + std::cos(_la*(real_type)+1.0) ; + _ppos[2] = this->_radC * + std::sin(_la*(real_type)+1.0) ; _mesh. _tria.push_node(_ppos, _inod) ; _mesh. - _tria.node(_inod)->fdim() = +2; + _tria.node(_inod)->fdim() = +0; _mesh. _tria.node(_inod)->feat() = +0; _mesh. - _tria.node(_inod)->topo() = +2; + _tria.node(_inod)->topo() = +2; + + _ppos[0] = this->_radA * + std::cos(_lo*(real_type)+9.0) * + std::cos(_la*(real_type)-1.0) ; + _ppos[1] = this->_radB * + std::sin(_lo*(real_type)+9.0) * + std::cos(_la*(real_type)-1.0) ; + _ppos[2] = this->_radC * + std::sin(_la*(real_type)-1.0) ; + _mesh. + _tria.push_node(_ppos, _inod) ; + _mesh. + _tria.node(_inod)->fdim() = +0; + _mesh. + _tria.node(_inod)->feat() = +0; + _mesh. + _tria.node(_inod)->topo() = +2; + } } /* -------------------------------------------------------- - * LINE-SURF: helper to compute intersections. + * HELPERS: predicates for intersection tests. + -------------------------------------------------------- + */ + + __normal_call void_type make_aabb ( + real_type *_apos, + real_type *_bpos, + real_type *_rmin, + real_type *_rmax + ) + { + /*- build an AABB that encloses a spheroidal arc-seg. */ + real_type _rEPS = this->_rEPS ; + + _rmin[0] = std::min( + _apos[0], _bpos[0]) ; + _rmin[1] = std::min( + _apos[1], _bpos[1]) ; + _rmin[2] = std::min( + _apos[2], _bpos[2]) ; + + _rmax[0] = std::max( + _apos[0], _bpos[0]) ; + _rmax[1] = std::max( + _apos[1], _bpos[1]) ; + _rmax[2] = std::max( + _apos[2], _bpos[2]) ; + + real_type _rmid[3] = { + (real_type) +.5 * _rmin[0] + + (real_type) +.5 * _rmax[0] , + (real_type) +.5 * _rmin[1] + + (real_type) +.5 * _rmax[1] , + (real_type) +.5 * _rmin[2] + + (real_type) +.5 * _rmax[2] , + } ; + + real_type _rlen; + _rlen = (real_type) +0. ; + _rlen = std::max ( + _rlen , _rmax[0]-_rmin[0]); + _rlen = std::max ( + _rlen , _rmax[1]-_rmin[1]); + _rlen = std::max ( + _rlen , _rmax[2]-_rmin[2]); + + _rlen*= (real_type) +.5 ; + _rlen+= _rEPS ; + + _rmin[0] = std::min( + _rmin[0], _rmid[0]-_rlen) ; + _rmin[1] = std::min( + _rmin[1], _rmid[1]-_rlen) ; + _rmin[2] = std::min( + _rmin[2], _rmid[2]-_rlen) ; + + _rmax[0] = std::max( + _rmax[0], _rmid[0]+_rlen) ; + _rmax[1] = std::max( + _rmax[1], _rmid[1]+_rlen) ; + _rmax[2] = std::max( + _rmax[2], _rmid[2]+_rlen) ; + } + + /* + -------------------------------------------------------- + * HELPERS: predicates for intersection tests. -------------------------------------------------------- */ @@ -316,8 +683,7 @@ ) { - /*- calc. intersection of a line & spheroidal surface */ - + /*- calc. intersection of a line & spheroidal surface */ bool_type _find = false ; real_type _mm[3] = { @@ -394,6 +760,638 @@ return ( _find ) ; } + /* + -------------------------------------------------------- + * HELPERS: predicates for projection to geom. + -------------------------------------------------------- + */ + + __normal_call void_type proj_surf ( + real_type *_psrc , + real_type *_pprj + ) + { + /*--------------------------- great-ellipse projector */ + real_type _zero[3] ; + _zero[0] = (real_type) +.0 ; + _zero[1] = (real_type) +.0 ; + _zero[2] = (real_type) +.0 ; + + real_type _ttaa, _ttbb ; + if (line_surf( + _zero, _psrc, _ttaa, _ttbb) ) + { + real_type _pmid[3] = { + _psrc[0] * (real_type) +.5 + + _zero[0] * (real_type) +.5 , + _psrc[1] * (real_type) +.5 + + _zero[1] * (real_type) +.5 , + _psrc[2] * (real_type) +.5 + + _zero[2] * (real_type) +.5 + } ; + real_type _pdel[3] = { + _psrc[0] * (real_type) +.5 - + _zero[0] * (real_type) +.5 , + _psrc[1] * (real_type) +.5 - + _zero[1] * (real_type) +.5 , + _psrc[2] * (real_type) +.5 - + _zero[2] * (real_type) +.5 + } ; + + if (_ttaa > (real_type)-1.) + { + _pprj[0] = + _pmid[0] + _ttaa*_pdel[0] ; + _pprj[1] = + _pmid[1] + _ttaa*_pdel[1] ; + _pprj[2] = + _pmid[2] + _ttaa*_pdel[2] ; + } + else + if (_ttbb > (real_type)-1.) + { + _pprj[0] = + _pmid[0] + _ttbb*_pdel[0] ; + _pprj[1] = + _pmid[1] + _ttbb*_pdel[1] ; + _pprj[2] = + _pmid[2] + _ttbb*_pdel[2] ; + } + } + } + + __normal_call void_type proj_curv ( + real_type *_psrc , + real_type *_pprj , + real_type *_apos , + real_type *_bpos + ) + { + /*--------------------------- great-ellipse projector */ + real_type _zero[3] ; + _zero[0] = (real_type) +.0 ; + _zero[1] = (real_type) +.0 ; + _zero[2] = (real_type) +.0 ; + + real_type _vnrm[3] ; + geometry::tria_norm_3d ( + _zero , _apos, + _bpos , _vnrm) ; + + real_type _ptmp[3] ; + geometry::proj_flat_3d ( + _psrc , _zero, + _vnrm , _ptmp) ; + + proj_surf(_ptmp, _pprj); + + real_type _anrm[3] ; + real_type _bnrm[3] ; + geometry::tria_norm_3d ( + _zero , _apos, + _pprj , _anrm) ; + + geometry::tria_norm_3d ( + _zero , _pprj, + _bpos , _bnrm) ; + + real_type _asgn = + geometry::dot_3d(_vnrm, _anrm); + + real_type _bsgn = + geometry::dot_3d(_vnrm, _bnrm); + + if (_asgn < (real_type)+0.|| + _bsgn < (real_type)+0.) + { + + real_type _alen = + geometry:: + lensqr_3d(_apos, _pprj); + + real_type _blen = + geometry:: + lensqr_3d(_bpos, _pprj); + + if (_alen <= _blen) + { + _pprj[0] = _apos[0]; + _pprj[1] = _apos[1]; + _pprj[2] = _apos[2]; + } + else + { + _pprj[0] = _bpos[0]; + _pprj[1] = _bpos[1]; + _pprj[2] = _bpos[2]; + } + } + } + + /* + -------------------------------------------------------- + * HELPERS: predicates for intersection tests. + -------------------------------------------------------- + */ + + template < + typename hits_func + > + __normal_call bool_type ball_test ( + ball_type &_ball , + edge_type &_edge , + real_type *_apos , + real_type *_bpos , + hits_func &_hfun , + iptr_type &_hnum + ) + { + /*--------------------------- bisect curve about ball */ + bool_type _okay = false ; + + real_type _pprj[3] ; + proj_curv( _ball. _pmid , + _pprj, + _apos, _bpos) ; + + _okay = _okay | + ball_kern( _ball, _edge , + _apos, _pprj, 0 , + _hfun, _hnum) ; + + _okay = _okay | + ball_kern( _ball, _edge , + _pprj, _bpos, 0 , + _hfun, _hnum) ; + + return _okay ; + } + + template < + typename hits_func + > + __normal_call bool_type ball_kern ( + ball_type &_ball , + edge_type &_edge , + real_type *_apos , + real_type *_bpos , + iptr_type _call , + hits_func &_hfun , + iptr_type &_hnum + ) + { + /*- calc. intersection of a ball & spheroidal arc-seg */ + geometry::hits_type + _htmp = geometry::face_hits ; + + if (_call++ > +32) return false ; + + /*--------------------------- call linear intersector */ + real_type _ppos[3] ; + real_type _qpos[3] ; + iptr_type _inum = + (iptr_type)geometry::ball_line_3d ( + _ball._pmid, + _ball._rrad, + _apos,_bpos, + _ppos,_qpos ) ; + + if (_inum >= (iptr_type) +1) + { + + if (_inum == (iptr_type) +1) + { + /*--------------------------- call hit output functor */ + real_type _pprj[3] ; + proj_surf(_ppos, _pprj) ; + + real_type _plen = + geometry:: + lensqr_3d(_ppos, _pprj) ; + + if (_plen < this->_rEPS * + this->_rEPS ) + { + _hfun(&_pprj[0], _htmp , + _edge.feat() , + _edge.topo() , + _edge.itag() ) ; + + _hnum += +1 ; + + return true ; + } + } + else + if (_inum == (iptr_type) +2) + { + /*--------------------------- call hit output functor */ + real_type _pprj[3] ; + real_type _qprj[3] ; + proj_surf(_ppos, _pprj) ; + proj_surf(_qpos, _qprj) ; + + real_type _plen = + geometry:: + lensqr_3d(_ppos, _pprj) ; + + real_type _qlen = + geometry:: + lensqr_3d(_qpos, _qprj) ; + + real_type _xlen = + std::max (_plen, _qlen) ; + + if (_xlen < this->_rEPS * + this->_rEPS ) + { + _hfun(&_pprj[0], _htmp , + _edge.feat() , + _edge.topo() , + _edge.itag() ) ; + + _hfun(&_qprj[0], _htmp , + _edge.feat() , + _edge.topo() , + _edge.itag() ) ; + + _hnum += +2 ; + + return true ; + } + } + + /*--------------------------- recursive arc bisection */ + bool_type _okay = false ; + + real_type _cpos[3] = { + _apos[0] * (real_type) +.5 + + _bpos[0] * (real_type) +.5 , + _apos[1] * (real_type) +.5 + + _bpos[1] * (real_type) +.5 , + _apos[2] * (real_type) +.5 + + _bpos[2] * (real_type) +.5 , + } ; + + real_type _cprj[3] ; + proj_surf(_cpos, _cprj) ; + + _okay = _okay | + ball_kern( _ball, _edge, + _apos, _cprj, _call, + _hfun, _hnum) ; + + _okay = _okay | + ball_kern( _ball, _edge, + _cprj, _bpos, _call, + _hfun, _hnum) ; + + return _okay ; + + } + + return false ; + } + + /* + -------------------------------------------------------- + * HELPERS: predicates for intersection tests. + -------------------------------------------------------- + */ + + template < + typename hits_func + > + __normal_call bool_type flat_kern ( + flat_type &_flat , + edge_type &_edge , + real_type *_apos , + real_type *_bpos , + iptr_type _call , + hits_func &_hfun , + iptr_type &_hnum + ) + { + /*- calc. intersection of a flat & spheroidal arc-seg */ + geometry::hits_type + _htmp = geometry::face_hits ; + + if (_call++ > +32) return false ; + + /*--------------------------- call linear intersector */ + real_type _xpos[3] ; + if (geometry::line_flat_3d ( + _flat._ppos, + _flat._nvec, + _apos,_bpos, + _xpos, true) ) + { + + /*--------------------------- call hit output functor */ + real_type _xprj[3] ; + proj_surf(_xpos, _xprj) ; + + real_type _xlen = + geometry:: + lensqr_3d(_xpos, _xprj) ; + + if (_xlen < this->_rEPS * + this->_rEPS ) + { + _hfun(&_xprj[0], _htmp , + _edge.feat() , + _edge.topo() , + _edge.itag() ) ; + + _hnum += +1 ; + + return true ; + } + + /*--------------------------- recursive arc bisection */ + bool_type _okay = false ; + + real_type _cpos[3] = { + _apos[0] * (real_type) +.5 + + _bpos[0] * (real_type) +.5 , + _apos[1] * (real_type) +.5 + + _bpos[1] * (real_type) +.5 , + _apos[2] * (real_type) +.5 + + _bpos[2] * (real_type) +.5 , + } ; + + real_type _cprj[3] ; + proj_surf(_cpos, _cprj) ; + + _okay = _okay | + flat_kern( _flat, _edge, + _apos, _cprj, _call, + _hfun, _hnum) ; + + _okay = _okay | + flat_kern( _flat, _edge, + _cprj, _bpos, _call, + _hfun, _hnum) ; + + return _okay ; + + } + + return false ; + } + + /* + -------------------------------------------------------- + * HELPERS: predicates for intersection tests. + -------------------------------------------------------- + */ + + template < + typename hits_func + > + __normal_call bool_type disc_kern ( + disc_type &_disc , + real_type *_apos , + real_type *_bpos , + iptr_type _call , + hits_func &_hfun + ) + { + /*- calc. intersection of a disc & spheroidal surface */ + if (is_inside(_apos) != + is_inside(_bpos) ) + { + if (_call++ > +3 ) + { + /*--------------------------- call linear intersector */ + line_type _ldat; + _ldat._ipos[0] = _apos[0] ; + _ldat._ipos[1] = _apos[1] ; + _ldat._ipos[2] = _apos[2] ; + _ldat._jpos[0] = _bpos[0] ; + _ldat._jpos[1] = _bpos[1] ; + _ldat._jpos[2] = _bpos[2] ; + + return intersect(_ldat, _hfun); + + } + else + { + /*--------------------------- recursive arc bisection */ + bool_type _okay = false ; + + real_type _cpos[4] = { + (real_type) +.5 * _apos[0] + + (real_type) +.5 * _bpos[0] , + (real_type) +.5 * _apos[1] + + (real_type) +.5 * _bpos[1] , + (real_type) +.5 * _apos[2] + + (real_type) +.5 * _bpos[2] + } ; + + real_type _cdir[4] = { + _cpos[0] - _disc._pmid[0] , + _cpos[1] - _disc._pmid[1] , + _cpos[2] - _disc._pmid[2] , + } ; + + _cdir[3] = + geometry::length_3d(_cdir); + _cdir[0]/= _cdir[3] ; + _cdir[1]/= _cdir[3] ; + _cdir[2]/= _cdir[3] ; + + _cpos[0] = _disc._pmid[0] + + _cdir[0] * _disc._rrad ; + _cpos[1] = _disc._pmid[1] + + _cdir[1] * _disc._rrad ; + _cpos[2] = _disc._pmid[2] + + _cdir[2] * _disc._rrad ; + + _okay = _okay | + disc_kern( _disc, + _apos , _cpos, _call, _hfun) ; + + _okay = _okay | + disc_kern( _disc, + _cpos , _bpos, _call, _hfun) ; + + return _okay ; + + } + } + + return false ; + } + + /* + -------------------------------------------------------- + * HELPERS: predicates for intersection tests. + -------------------------------------------------------- + */ + + template < + typename hits_func + > + class flat_intersect + { + /*------------------ flat-geom intersection predicate */ + public : + geom_type &_geom ; + flat_type &_flat ; + + hits_func &_hfun ; + + bool_type _find ; + iptr_type _hnum ; + + public : + flat_intersect operator = ( + flat_intersect & + ) = delete ; + flat_intersect operator = ( + flat_intersect&& + ) = delete ; + + public : + /*----------------------- construct using _src. obj. */ + __normal_call flat_intersect ( + flat_type &_fsrc , + geom_type &_gsrc , + hits_func &_hsrc + ) : _geom( _gsrc ) , + _flat( _fsrc ) , + _hfun( _hsrc ) + { + this->_hnum = + 0 ; + this->_find = false ; + } + /*----------------------- all intersection about node */ + __normal_call void_type operator() ( + typename + tree_type::item_data *_iptr + ) + { + for ( ; _iptr != nullptr; + _iptr = _iptr->_next ) + { + /*--------------- flat/edge intersection test */ + iptr_type _epos = + _iptr->_data.ipos() ; + + iptr_type _enod[2]; + _enod[0] =_geom. + _mesh._set2[_epos].node(0) ; + _enod[1] =_geom. + _mesh._set2[_epos].node(1) ; + + /*--------------- call output function on hit */ + _geom .flat_kern ( + this->_flat, + this->_geom. + _mesh ._set2[_epos] , + &this->_geom. + _mesh._set1[_enod[ 0]].pval(0) , + &this->_geom. + _mesh._set1[_enod[ 1]].pval(0) , + + 0 , + this->_hfun, + this->_hnum) ; + + this->_find = + this->_find | (this->_hnum!=0) ; + } + } + + } ; + + template < + typename hits_func + > + class ball_intersect + { + /*------------------ ball-geom intersection predicate */ + public : + geom_type &_geom ; + ball_type &_ball ; + + real_type _rmin[3] ; + real_type _rmax[3] ; + + hits_func &_hfun ; + + bool_type _find ; + iptr_type _hnum ; + + public : + ball_intersect operator = ( + ball_intersect & + ) = delete ; + ball_intersect operator = ( + ball_intersect&& + ) = delete ; + + public : + /*----------------------- construct using _src. obj. */ + __normal_call ball_intersect ( + ball_type &_bsrc , + real_type *_psrc , + real_type *_qsrc , + geom_type &_gsrc , + hits_func &_hsrc + ) : _geom( _gsrc ) , + _ball( _bsrc ) , + _hfun( _hsrc ) + { + this->_rmin[0] = _psrc[0] ; + this->_rmin[1] = _psrc[1] ; + this->_rmin[2] = _psrc[2] ; + + this->_rmax[0] = _qsrc[0] ; + this->_rmax[1] = _qsrc[1] ; + this->_rmax[2] = _qsrc[2] ; + + this->_hnum = + 0 ; + this->_find = false ; + } + /*----------------------- all intersection about node */ + __normal_call void_type operator() ( + typename + tree_type::item_data *_iptr + ) + { + for ( ; _iptr != nullptr; + _iptr = _iptr->_next ) + { + /*--------------- ball/edge intersection test */ + iptr_type _epos = + _iptr->_data.ipos() ; + + iptr_type _enod[2]; + _enod[0] =_geom. + _mesh._set2[_epos].node(0) ; + _enod[1] =_geom. + _mesh._set2[_epos].node(1) ; + + /*--------------- call output function on hit */ + _geom. ball_test ( + this->_ball, + this->_geom. + _mesh ._set2[_epos] , + &this->_geom. + _mesh._set1[_enod[ 0]].pval(0) , + &this->_geom. + _mesh._set1[_enod[ 1]].pval(0) , + this->_hfun, + this->_hnum) ; + + this->_find = + this->_find | (this->_hnum!=0) ; + } + } + + } ; + /* -------------------------------------------------------- * INTERSECT: find FLAT/1-GEOM. intersections. @@ -408,10 +1406,29 @@ hits_func &_hfun ) { - __unreferenced(_flat) ; - __unreferenced(_hfun) ; - - return ( false ) ; + /*------------------ tree-flat intersection predicate */ + typedef + geom_tree::aabb_pred_flat_3 < + real_type, + iptr_type > tree_pred ; + + /*------------------ dual-face intersection predicate */ + typedef + flat_intersect < + hits_func > hits_pred ; + + /*------------------ call actual intersection testing */ + tree_pred _pred(_flat. _ppos, + _flat. _nvec, + _flat. _rmin, + _flat. _rmax) ; + hits_pred _func(_flat, + *this, _hfun) ; + + this->_ebox.find(_pred,_func) ; + + /*------------------ _TRUE if any intersections found */ + return ( _func._find ) ; } /* @@ -428,10 +1445,39 @@ hits_func &_hfun ) { - __unreferenced(_ball) ; - __unreferenced(_hfun) ; - - return ( false ) ; + /*------------------ tree-ball intersection predicate */ + typedef + geom_tree::aabb_pred_ball_3 < + real_type, + iptr_type > tree_pred ; + + /*------------------ ball-line intersection predicate */ + typedef + ball_intersect < + hits_func > hits_pred ; + + /*------------------ call actual intersection testing */ + real_type _rmin[3] = { + _ball._pmid[0] -_ball. _rrad, + _ball._pmid[1] -_ball. _rrad, + _ball._pmid[2] -_ball. _rrad + } ; + real_type _rmax[3] = { + _ball._pmid[0] +_ball. _rrad, + _ball._pmid[1] +_ball. _rrad, + _ball._pmid[2] +_ball. _rrad + } ; + + tree_pred _pred(_ball. _pmid, + _ball. _rrad) ; + hits_pred _func(_ball, + _rmin, _rmax, + *this, _hfun) ; + + this->_ebox.find(_pred,_func) ; + + /*------------------ _TRUE if any intersections found */ + return ( _func._find ) ; } /* @@ -450,70 +1496,73 @@ { bool_type _find = false ; - real_type _ipos[3] ; - _ipos[0] =_line._ipos[0]; - _ipos[1] =_line._ipos[1]; - _ipos[2] =_line._ipos[2]; + real_type _ipos[3] ; + _ipos[0] = _line._ipos[ 0]; + _ipos[1] = _line._ipos[ 1]; + _ipos[2] = _line._ipos[ 2]; - real_type _jpos[3] ; - _jpos[0] =_line._jpos[0]; - _jpos[1] =_line._jpos[1]; - _jpos[2] =_line._jpos[2]; + real_type _jpos[3] ; + _jpos[0] = _line._jpos[ 0]; + _jpos[1] = _line._jpos[ 1]; + _jpos[2] = _line._jpos[ 2]; real_type _ttaa, _ttbb; if (line_surf(_ipos, _jpos, _ttaa, _ttbb)) { - real_type _pmid[3] = { - _jpos[0] * (real_type)+.5 + - _ipos[0] * (real_type)+.5 , - _jpos[1] * (real_type)+.5 + - _ipos[1] * (real_type)+.5 , - _jpos[2] * (real_type)+.5 + - _ipos[2] * (real_type)+.5 - } ; - real_type _pdel[3] = { - _jpos[0] * (real_type)+.5 - - _ipos[0] * (real_type)+.5 , - _jpos[1] * (real_type)+.5 - - _ipos[1] * (real_type)+.5 , - _jpos[2] * (real_type)+.5 - - _ipos[2] * (real_type)+.5 - } ; - - real_type _apos[3] = { - _pmid[0] + _ttaa * _pdel[0] , - _pmid[1] + _ttaa * _pdel[1] , - _pmid[2] + _ttaa * _pdel[2] - } ; - real_type _bpos[3] = { - _pmid[0] + _ttbb * _pdel[0] , - _pmid[1] + _ttbb * _pdel[1] , - _pmid[2] + _ttbb * _pdel[2] - } ; - - char_type _hits = - geometry::null_hits ; - char_type _feat = +2; - char_type _topo = +2; - iptr_type _itag = +0; - - if (_ttaa >= (real_type)-1. && - _ttaa <= (real_type)+1. ) - { - _find = true ; - _hfun ( _apos, _hits , - _feat, _topo , - _itag) ; - } - - if (_ttbb >= (real_type)-1. && - _ttbb <= (real_type)+1. ) - { - _find = true ; - _hfun ( _bpos, _hits , - _feat, _topo , - _itag) ; - } + + real_type _pmid[3] = { + _jpos[0] * (real_type)+.5 + + _ipos[0] * (real_type)+.5 , + _jpos[1] * (real_type)+.5 + + _ipos[1] * (real_type)+.5 , + _jpos[2] * (real_type)+.5 + + _ipos[2] * (real_type)+.5 + } ; + real_type _pdel[3] = { + _jpos[0] * (real_type)+.5 - + _ipos[0] * (real_type)+.5 , + _jpos[1] * (real_type)+.5 - + _ipos[1] * (real_type)+.5 , + _jpos[2] * (real_type)+.5 - + _ipos[2] * (real_type)+.5 + } ; + + real_type _apos[3] = { + _pmid[0] + _ttaa * _pdel[0] , + _pmid[1] + _ttaa * _pdel[1] , + _pmid[2] + _ttaa * _pdel[2] + } ; + + real_type _bpos[3] = { + _pmid[0] + _ttbb * _pdel[0] , + _pmid[1] + _ttbb * _pdel[1] , + _pmid[2] + _ttbb * _pdel[2] + } ; + + char_type _hits = + geometry::null_hits ; + char_type _feat = +2; + char_type _topo = +2; + iptr_type _itag = +0; + + if (_ttaa >= (real_type)-1.) + if (_ttaa <= (real_type)+1.) + { + _find = true ; + _hfun ( _apos, _hits , + _feat, _topo , + _itag) ; + } + + if (_ttbb >= (real_type)-1.) + if (_ttbb <= (real_type)+1.) + { + _find = true ; + _hfun ( _bpos, _hits , + _feat, _topo , + _itag) ; + } + } return ( _find ) ; @@ -536,88 +1585,88 @@ { bool_type _find = false ; - real_type _bvec[3] = { - _sbal[0] - _disc._pmid[0] , - _sbal[1] - _disc._pmid[1] , - _sbal[2] - _disc._pmid[2] - } ; - real_type _blen = - geometry::length_3d(_bvec) ; - _bvec[0]/= _blen ; - _bvec[1]/= _blen ; - _bvec[2]/= _blen ; - - _bvec[0]*= _disc._rrad; - _bvec[1]*= _disc._rrad; - _bvec[2]*= _disc._rrad; - - real_type _tpos[3] = { - _disc._pmid[0] + _bvec[0] , - _disc._pmid[1] + _bvec[1] , - _disc._pmid[2] + _bvec[2] , - } ; - + __unreferenced(_sbal) ; + real_type _circ[3] = { (real_type) +0. , (real_type) +0. , (real_type) +0. } ; - real_type _ttaa; - real_type _ttbb; - if (line_surf(_circ, _tpos, _ttaa, _ttbb)) - { - real_type _pmid[3] = { - _tpos[0] * (real_type)+.5 + - _circ[0] * (real_type)+.5 , - _tpos[1] * (real_type)+.5 + - _circ[1] * (real_type)+.5 , - _tpos[2] * (real_type)+.5 + - _circ[2] * (real_type)+.5 - } ; - real_type _pdel[3] = { - _tpos[0] * (real_type)+.5 - - _circ[0] * (real_type)+.5 , - _tpos[1] * (real_type)+.5 - - _circ[1] * (real_type)+.5 , - _tpos[2] * (real_type)+.5 - - _circ[2] * (real_type)+.5 - } ; - - real_type _apos[3] = { - _pmid[0] + _ttaa * _pdel[0] , - _pmid[1] + _ttaa * _pdel[1] , - _pmid[2] + _ttaa * _pdel[2] - } ; - real_type _bpos[3] = { - _pmid[0] + _ttbb * _pdel[0] , - _pmid[1] + _ttbb * _pdel[1] , - _pmid[2] + _ttbb * _pdel[2] - } ; - - char_type _hits = - geometry::null_hits ; - char_type _feat = +2; - char_type _topo = +2; - iptr_type _itag = +0; - - if (_ttaa > (real_type)-1.) - { - _find = true ; - _hfun ( _apos, _hits , - _feat, _topo , - _itag) ; - } + real_type _pdir[4] = { + _disc._pmid[0] - _circ[0] , + _disc._pmid[1] - _circ[1] , + _disc._pmid[2] - _circ[2] , + (real_type) +0. } ; - if (_ttbb > (real_type)-1.) - { - _find = true ; - _hfun ( _bpos, _hits , - _feat, _topo , - _itag) ; - } - } - - return ( _find ) ; + real_type _ddir[4] ; + geometry::cross_3d( + _pdir, _disc._nvec, _ddir) ; + + _pdir[3] = + geometry::length_3d(_pdir) ; + _pdir[0]/= _pdir[3] ; + _pdir[1]/= _pdir[3] ; + _pdir[2]/= _pdir[3] ; + + _ddir[3] = + geometry::length_3d(_ddir) ; + _ddir[0]/= _ddir[3] ; + _ddir[1]/= _ddir[3] ; + _ddir[2]/= _ddir[3] ; + + real_type _apos[3] = { + _disc._pmid[0] - + _disc._rrad *_pdir[0] , + _disc._pmid[1] - + _disc._rrad *_pdir[1] , + _disc._pmid[2] - + _disc._rrad *_pdir[2] + } ; + + real_type _bpos[3] = { + _disc._pmid[0] + + _disc._rrad *_ddir[0] , + _disc._pmid[1] + + _disc._rrad *_ddir[1] , + _disc._pmid[2] + + _disc._rrad *_ddir[2] + } ; + + real_type _cpos[3] = { + _disc._pmid[0] + + _disc._rrad *_pdir[0] , + _disc._pmid[1] + + _disc._rrad *_pdir[1] , + _disc._pmid[2] + + _disc._rrad *_pdir[2] + } ; + + real_type _dpos[3] = { + _disc._pmid[0] - + _disc._rrad *_ddir[0] , + _disc._pmid[1] - + _disc._rrad *_ddir[1] , + _disc._pmid[2] - + _disc._rrad *_ddir[2] + } ; + + _find = _find | + disc_kern(_disc, + _apos , _bpos, +0, _hfun) ; + + _find = _find | + disc_kern(_disc, + _bpos , _cpos, +0, _hfun) ; + + _find = _find | + disc_kern(_disc, + _cpos , _dpos, +0, _hfun) ; + + _find = _find | + disc_kern(_disc, + _dpos , _apos, +0, _hfun) ; + + return _find ; } /* @@ -642,7 +1691,7 @@ + _zz * _zz <= (real_type) +1. ) ; } - + } ; diff --git a/src/libcpp/geom_type/geom_mesh_euclidean_2.hpp b/src/libcpp/geom_type/geom_mesh_euclidean_2.hpp index 9630c0f..05154a0 100644 --- a/src/libcpp/geom_type/geom_mesh_euclidean_2.hpp +++ b/src/libcpp/geom_type/geom_mesh_euclidean_2.hpp @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 19 January, 2019 + * Last updated: 02 March, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -359,18 +359,20 @@ typename list_type , typename geom_opts > - __normal_call char_type node_feat ( + __normal_call void_type node_feat ( iptr_type *_node , list_type &_aset , + char_type &_feat , + char_type &_topo , geom_opts &_opts ) { - char_type _feat = - (_aset.count() == +2) - ? null_feat : soft_feat ; - real_type _DtoR = (real_type) +3.1415926536 / 180. ; + + real_type _ZERO = -1. + + std::numeric_limits + ::epsilon(); real_type _phi1 = (real_type)+180. - _opts.phi1(); @@ -384,9 +386,14 @@ __unreferenced(_node) ; + _feat = null_feat ; + _topo = (char_type)_aset.count(); + for (auto _ipos = _aset.head() ; _ipos != _aset.tend() ; ++_ipos ) + { + char_type _tbad = +1 ; for (auto _jpos = _ipos+1 ; _jpos != _aset.tend() ; ++_jpos ) @@ -432,6 +439,8 @@ else _acos *= (real_type)-1.; + if (_acos >= _ZERO) + { if (_acos <= _hard) { _feat = @@ -443,9 +452,22 @@ _feat = std::max(_feat, soft_feat) ; } - } - - return _feat ; + } + else + { + if (_tbad >= + 1 ) + { + _topo -= _tbad--; + } + } + } + } + { + if (_topo != + 0 ) + if (_topo != + 2 ) + _feat = + std::max(_feat, soft_feat) ; + } } /* @@ -491,8 +513,6 @@ } } - if (_opts.feat() ) - { /*---------------------------------- find sharp feat. */ for (auto _iter = this->_tria._set1.head() ; @@ -501,22 +521,33 @@ ++_iter ) { /*---------------------------------- find sharp 0-dim */ - if (_iter->mark() >= +0) + if (_iter->mark() >= +0 ) { + if (_iter->itag() <= -1 || + _opts .feat() ) + { + /*---------------------------------- set geo.-defined */ _eadj.set_count (0); - this->_tria.node_edge ( - &_iter->node(0), _eadj) ; + this->_tria.node_edge ( + &_iter->node (0), _eadj) ; - _iter->topo () = - (char_type)_eadj.count() ; - - _iter->feat () = node_feat( - &_iter->node(0), _eadj, - _opts) ; + node_feat ( + &_iter->node (0), + _eadj , + _iter->feat () , + _iter->topo () , + _opts ) ; + + if (_iter->itag() <= -1) + { + /*---------------------------------- set user-defined */ + _iter->feat () = + std::max(_iter->feat () , + soft_feat) ; + } + } } - } - } for (auto _iter = @@ -605,7 +636,7 @@ if (_imin < (iptr_type) +0 && _imax < (iptr_type) +0 ) __assert( _imin >= +0 && - "INIT-PART: -ve part index"); + "GEOM::INIT-PART: -ve part index!") ; /*----------------------------- push unique PART id's */ auto _flag = @@ -746,24 +777,6 @@ } } ; - /*----------------------------- setup user-feat. face */ - for (auto _iter = - this->_tria._set2.head() ; - _iter != - this->_tria._set2.tend() ; - ++_iter ) - { - if (_iter->mark() >= +0 && - _iter->itag() <= -1 ) - { - this->_tria._set1[ - _iter->node(0)].itag() = -1 ; - - this->_tria._set1[ - _iter->node(1)].itag() = -1 ; - } - } - /*----------------------------- init. aabb at -+ inf. */ for(auto _idim = 2; _idim-- != 0 ; ) { @@ -854,8 +867,7 @@ { if (_iter->mark() >= +0 ) { - if (_opts .feat() && - _iter->feat() != null_feat) + if (_iter->feat() != null_feat) { /*----------------------------- push any 'real' feat. */ iptr_type _node = -1 ; @@ -873,6 +885,10 @@ _mesh._tria.node (_node)->topo() = _iter->topo() ; + + _mesh._tria.node + (_node)->part() + = _iter->itag() ; } } else @@ -894,6 +910,10 @@ _mesh._tria.node (_node)->topo() = _iter->topo() ; + + _mesh._tria.node + (_node)->part() + = _iter->itag() ; } } } @@ -1268,12 +1288,11 @@ hits_func &_hfun ) { - /*------------------ tree-rect intersection predicate */ + /*------------------ tree-circ intersection predicate */ typedef - geom_tree::aabb_pred_rect_k < + geom_tree::aabb_pred_ball_2 < real_type, - iptr_type, - +2 > tree_pred ; + iptr_type > tree_pred ; /*------------------ ball-line intersection predicate */ typedef @@ -1281,16 +1300,8 @@ hits_func > hits_pred ; /*------------------ call actual intersection testing */ - real_type _rmin[ 2] = { - _ball._pmid[0] - _ball._rrad , - _ball._pmid[1] - _ball._rrad - } ; - real_type _rmax[ 2] = { - _ball._pmid[0] + _ball._rrad , - _ball._pmid[1] + _ball._rrad - } ; - - tree_pred _pred(_rmin, _rmax) ; + tree_pred _pred(_ball. _pmid, + _ball. _rrad) ; hits_pred _func(_ball. _pmid, _ball. _rrad, *this, _hfun) ; @@ -1317,10 +1328,9 @@ { /*------------------ tree-line intersection predicate */ typedef - geom_tree::aabb_pred_line_k < + geom_tree::aabb_pred_line_2 < real_type, - iptr_type, - +2 > tree_pred ; + iptr_type > tree_pred ; /*------------------ tria-line intersection predicate */ typedef @@ -1467,10 +1477,9 @@ /*-------------- tree-line intersection predicate */ typedef - geom_tree::aabb_pred_line_k < + geom_tree::aabb_pred_line_2 < real_type, - iptr_type, - +2 > tree_pred ; + iptr_type > tree_pred ; /*-------------- tria-line intersection predicate */ typedef @@ -1551,10 +1560,9 @@ /*-------------- tree-line intersection predicate */ typedef - geom_tree::aabb_pred_line_k < + geom_tree::aabb_pred_line_2 < real_type, - iptr_type, - +2 > tree_pred ; + iptr_type > tree_pred ; /*-------------- tria-line intersection predicate */ typedef diff --git a/src/libcpp/geom_type/geom_mesh_euclidean_3.hpp b/src/libcpp/geom_type/geom_mesh_euclidean_3.hpp index a46783b..a4fcb77 100644 --- a/src/libcpp/geom_type/geom_mesh_euclidean_3.hpp +++ b/src/libcpp/geom_type/geom_mesh_euclidean_3.hpp @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 19 January, 2019 + * Last updated: 02 March, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -415,12 +415,14 @@ _iter != _lsrc.tend() ; ++_iter ) { - if (this->_tria . - _set2[*_iter]. - feat() != null_feat) - { - _ldst.push_tail (*_iter) ; - } + auto _feat = + this-> + _tria._set2[ *_iter].feat() ; + + if (_feat != null_feat) + { + _ldst.push_tail(*_iter) ; + } } } @@ -428,19 +430,20 @@ typename list_type , typename geom_opts > - __normal_call char_type node_feat ( + __normal_call void_type node_feat ( iptr_type *_node , list_type &_aset , + char_type &_feat , + char_type &_topo , geom_opts &_opts ) { - char_type _feat = - (_aset.count() == +0) || - (_aset.count() == +2) - ? null_feat : soft_feat ; - real_type _DtoR = (real_type) +3.1415926536 / 180.0; + + real_type _ZERO = -1. + + std::numeric_limits + ::epsilon(); real_type _phi1 = (real_type)+180. - _opts.phi1(); @@ -454,9 +457,14 @@ __unreferenced(_node) ; + _feat = null_feat ; + _topo = (char_type)_aset.count(); + for (auto _ipos = _aset.head() ; _ipos != _aset.tend() ; ++_ipos ) + { + char_type _tbad = +1 ; for (auto _jpos = _ipos+1 ; _jpos != _aset.tend() ; ++_jpos ) @@ -502,6 +510,8 @@ else _acos *= (real_type)-1.; + if (_acos >= _ZERO) + { if (_acos <= _hard) { _feat = @@ -512,10 +522,23 @@ { _feat = std::max(_feat, soft_feat) ; + } + } + else + { + if (_tbad >= + 1 ) + { + _topo -= _tbad--; + } } - } - - return _feat ; + } + } + { + if (_topo != + 0 ) + if (_topo != + 2 ) + _feat = + std::max(_feat, soft_feat) ; + } } template < @@ -631,18 +654,20 @@ typename list_type , typename geom_opts > - __normal_call char_type edge_feat ( + __normal_call void_type edge_feat ( iptr_type *_enod , list_type &_aset , + char_type &_feat , + char_type &_topo , geom_opts &_opts ) { - char_type _feat = - (_aset.count() == +2) - ? null_feat : soft_feat ; - real_type _DtoR = (real_type) +3.1415926536 / 180.0; + + real_type _ZERO = -1. + + std::numeric_limits + ::epsilon(); real_type _phi2 = (real_type)+180. - _opts.phi2(); @@ -653,10 +678,15 @@ std::cos( _phi2 * _DtoR) ; real_type _soft = std::cos( _eta2 * _DtoR) ; + + _feat = null_feat ; + _topo = (char_type)_aset.count(); for (auto _ipos = _aset.head() ; _ipos != _aset.tend() ; ++_ipos ) + { + char_type _tbad = +1 ; for (auto _jpos = _ipos+1 ; _jpos != _aset.tend() ; ++_jpos ) @@ -755,6 +785,8 @@ else _acos *= (real_type)-1.; + if (_acos >= _ZERO) + { if (_acos <= _hard) { _feat = @@ -766,9 +798,22 @@ _feat = std::max(_feat, soft_feat) ; } + } + else + { + if (_tbad >= + 1 ) + { + _topo -= _tbad--; + } + } + } + } + { + if (_topo != + 0 ) + if (_topo != + 2 ) + _feat = + std::max(_feat, soft_feat) ; } - - return _feat ; } /* @@ -784,9 +829,17 @@ geom_opts &_opts ) { - containers::array _eadj ; - containers::array _fadj ; - containers::array _set2 ; + containers:: + array _eadj ; + containers:: + array _fadj ; + containers:: + array _ebnd ; + + containers:: + array _nmrk ; + containers:: + array _emrk ; /*---------------------------------- init. geom feat. */ for (auto _iter = @@ -827,29 +880,55 @@ } } - if (_opts.feat() ) - { /*---------------------------------- find sharp feat. */ + _nmrk.set_count ( + this->_tria._set1.count() , + containers::loose_alloc,-1) ; + + _emrk.set_count ( + this->_tria._set2.count() , + containers::loose_alloc,-1) ; + + iptr_type _nnum = +0 ; + iptr_type _enum = +0 ; + for (auto _epos = this->_tria._set2.head() ; _epos != this->_tria._set2.tend() ; - ++_epos ) + ++_epos, ++_enum) { /*---------------------------------- find sharp 1-dim */ if (_epos->mark() >= +0) { + if (_epos->itag() <= -1 || + _emrk[_enum] >= +1 || + _opts .feat() ) + { + /*---------------------------------- set geo.-defined */ _fadj.set_count (0); - this->_tria.edge_tri3 ( - &_epos->node(0), _fadj) ; - - _epos->topo () = - (char_type)_fadj.count() ; - - _epos->feat () = edge_feat( - &_epos->node(0), _fadj, - _opts) ; + this->_tria.edge_tri3 ( + &_epos->node (0), _fadj) ; + + edge_feat ( + &_epos->node (0), + _fadj , + _epos->feat () , + _epos->topo () , + _opts ) ; + + if (_epos->itag() <= -1) + { + /*---------------------------------- set user-defined */ + _epos->feat () = + std::max(_epos->feat () , + soft_feat) ; + + _nmrk[_epos->node(0)] = +1; + _nmrk[_epos->node(1)] = +1; + } + } } } @@ -857,33 +936,45 @@ this->_tria._set1.head() ; _npos != this->_tria._set1.tend() ; - ++_npos ) + ++_npos, ++_nnum) { /*---------------------------------- find sharp 0-dim */ - if (_npos->mark() >= +0) + if (_npos->mark() >= +0 ) + { + if (_npos->itag() <= -1 || + _nmrk[_nnum] >= +1 || + _opts .feat() ) { + /*---------------------------------- set geo.-defined */ _eadj.set_count (0); _fadj.set_count (0); - _set2.set_count (0); - - this->_tria.node_edge ( - &_npos->node(0), _eadj) ; + this->_tria.node_edge ( + &_npos->node (0), _eadj) ; - this->_tria.node_tri3 ( - &_npos->node(0), _fadj) ; + this->_tria.node_tri3 ( + &_npos->node (0), _fadj) ; + + _ebnd.set_count (0); - feat_list(_eadj, _set2); + feat_list (_eadj, _ebnd); - _npos->topo () = - (char_type)_set2.count() ; + node_feat ( + &_npos->node (0), + _ebnd , + _npos->feat () , + _npos->topo () , + _opts ) ; - _npos->feat () = node_feat( - &_npos->node(0), _set2, - _opts) ; + if (_npos->itag () <= -1) + { + /*---------------------------------- set user-defined */ + _npos->feat () = + std::max(_npos->feat () , + soft_feat) ; + } + } } - } - } for (auto _iter = @@ -992,7 +1083,7 @@ if (_imin < (iptr_type) +0 && _imax < (iptr_type) +0 ) __assert( _imin >= +0 && - "INIT-PART:-ve part index") ; + "GEOM::INIT-PART: -ve part index!") ; /*----------------------------- push unique PART id's */ auto _flag = @@ -1170,44 +1261,6 @@ } } ; - /*----------------------------- setup user-feat. face */ - for (auto _iter = - this->_tria._set2.head() ; - _iter != - this->_tria._set2.tend() ; - ++_iter ) - { - if (_iter->mark() >= +0 && - _iter->itag() <= -1 ) - { - this->_tria._set1[ - _iter->node(0)].itag() = -1 ; - - this->_tria._set1[ - _iter->node(1)].itag() = -1 ; - } - } - - for (auto _iter = - this->_tria._set3.head() ; - _iter != - this->_tria._set3.tend() ; - ++_iter ) - { - if (_iter->mark() >= +0 && - _iter->itag() <= -1 ) - { - this->_tria._set1[ - _iter->node(0)].itag() = -1 ; - - this->_tria._set1[ - _iter->node(1)].itag() = -1 ; - - this->_tria._set1[ - _iter->node(2)].itag() = -1 ; - } - } - /*----------------------------- init. aabb at -+ inf. */ for(auto _idim = 3; _idim-- != 0 ; ) { @@ -1313,8 +1366,7 @@ { if (_iter->mark() >= +0 ) { - if (_opts .feat() && - _iter->feat() != null_feat) + if (_iter->feat() != null_feat) { /*----------------------------- push any 'real' feat. */ iptr_type _node = -1 ; @@ -1332,6 +1384,10 @@ _mesh._tria.node (_node)->topo() = _iter->topo() ; + + _mesh._tria.node + (_node)->part() + = _iter->itag() ; } } else @@ -1353,6 +1409,10 @@ _mesh._tria.node (_node)->topo() = _iter->topo() ; + + _mesh._tria.node + (_node)->part() + = _iter->itag() ; } } } @@ -1965,12 +2025,11 @@ hits_func &_hfun ) { - /*------------------ tree-rect intersection predicate */ + /*------------------ tree-flat intersection predicate */ typedef - geom_tree::aabb_pred_rect_k < + geom_tree::aabb_pred_flat_3 < real_type, - iptr_type, - +3 > tree_pred ; + iptr_type > tree_pred ; /*------------------ dual-face intersection predicate */ typedef @@ -1978,7 +2037,9 @@ hits_func > hits_pred ; /*------------------ call actual intersection testing */ - tree_pred _pred(_flat. _rmin, + tree_pred _pred(_flat. _ppos, + _flat. _nvec, + _flat. _rmin, _flat. _rmax) ; hits_pred _func(_flat. _ppos, _flat. _nvec, @@ -2006,10 +2067,9 @@ { /*------------------ tree-line intersection predicate */ typedef - geom_tree::aabb_pred_line_k < + geom_tree::aabb_pred_line_3 < real_type, - iptr_type, - +3 > tree_pred ; + iptr_type > tree_pred ; /*------------------ tria-line intersection predicate */ typedef @@ -2043,31 +2103,20 @@ hits_func &_hfun ) { - /*------------------ tree-rect intersection predicate */ + /*------------------ tree-ball intersection predicate */ typedef - geom_tree::aabb_pred_rect_k < + geom_tree::aabb_pred_ball_3 < real_type, - iptr_type, - +3 > tree_pred ; + iptr_type > tree_pred ; /*------------------ ball-line intersection predicate */ typedef ball_line_pred < hits_func > hits_pred ; - /*------------------ call actual intersection testing */ - real_type _rmin[ 3] = { - _ball._pmid[0] - _ball._rrad , - _ball._pmid[1] - _ball._rrad , - _ball._pmid[2] - _ball._rrad - } ; - real_type _rmax[ 3] = { - _ball._pmid[0] + _ball._rrad , - _ball._pmid[1] + _ball._rrad , - _ball._pmid[2] + _ball._rrad - } ; - - tree_pred _pred(_rmin, _rmax) ; + /*------------------ call actual intersection testing */ + tree_pred _pred(_ball. _pmid, + _ball. _rrad) ; hits_pred _func(_ball. _pmid, _ball. _rrad, *this, _hfun) ; @@ -2093,12 +2142,11 @@ hits_func &_hfun ) { - /*------------------ tree-rect intersection predicate */ + /*------------------ tree-ball intersection predicate */ typedef - geom_tree::aabb_pred_rect_k < + geom_tree::aabb_pred_ball_3 < real_type, - iptr_type, - +3 > tree_pred ; + iptr_type > tree_pred ; /*------------------ circ_tria intersection predicate */ typedef @@ -2108,18 +2156,8 @@ __unreferenced(_sbal) ; /*------------------ call actual intersection testing */ - real_type _rmin[ 3] = { - _disc._pmid[0] - _disc._rrad , - _disc._pmid[1] - _disc._rrad , - _disc._pmid[2] - _disc._rrad - } ; - real_type _rmax[ 3] = { - _disc._pmid[0] + _disc._rrad , - _disc._pmid[1] + _disc._rrad , - _disc._pmid[2] + _disc._rrad - } ; - - tree_pred _pred(_rmin, _rmax) ; + tree_pred _pred(_disc. _pmid, + _disc. _rrad) ; hits_pred _func(_disc. _pmid, _disc. _nvec, _disc. _rrad, @@ -2266,10 +2304,9 @@ /*-------------- tree-line intersection predicate */ typedef - geom_tree::aabb_pred_line_k < + geom_tree::aabb_pred_line_3 < real_type, - iptr_type, - +3 > tree_pred ; + iptr_type > tree_pred ; /*-------------- tria-line intersection predicate */ typedef @@ -2361,10 +2398,9 @@ /*-------------- tree-line intersection predicate */ typedef - geom_tree::aabb_pred_line_k < + geom_tree::aabb_pred_line_3 < real_type, - iptr_type, - +3 > tree_pred ; + iptr_type > tree_pred ; /*-------------- tria-line intersection predicate */ typedef diff --git a/src/libcpp/iter_mesh/iter_divs_2.inc b/src/libcpp/iter_mesh/iter_divs_2.inc index f1be73a..c6493aa 100644 --- a/src/libcpp/iter_mesh/iter_divs_2.inc +++ b/src/libcpp/iter_mesh/iter_divs_2.inc @@ -31,9 +31,9 @@ * -------------------------------------------------------- * - * Last updated: 27 December, 2018 + * Last updated: 02 March, 2019 * - * Copyright 2013-2018 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -60,6 +60,7 @@ iter_opts &_opts , iptr_type _enum , bool_type &_okay , + iptr_type &_inew , iptr_list &_tset , iptr_list &_tnew , real_list &_tsrc , @@ -72,7 +73,7 @@ real_type _lmin = (real_type) +1.00E+00 , real_type _qinc - = (real_type) +1.00E-04 + = (real_type) +1.00E-02 ) { __unreferenced(_pred) ; // for MSVC... @@ -161,7 +162,9 @@ _tsrc[0] , _tsrc[1]) ; /*--------------------------------- exit if too good! */ - if (_TMIN>=_TLIM) return ; + real_type _TSQR = std::sqrt(_TLIM) ; + + if (_TMIN>_TSQR ) return ; /*--------------------------------- get adjacent ball */ real_type _ibal[_dims+1] ; @@ -267,6 +270,8 @@ iptr_type _nnew = _mesh.push_node(_ndat) ; + _inew = _nnew ; + auto _nptr = _mesh._set1.head() + _nnew ; @@ -337,27 +342,15 @@ move_node( _geom, _mesh , _hfun, _pred, _hval , - _opts, _nptr, +1 , - _move, _tnew, - _ttmp, _tdst, - _dtmp, _ddst, - _minC, _TLIM, - _minD, _DLIM ) ; - - if (_move > 0) continue ; - - move_node( _geom, _mesh , - _hfun, _pred, _hval , - _opts, _nptr, +2 , + _opts, _nptr, + _odt_kind, _move, _tnew, _ttmp, _tdst, _dtmp, _ddst, _minC, _TLIM, _minD, _DLIM ) ; - if (_move > 0) continue ; - - break ; + if (_move <= +0) break ; } /*--------------------------------- compare cost scr. */ @@ -366,11 +359,8 @@ loop_tscr(_mesh, _pred, _tnew, _tdst) ; - real_type constexpr - _GOOD = (real_type) +1.00; - move_okay(_tdst, _tsrc, _okay, - std::sqrt(_GOOD) , + _TSQR, _qinc) ; if (_okay) diff --git a/src/libcpp/iter_mesh/iter_dual_2.inc b/src/libcpp/iter_mesh/iter_dual_2.inc index f59ad80..e2d35d3 100644 --- a/src/libcpp/iter_mesh/iter_dual_2.inc +++ b/src/libcpp/iter_mesh/iter_dual_2.inc @@ -133,7 +133,7 @@ _node->pval (_dims); real_type _qlim = _qmin + - (real_type) +1.0E-004 ; + (real_type) +1.0E-003 ; _tnum = (iptr_type) +0 ; for (auto _tria = _tset.head(), diff --git a/src/libcpp/iter_mesh/iter_mesh_2.hpp b/src/libcpp/iter_mesh/iter_mesh_2.hpp index 11a66a4..f8604c1 100644 --- a/src/libcpp/iter_mesh/iter_mesh_2.hpp +++ b/src/libcpp/iter_mesh/iter_mesh_2.hpp @@ -31,9 +31,9 @@ * -------------------------------------------------------- * - * Last updated: 27 December, 2018 + * Last updated: 02 March, 2019 * - * Copyright 2013-2018 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -75,7 +75,14 @@ iptr_type static constexpr _dims = pred_type::_dims ; - + + char_type static + constexpr _odt_kind = +1 ; // optimal delaunay + char_type static + constexpr _cvt_kind = +2 ; // centroidal voro. + char_type static + constexpr _sdQ_kind = +3 ; // local cost grad. + typedef mesh::iter_params < real_type , iptr_type > iter_opts ; @@ -448,7 +455,7 @@ real_list &_hval , iter_opts &_opts , node_iter _node , - iptr_type _kind , + char_type _kind , bool_type &_okay , iptr_list &_tset , real_list &_told , @@ -484,7 +491,7 @@ real_list &_hval , iter_opts &_opts , node_iter _node , - iptr_type _kind , + char_type _kind , iptr_type &_move , iptr_list &_tset , real_list &_told , @@ -511,37 +518,40 @@ real_type _proj [_dims] = { (real_type) +0.0 } ; - real_type _ladj = (real_type) + 0.0 ; - - /*---------------- calc. line search direction vector */ - if (_kind == +1 ) - { - _odt_move_2 ( - _mesh, _hfun, _pred, - _hval, _tset, _node, - _line, _ladj) ; - } - else - if (_TMIN<=_TLIM) - { - grad_move_2 ( - _mesh, _hfun, _pred, - _tset, _node, _told, - _line, _ladj) ; - } - else { return ; } - - /*---------------- scale line search direction vector */ - real_type _llen = std:: - sqrt(_pred.length_sq(_line)) ; - - real_type _xtol = - (real_type)+.1 * _opts.qtol() ; - - if (_llen<= - _ladj * _xtol) return; - - real_type _scal = // overrelaxation + real_type _ladj = (real_type) + 0.0 ; + + /*---------------- calc. line search direction vector */ + if (_kind == _odt_kind) + { + _odt_move_2 ( + _mesh, _hfun, _pred, + _hval, _tset, _node, + _line, _ladj) ; + } + else + if (_kind == _sdQ_kind) + { + if (_TMIN<=_TLIM) + { + grad_move_2 ( + _mesh, _hfun, _pred, + _tset, _node, _told, + _line, _ladj) ; + } + else { return ; } + } + + /*---------------- scale line search direction vector */ + real_type _llen = std:: + sqrt(_pred.length_sq(_line)) ; + + real_type _xtol = + (real_type)+.1 * _opts.qtol() ; + + if (_llen<= + _ladj * _xtol) return; + + real_type _scal = // overrelaxation _llen * (real_type)5./3. ; /*---------------- do backtracking line search iter's */ @@ -783,13 +793,13 @@ cost_pair const&_idat , cost_pair const&_jdat ) const - { return _idat._cost < - _jdat._cost ; + { return + _idat._cost < _jdat._cost ; } } ; typedef containers:: - array cost_list ; + arraycost_list ; iptr_list _eset ; cost_list _sset ; @@ -936,14 +946,15 @@ ++_iter ) { auto _sift = std::min ( - (size_t) +8 , + (size_t) + 8, (size_t) (_aset.tend()-_iter) ) ; auto _next = _iter + std::rand() % _sift ; std::swap(*_iter,*_next); - } + } + } } @@ -1024,7 +1035,7 @@ /*---------------- attempt a GRAD-based smoothing */ move_dual( _geom, _mesh , _hfun, _pred, _hval , - _opts, _node, + _opts, _node, _okay, _tset, _dold, _dnew, _DMIN, _DLIM ) ; @@ -1090,7 +1101,8 @@ /*---------------- attempt a CCVT-style smoothing */ move_node( _geom, _mesh , _hfun, _pred, _hval , - _opts, _node, +1 , + _opts, _node, + _odt_kind, _okay, _tset, _told, _tnew, _dold, _dnew, @@ -1102,7 +1114,8 @@ /*---------------- attempt a GRAD-based smoothing */ move_node( _geom, _mesh , _hfun, _pred, _hval , - _opts, _node, +2 , + _opts, _node, + _sdQ_kind, _okay, _tset, _told, _tnew, _dold, _dnew, @@ -1348,48 +1361,38 @@ iptr_type static constexpr _DEG_MAX = (iptr_type) +8 ; - # define __markedge \ + # define __marknode \ init_mark( _mesh, _nmrk, _emrk, \ _tmrk, std::max(_imrk - 0, +0) ) ; \ if (std::abs( \ - _nmrk[_enod[0]])!= _imrk) \ - { \ - if (_nmrk[_enod[0]] >= +0) \ - { \ - _nmrk[_enod[0]] = +_imrk; \ - } \ - else \ - { \ - _nmrk[_enod[0]] = -_imrk; \ - } \ - _nset.push_tail(_enod[0]) ; \ - } \ - if (std::abs( \ - _nmrk[_enod[1]])!= _imrk) \ + _nmrk[_nnew])!= _imrk) \ { \ - if (_nmrk[_enod[1]] >= +0) \ + if (_nmrk[_nnew] >= +0) \ { \ - _nmrk[_enod[1]] = +_imrk; \ + _nmrk[_nnew] = +_imrk; \ } \ else \ { \ - _nmrk[_enod[1]] = -_imrk; \ + _nmrk[_nnew] = -_imrk; \ } \ - _nset.push_tail(_enod[1]) ; \ + _nset.push_tail(_nnew) ; \ } \ + + + _nzip = +0 ; _ndiv = +0 ; - _nzip = +0; _ndiv = +0 ; - - iptr_list _iset, _jset , _eset ; - iptr_list _aset, _bset , _cset ; - real_list _told, _tnew , _ttmp ; - real_list _dold, _dnew ; + iptr_list _aset, _bset, _cset ; + iptr_list _eset, _done; + iptr_list _iset, _jset; + real_list _told, _tnew, _ttmp ; + real_list _dold, _dnew; for (auto _node = _mesh._set1.count() ; _node-- != +0; ) { /*--------------------- scan nodes and zip//div edges */ - if (_mesh._set1[_node].mark() >= +0 && + if ( _mesh. + _set1[_node].mark () >= +0 && std::abs ( _nmrk[_node]) >= _imrk - 2 ) { @@ -1423,7 +1426,7 @@ _eadj != _tend ; _eadj += _einc ) { - auto _eptr = + auto _eptr = _mesh._set2.head() + *_eadj; iptr_type _enod[2] ; @@ -1441,12 +1444,15 @@ (real_type) -0.500 ; real_type _ltol = (real_type) +0.500 ; + + iptr_type _nnew = -1; if (_opts.div_()) _div_edge( _geom, _mesh, _hfun, _pred, _hval, _opts,*_eadj, - _move, _iset, _jset, + _move, _nnew, + _iset, _jset, _told, _tnew, _ttmp, _dold, _dnew, _TLIM, _DLIM, @@ -1454,23 +1460,26 @@ if (_move) { - __markedge; _ndiv += +1; break ; + __marknode; _ndiv += +1; break ; } } else { + iptr_type _nnew = -1; + if (_opts.div_()) _div_edge( _geom, _mesh, _hfun, _pred, _hval, _opts,*_eadj, - _move, _iset, _jset, + _move, _nnew, + _iset, _jset, _told, _tnew, _ttmp, _dold, _dnew, _TLIM, _DLIM) ; if (_move) { - __markedge; _ndiv += +1; break ; + __marknode; _ndiv += +1; break ; } } } @@ -1486,11 +1495,14 @@ real_type _ltol = (real_type) +2.000 ; + iptr_type _nnew = -1; + if (_opts.zip_()) _zip_edge( _geom, _mesh, _hfun, _pred, _hval, _opts,*_eadj, - _move, _iset, _jset, + _move, _nnew, + _iset, _jset, _aset, _bset, _cset, _told, _tnew, _ttmp, _dold, _dnew, @@ -1499,16 +1511,19 @@ if (_move) { - __markedge; _nzip += +1; break ; + __marknode; _nzip += +1; break ; } } else { + iptr_type _nnew = -1; + if (_opts.zip_()) _zip_edge( _geom, _mesh, _hfun, _pred, _hval, _opts,*_eadj, - _move, _iset, _jset, + _move, _nnew, + _iset, _jset, _aset, _bset, _cset, _told, _tnew, _ttmp, _dold, _dnew, @@ -1516,7 +1531,7 @@ if (_move) { - __markedge; _nzip += +1; break ; + __marknode; _nzip += +1; break ; } } } @@ -1526,7 +1541,7 @@ } } - # undef __markedge + # undef __marknode } @@ -1595,9 +1610,11 @@ # ifdef __use_timers typename std ::chrono:: - high_resolution_clock::time_point _ttic ; + high_resolution_clock:: + time_point _ttic ; typename std ::chrono:: - high_resolution_clock::time_point _ttoc ; + high_resolution_clock:: + time_point _ttoc ; typename std ::chrono:: high_resolution_clock _time; diff --git a/src/libcpp/iter_mesh/iter_node_2.inc b/src/libcpp/iter_mesh/iter_node_2.inc index 5717055..0d51696 100644 --- a/src/libcpp/iter_mesh/iter_node_2.inc +++ b/src/libcpp/iter_mesh/iter_node_2.inc @@ -423,6 +423,8 @@ _ladj += _lsqr ; } + _ladj = std::sqrt (_ladj/_tnum) ; + /*------------------ adj. gradients w.r.t. W: dQ / dw */ real_type _qbar , _qlow ; _qbar = (real_type) +1. ; @@ -432,7 +434,7 @@ iptr_type _hnum = +1 ; real_type _qlim = _qmin + - (real_type) +1.0E-004 ; + (real_type) +1.0E-003 ; _tnum = (iptr_type) +0 ; for (auto _tria = _tset.head() , @@ -545,7 +547,7 @@ } } - if (_tnum > +0) + if (_lnum > +0) { for (auto _idim = _dims; _idim-- != +0; ) { @@ -553,8 +555,6 @@ } _qlow /= _lnum ; - _ladj /= _tnum ; - _ladj = std::sqrt (_ladj) ; } if (_hnum > +0) diff --git a/src/libcpp/iter_mesh/iter_params.hpp b/src/libcpp/iter_mesh/iter_params.hpp index 575727b..6212e13 100644 --- a/src/libcpp/iter_mesh/iter_params.hpp +++ b/src/libcpp/iter_mesh/iter_params.hpp @@ -31,9 +31,9 @@ * -------------------------------------------------------- * - * Last updated: 14 December, 2017 + * Last updated: 02 March, 2019 * - * Copyright 2013-2017 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -174,57 +174,6 @@ } ; - - /* - -------------------------------------------------------- - * ITER-TIMERS: cpu timers for ITER-MESH-K - -------------------------------------------------------- - */ - - template < - typename R , - typename I - > - class iter_timers - { - public : - - typedef R real_type ; - typedef I iptr_type ; - - typedef iter_timers self_type ; - - real_type _iter_full = (real_type) +0. ; - - real_type _move_full = (real_type) +0. ; - real_type _topo_full = (real_type) +0. ; - real_type _zips_full = (real_type) +0. ; - - public : - - /*-------------------------------------- elapsed time */ - - # ifdef __use_timers - - __inline_call double time_span ( - typename std:: - chrono::high_resolution_clock - ::time_point const& _ttic, - typename std:: - chrono::high_resolution_clock - ::time_point const& _ttoc - ) - { - return (double)( - std::chrono::duration_cast< - std::chrono::microseconds > - (_ttoc-_ttic).count()) / +1.0E+06 ; - } - - # endif//__use_timers - - } ; - } diff --git a/src/libcpp/iter_mesh/iter_timers.hpp b/src/libcpp/iter_mesh/iter_timers.hpp new file mode 100644 index 0000000..72865ea --- /dev/null +++ b/src/libcpp/iter_mesh/iter_timers.hpp @@ -0,0 +1,107 @@ + + /* + -------------------------------------------------------- + * ITER-TIMERS: timers for ITER-MESH-K. + -------------------------------------------------------- + * + * This program may be freely redistributed under the + * condition that the copyright notices (including this + * entire header) are not removed, and no compensation + * is received through use of the software. Private, + * research, and institutional use is free. You may + * distribute modified versions of this code UNDER THE + * CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE + * TO IT IN THE SAME FILE REMAIN UNDER COPYRIGHT OF THE + * ORIGINAL AUTHOR, BOTH SOURCE AND OBJECT CODE ARE + * MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR + * NOTICE IS GIVEN OF THE MODIFICATIONS. Distribution + * of this code as part of a commercial system is + * permissible ONLY BY DIRECT ARRANGEMENT WITH THE + * AUTHOR. (If you are not directly supplying this + * code to a customer, and you are instead telling them + * how they can obtain it for free, then you are not + * required to make any arrangement with me.) + * + * Disclaimer: Neither I nor: Columbia University, The + * Massachusetts Institute of Technology, The + * University of Sydney, nor The National Aeronautics + * and Space Administration warrant this code in any + * way whatsoever. This code is provided "as-is" to be + * used at your own risk. + * + -------------------------------------------------------- + * + * Last updated: 02 March, 2019 + * + * Copyright 2013-2019 + * Darren Engwirda + * de2363@columbia.edu + * https://github.com/dengwirda/ + * + -------------------------------------------------------- + */ + +# pragma once + +# ifndef __ITER_TIMERS__ +# define __ITER_TIMERS__ + + namespace mesh { + + /* + -------------------------------------------------------- + * ITER-TIMERS: cpu timers for ITER-MESH-K + -------------------------------------------------------- + */ + + template < + typename R , + typename I + > + class iter_timers + { + public : + + typedef R real_type ; + typedef I iptr_type ; + + typedef iter_timers self_type ; + + real_type _iter_full = (real_type) +0. ; + + real_type _move_full = (real_type) +0. ; + real_type _topo_full = (real_type) +0. ; + real_type _zips_full = (real_type) +0. ; + + public : + + /*-------------------------------------- elapsed time */ + + # ifdef __use_timers + + __inline_call double time_span ( + typename std:: + chrono::high_resolution_clock + ::time_point const& _ttic, + typename std:: + chrono::high_resolution_clock + ::time_point const& _ttoc + ) + { + return (double)( + std::chrono::duration_cast< + std::chrono::microseconds > + (_ttoc-_ttic).count()) / +1.0E+06 ; + } + + # endif//__use_timers + + } ; + + } + + +# endif //__ITER_TIMERS__ + + + diff --git a/src/libcpp/iter_mesh/iter_zips_2.inc b/src/libcpp/iter_mesh/iter_zips_2.inc index 09f3b41..f619d66 100644 --- a/src/libcpp/iter_mesh/iter_zips_2.inc +++ b/src/libcpp/iter_mesh/iter_zips_2.inc @@ -31,9 +31,9 @@ * -------------------------------------------------------- * - * Last updated: 27 December, 2018 + * Last updated: 02 March, 2019 * - * Copyright 2013-2018 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -60,6 +60,7 @@ iter_opts &_opts , iptr_type _edge , bool_type &_okay , + iptr_type &_nnew , iptr_list &_iset , iptr_list &_jset , iptr_list &_aset , @@ -75,7 +76,7 @@ real_type _lmax = (real_type) +8.00E-01 , real_type _qinc - = (real_type) +1.00E-04 + = (real_type) +1.00E-02 ) { __unreferenced(_pred) ; // for MSVC... @@ -237,7 +238,9 @@ real_type _TMIN = std::min( _minA, _minB); - if (_TMIN>_TLIM) return; + real_type _TSQR = std::sqrt(_TLIM) ; + + if (_TMIN>_TSQR) return; /*--------------------------------- get adjacent ball */ real_type _pbal[_dims+1] = { @@ -339,6 +342,8 @@ iptr_type _inew = _mesh.push_node(_ndat) ; + + _nnew = _inew ; auto _nptr = _mesh._set1.head() + _inew ; @@ -440,27 +445,15 @@ move_node( _geom, _mesh , _hfun, _pred, _hval , - _opts, _nptr, +1 , + _opts, _nptr, + _odt_kind, _move, _cset, _ttmp, _tdst, _dtmp, _ddst, _minC, _TLIM, _minD, _DLIM ) ; - if (_move > 0) continue ; - - move_node( _geom, _mesh , - _hfun, _pred, _hval , - _opts, _nptr, +2 , - _move, _cset, - _ttmp, _tdst, - _dtmp, _ddst, - _minC, _TLIM, - _minD, _DLIM ) ; - - if (_move > 0) continue ; - - break ; + if (_move <= +0) break ; } /*--------------------------------- compare cost scr. */ @@ -469,11 +462,8 @@ loop_tscr(_mesh, _pred, _cset, _tdst) ; - real_type constexpr - _GOOD = (real_type) +1.00; - move_okay(_tdst, _tsrc, _okay, - std::sqrt(_GOOD) , + _TSQR, _qinc) ; if (_okay) diff --git a/src/libcpp/itermesh.hpp b/src/libcpp/itermesh.hpp index fb6153e..0c08cc4 100644 --- a/src/libcpp/itermesh.hpp +++ b/src/libcpp/itermesh.hpp @@ -31,9 +31,9 @@ * ------------------------------------------------------------ * - * Last updated: 13 August, 2018 + * Last updated: 02 March, 2019 * - * Copyright 2013-2018 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -56,6 +56,7 @@ # include "meshtype.hpp" # include "iter_mesh/iter_params.hpp" +# include "iter_mesh/iter_timers.hpp" # include "iter_mesh/iter_mesh_euclidean_2.hpp" # include "iter_mesh/iter_pred_euclidean_2.hpp" diff --git a/src/libcpp/libbasic.hpp b/src/libcpp/libbasic.hpp index 6fef1aa..fb48cd9 100644 --- a/src/libcpp/libbasic.hpp +++ b/src/libcpp/libbasic.hpp @@ -31,9 +31,9 @@ * ------------------------------------------------------------ * - * Last updated: 10 September, 2017 + * Last updated: 20 February, 2019 * - * Copyright 2013-2017 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -152,14 +152,28 @@ ------------------------------------------------------------ */ -# define __isflip(__i) ( (__i)<0) +# define __isflip(__i) ( (__i) < 0) -# define __doflip(__i) (-(__i)-2) - -# define __unflip(__i) (((__i)<0) \ - ? __doflip(__i) : (__i) ) +# define __doflip(__i) (-(__i) - 2) +# define __unflip(__i) (((__i) < 0) \ + ? __doflip(__i) : (__i) ) + +/* +------------------------------------------------------------ + * integer "flip" routines +------------------------------------------------------------ + */ +# define __setbit(x,b) ((x)|= (1ULL<<(b))) + +# define __popbit(x,b) ((x)&= ~(1ULL<<(b))) + +# define __flpbit(x,b) ((x)^= (1ULL<<(b))) + +# define __chkbit(x,b) (!!((x)&(1ULL<<(b))) ) + + # endif //__LIBBASIC__ diff --git a/src/libcpp/mesh_func/hfun_mesh_euclidean_2.hpp b/src/libcpp/mesh_func/hfun_mesh_euclidean_2.hpp index ab8f8f5..f0ababa 100644 --- a/src/libcpp/mesh_func/hfun_mesh_euclidean_2.hpp +++ b/src/libcpp/mesh_func/hfun_mesh_euclidean_2.hpp @@ -115,6 +115,10 @@ + 2 , tree_node, allocator > tree_type ; + + typedef geom_tree::aabb_pred_node_2 < + real_type, + iptr_type > tree_pred ; public : @@ -299,12 +303,9 @@ } } } + } ; - typedef geom_tree::aabb_pred_node_k < - real_type , - iptr_type , 2 > tree_pred ; - /* -------------------------------------------------------- * FIND-PRED: TRUE if PPOS is in TPOS. diff --git a/src/libcpp/mesh_func/hfun_mesh_euclidean_3.hpp b/src/libcpp/mesh_func/hfun_mesh_euclidean_3.hpp index ef0e3fa..45ef6c3 100644 --- a/src/libcpp/mesh_func/hfun_mesh_euclidean_3.hpp +++ b/src/libcpp/mesh_func/hfun_mesh_euclidean_3.hpp @@ -119,6 +119,10 @@ + 3 , tree_node, allocator > tree_type ; + + typedef geom_tree::aabb_pred_node_3 < + real_type, + iptr_type > tree_pred ; public : @@ -318,10 +322,6 @@ } ; - typedef geom_tree::aabb_pred_node_k < - real_type , - iptr_type , 3 > tree_pred ; - /* -------------------------------------------------------- * TRIA-PRED: TRUE if PPOS is in TPOS. diff --git a/src/libcpp/rdel_mesh/rdel_base_3.hpp b/src/libcpp/rdel_mesh/rdel_base_3.hpp index f40f5c6..f552b88 100644 --- a/src/libcpp/rdel_mesh/rdel_base_3.hpp +++ b/src/libcpp/rdel_mesh/rdel_base_3.hpp @@ -547,7 +547,7 @@ } /*--------------------------- test loc. intersections */ - auto _imin = _pred._list.tend() ; + auto _iful = _pred._list.tend() ; auto _imax = _pred._list.tend() ; real_type _RTOL = _rEPS*_ebal[3]; @@ -555,11 +555,11 @@ real_type _dmax = -std::numeric_limits ::infinity() ; - real_type _dmin = - +std::numeric_limits + real_type _dful = + -std::numeric_limits ::infinity() ; - bool_type _safe, _have = false ; + bool_type _safe ; for (auto _iter = _pred._list.head() ; _iter != @@ -586,20 +586,50 @@ /*--------------------------- keep furthest from ball */ if (_dsqr > _dmax ) { - _have = true ; _dmax = _dsqr ; _imax = _iter ; } - if (_dsqr < _dmin ) + if (_dsqr > _dful && + _safe ) { - _have = true ; - _dmin = _dsqr ; - _imin = _iter ; + _dful = _dsqr ; + _iful = _iter ; } } } - if ( _have ) + if (_iful != + _pred._list.tend() ) + { + /*--------------------------- keep best intersections */ + _sbal[ 0] = _iful->pval(0); + _sbal[ 1] = _iful->pval(1); + _sbal[ 2] = _iful->pval(2); + + _part = _iful->itag (); + + _hits = _iful->hits (); + _feat = _iful->feat (); + _topo = _iful->topo (); + + /*--------------------------- eval. surf. ball radius */ + _sbal[ 3]+= + geometry::lensqr_3d(_sbal , + &_mesh._tria. + node(_enod[0])->pval(0)) ; + _sbal[ 3]+= + geometry::lensqr_3d(_sbal , + &_mesh._tria. + node(_enod[1])->pval(0)) ; + + _sbal[ 3]/= (real_type) +2.; + + /*--------------------------- return restricted state */ + return ( true ) ; + } + else + if (_imax != + _pred._list.tend() ) { /*--------------------------- keep best intersections */ _sbal[ 0] = _imax->pval(0); @@ -807,7 +837,7 @@ _hset[5][1] = _fnod[2] ; /*--------------------------- test loc. intersections */ - auto _imin = _pred._list.tend() ; + auto _iful = _pred._list.tend() ; auto _imax = _pred._list.tend() ; real_type _RTOL = _rEPS*_fbal[3]; @@ -815,11 +845,11 @@ real_type _dmax = -std::numeric_limits ::infinity() ; - real_type _dmin = - +std::numeric_limits + real_type _dful = + -std::numeric_limits ::infinity() ; - bool_type _safe, _have = false ; + bool_type _safe ; for (auto _iter = _pred._list.head() ; _iter != @@ -846,20 +876,52 @@ /*--------------------------- keep furthest from ball */ if (_dsqr > _dmax ) { - _have = true ; _dmax = _dsqr ; _imax = _iter ; } - if (_dsqr < _dmin ) + if (_dsqr < _dful && + _safe ) { - _have = true ; - _dmin = _dsqr ; - _imin = _iter ; + _dful = _dsqr ; + _iful = _iter ; } } } - if ( _have ) + if (_iful != + _pred._list.tend() ) + { + /*--------------------------- keep best intersections */ + _sbal[ 0] = _iful->pval(0); + _sbal[ 1] = _iful->pval(1); + _sbal[ 2] = _iful->pval(2); + + _part = _iful->itag (); + _feat = _iful->feat (); + _topo = _iful->topo (); + + /*--------------------------- eval. surf. ball radius */ + _sbal[ 3]+= + geometry::lensqr_3d(_sbal , + &_mesh._tria. + node(_fnod[0])->pval(0)) ; + _sbal[ 3]+= + geometry::lensqr_3d(_sbal , + &_mesh._tria. + node(_fnod[1])->pval(0)) ; + _sbal[ 3]+= + geometry::lensqr_3d(_sbal , + &_mesh._tria. + node(_fnod[2])->pval(0)) ; + + _sbal[ 3]/= (real_type) +3.; + + /*--------------------------- return restricted state */ + return ( true ) ; + } + else + if (_imax != + _pred._list.tend() ) { /*--------------------------- keep best intersections */ _sbal[ 0] = _imax->pval(0); diff --git a/src/libcpp/rdel_mesh/rdel_cost_delfront_2.inc b/src/libcpp/rdel_mesh/rdel_cost_delfront_2.inc index 79b50bd..01d34a1 100644 --- a/src/libcpp/rdel_mesh/rdel_cost_delfront_2.inc +++ b/src/libcpp/rdel_mesh/rdel_cost_delfront_2.inc @@ -208,7 +208,7 @@ } /*------------------------- calc. refinement priority */ - _tdat._mark = +3 ; + _tdat._mark = +0 ; # ifdef __bias_BNDS for(_enum = +3; _enum-- != +0; ) diff --git a/src/libcpp/rdel_mesh/rdel_cost_delfront_3.inc b/src/libcpp/rdel_mesh/rdel_cost_delfront_3.inc index f31f738..24f1b9f 100644 --- a/src/libcpp/rdel_mesh/rdel_cost_delfront_3.inc +++ b/src/libcpp/rdel_mesh/rdel_cost_delfront_3.inc @@ -216,7 +216,7 @@ } /*------------------------- calc. refinement priority */ - _fdat._mark = +3 ; + _fdat._mark = +0 ; # ifdef __bias_BNDS for(_enum = +3; _enum-- != +0; ) @@ -459,7 +459,7 @@ } /*------------------------- calc. refinement priority */ - _tdat._mark = +4 ; + _tdat._mark = +0 ; # ifdef __bias_BNDS for(_fnum = +4; _fnum-- != +0; ) diff --git a/src/libcpp/rdel_mesh/rdel_mesh_2.hpp b/src/libcpp/rdel_mesh/rdel_mesh_2.hpp index 8f96eee..22d6d9f 100644 --- a/src/libcpp/rdel_mesh/rdel_mesh_2.hpp +++ b/src/libcpp/rdel_mesh/rdel_mesh_2.hpp @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 24 January, 2019 + * Last updated: 15 February, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -222,7 +222,8 @@ node_data const& _idat, node_data const& _jdat ) const - { return _idat._pass < _jdat._pass ; + { return _idat._node[0] < + _jdat._node[0] ; } } ; class ball_pred @@ -1324,35 +1325,35 @@ _dump.push("\n") ; _dump.push("\n") ; - _dump.push(" |TYPE-1| (edge) = ") ; + _dump.push(" EDGE-CIRC = ") ; _dump.push(std::to_string( - _enod[rdel_opts::circ_kind])) ; + _enod[rdel_opts::circ_kind])); _dump.push("\n") ; - _dump.push(" |TYPE-2| (edge) = ") ; + _dump.push(" EDGE-OFFH = ") ; _dump.push(std::to_string( - _enod[rdel_opts::offh_kind])) ; + _enod[rdel_opts::offH_kind])); _dump.push("\n") ; - _dump.push(" |TYPE-D| (edge) = ") ; + _dump.push(" EDGE-OFFT = ") ; _dump.push(std::to_string( - _enod[rdel_opts::disk_kind])) ; + _enod[rdel_opts::offT_kind])); _dump.push("\n") ; _dump.push("\n") ; - _dump.push(" |TYPE-1| (tria) = ") ; + _dump.push(" TRIA-CIRC = ") ; _dump.push(std::to_string( - _tnod[rdel_opts::circ_kind])) ; + _tnod[rdel_opts::circ_kind])); _dump.push("\n") ; - _dump.push(" |TYPE-2| (tria) = ") ; + _dump.push(" TRIA-SINK = ") ; _dump.push(std::to_string( - _tnod[rdel_opts::offh_kind])) ; + _tnod[rdel_opts::sink_kind])); _dump.push("\n") ; - _dump.push(" |TYPE-3| (tria) = ") ; + _dump.push(" TRIA-OFFH = ") ; _dump.push(std::to_string( - _tnod[rdel_opts::offc_kind])) ; + _tnod[rdel_opts::offH_kind])); _dump.push("\n") ; - _dump.push(" |TYPE-4| (tria) = ") ; + _dump.push(" TRIA-OFFC = ") ; _dump.push(std::to_string( - _tnod[rdel_opts::sink_kind])) ; + _tnod[rdel_opts::offC_kind])); _dump.push("\n") ; } diff --git a/src/libcpp/rdel_mesh/rdel_mesh_3.hpp b/src/libcpp/rdel_mesh/rdel_mesh_3.hpp index b4e7043..dca4ceb 100644 --- a/src/libcpp/rdel_mesh/rdel_mesh_3.hpp +++ b/src/libcpp/rdel_mesh/rdel_mesh_3.hpp @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 24 January, 2019 + * Last updated: 15 February, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -243,7 +243,8 @@ node_data const& _idat, node_data const& _jdat ) const - { return _idat._pass < _jdat._pass ; + { return _idat._node[0] < + _jdat._node[0] ; } } ; class ball_pred @@ -1168,9 +1169,16 @@ _mode = face_mode; _irDT = true; // init. new face in rDT + + /* + init_disc( _geom, _hfun, + _mesh, + _epro, _fpro, _pass, + _mode, _args) ; + */ init_rdel( _geom, _hfun, - _mesh, false, + _mesh, true, _nnew, _tnew, _edat, _escr, _fdat, _fscr, @@ -1606,57 +1614,57 @@ _dump.push("\n") ; _dump.push("\n") ; - _dump.push(" |TYPE-1| (edge) = "); + _dump.push(" EDGE-CIRC = ") ; _dump.push(std::to_string( _enod[rdel_opts::circ_kind])); _dump.push("\n") ; - _dump.push(" |TYPE-2| (edge) = "); + _dump.push(" EDGE-OFFH = ") ; _dump.push(std::to_string( - _enod[rdel_opts::offh_kind])); + _enod[rdel_opts::offH_kind])); _dump.push("\n") ; - _dump.push(" |TYPE-D| (edge) = "); + _dump.push(" EDGE-OFFT = ") ; _dump.push(std::to_string( - _enod[rdel_opts::disk_kind])); + _enod[rdel_opts::offT_kind])); _dump.push("\n") ; _dump.push("\n") ; - _dump.push(" |TYPE-1| (face) = "); + _dump.push(" FACE-CIRC = ") ; _dump.push(std::to_string( _fnod[rdel_opts::circ_kind])); _dump.push("\n") ; - _dump.push(" |TYPE-2| (face) = "); + _dump.push(" FACE-SINK = ") ; _dump.push(std::to_string( - _fnod[rdel_opts::offh_kind])); + _fnod[rdel_opts::sink_kind])); _dump.push("\n") ; - _dump.push(" |TYPE-3| (face) = "); + _dump.push(" FACE-OFFH = ") ; _dump.push(std::to_string( - _fnod[rdel_opts::offc_kind])); + _fnod[rdel_opts::offH_kind])); _dump.push("\n") ; - _dump.push(" |TYPE-4| (face) = "); + _dump.push(" FACE-OFFC = ") ; _dump.push(std::to_string( - _fnod[rdel_opts::sink_kind])); + _fnod[rdel_opts::offC_kind])); _dump.push("\n") ; - _dump.push(" |TYPE-D| (face) = "); + _dump.push(" FACE-OFFT = ") ; _dump.push(std::to_string( - _fnod[rdel_opts::disk_kind])); + _fnod[rdel_opts::offT_kind])); _dump.push("\n") ; _dump.push("\n") ; - _dump.push(" |TYPE-1| (tria) = "); + _dump.push(" TRIA-CIRC = ") ; _dump.push(std::to_string( _tnod[rdel_opts::circ_kind])); _dump.push("\n") ; - _dump.push(" |TYPE-2| (tria) = "); + _dump.push(" TRIA-SINK = ") ; _dump.push(std::to_string( - _tnod[rdel_opts::offh_kind])); + _tnod[rdel_opts::sink_kind])); _dump.push("\n") ; - _dump.push(" |TYPE-3| (tria) = "); + _dump.push(" TRIA-OFFH = ") ; _dump.push(std::to_string( - _tnod[rdel_opts::offc_kind])); + _tnod[rdel_opts::offH_kind])); _dump.push("\n") ; - _dump.push(" |TYPE-4| (tria) = "); + _dump.push(" TRIA-OFFC = ") ; _dump.push(std::to_string( - _tnod[rdel_opts::sink_kind])); + _tnod[rdel_opts::offC_kind])); _dump.push("\n") ; } diff --git a/src/libcpp/rdel_mesh/rdel_offh_delfront_2.inc b/src/libcpp/rdel_mesh/rdel_offh_delfront_2.inc index 0671dbe..2842bd7 100644 --- a/src/libcpp/rdel_mesh/rdel_offh_delfront_2.inc +++ b/src/libcpp/rdel_mesh/rdel_offh_delfront_2.inc @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 22 January, 2019 + * Last updated: 19 February, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -71,6 +71,9 @@ real_type static const _rtol = (real_type)+1.00E-03 ; + real_type static const _dtol = + (real_type)+8.75E-01 ; + real_type static const _epsb = (real_type)+9.50E-01 ; @@ -155,7 +158,7 @@ else // adv.-front limiter { _kind = - rdel_opts::offh_kind ; + rdel_opts::offH_kind ; } _dist -= _pert * _dist ; @@ -179,12 +182,18 @@ _mesh._tria.find_node ( // test voro. limiter &_ppos[0], _dual, _enod[ 0]) ; - if (_dual != _enod[0] ) + real_type _dlen = + geometry::length_2d( + _ppos, &_mesh._tria. + node(_dual)->pval(0)) ; + + if (_dual != _enod [ +0] ) + if (_dlen <= _dtol*_hval ) { return rdel_opts::null_kind ; } - - if (_move < _rtol * _hval) break ; + + if (_move <= _rtol*_hval ) break ; } return ( _kind ) ; @@ -259,14 +268,14 @@ for(iptr_type _iter = +0; _iter++ != +4 ; ) { - _kind = rdel_opts::offh_kind ; - _nsiz[2] = _hfun.eval(_ppos, _hint) ; real_type _hmid = _nsiz[0] * (real_type)1./4. + _nsiz[1] * (real_type)1./4. + _nsiz[2] * (real_type)2./4. ; + + _kind = rdel_opts::offH_kind ; real_type _dist ; real_type _dsqr = @@ -290,13 +299,13 @@ { // adv.-front limiter _dist = _alth * _hmid; _kind = - rdel_opts::offh_kind; + rdel_opts::offH_kind; } if (_dist >= _dvec[3]) // off-centre limiter { _dist = _dvec[3]; _kind = - rdel_opts::offc_kind; + rdel_opts::offC_kind; } if (_dist >= _dvec[2] * _epsb) { // circumball limiter diff --git a/src/libcpp/rdel_mesh/rdel_offh_delfront_3.inc b/src/libcpp/rdel_mesh/rdel_offh_delfront_3.inc index a79d96c..d2aef15 100644 --- a/src/libcpp/rdel_mesh/rdel_offh_delfront_3.inc +++ b/src/libcpp/rdel_mesh/rdel_offh_delfront_3.inc @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 22 January, 2019 + * Last updated: 19 February, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -71,6 +71,9 @@ real_type static const _rtol = (real_type)+1.00E-03 ; + real_type static const _dtol = + (real_type)+8.75E-01 ; + real_type static const _epsb = (real_type)+9.50E-01 ; @@ -164,7 +167,7 @@ else // adv.-front limiter { _kind = - rdel_opts::offh_kind ; + rdel_opts::offH_kind ; } _dist -= _pert * _dist ; @@ -189,12 +192,18 @@ _mesh._tria.find_node ( // test voro. limiter &_ppos[0], _dual, _enod[ 0]) ; - if (_dual != _enod[0] ) + real_type _dlen = + geometry::length_3d( + _ppos, &_mesh._tria. + node(_dual)->pval(0)) ; + + if (_dual != _enod [ +0] ) + if (_dlen <= _dtol*_hval ) { return rdel_opts::null_kind ; } - - if (_move < _rtol * _hval) break ; + + if (_move <= _rtol*_hval ) break ; } return (_kind ) ; @@ -230,6 +239,9 @@ real_type static const _rtol = (real_type)+1.00E-03 ; + real_type static const _dtol = + (real_type)+8.75E-01 ; + real_type static const _epsb = (real_type)+9.50E-01 ; @@ -308,7 +320,7 @@ real_type _near = (real_type)+.33 * _frad ; - _kind = rdel_opts::offh_kind ; + _kind = rdel_opts::offH_kind ; if (_dsqr < _near) { // min.-space limiter @@ -325,13 +337,13 @@ { // adv.-front limiter _dist = _alth * _hval; _kind = - rdel_opts::offh_kind; + rdel_opts::offH_kind; } if (_dist >= _dvec[4]) // off-centre limiter { _dist = _dvec[4]; _kind = - rdel_opts::offc_kind; + rdel_opts::offC_kind; } if (_dist >= _dvec[3] * _epsb) { // circumball limiter @@ -367,14 +379,20 @@ _mesh._tria.find_node ( // test voro. limiter &_ppos[0], _dual, _enod[ 0]) ; - - if (_dual != _enod[0] && - _dual != _enod[1] ) + + real_type _dlen = + geometry::length_3d( + _ppos, &_mesh._tria. + node(_dual)->pval(0)) ; + + if (_dual != _enod [ +0] && + _dual != _enod [ +1] ) + if (_dlen <= _dtol*_hmin ) { return rdel_opts::null_kind ; } - - if (_move < _rtol * _hmin) break ; + + if (_move <= _rtol*_hmin ) break ; } return (_kind ) ; @@ -465,6 +483,8 @@ _nsiz[2] * (real_type)1./6. + _nsiz[3] * (real_type)3./6. ; + _kind = rdel_opts::offH_kind ; + real_type _dist ; real_type _dsqr = _hmid*_hmid - _fbal[ 3] ; @@ -472,8 +492,6 @@ real_type _near = _fbal[3] * (real_type)1./3. ; - _kind = rdel_opts::offh_kind ; - if (_dsqr < _near) { // min.-space limiter _dist = std::numeric_limits @@ -489,14 +507,14 @@ { // adv.-front limiter _dist = _alth * _hmid; _kind = - rdel_opts::offh_kind; + rdel_opts::offH_kind; } if (_dist >= _dvec[4]) // off-centre limiter { _dist = _dvec[4]; _kind = - rdel_opts::offc_kind; + rdel_opts::offC_kind; } if (_dist >= _dvec[3] * _epsb) diff --git a/src/libcpp/rdel_mesh/rdel_params.hpp b/src/libcpp/rdel_mesh/rdel_params.hpp index edab000..faf08ab 100644 --- a/src/libcpp/rdel_mesh/rdel_params.hpp +++ b/src/libcpp/rdel_mesh/rdel_params.hpp @@ -31,9 +31,9 @@ * -------------------------------------------------------- * - * Last updated: 04 September, 2017 + * Last updated: 20 February, 2019 * - * Copyright 2013-2017 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -67,11 +67,16 @@ typedef mesh::rdel_params self_type ; - enum node_kind{ null_kind , - fail_kind, circ_kind , - offh_kind, offc_kind , - sink_kind, disk_kind , - last_kind} ; + enum node_kind { + null_kind = +0 , + fail_kind , + circ_kind , // "circ"-type refinement + sink_kind , // "sink"-type off-centre + offH_kind , // "size"-type off-centre + offC_kind , // "circ"-type off-centre + offE_kind , // "err."-type off-centre + offT_kind , // "topo"-type off-centre + last_kind } ; iptr_type _verb ; // logfile output verbosity @@ -87,6 +92,10 @@ iptr_type _dims ; // topo. dimensions to mesh + iptr_type _iter ; // max. num. refine iter. + + iptr_type _rule ; // rule for cell refinement + real_type _siz1 ; // 1-dim. element size mul. real_type _siz2 ; // 2-dim. element size mul. real_type _siz3 ; // 3-dim. element size mul. @@ -105,11 +114,9 @@ real_type _vol3 ; // volume-length ratio - bool_type _top1 ; // impose "1-manifold-ness" - bool_type _top2 ; // impose "2-manifold-ness" + bool_type _top1 ; // impose 1-"manifold-ness" + bool_type _top2 ; // impose 2-"manifold-ness" - iptr_type _iter ; // max. no. refinement iter. - public : __static_call @@ -121,20 +128,34 @@ __inline_call real_type init_siz2 ( ) { return .5 *(4./3. + - 2./(1.+std::sqrt(1./3.))) ; + 2. / (1.+std::sqrt(1./3.))) ; } __static_call __inline_call real_type init_siz3 ( ) { return .5 *(4./3. + - 2./(1.+std::sqrt(3./8.))) ; + 2. / (1.+std::sqrt(3./8.))) ; + } + + __static_call + __inline_call iptr_type init_rule ( + ) + { + iptr_type _rule = +0 ; + __setbit( _rule, offH_kind) ; + __setbit( _rule, offC_kind) ; + //__setbit( _rule, offT_kind) ; + __setbit( _rule, sink_kind) ; + + return _rule ; } __static_call __inline_call iptr_type init_iter ( ) { return iptr_type ( - std::numeric_limits::max()) ; + std::numeric_limits + ::max()) ; } public : @@ -156,6 +177,10 @@ _dims(iptr_type(+ 3)) , + _iter(init_iter()) , + + _rule(init_rule()) , + _siz1(init_siz1()) , _siz2(init_siz2()) , _siz3(init_siz3()) , @@ -175,9 +200,7 @@ _vol3(real_type(+.000)) , _top1(bool_type(false)) , - _top2(bool_type(false)) , - - _iter(init_iter()) + _top2(bool_type(false)) { // load default values } @@ -187,6 +210,11 @@ { return this->_verb ; } + __inline_call iptr_type & rule ( + ) + { return this->_rule ; + } + __inline_call iptr_type & iter ( ) { return this->_iter ; @@ -294,6 +322,11 @@ { return this->_verb ; } + __inline_call iptr_type const& rule ( + ) const + { return this->_rule ; + } + __inline_call iptr_type const& iter ( ) const { return this->_iter ; @@ -397,60 +430,6 @@ } ; - /* - -------------------------------------------------------- - * RDEL-TIMERS: cpu timers for RDEL-MESH-K - -------------------------------------------------------- - */ - - template < - typename R , - typename I - > - class rdel_timers - { - public : - - typedef R real_type ; - typedef I iptr_type ; - - typedef rdel_timers self_type ; - - real_type _mesh_seed = (real_type) +0. ; - real_type _node_init = (real_type) +0. ; - real_type _node_rule = (real_type) +0. ; - real_type _edge_init = (real_type) +0. ; - real_type _edge_rule = (real_type) +0. ; - real_type _face_init = (real_type) +0. ; - real_type _face_rule = (real_type) +0. ; - real_type _tria_init = (real_type) +0. ; - real_type _tria_rule = (real_type) +0. ; - - public : - - /*-------------------------------------- elapsed time */ - - # ifdef __use_timers - - __inline_call double time_span ( - typename std:: - chrono::high_resolution_clock - ::time_point const& _ttic, - typename std:: - chrono::high_resolution_clock - ::time_point const& _ttoc - ) - { - return (double)( - std::chrono::duration_cast< - std::chrono::microseconds > - (_ttoc-_ttic).count()) / +1.0E+06 ; - } - - # endif//__use_timers - - } ; - } # endif // __RDEL_PARAMS__ diff --git a/src/libcpp/rdel_mesh/rdel_pred_delaunay_2.hpp b/src/libcpp/rdel_mesh/rdel_pred_delaunay_2.hpp index cfedccd..71589be 100644 --- a/src/libcpp/rdel_mesh/rdel_pred_delaunay_2.hpp +++ b/src/libcpp/rdel_mesh/rdel_pred_delaunay_2.hpp @@ -58,10 +58,10 @@ * Procedia Engineering, 163, pp. 84-96, * http://dx.doi.org/10.1016/j.proeng.2016.11.024 * - * D. Engwirda, (2014): "Locally-optimal Delaunay- - * refinement and optimisation-based mesh generation", - * Ph.D. Thesis, School of Mathematics and Statistics, - * Univ. of Sydney. + * D. Engwirda, (2014): "Locally-optimal Delaunay- + * refinement and optimisation-based mesh generation", + * Ph.D. Thesis, School of Mathematics and Statistics, + * Univ. of Sydney. * http://hdl.handle.net/2123/13148 * * which is based on various previous works, including diff --git a/src/libcpp/rdel_mesh/rdel_pred_delaunay_3.hpp b/src/libcpp/rdel_mesh/rdel_pred_delaunay_3.hpp index 6036130..3a2a915 100644 --- a/src/libcpp/rdel_mesh/rdel_pred_delaunay_3.hpp +++ b/src/libcpp/rdel_mesh/rdel_pred_delaunay_3.hpp @@ -58,10 +58,10 @@ * Procedia Engineering, 163, pp. 84-96, * http://dx.doi.org/10.1016/j.proeng.2016.11.024 * - * D. Engwirda, (2014): "Locally-optimal Delaunay- - * refinement and optimisation-based mesh generation", - * Ph.D. Thesis, School of Mathematics and Statistics, - * Univ. of Sydney. + * D. Engwirda, (2014): "Locally-optimal Delaunay- + * refinement and optimisation-based mesh generation", + * Ph.D. Thesis, School of Mathematics and Statistics, + * Univ. of Sydney. * http://hdl.handle.net/2123/13148 * * which is based on various previous works, including diff --git a/src/libcpp/rdel_mesh/rdel_pred_delfront_2.hpp b/src/libcpp/rdel_mesh/rdel_pred_delfront_2.hpp index 01a41db..b357964 100644 --- a/src/libcpp/rdel_mesh/rdel_pred_delfront_2.hpp +++ b/src/libcpp/rdel_mesh/rdel_pred_delfront_2.hpp @@ -152,16 +152,16 @@ class edge_data { public : - iptr_type _mark ; - iptr_type _imax ; + iptr_type _mark = +0; + iptr_type _imax = +0; real_type _cost ; } ; class tria_data { public : - iptr_type _mark ; - iptr_type _imax ; + iptr_type _mark = +0; + iptr_type _imax = +0; real_type _cost ; } ; @@ -375,15 +375,18 @@ __assert( false && "EDGE-NODE: interior facet!") ; - return rdel_opts:: null_kind ; + return ( _kind ) ; } if (_kind == rdel_opts::null_kind ) { /*----------------------- attempt offcentre placement */ - _kind = edge_offh(_geom , - _size, _mesh , _enod , - _pmax, _ppos , _args ) ; + if(__chkbit(_args.rule(), + rdel_opts + ::offH_kind) ) + _kind = edge_offh( _geom, + _size, _mesh , _enod, + _pmax, _ppos , _args) ; } if (_kind == rdel_opts::null_kind || _kind == rdel_opts::circ_kind ) @@ -441,19 +444,19 @@ _tbal[1] = _mesh. _tria.tria(_tpos)->circ(1); - _tbal[2] = (real_type)+.0 ; + _tbal[2] = (real_type)+0. ; _tbal[2]+= - geometry::lensqr_2d (_tbal, - &_mesh._tria.node( - _tnod[0])->pval(0)) ; + geometry::lensqr_2d (_tbal, + &_mesh._tria.node( + _tnod[0])->pval( 0)) ; _tbal[2]+= - geometry::lensqr_2d (_tbal, - &_mesh._tria.node( - _tnod[1])->pval(0)) ; + geometry::lensqr_2d (_tbal, + &_mesh._tria.node( + _tnod[1])->pval( 0)) ; _tbal[2]+= - geometry::lensqr_2d (_tbal, - &_mesh._tria.node( - _tnod[2])->pval(0)) ; + geometry::lensqr_2d (_tbal, + &_mesh._tria.node( + _tnod[2])->pval( 0)) ; _tbal[2]/= (real_type)+3. ; @@ -552,19 +555,20 @@ /*--------------------------------- get face indexing */ iptr_type _enod [ +3]; mesh_type::tria_type::tria_type:: - face_node(_enod, _emin, +2, +1); + face_node(_enod, _emin, 2, 1) ; _enod[0] = _mesh._tria. - tria(_tpos)->node(_enod[0]); + tria(_tpos)->node( _enod[0]) ; _enod[1] = _mesh._tria. - tria(_tpos)->node(_enod[1]); + tria(_tpos)->node( _enod[1]) ; /*----------------------------------- calc. edge-ball */ real_type _ebal [ +3]; geometry::circ_ball_2d(_ebal, - &_mesh._tria. - node(_enod[0])->pval(0), - &_mesh._tria. - node(_enod[1])->pval(0)) ; + &_mesh. + _tria.node(_enod[0])->pval(0) , + &_mesh. + _tria.node(_enod[1])->pval(0) + ) ; /*------------------------- assemble "frontal" vector */ real_type _dvec[4] ; @@ -598,19 +602,25 @@ if (_kind == rdel_opts::null_kind ) { /*----------------------- attempt offcentre placement */ - _kind = tria_offh(_geom, - _size, _mesh , _enod, - _ebal, _tbal , _dvec, - _ppos, _args ) ; + if(__chkbit(_args.rule(), + rdel_opts + ::offH_kind) ) + _kind = tria_offh( _geom, + _size , _mesh, _enod, + _ebal , _tbal, _dvec, + _ppos , _args) ; } if (_kind == rdel_opts::null_kind || _kind == rdel_opts::circ_kind ) { /*----------------------- attempt sink-type placement */ - _kind = tria_sink(_geom, - _size, _mesh , - _tpos, _tnod , _tbal, - _ppos, _args ) ; + if(__chkbit(_args.rule(), + rdel_opts + ::sink_kind) ) + _kind = tria_sink( _geom, + _size , _mesh, + _tpos , _tnod, _tbal, + _ppos , _args) ; } if (_kind == rdel_opts::null_kind || _kind == rdel_opts::circ_kind ) diff --git a/src/libcpp/rdel_mesh/rdel_pred_delfront_3.hpp b/src/libcpp/rdel_mesh/rdel_pred_delfront_3.hpp index e99eecd..80c36c0 100644 --- a/src/libcpp/rdel_mesh/rdel_pred_delfront_3.hpp +++ b/src/libcpp/rdel_mesh/rdel_pred_delfront_3.hpp @@ -58,10 +58,10 @@ * Procedia Engineering, 163, pp. 84-96, * http://dx.doi.org/10.1016/j.proeng.2016.11.024 * - * D. Engwirda, (2014): "Locally-optimal Delaunay- - * refinement and optimisation-based mesh generation", - * Ph.D. Thesis, School of Mathematics and Statistics, - * Univ. of Sydney. + * D. Engwirda, (2014): "Locally-optimal Delaunay- + * refinement and optimisation-based mesh generation", + * Ph.D. Thesis, School of Mathematics and Statistics, + * Univ. of Sydney. * http://hdl.handle.net/2123/13148 * * which is based on various previous works, including @@ -87,13 +87,13 @@ * of "off-centre" refinement rules - a generalisation * of ideas introduced in: * - * S. Rebay, (1993): "Efficient Unstructured Mesh - * Generation by Means of Delaunay Triangulation and + * S. Rebay, (1993): "Efficient Unstructured Mesh + * Generation by Means of Delaunay Triangulation and * the Bowyer-Watson Algorithm", J. Comp. Phys., 106, * pp. 125-138 * https://doi.org/10.1006/jcph.1993.1097 * - * H. Erten, A. Üngör, (2009): "Quality Triangulations + * H. Erten, A. Üngör, (2009): "Quality Triangulations * with Locally Optimal Steiner Points", SIAM J. Sci. * Comp., 31, pp. 2103-2130, * https://doi.org/10.1137/080716748 @@ -151,24 +151,24 @@ class edge_data { public : - iptr_type _mark ; - iptr_type _imax ; + iptr_type _mark = +0; + iptr_type _imax = +0; real_type _cost ; } ; class face_data { public : - iptr_type _mark ; - iptr_type _imax ; + iptr_type _mark = +0; + iptr_type _imax = +0; real_type _cost ; } ; class tria_data { public : - iptr_type _mark ; - iptr_type _imax ; + iptr_type _mark = +0; + iptr_type _imax = +0; real_type _cost ; } ; @@ -508,15 +508,18 @@ __assert( false && "EDGE-NODE: interior facet!") ; - return rdel_opts:: null_kind ; + return ( _kind ) ; } if (_kind == rdel_opts::null_kind ) { /*----------------------- attempt offcentre placement */ - _kind = edge_offh(_geom , - _size, _mesh , _enod , - _pmax, _ppos , _args ) ; + if(__chkbit(_args.rule(), + rdel_opts + ::offH_kind) ) + _kind = edge_offh( _geom, + _size , _mesh, _enod, + _pmax , _ppos, _args) ; } if (_kind == rdel_opts::null_kind || _kind == rdel_opts::circ_kind ) @@ -580,7 +583,7 @@ __assert( false && "FACE-NODE: interior facet!") ; - return rdel_opts:: null_kind ; + return ( _kind ) ; } /*--------------------------------- find edge lengths */ @@ -759,20 +762,26 @@ if (_kind == rdel_opts::null_kind ) { /*----------------------- attempt offcentre placement */ - _kind = face_offh(_geom, - _size, _mesh , _enod, - _pmid, _evec , _dvec, - _pmax, - _ppos, _args ) ; + if(__chkbit(_args.rule(), + rdel_opts + ::offH_kind) ) + _kind = face_offh( _geom, + _size , _mesh, _enod, + _pmid , _evec, _dvec, + _pmax , + _ppos , _args) ; } if (_kind == rdel_opts::null_kind || _kind == rdel_opts::circ_kind ) { /*----------------------- attempt sink-type placement */ - _kind = face_sink(_geom, - _size, _mesh , _tadj, - _enod, _pmax , - _ppos, _args ) ; + if(__chkbit(_args.rule(), + rdel_opts + ::sink_kind) ) + _kind = face_sink( _geom, + _size , _mesh, _tadj, + _enod , _pmax, + _ppos , _args) ; } if (_kind == rdel_opts::null_kind || _kind == rdel_opts::circ_kind ) @@ -835,23 +844,23 @@ _tbal[2] = _mesh. _tria.tria(_tpos)->circ(2); - _tbal[3] = (real_type)+.0 ; + _tbal[3] = (real_type)+0. ; _tbal[3]+= - geometry::lensqr_3d (_tbal, - &_mesh._tria.node( - _tnod[0])->pval(0)) ; + geometry::lensqr_3d (_tbal, + &_mesh._tria.node( + _tnod[0])->pval( 0)) ; _tbal[3]+= - geometry::lensqr_3d (_tbal, - &_mesh._tria.node( - _tnod[1])->pval(0)) ; + geometry::lensqr_3d (_tbal, + &_mesh._tria.node( + _tnod[1])->pval( 0)) ; _tbal[3]+= - geometry::lensqr_3d (_tbal, - &_mesh._tria.node( - _tnod[2])->pval(0)) ; + geometry::lensqr_3d (_tbal, + &_mesh._tria.node( + _tnod[2])->pval( 0)) ; _tbal[3]+= - geometry::lensqr_3d (_tbal, - &_mesh._tria.node( - _tnod[3])->pval(0)) ; + geometry::lensqr_3d (_tbal, + &_mesh._tria.node( + _tnod[3])->pval( 0)) ; _tbal[3]/= (real_type)+4. ; @@ -1014,20 +1023,22 @@ /*----------------------------------- calc. edge-ball */ real_type _elen ; _elen = geometry::lensqr_3d ( - &_mesh._tria. - node(_enod[0])->pval(0), - &_mesh._tria. - node(_enod[1])->pval(0)) ; + &_mesh. + _tria.node(_enod[0])->pval(0) , + &_mesh. + _tria.node(_enod[1])->pval(0) + ) ; /*----------------------------------- calc. face-ball */ real_type _fbal [ +4]; geometry::circ_ball_3d(_fbal, - &_mesh._tria. - node(_fnod[0])->pval(0), - &_mesh._tria. - node(_fnod[1])->pval(0), - &_mesh._tria. - node(_fnod[2])->pval(0)) ; + &_mesh. + _tria.node(_fnod[0])->pval(0) , + &_mesh. + _tria.node(_fnod[1])->pval(0) , + &_mesh. + _tria.node(_fnod[2])->pval(0) + ) ; /*------------------------- assemble "frontal" vector */ real_type _dvec[5] ; @@ -1062,19 +1073,25 @@ if (_kind == rdel_opts::null_kind ) { /*----------------------- attempt offcentre placement */ - _kind = tria_offh(_geom, - _size, _mesh , _fnod, - _fbal, _tbal , _dvec, - _ppos, _args ) ; + if(__chkbit(_args.rule(), + rdel_opts + ::offH_kind) ) + _kind = tria_offh( _geom, + _size , _mesh, _fnod, + _fbal , _tbal, _dvec, + _ppos , _args) ; } if (_kind == rdel_opts::null_kind || _kind == rdel_opts::circ_kind ) { /*----------------------- attempt sink-type placement */ - _kind = tria_sink(_geom, - _size, _mesh , - _tpos, _tnod , _tbal, - _ppos, _args ) ; + if(__chkbit(_args.rule(), + rdel_opts + ::sink_kind) ) + _kind = tria_sink( _geom, + _size , _mesh, + _tpos , _tnod, _tbal, + _ppos , _args) ; } if (_kind == rdel_opts::null_kind || _kind == rdel_opts::circ_kind ) diff --git a/src/libcpp/rdel_mesh/rdel_refine_ball_3.inc b/src/libcpp/rdel_mesh/rdel_refine_ball_3.inc index 42a6a5a..c3e49f7 100644 --- a/src/libcpp/rdel_mesh/rdel_refine_ball_3.inc +++ b/src/libcpp/rdel_mesh/rdel_refine_ball_3.inc @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 25 January, 2019 + * Last updated: 15 February, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -412,5 +412,158 @@ } } + /* + -------------------------------------------------------- + * INIT-DISC: setup disc "collar". + -------------------------------------------------------- + */ + + __static_call + __normal_call void_type init_disc ( + geom_type &_geom , + hfun_type &_hfun , + mesh_type &_mesh , + typename + mesh_type::edge_list &_epro , + typename + mesh_type::face_list &_fpro , + iptr_type _mode , + iptr_type _pass , + rdel_opts &_args + ) + { + real_type static constexpr + _DTOL = (real_type) 7./4. ; + + /*--------------------- push finalised "collar" nodes */ + __unreferenced(_hfun) ; + __unreferenced(_mode) ; + __unreferenced(_pass) ; + __unreferenced(_args) ; + __unreferenced(_epro) ; + + return; + + for (auto _epos = + _mesh._eset._lptr.head () ; + _epos != + _mesh._eset._lptr.tend () ; + ++_epos ) + { + if ( *_epos == nullptr) continue ; + + for (auto _edge = *_epos ; + _edge != nullptr; + _edge = _edge->_next ) + { + if (_edge->_data. + _feat == hard_feat ) + { + + real_type _ebal[ 4]; + geometry + ::circ_ball_3d (_ebal , + &_mesh._tria.node( + _edge->_data. + _node[0])->pval(0), + &_mesh._tria.node( + _edge->_data. + _node[1])->pval(0) + ) ; + + real_type _evec[ 3]; + geometry::vector_3d( + &_mesh._tria.node( + _edge->_data. + _node[0])->pval(0), + &_mesh._tria.node( + _edge->_data. + _node[1])->pval(0), + _evec ) ; + + typename + geom_type::disc_type _hdat ; + _hdat._pmid[0] = _ebal[0] ; + _hdat._pmid[1] = _ebal[1] ; + _hdat._pmid[2] = _ebal[2] ; + + _hdat._nvec[0] = _evec[0] ; + _hdat._nvec[1] = _evec[1] ; + _hdat._nvec[2] = _evec[2] ; + + _hdat._rrad = + _DTOL* std::sqrt(_ebal[3]) ; + + mesh::keep_all_3d < + real_type , + iptr_type > _halo ; + + + //!! + real_type _junk[3] = {+0.} ; + + _geom. + intersect(_hdat, _junk, _halo) ; + + + for (auto _iter = + _halo._list.head() ; + _iter != + _halo._list.tend() ; + ++_iter ) + { + iptr_type _hint = + _mesh._tria.node(_edge-> + _data._node[0])->next(); + + /* + iptr_type _near = -1; + _mesh._tria.find_node( + &_iter->pval(0), _near, + _edge-> + _data._node[0]) ; + + if (_near != + _edge->_data._node[0] && + _near != + _edge->_data._node[1] ) + continue ; + */ + + iptr_type _node = -1; + if (_mesh._tria.push_node( + &_iter->pval(0), _node, + _hint ) ) + { + _mesh._tria.node( + _node)->idxh() = + hfun_type::null_hint(); + + _mesh._tria.node( + _node)->fdim() = 2 ; + + _mesh._tria.node( + _node)->feat() = 0 ; + + _mesh._tria.node( + _node)->topo() = 2 ; + + face_data _fdat; + _fdat._node[0] = + _edge ->_data._node [0]; + _fdat._node[1] = + _edge ->_data._node [1]; + _fdat._node[2] = _node ; + + _fpro.push(_fdat) ; + } + } + + } + + } + } + } + diff --git a/src/libcpp/rdel_mesh/rdel_refine_base_2.inc b/src/libcpp/rdel_mesh/rdel_refine_base_2.inc index 857ba0c..4ee8218 100644 --- a/src/libcpp/rdel_mesh/rdel_refine_base_2.inc +++ b/src/libcpp/rdel_mesh/rdel_refine_base_2.inc @@ -59,12 +59,12 @@ iptr_type _sign = -1 ; tria_data _tdat; - _tdat._node[0] = - _mesh._tria.tria(_tpos)->node(+0) ; - _tdat._node[1] = - _mesh._tria.tria(_tpos)->node(+1) ; - _tdat._node[2] = - _mesh._tria.tria(_tpos)->node(+2) ; + _tdat._node[0] = _mesh. + _tria.tria(_tpos)->node(+0) ; + _tdat._node[1] = _mesh. + _tria.tria(_tpos)->node(+1) ; + _tdat._node[2] = _mesh. + _tria.tria(_tpos)->node(+2) ; typename mesh_type:: tria_list:: diff --git a/src/libcpp/rdel_mesh/rdel_refine_base_3.inc b/src/libcpp/rdel_mesh/rdel_refine_base_3.inc index 8c0b305..1c80ccf 100644 --- a/src/libcpp/rdel_mesh/rdel_refine_base_3.inc +++ b/src/libcpp/rdel_mesh/rdel_refine_base_3.inc @@ -59,14 +59,14 @@ iptr_type _sign = -1; tria_data _tdat; - _tdat._node[0] = - _mesh._tria.tria(_tpos)->node(+0) ; - _tdat._node[1] = - _mesh._tria.tria(_tpos)->node(+1) ; - _tdat._node[2] = - _mesh._tria.tria(_tpos)->node(+2) ; - _tdat._node[3] = - _mesh._tria.tria(_tpos)->node(+3) ; + _tdat._node[0] = _mesh. + _tria.tria(_tpos)->node(+0) ; + _tdat._node[1] = _mesh. + _tria.tria(_tpos)->node(+1) ; + _tdat._node[2] = _mesh. + _tria.tria(_tpos)->node(+2) ; + _tdat._node[3] = _mesh. + _tria.tria(_tpos)->node(+3) ; typename mesh_type:: tria_list:: diff --git a/src/libcpp/rdel_mesh/rdel_refine_face_3.inc b/src/libcpp/rdel_mesh/rdel_refine_face_3.inc index b3bb5a4..3afc20b 100644 --- a/src/libcpp/rdel_mesh/rdel_refine_face_3.inc +++ b/src/libcpp/rdel_mesh/rdel_refine_face_3.inc @@ -449,7 +449,7 @@ if (_kind != rdel_opts::null_kind) { if (_tdim != +3) - { + { /*------------------------- re-queue face if un-split */ _tscr.push_tail(_qdat) ; } diff --git a/src/libcpp/rdel_mesh/rdel_refine_topo_2.inc b/src/libcpp/rdel_mesh/rdel_refine_topo_2.inc index 21f5887..b30661c 100644 --- a/src/libcpp/rdel_mesh/rdel_refine_topo_2.inc +++ b/src/libcpp/rdel_mesh/rdel_refine_topo_2.inc @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 25 January, 2019 + * Last updated: 18 February, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -117,11 +117,11 @@ continue ; if( _edup.find(_edat, _eptr)) continue ; - else - { _edup.push(_edat) ; - } + + _edup.push (_edat ) ; - _eset.push_tail(_mptr->_data ) ; + _eset.push_tail( + _mptr->_data ) ; } } @@ -140,21 +140,44 @@ geom_type &_geom, mesh_type &_mesh, edat_list &_eset, + iptr_type _npos, + real_type *_ppos, real_type *_pmax, char_type &_tdim, - iptr_type &_tadj + iptr_type &_tadj, + rdel_opts &_opts ) { typename rdel_opts::node_kind _kind = rdel_opts::null_kind ; - _pmax[2] = (real_type)+0.; + real_type _rEPS = std::pow ( + std::numeric_limits + ::epsilon (), +.75) ; + + real_type _vnrm[3] ; + _vnrm[0] = (real_type) +0. ; + _vnrm[1] = (real_type) +0. ; + _vnrm[2] = (real_type) +0. ; + + _pmax[2] = (real_type) +0. ; /*--------------------- find max. SDB in edge-set */ for (auto _edge = _eset.head() ; _edge != _eset.tend() ; ++_edge ) { + iptr_type _enod[ +3] ; + mesh_type::tria_type:: + tria_type::face_node( + _enod, _edge->_eadj, 2, 1) ; + _enod[0] =_mesh._tria. + tria(_edge->_tadj) + ->node( _enod[ 0]) ; + _enod[1] =_mesh._tria. + tria(_edge->_tadj) + ->node( _enod[ 1]) ; + char_type _feat, _topo ; iptr_type _part; real_type _fbal[ +3] ; @@ -167,7 +190,34 @@ _feat, _topo , _part ) ) __assert( false && "TOPO-NODE: edge-ball" ) ; - + + real_type _enrm[ +3] ; + geometry::line_norm_2d ( + &_mesh._tria. + node(_enod[0])->pval(0) , + &_mesh._tria. + node(_enod[1])->pval(0) , + _enrm) ; + + _enrm[2] = + geometry::length_2d(_enrm) ; + _enrm[0]/= _enrm[2] ; + _enrm[1]/= _enrm[2] ; + + real_type _sign = + geometry::dot_2d(_enrm, _vnrm) ; + + if (_sign >= (real_type) +0. ) + { + _vnrm[0]+= _enrm[0]; + _vnrm[1]+= _enrm[1]; + } + else + { + _vnrm[0]-= _enrm[0]; + _vnrm[1]-= _enrm[1]; + } + if (_pmax[2] < _sbal[2]) { _pmax[2] = _sbal[2]; @@ -183,9 +233,104 @@ _pmax[1] = _sbal[1]; } } + + if(__chkbit(_opts.rule(), + rdel_opts + ::offT_kind) ) + { + /*--------------------- calc. "skinny" off-centre */ + _vnrm[2] = + geometry::length_2d(_vnrm) ; + _vnrm[0]/= _vnrm[2] ; + _vnrm[1]/= _vnrm[2] ; + + real_type _rrad = + std::sqrt(_pmax[ 2]); + + if (_vnrm[2] >= _rEPS * _rrad) + { + /*--------------------- project on local "normal" */ + _vnrm[0]*= _rrad * + (real_type) +2. ; + _vnrm[1]*= _rrad * + (real_type) +2. ; + + typename + geom_type::line_type _ldat ; + _ldat._ipos[0] = _ppos[0] + - _vnrm[0]; + _ldat._ipos[1] = _ppos[1] + - _vnrm[1]; - /*--------------------- return steiner point kind */ - return ( _kind ) ; + _ldat._jpos[0] = _ppos[0] + + _vnrm[0]; + _ldat._jpos[1] = _ppos[1] + + _vnrm[1]; + + mesh::keep_all_2d < + real_type , + iptr_type > _halo ; + + _geom.intersect( _ldat , + _halo ) ; + + auto _ioff = _halo._list.tend(); + + real_type _dmin = + _rEPS * std::pow(_rrad, +2); + + real_type static constexpr + _NEAR = (real_type) +.75 ; + + real_type _best = + +std::numeric_limits + ::infinity(); + + for (auto _iter = + _halo._list.head() ; + _iter != + _halo._list.tend() ; + ++_iter ) + { + iptr_type _near = -1; + _mesh._tria.find_node ( + &_iter->pval(0), _near , _npos) ; + + auto _nptr = &_mesh. + _tria.node(_near)->pval(0) ; + + real_type _nsqr = + geometry::lensqr_2d ( + &_iter->pval(0), _nptr) ; + + real_type _dsqr = + geometry::lensqr_2d ( + &_iter->pval(0), _ppos) ; + + if (_nsqr > _dsqr * + _NEAR ) + if (_dsqr < _best && + _dsqr > _dmin ) + { + _kind = + rdel_opts::offT_kind ; + + _best = _dsqr ; + + _ioff = _iter ; + } + } + + if (_ioff != _halo._list.tend()) + { + _pmax[0] = _ioff->pval (0) ; + _pmax[1] = _ioff->pval (1) ; + } + + } + } + + return _kind ; } /* @@ -320,9 +465,14 @@ char_type _tmax = -1; iptr_type _hint = -1; + auto _ppos = + &_mesh . + _tria.node (_npos)->pval(0) ; + _kind = topo_node( _geom, _mesh , _eset, - _pmax , _tmax, _hint) ; + _npos , _ppos, _pmax, + _tmax , _hint, _args) ; /*------------------------- push node via constraints */ _kind = push_node( _geom, diff --git a/src/libcpp/rdel_mesh/rdel_refine_topo_3.inc b/src/libcpp/rdel_mesh/rdel_refine_topo_3.inc index d6aa98d..60286a2 100644 --- a/src/libcpp/rdel_mesh/rdel_refine_topo_3.inc +++ b/src/libcpp/rdel_mesh/rdel_refine_topo_3.inc @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 25 January, 2019 + * Last updated: 18 February, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -231,7 +231,7 @@ } } - return ( _eset.count() == _topo ) ; + return ( _eset.count() <= _topo ) ; //!! what about missing edges? } /* @@ -429,26 +429,49 @@ __static_call __normal_call - typename rdel_opts::node_kind topo_node ( + typename rdel_opts::node_kind etop_node ( geom_type &_geom, mesh_type &_mesh, edat_list &_eset, - fdat_list &_fset, + iptr_type _npos, + real_type *_ppos, real_type *_pmax, char_type &_tdim, - iptr_type &_tadj + iptr_type &_tadj, + rdel_opts &_opts ) { typename rdel_opts::node_kind _kind = rdel_opts::null_kind ; - _pmax[3] = (real_type)+0.; + real_type _rEPS = std::pow ( + std::numeric_limits + ::epsilon (), +.75) ; + + real_type _vdir[4] ; + _vdir[0] = (real_type) +0. ; + _vdir[1] = (real_type) +0. ; + _vdir[2] = (real_type) +0. ; + _vdir[3] = (real_type) +0. ; + + _pmax[3] = (real_type) +0. ; /*--------------------- find max. SDB in edge-set */ for (auto _edge = _eset.head() ; _edge != _eset.tend() ; ++_edge ) { + iptr_type _enod[ +3] ; + mesh_type::tria_type:: + tria_type::face_node( + _enod, _edge->_eadj, 3, 1) ; + _enod[0] =_mesh._tria. + tria(_edge->_tadj) + ->node( _enod[ 0]) ; + _enod[1] =_mesh._tria. + tria(_edge->_tadj) + ->node( _enod[ 1]) ; + char_type _hits; char_type _feat, _topo ; iptr_type _part; @@ -463,7 +486,37 @@ _topo, _part ) ) __assert( false && "TOPO-NODE: edge-ball" ) ; - + + real_type _edir[ +4] ; + geometry::vector_3d ( + &_mesh._tria. + node(_enod[0])->pval(0) , + &_mesh._tria. + node(_enod[1])->pval(0) , + _edir) ; + + _edir[3] = + geometry::length_3d(_edir) ; + _edir[0]/= _edir[3] ; + _edir[1]/= _edir[3] ; + _edir[2]/= _edir[3] ; + + real_type _sign = + geometry::dot_3d(_edir, _vdir) ; + + if (_sign >= (real_type) +0. ) + { + _vdir[0]+= _edir[0]; + _vdir[1]+= _edir[1]; + _vdir[2]+= _edir[2]; + } + else + { + _vdir[0]-= _edir[0]; + _vdir[1]-= _edir[1]; + _vdir[2]-= _edir[2]; + } + if (_pmax[3] < _sbal[3]) { _pmax[3] = _sbal[3]; @@ -480,11 +533,162 @@ _pmax[2] = _sbal[2]; } } + + if(__chkbit(_opts.rule(), + rdel_opts + ::offT_kind) ) + { + /*--------------------- calc. "skinny" off-centre */ + _vdir[3] = + geometry::length_3d(_vdir) ; + _vdir[0]/= _vdir[3] ; + _vdir[1]/= _vdir[3] ; + _vdir[2]/= _vdir[3] ; + + real_type _rrad = + std::sqrt(_pmax[ 3]); + + if (_vdir[3] >= _rEPS * _rrad) + { + /*--------------------- project on local "normal" */ + typename + geom_type::flat_type _fdat ; + _fdat._ppos[0] = _ppos[0]; + _fdat._ppos[1] = _ppos[1]; + _fdat._ppos[2] = _ppos[2]; + _fdat._nvec[0] = _vdir[0]; + _fdat._nvec[1] = _vdir[1]; + _fdat._nvec[2] = _vdir[2]; + + _fdat._rmin[0] = _ppos[0] + - (real_type)+2.*_rrad ; + _fdat._rmin[1] = _ppos[1] + - (real_type)+2.*_rrad ; + _fdat._rmin[2] = _ppos[2] + - (real_type)+2.*_rrad ; + + _fdat._rmax[0] = _ppos[0] + + (real_type)+2.*_rrad ; + _fdat._rmax[1] = _ppos[1] + + (real_type)+2.*_rrad ; + _fdat._rmax[2] = _ppos[2] + + (real_type)+2.*_rrad ; + + mesh::keep_all_3d < + real_type , + iptr_type > _halo ; + + _geom.intersect( _fdat , + _halo ) ; + + auto _ioff = _halo._list.tend(); + + real_type _dmin = + _rEPS * std::pow(_rrad, +2); + + real_type static constexpr + _NEAR = (real_type) +.75 ; + + real_type _best = + +std::numeric_limits + ::infinity(); + + for (auto _iter = + _halo._list.head() ; + _iter != + _halo._list.tend() ; + ++_iter ) + { + iptr_type _near = -1; + _mesh._tria.find_node ( + &_iter->pval(0), _near , _npos) ; + + auto _nptr = &_mesh. + _tria.node(_near)->pval(0) ; + + real_type _nsqr = + geometry::lensqr_3d ( + &_iter->pval(0), _nptr) ; + + real_type _dsqr = + geometry::lensqr_3d ( + &_iter->pval(0), _ppos) ; + + if (_nsqr > _dsqr * + _NEAR ) + if (_dsqr < _best && + _dsqr > _dmin ) + { + _kind = + rdel_opts::offT_kind ; + + _best = _dsqr ; + + _ioff = _iter ; + } + } + + if (_ioff != _halo._list.tend()) + { + _pmax[0] = _ioff->pval (0) ; + _pmax[1] = _ioff->pval (1) ; + _pmax[2] = _ioff->pval (2) ; + } + + } + } + + return _kind ; + } + + __static_call + __normal_call + typename rdel_opts::node_kind ftop_node ( + geom_type &_geom, + mesh_type &_mesh, + fdat_list &_fset, + iptr_type _npos, + real_type *_ppos, + real_type *_pmax, + char_type &_tdim, + iptr_type &_tadj, + rdel_opts &_opts + ) + { + typename rdel_opts::node_kind + _kind = rdel_opts::null_kind ; + + real_type _rEPS = std::pow ( + std::numeric_limits + ::epsilon (), +.75) ; + + real_type _vnrm[4] ; + _vnrm[0] = (real_type) +0. ; + _vnrm[1] = (real_type) +0. ; + _vnrm[2] = (real_type) +0. ; + _vnrm[3] = (real_type) +0. ; + + _pmax[3] = (real_type) +0. ; + /*--------------------- find max. SDB in face-set */ for (auto _face = _fset.head() ; _face != _fset.tend() ; ++_face ) { + iptr_type _fnod[ +4] ; + mesh_type::tria_type:: + tria_type::face_node( + _fnod, _face->_fadj, 3, 2) ; + _fnod[0] =_mesh._tria. + tria(_face->_tadj) + ->node( _fnod[ 0]) ; + _fnod[1] =_mesh._tria. + tria(_face->_tadj) + ->node( _fnod[ 1]) ; + _fnod[2] =_mesh._tria. + tria(_face->_tadj) + ->node( _fnod[ 2]) ; + char_type _feat, _topo ; iptr_type _part; real_type _fbal[ +4] ; @@ -497,7 +701,39 @@ _feat, _topo , _part) ) __assert( false && "TOPO-NODE: face-ball" ) ; - + + real_type _fnrm[ +4] ; + geometry::tria_norm_3d ( + &_mesh._tria. + node(_fnod[0])->pval(0) , + &_mesh._tria. + node(_fnod[1])->pval(0) , + &_mesh._tria. + node(_fnod[2])->pval(0) , + _fnrm) ; + + _fnrm[3] = + geometry::length_3d(_fnrm) ; + _fnrm[0]/= _fnrm[3] ; + _fnrm[1]/= _fnrm[3] ; + _fnrm[2]/= _fnrm[3] ; + + real_type _sign = + geometry::dot_3d(_fnrm, _vnrm) ; + + if (_sign >= (real_type) +0. ) + { + _vnrm[0]+= _fnrm[0]; + _vnrm[1]+= _fnrm[1]; + _vnrm[2]+= _fnrm[2]; + } + else + { + _vnrm[0]-= _fnrm[0]; + _vnrm[1]-= _fnrm[1]; + _vnrm[2]-= _fnrm[2]; + } + if (_pmax[3] < _sbal[3]) { _pmax[3] = _sbal[3]; @@ -514,8 +750,112 @@ _pmax[2] = _sbal[2]; } } - /*--------------------- return steiner point kind */ - return ( _kind ) ; + + if(__chkbit(_opts.rule(), + rdel_opts + ::offT_kind) ) + { + /*--------------------- calc. "skinny" off-centre */ + _vnrm[3] = + geometry::length_3d(_vnrm) ; + _vnrm[0]/= _vnrm[3] ; + _vnrm[1]/= _vnrm[3] ; + _vnrm[2]/= _vnrm[3] ; + + real_type _rrad = + std::sqrt(_pmax[ 3]); + + if (_vnrm[3] >= _rEPS * _rrad) + { + /*--------------------- project on local "normal" */ + _vnrm[0]*= _rrad * + (real_type) +2. ; + _vnrm[1]*= _rrad * + (real_type) +2. ; + _vnrm[2]*= _rrad * + (real_type) +2. ; + + typename + geom_type::line_type _ldat ; + _ldat._ipos[0] = _ppos[0] + - _vnrm[0]; + _ldat._ipos[1] = _ppos[1] + - _vnrm[1]; + _ldat._ipos[2] = _ppos[2] + - _vnrm[2]; + + _ldat._jpos[0] = _ppos[0] + + _vnrm[0]; + _ldat._jpos[1] = _ppos[1] + + _vnrm[1]; + _ldat._jpos[2] = _ppos[2] + + _vnrm[2]; + + mesh::keep_all_3d < + real_type , + iptr_type > _halo ; + + _geom.intersect( _ldat , + _halo ) ; + + auto _ioff = _halo._list.tend(); + + real_type _dmin = + _rEPS * std::pow(_rrad, +2); + + real_type static constexpr + _NEAR = (real_type) +.75 ; + + real_type _best = + +std::numeric_limits + ::infinity(); + + for (auto _iter = + _halo._list.head() ; + _iter != + _halo._list.tend() ; + ++_iter ) + { + iptr_type _near = -1; + _mesh._tria.find_node ( + &_iter->pval(0), _near , _npos) ; + + auto _nptr = &_mesh. + _tria.node(_near)->pval(0) ; + + real_type _nsqr = + geometry::lensqr_3d ( + &_iter->pval(0), _nptr) ; + + real_type _dsqr = + geometry::lensqr_3d ( + &_iter->pval(0), _ppos) ; + + if (_nsqr > _dsqr * + _NEAR ) + if (_dsqr < _best && + _dsqr > _dmin ) + { + _kind = + rdel_opts::offT_kind ; + + _best = _dsqr ; + + _ioff = _iter ; + } + } + + if (_ioff != _halo._list.tend()) + { + _pmax[0] = _ioff->pval (0) ; + _pmax[1] = _ioff->pval (1) ; + _pmax[2] = _ioff->pval (2) ; + } + + } + } + + return _kind ; } /* @@ -662,10 +1002,15 @@ real_type _pmax [ +4] ; char_type _tmax = -1; iptr_type _hint = -1; + + auto _ppos = + &_mesh . + _tria.node (_npos)->pval(0) ; - _kind = topo_node( _geom, - _mesh , _eset, _fset, - _pmax , _tmax, _hint) ; + _kind = etop_node( _geom, + _mesh , _eset, + _npos , _ppos, _pmax, + _tmax , _hint, _args) ; /*------------------------- push node via constraints */ _kind = push_node( _geom, @@ -828,12 +1173,15 @@ real_type _pmax [ +4] ; char_type _tmax = -1; iptr_type _hint = -1; + + auto _ppos = + &_mesh . + _tria.node (_npos)->pval(0) ; - _eset.set_count( +0) ; //!! - - _kind = topo_node( _geom, - _mesh , _eset, _fset, - _pmax , _tmax, _hint) ; + _kind = ftop_node( _geom, + _mesh , _fset, + _npos , _ppos, _pmax, + _tmax , _hint, _args) ; /*------------------------- push node via constraints */ _kind = push_node( _geom, diff --git a/src/libcpp/rdel_mesh/rdel_timers.hpp b/src/libcpp/rdel_mesh/rdel_timers.hpp new file mode 100644 index 0000000..1369bf8 --- /dev/null +++ b/src/libcpp/rdel_mesh/rdel_timers.hpp @@ -0,0 +1,110 @@ + + /* + -------------------------------------------------------- + * RDEL-TIMERS: CPU timers for RDEL-MESH-K. + -------------------------------------------------------- + * + * This program may be freely redistributed under the + * condition that the copyright notices (including this + * entire header) are not removed, and no compensation + * is received through use of the software. Private, + * research, and institutional use is free. You may + * distribute modified versions of this code UNDER THE + * CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE + * TO IT IN THE SAME FILE REMAIN UNDER COPYRIGHT OF THE + * ORIGINAL AUTHOR, BOTH SOURCE AND OBJECT CODE ARE + * MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR + * NOTICE IS GIVEN OF THE MODIFICATIONS. Distribution + * of this code as part of a commercial system is + * permissible ONLY BY DIRECT ARRANGEMENT WITH THE + * AUTHOR. (If you are not directly supplying this + * code to a customer, and you are instead telling them + * how they can obtain it for free, then you are not + * required to make any arrangement with me.) + * + * Disclaimer: Neither I nor: Columbia University, The + * Massachusetts Institute of Technology, The + * University of Sydney, nor The National Aeronautics + * and Space Administration warrant this code in any + * way whatsoever. This code is provided "as-is" to be + * used at your own risk. + * + -------------------------------------------------------- + * + * Last updated: 20 February, 2019 + * + * Copyright 2013-2019 + * Darren Engwirda + * de2363@columbia.edu + * https://github.com/dengwirda/ + * + -------------------------------------------------------- + */ + +# pragma once + +# ifndef __RDEL_TIMERS__ +# define __RDEL_TIMERS__ + + namespace mesh { + + /* + -------------------------------------------------------- + * RDEL-TIMERS: cpu timers for RDEL-MESH-K + -------------------------------------------------------- + */ + + template < + typename R , + typename I + > + class rdel_timers + { + public : + + typedef R real_type ; + typedef I iptr_type ; + + typedef rdel_timers self_type ; + + real_type _mesh_seed = (real_type) +0. ; + real_type _node_init = (real_type) +0. ; + real_type _node_rule = (real_type) +0. ; + real_type _edge_init = (real_type) +0. ; + real_type _edge_rule = (real_type) +0. ; + real_type _face_init = (real_type) +0. ; + real_type _face_rule = (real_type) +0. ; + real_type _tria_init = (real_type) +0. ; + real_type _tria_rule = (real_type) +0. ; + + public : + + /*-------------------------------------- elapsed time */ + + # ifdef __use_timers + + __inline_call double time_span ( + typename std:: + chrono::high_resolution_clock + ::time_point const& _ttic, + typename std:: + chrono::high_resolution_clock + ::time_point const& _ttoc + ) + { + return (double)( + std::chrono::duration_cast< + std::chrono::microseconds > + (_ttoc-_ttic).count()) / +1.0E+06 ; + } + + # endif//__use_timers + + } ; + + } + +# endif // __RDEL_TIMERS__ + + + diff --git a/src/libcpp/rdel_mesh/rdel_update_face_2.inc b/src/libcpp/rdel_mesh/rdel_update_face_2.inc index f1f6ad7..1dd48d6 100644 --- a/src/libcpp/rdel_mesh/rdel_update_face_2.inc +++ b/src/libcpp/rdel_mesh/rdel_update_face_2.inc @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 17 January, 2019 + * Last updated: 02 February, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -180,6 +180,8 @@ rdel_opts &_opts ) { + if(!_geom.have_feat(1))return ; + /*-------------------------------- correct node dims? */ iptr_type _fdim =+0; for (auto _node =+3; _node-- != +0; ) @@ -310,6 +312,8 @@ { /*-------------------------------- check "restricted" */ { + if(!_geom.have_feat(1)) return; + iptr_type _tnod[ +3] ; _tnod[0] = _mesh. _tria.tria(_tpos)->node(0); diff --git a/src/libcpp/rdel_mesh/rdel_update_face_3.inc b/src/libcpp/rdel_mesh/rdel_update_face_3.inc index 7d286a6..d6ac0cf 100644 --- a/src/libcpp/rdel_mesh/rdel_update_face_3.inc +++ b/src/libcpp/rdel_mesh/rdel_update_face_3.inc @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 17 January, 2019 + * Last updated: 02 February, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -223,6 +223,8 @@ rdel_opts &_opts ) { + if(!_geom.have_feat(1))return ; + /*-------------------------------- correct node dims? */ iptr_type _fdim =+0; for (auto _node =+4; _node-- != +0; ) @@ -342,6 +344,8 @@ rdel_opts &_opts ) { + if(!_geom.have_feat(2))return ; + /*-------------------------------- correct node dims? */ iptr_type _fdim =+0; for (auto _node =+4; _node-- != +0; ) @@ -477,6 +481,8 @@ ) { { + if(!_geom.have_feat(2)) return; + /*---------------------------- check "restricted" */ iptr_type _tnod[ +4] ; _tnod[0] = _mesh. diff --git a/src/libcpp/rdelmesh.hpp b/src/libcpp/rdelmesh.hpp index 568404d..1e1c0cf 100644 --- a/src/libcpp/rdelmesh.hpp +++ b/src/libcpp/rdelmesh.hpp @@ -31,9 +31,9 @@ * ------------------------------------------------------------ * - * Last updated: 29 December, 2018 + * Last updated: 20 February, 2019 * - * Copyright 2013-2018 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -74,6 +74,7 @@ # include "rdel_mesh/rdel_complex_3.hpp" # include "rdel_mesh/rdel_params.hpp" +# include "rdel_mesh/rdel_timers.hpp" # include "rdel_mesh/rdel_filt_k.hpp" diff --git a/src/libcpp/tessellate/del_tri_euclidean_2.hpp b/src/libcpp/tessellate/del_tri_euclidean_2.hpp index 0868a69..a596cad 100644 --- a/src/libcpp/tessellate/del_tri_euclidean_2.hpp +++ b/src/libcpp/tessellate/del_tri_euclidean_2.hpp @@ -31,7 +31,7 @@ * ------------------------------------------------------------ * - * Last updated: 22 January, 2019 + * Last updated: 18 February, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -96,29 +96,41 @@ __inline_call bool_type operator() ( mesh_type &_mesh, __const_ptr ( real_type) _ppos, - real_type _dmin, + real_type &_dmin, iptr_type &_nmin, iptr_type _tpos, - iptr_type &_fadj + iptr_type &_fmin ) const { bool_type _done = true ; iptr_type _fpos ; + iptr_type _fadj ; + iptr_type _tadj ; for(_fpos = 3; _fpos-- != 0; ) { /*--------------- check dist. wrt. k-th face apex */ + _tadj = _mesh. + tria(_tpos)->next(_fpos) ; + _fadj = _mesh. + tria(_tpos)->fpos(_fpos) ; + + _tadj = __unflip(_tadj) ; + + if(_tadj == _mesh.null_flag()) + continue ; + iptr_type _fnod[ 3] ; mesh_type::tria_type:: - face_node(_fnod, _fpos, 2, 1) ; + face_node(_fnod, _fadj, 2, 1) ; /* _fnod[0] = _mesh. - tria(_tpos)->node(_fnod[0]); + tria(_tadj)->node(_fnod[0]); _fnod[1] = _mesh. - tria(_tpos)->node(_fnod[1]); + tria(_tadj)->node(_fnod[1]); */ _fnod[2] = _mesh. - tria(_tpos)->node(_fnod[2]); + tria(_tadj)->node(_fnod[2]); iptr_type _apex = _fnod[2] ; @@ -131,7 +143,7 @@ _done = false; _dmin = _dsqr; _nmin = _apex; - _fadj = _fpos; + _fmin = _fpos; } } diff --git a/src/libcpp/tessellate/del_tri_euclidean_3.hpp b/src/libcpp/tessellate/del_tri_euclidean_3.hpp index 281aefb..88a2ac6 100644 --- a/src/libcpp/tessellate/del_tri_euclidean_3.hpp +++ b/src/libcpp/tessellate/del_tri_euclidean_3.hpp @@ -31,7 +31,7 @@ * ------------------------------------------------------------ * - * Last updated: 22 January, 2019 + * Last updated: 18 February, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -96,21 +96,33 @@ __inline_call bool_type operator() ( mesh_type &_mesh, __const_ptr ( real_type) _ppos, - real_type _dmin, + real_type &_dmin, iptr_type &_nmin, iptr_type _tpos, - iptr_type &_fadj + iptr_type &_fmin ) const { bool_type _done = true ; iptr_type _fpos ; + iptr_type _tadj ; + iptr_type _fadj ; for(_fpos = 4; _fpos-- != 0; ) { /*--------------- check dist. wrt. k-th face apex */ + _tadj = _mesh. + tria(_tpos)->next(_fpos) ; + _fadj = _mesh. + tria(_tpos)->fpos(_fpos) ; + + _tadj = __unflip(_tadj) ; + + if(_tadj == _mesh.null_flag()) + continue ; + iptr_type _fnod[ 4] ; mesh_type::tria_type:: - face_node(_fnod, _fpos, 3, 2) ; + face_node(_fnod, _fadj, 3, 2) ; /* _fnod[0] = _mesh. tria(_tpos)->node(_fnod[0]); @@ -120,7 +132,7 @@ tria(_tpos)->node(_fnod[2]); */ _fnod[3] = _mesh. - tria(_tpos)->node(_fnod[3]); + tria(_tadj)->node(_fnod[3]); iptr_type _apex = _fnod[3] ; @@ -133,7 +145,7 @@ _done = false; _dmin = _dsqr; _nmin = _apex; - _fadj = _fpos; + _fmin = _fpos; } } diff --git a/src/libcpp/tessellate/delaunay_tri_k.hpp b/src/libcpp/tessellate/delaunay_tri_k.hpp index 92815e1..fb86c72 100644 --- a/src/libcpp/tessellate/delaunay_tri_k.hpp +++ b/src/libcpp/tessellate/delaunay_tri_k.hpp @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 22 January, 2019 + * Last updated: 18 February, 2019 * * Copyright 2013-2019 * Darren Engwirda diff --git a/src/libcpp/tessellate/delaunay_walk_mesh.inc b/src/libcpp/tessellate/delaunay_walk_mesh.inc index 5fb9bc2..1a3ccbd 100644 --- a/src/libcpp/tessellate/delaunay_walk_mesh.inc +++ b/src/libcpp/tessellate/delaunay_walk_mesh.inc @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 22 January, 2019 + * Last updated: 18 February, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -136,8 +136,9 @@ iptr_type _next = +0 ; iptr_type _face = +0 ; - for (auto _iter = _ntri; - _iter-- != +0; _elem = _next) + for (auto _iter = this->_tset.count () ; + _iter-- != +0; + _elem = _next) { /*--------------- test node-in-tria condition */ _done = typename tria_pred:: @@ -178,12 +179,9 @@ { /*------------- walk mesh until we find nearest vert. */ bool_type _done = false; - bool_type _scan = false; - + iptr_type _elem ; iptr_type _nnod = (iptr_type)this->_nset.count () ; - iptr_type _ntri = - (iptr_type)this->_tset.count () ; if (this->_nset.count() == +0 || this->_tset.count() == +0 ) @@ -194,83 +192,54 @@ if (_hint < +0 || _hint>=_nnod) { /*----------------------------- not a valid index */ - _node = +0 ; - _scan = true ; + _elem = +0 ; } else if (node(_hint)->mark() <= -1 ) { /*----------------------------- not a valid _node */ - _node = +0 ; - _scan = true ; + _elem = +0 ; } else { - _node = _hint ; + _elem = node(_hint)->next(); } - - /*------------- randomised scan: locate initial _node */ - if (_scan) - { - - real_type _best = - +std::numeric_limits - ::infinity () ; - - /*--------------------- scale no. iter with dims. */ - iptr_type _imax = (iptr_type) - std::pow ((double) _nnod, - +1.0 / tria_pred::_dims) / 4 + 1 ; - for (auto _iter = _imax; _iter-- != +0; ) - { - /*----------------------- randomise selection */ - iptr_type _next = - (iptr_type)std::rand() - % this->_nset.count(); + walk_tria_near ( + _ppos, _elem, _elem); - /*----------------------- reject "null" nodes */ - if (node(_next)->mark() < +0) - continue ; + /*----------------------------- nearest from cell */ + real_type _dmin = +std :: + numeric_limits::infinity () ; - /*----------------------- calc. dist. to node */ - real_type _dist = - tria_pred::lensqr_kd(_ppos, - &node(_next)->pval(0)); + for (auto _inod = tria_pred::_dims + 1 ; + _inod-- != +0 ; ) + { + auto _next = tria(_elem) + ->node(_inod) ; - if (_best > _dist) + auto _dsqr = + tria_pred::lensqr_kd( + _ppos , &node(_next)->pval(0)) ; + + if (_dsqr < _dmin) { - /*----------------------- keep closest so far */ - _best = _dist; + _dmin = _dsqr; _node = _next; } } - } - - /*------------- ensure that we have a non--null _node */ - for ( ; _node != _nnod; ++_node) - { - if (node(_node)->mark()>=+0) break ; - } - - /*------------- walk elements toward _ppos from _node */ - iptr_type _elem = - node(_node)->next () ; - - real_type _dmin = - tria_pred::lensqr_kd(_ppos, - &node(_node)->pval(0)); - + /*------------- walk elements toward _ppos from _elem */ iptr_type _next = +0 ; iptr_type _face = +0 ; - for (auto _iter = _ntri; - _iter-- != +0; _elem = _next) + for (auto _iter = this->_tset.count (); + _iter-- != +0 ; + _elem = _next) { /*--------------- test nearest-node condition */ _done = typename tria_pred:: - template near_pred() ( + template near_pred()( *this, _ppos, _dmin, _node, _elem, _face) ;