From b312e2d4c31b83d2eb2701e943a90d8c9c6e4587 Mon Sep 17 00:00:00 2001 From: dengwirda Date: Fri, 25 Jan 2019 21:41:58 -0500 Subject: [PATCH] 0.9.8 updates: "Off-centre" fix -- dual cell checks --- example.jig | 13 + inc/lib_jigsaw.h | 25 +- src/jigsaw.cpp | 6 +- src/libcpp/hj_solver/hj_solver_2.hpp | 41 ++ src/libcpp/libbasic.hpp | 2 - src/libcpp/rdel_mesh/rdel_base_2.hpp | 69 ++- src/libcpp/rdel_mesh/rdel_base_3.hpp | 213 ++++----- src/libcpp/rdel_mesh/rdel_make_2.hpp | 105 +---- src/libcpp/rdel_mesh/rdel_make_3.hpp | 109 +---- src/libcpp/rdel_mesh/rdel_mesh_2.hpp | 175 ++++---- src/libcpp/rdel_mesh/rdel_mesh_3.hpp | 275 +++++------- src/libcpp/rdel_mesh/rdel_offh_delfront_2.inc | 21 +- src/libcpp/rdel_mesh/rdel_offh_delfront_3.inc | 39 +- src/libcpp/rdel_mesh/rdel_refine_ball_2.inc | 223 ++++++---- src/libcpp/rdel_mesh/rdel_refine_ball_3.inc | 228 ++++++---- src/libcpp/rdel_mesh/rdel_refine_base_2.inc | 36 +- src/libcpp/rdel_mesh/rdel_refine_base_3.inc | 44 +- src/libcpp/rdel_mesh/rdel_refine_face_2.inc | 8 +- src/libcpp/rdel_mesh/rdel_refine_face_3.inc | 24 +- src/libcpp/rdel_mesh/rdel_refine_topo_2.inc | 144 +++--- src/libcpp/rdel_mesh/rdel_refine_topo_3.inc | 407 +++++++++-------- src/libcpp/rdel_mesh/rdel_test_bounds_2.inc | 69 ++- src/libcpp/rdel_mesh/rdel_test_bounds_3.inc | 130 +++--- src/libcpp/rdel_mesh/rdel_update_topo_2.inc | 223 ++++++++++ src/libcpp/rdel_mesh/rdel_update_topo_3.inc | 416 ++++++++++++++++++ src/libcpp/tessellate/del_tri_euclidean_2.hpp | 64 ++- src/libcpp/tessellate/del_tri_euclidean_3.hpp | 66 ++- ...y_walk_tria.inc => delaunay_scan_tria.inc} | 6 +- src/libcpp/tessellate/delaunay_tri_k.hpp | 57 ++- ...y_walk_node.inc => delaunay_walk_mesh.inc} | 145 +++++- src/marche.hpp | 279 ++++++++++++ 31 files changed, 2373 insertions(+), 1289 deletions(-) create mode 100644 src/libcpp/hj_solver/hj_solver_2.hpp create mode 100644 src/libcpp/rdel_mesh/rdel_update_topo_2.inc create mode 100644 src/libcpp/rdel_mesh/rdel_update_topo_3.inc rename src/libcpp/tessellate/{delaunay_walk_tria.inc => delaunay_scan_tria.inc} (95%) rename src/libcpp/tessellate/{delaunay_walk_node.inc => delaunay_walk_mesh.inc} (55%) create mode 100644 src/marche.hpp diff --git a/example.jig b/example.jig index 7b4fbcf..6614696 100644 --- a/example.jig +++ b/example.jig @@ -49,6 +49,19 @@ MESH_FILE = out/bunny.msh +# +# OPTIONAL fields (INIT): +# ---------------------- +# + +# ---> INIT_FILE - 'INITNAME.MSH', a string containing the +# name of the initial distribution file (is required +# at input). +# + +# INIT_FILE = *.msh + + # # OPTIONAL fields (GEOM): # ---------------------- diff --git a/inc/lib_jigsaw.h b/inc/lib_jigsaw.h index b959a6a..b372aee 100644 --- a/inc/lib_jigsaw.h +++ b/inc/lib_jigsaw.h @@ -14,7 +14,7 @@ * JIGSAW: Interface to the JIGSAW meshing library. -------------------------------------------------------- * - * Last updated: 19 January, 2019 + * Last updated: 22 January, 2019 * * Copyright 2013 -- 2019 * Darren Engwirda @@ -91,8 +91,8 @@ */ jigsaw_msh_t *_init, - /* HFUN (OPTIONAL): mesh-spacing function H(x). - * => NULL for empty H(x). + /* HFUN (OPTIONAL): mesh-spacing function h(x). + * => NULL for empty h(x). */ jigsaw_msh_t *_hfun, @@ -103,7 +103,7 @@ /* -------------------------------------------------------- - * generate rDEL via JIGSAW. + * compute rDT's via TRIPOD. -------------------------------------------------------- */ @@ -126,6 +126,23 @@ */ jigsaw_msh_t *_mesh ) ; + + /* + -------------------------------------------------------- + * limit |df/dx| via MARCHE. + -------------------------------------------------------- + */ + + extern indx_t marche ( + + /* JCFG (REQUIRED): settings obj. definition. + */ + jigsaw_jig_t *_jcfg, + + /* HFUN (REQUIRED): apply limiter to |df/dx|. + */ + jigsaw_msh_t *_ffun + ) ; /* -------------------------------------------------------- diff --git a/src/jigsaw.cpp b/src/jigsaw.cpp index 90108e6..a2a78d8 100644 --- a/src/jigsaw.cpp +++ b/src/jigsaw.cpp @@ -34,7 +34,7 @@ -------------------------------------------------------- * * JIGSAW release 0.9.8.x - * Last updated: 19 January, 2019 + * Last updated: 22 January, 2019 * * Copyright 2013 -- 2019 * Darren Engwirda @@ -148,7 +148,7 @@ * gradient-based optimisation: * * L. Chen, J.C. Xu, (2004): "Optimal Delaunay - * triangulations, J. Comp. Math., 22, pp. 299–308, + * triangulations, J. Comp. Math., 22, pp. 299-308, * https://www.jstor.org/stable/43693155 * * L.A. Freitag, C. Ollivier-Gooch, (1997): "Tetrahedral @@ -756,7 +756,7 @@ # include "jigsaw.hpp" # include "tripod.hpp" -// # include "marche.hpp" + # include "marche.hpp" // # include "stitch.hpp" diff --git a/src/libcpp/hj_solver/hj_solver_2.hpp b/src/libcpp/hj_solver/hj_solver_2.hpp new file mode 100644 index 0000000..47a5c3a --- /dev/null +++ b/src/libcpp/hj_solver/hj_solver_2.hpp @@ -0,0 +1,41 @@ + + template < + typename F + > + class hj_mesh_2 + { + + __static_call + __normal_call void_type limit_edge_2 ( + ) + { + } + + __static_call + __normal_call void_type limit_tria_3 ( + ) + { + } + + __static_call + __normal_call void_type limit_mesh ( + ffun_type &_ffun, + real_type _DFDX + ) + { + + containers::priority_map<> _sort; + + // push nodes onto pq + // pop min at each pass + // limit any cell neighbours + // if a cell is limited, update pq of its nodes + + // that's it! + + } + + } ; + + + diff --git a/src/libcpp/libbasic.hpp b/src/libcpp/libbasic.hpp index 7036966..6fef1aa 100644 --- a/src/libcpp/libbasic.hpp +++ b/src/libcpp/libbasic.hpp @@ -60,8 +60,6 @@ # pragma warning(disable:4503) // decorated name length # pragma warning(disable:4458) // shadowing -//!! also the damn MSVC non-standard "extensions" somehow - # elif defined(__LLVM__) # elif defined(__GNUC__) diff --git a/src/libcpp/rdel_mesh/rdel_base_2.hpp b/src/libcpp/rdel_mesh/rdel_base_2.hpp index c678952..d2ac909 100644 --- a/src/libcpp/rdel_mesh/rdel_base_2.hpp +++ b/src/libcpp/rdel_mesh/rdel_base_2.hpp @@ -29,6 +29,15 @@ * way whatsoever. This code is provided "as-is" to be * used at your own risk. * + -------------------------------------------------------- + * + * Last updated: 24 January, 2019 + * + * Copyright 2013-2019 + * Darren Engwirda + * de2363@columbia.edu + * https://github.com/dengwirda/ + * -------------------------------------------------------- * * This class defines the basic "restricted" delaunay @@ -44,10 +53,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 * * building on various previous works on rDT methods, @@ -128,8 +137,7 @@ __normal_call double half_sign ( __const_ptr (double) _pp, __const_ptr (double) _pa, - __const_ptr (double) _pb, - double _rt + __const_ptr (double) _pb ) { /*-------- helper: eval. sign w.r.t. half-plane [a,b] */ @@ -155,34 +163,7 @@ _ab[0] * _mp[0] + _ab[1] * _mp[1] ; - if (_dp < -2.0 * _rt || - _dp > +2.0 * _rt ) return ((double) _dp ); - - /*-------- fall-back to double-double if near to zero */ - dd_flt _PM[2]; - _PM[0] = _pa[0] * 0.5 ; - _PM[1] = _pa[1] * 0.5 ; - _PM[0]+= _pb[0] * 0.5 ; - _PM[1]+= _pb[1] * 0.5 ; - - dd_flt _AB[2]; - _AB[0] = _pb[0] ; - _AB[1] = _pb[1] ; - _AB[0]-= _pa[0] ; - _AB[1]-= _pa[1] ; - - dd_flt _MP[2]; - _MP[0] = _pp[0] ; - _MP[1] = _pp[1] ; - _MP[0]-= _pm[0] ; - _MP[1]-= _pm[1] ; - - dd_flt _DP = - _AB[0] * _MP[0] + - _AB[1] * _MP[1] ; - - return ((double) _DP ); } template < @@ -222,20 +203,19 @@ _BPOS[1] = _mesh. _tria.node(_bnod)->pval(1) ; - double _sign = +0.0; + double _sign = +0.0 ; if (_bnod > _anod) _sign = + half_sign ( (double*) _PPOS, (double*) _APOS, - (double*) _BPOS, - (double ) _rtol) ; + (double*) _BPOS) ; else _sign = - half_sign ( (double*) _PPOS, (double*) _BPOS, - (double*) _APOS, - (double ) _rtol) ; + (double*) _APOS) ; + /*-------- "fatten" dual cavity to filter imprecision */ if (_sign >= -_rtol && _sign <= +_rtol) _safe = false; @@ -265,7 +245,7 @@ iptr_type &_part ) { - real_type const _rEPS = + real_type static const _rEPS = std::pow(std::numeric_limits ::epsilon(),+.67); @@ -417,9 +397,9 @@ _pred._list.tend() ; ++_iter ) { - if (clip_dual(_mesh, _hset , - &_iter->pval( 0), _safe, - _RTOL) ) + if (clip_dual( _mesh, _hset , + &_iter->pval( 0), + _safe, _RTOL) ) { /*--------------------------- dist to face circumball */ real_type _dsqr = @@ -427,7 +407,10 @@ _ebal , &_iter->pval( 0)) ; - if(!_safe && _dsqr > _RTOL) + real_type _dtol = + (real_type) 1./3. * _ebal[2] ; + + if(!_safe && _dsqr > _dtol) /*--------------------------- prune near-degeneracies */ continue ; diff --git a/src/libcpp/rdel_mesh/rdel_base_3.hpp b/src/libcpp/rdel_mesh/rdel_base_3.hpp index d09c799..f40f5c6 100644 --- a/src/libcpp/rdel_mesh/rdel_base_3.hpp +++ b/src/libcpp/rdel_mesh/rdel_base_3.hpp @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 09 January, 2019 + * Last updated: 24 January, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -138,76 +138,7 @@ geom_type::disc_type disc_type ; typedef typename geom_type::ball_type ball_type ; - - /* - -------------------------------------------------------- - * EDGE-LOOP: assemble tria.'s adj. to edge. - -------------------------------------------------------- - */ - template < - typename list_type - > - __static_call - __normal_call void_type edge_loop ( - mesh_type &_mesh, - iptr_type *_enod, - iptr_type _tadj, - iptr_type _fadj, - list_type &_loop - ) - { - iptr_type _tcur = _tadj ; - iptr_type _fcur = _fadj ; - - iptr_type _null = - _mesh._tria.null_flag (); - - while (true) - { - _loop.push_tail ( _tcur); - - iptr_type _tpos = _tcur ; - iptr_type _fpos ; - iptr_type _tmrk ; - for(_fpos = +4; _fpos-- != +0; ) - { - if(_fpos == _fcur) continue; - - iptr_type _fnod[ +4] ; - mesh_type::tria_type:: - tria_type:: - face_node(_fnod, _fpos, 3, 2); - _fnod[ +0] = _mesh._tria. - tria(_tpos)->node(_fnod[ 0]); - _fnod[ +1] = _mesh._tria. - tria(_tpos)->node(_fnod[ 1]); - _fnod[ +2] = _mesh._tria. - tria(_tpos)->node(_fnod[ 2]); - - iptr_type _same = +0; - if (_fnod[0] == _enod[0] || - _fnod[0] == _enod[1] ) - _same += +1 ; - if (_fnod[1] == _enod[0] || - _fnod[1] == _enod[1] ) - _same += +1 ; - if (_fnod[2] == _enod[0] || - _fnod[2] == _enod[1] ) - _same += +1 ; - - if (_same == +2) { break ; } - } - - _mesh._tria.find_pair( - _tpos, _tcur, - _fpos, _fcur, _tmrk) ; - - if (_tcur == _tadj) { break ; } - if (_tcur == _null) { break ; } - } - } - /* -------------------------------------------------------- * CLIP-DUAL: test pt. wrt dual halfplanes. @@ -218,8 +149,7 @@ __normal_call double half_sign ( __const_ptr (double) _pp, __const_ptr (double) _pa, - __const_ptr (double) _pb, - double _rt + __const_ptr (double) _pb ) { /*-------- helper: eval. sign w.r.t. half-plane [a,b] */ @@ -252,41 +182,7 @@ _ab[1] * _mp[1] + _ab[2] * _mp[2] ; - if (_dp < -2.0 * _rt || - _dp > +2.0 * _rt ) return ((double) _dp ); - - /*-------- fall-back to double-double if near to zero */ - dd_flt _PM[3]; - _PM[0] = _pa[0] * 0.5 ; - _PM[1] = _pa[1] * 0.5 ; - _PM[2] = _pa[2] * 0.5 ; - _PM[0]+= _pb[0] * 0.5 ; - _PM[1]+= _pb[1] * 0.5 ; - _PM[2]+= _pb[2] * 0.5 ; - - dd_flt _AB[3]; - _AB[0] = _pb[0] ; - _AB[1] = _pb[1] ; - _AB[2] = _pb[2] ; - _AB[0]-= _pa[0] ; - _AB[1]-= _pa[1] ; - _AB[2]-= _pa[2] ; - - dd_flt _MP[3]; - _MP[0] = _pp[0] ; - _MP[1] = _pp[1] ; - _MP[2] = _pp[2] ; - _MP[0]-= _pm[0] ; - _MP[1]-= _pm[1] ; - _MP[2]-= _pm[2] ; - - dd_flt _DP = - _AB[0] * _MP[0] + - _AB[1] * _MP[1] + - _AB[2] * _MP[2] ; - - return ((double) _DP ); } template < @@ -336,15 +232,14 @@ _sign = + half_sign ( (double*) _PPOS, (double*) _APOS, - (double*) _BPOS, - (double ) _rtol) ; + (double*) _BPOS) ; else _sign = - half_sign ( (double*) _PPOS, (double*) _BPOS, - (double*) _APOS, - (double ) _rtol) ; + (double*) _APOS) ; + /*-------- "fatten" dual cavity to filter imprecision */ if (_sign >= -_rtol && _sign <= +_rtol) _safe = false; @@ -355,6 +250,76 @@ return true ; } + /* + -------------------------------------------------------- + * EDGE-LOOP: assemble tria.'s adj. to edge. + -------------------------------------------------------- + */ + + template < + typename list_type + > + __static_call + __normal_call void_type edge_loop ( + mesh_type &_mesh, + iptr_type *_enod, + iptr_type _tadj, + iptr_type _fadj, + list_type &_loop + ) + { + /*----------------- assemble edge-adj. cells in-order */ + iptr_type _tcur = _tadj ; + iptr_type _fcur = _fadj ; + + iptr_type _null = + _mesh._tria.null_flag (); + + while (true) + { + _loop.push_tail ( _tcur); + + iptr_type _tpos = _tcur ; + iptr_type _fpos ; + iptr_type _tmrk ; + for(_fpos = +4; _fpos-- != +0; ) + { + if(_fpos == _fcur) continue; + + iptr_type _fnod[ +4] ; + mesh_type::tria_type:: + tria_type:: + face_node(_fnod, _fpos, 3, 2); + _fnod[ +0] = _mesh._tria. + tria(_tpos)->node(_fnod[ 0]); + _fnod[ +1] = _mesh._tria. + tria(_tpos)->node(_fnod[ 1]); + _fnod[ +2] = _mesh._tria. + tria(_tpos)->node(_fnod[ 2]); + + iptr_type _same = +0; + if (_fnod[0] == _enod[0] || + _fnod[0] == _enod[1] ) + _same += +1 ; + if (_fnod[1] == _enod[0] || + _fnod[1] == _enod[1] ) + _same += +1 ; + if (_fnod[2] == _enod[0] || + _fnod[2] == _enod[1] ) + _same += +1 ; + + if (_same == +2) { break ; } + } + + _mesh._tria.find_pair( + _tpos, _tcur, + _fpos, _fcur, _tmrk) ; + + if (_tcur == _tadj) { break ; } + if (_tcur == _null) { break ; } + } + } + /* -------------------------------------------------------- * EDGE-BALL: calc. edge-based circumballs. @@ -375,7 +340,7 @@ iptr_type &_part ) { - real_type const _rEPS = + real_type static const _rEPS = std::pow(std::numeric_limits ::epsilon(),+.67); @@ -601,9 +566,9 @@ _pred._list.tend() ; ++_iter ) { - if (clip_dual(_mesh, _hset , - &_iter->pval( 0), _safe, - _RTOL) ) + if (clip_dual( _mesh, _hset , + &_iter->pval( 0), + _safe, _RTOL) ) { /*--------------------------- dist to face circumball */ real_type _dsqr = @@ -611,7 +576,10 @@ _ebal , &_iter->pval( 0)) ; - if(!_safe && _dsqr > _RTOL) + real_type _dtol = + (real_type) 1./3. * _ebal[3] ; + + if(!_safe && _dsqr > _dtol) /*--------------------------- prune near-degeneracies */ continue; @@ -682,7 +650,7 @@ iptr_type &_part ) { - real_type const _rEPS = + real_type static const _rEPS = std::pow(std::numeric_limits ::epsilon(),+.67); @@ -858,9 +826,9 @@ _pred._list.tend() ; ++_iter ) { - if (clip_dual(_mesh, _hset , - &_iter->pval( 0), _safe, - _RTOL) ) + if (clip_dual( _mesh, _hset , + &_iter->pval( 0), + _safe, _RTOL) ) { /*--------------------------- dist to face circumball */ real_type _dsqr = @@ -868,7 +836,10 @@ _fbal , &_iter->pval( 0)) ; - if(!_safe && _dsqr > _RTOL) + real_type _dtol = + (real_type) 1./3. * _fbal[3] ; + + if(!_safe && _dsqr > _dtol) /*--------------------------- prune near-degeneracies */ continue ; @@ -995,7 +966,7 @@ } -# endif //__RDEL_PRED_BASE_2__ +# endif //__RDEL_PRED_BASE_3__ diff --git a/src/libcpp/rdel_mesh/rdel_make_2.hpp b/src/libcpp/rdel_mesh/rdel_make_2.hpp index 663fef6..9dc356c 100644 --- a/src/libcpp/rdel_mesh/rdel_make_2.hpp +++ b/src/libcpp/rdel_mesh/rdel_make_2.hpp @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 17 January, 2019 + * Last updated: 24 January, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -108,8 +108,6 @@ typedef typename mesh_type::node_data node_data ; - typedef typename - mesh_type::ball_data ball_data ; typedef typename mesh_type::edge_data edge_data ; typedef typename @@ -130,49 +128,6 @@ real_type , iptr_type > rdel_stat ; - - /* - -------------------------------------------------------- - * INIT-BALL: add new ball to restricted-tria. - -------------------------------------------------------- - */ - - __static_call - __normal_call void_type init_ball ( - mesh_type &_mesh, - geom_type &_geom, - iptr_type _npos, - char_type _kind, - iptr_type &_nbal, - rdel_opts &_opts - ) - { - __unreferenced ( _geom ) ; - __unreferenced ( _opts ) ; - - if (_mesh._tria. - node(_npos)->feat()==hard_feat) - { - /*---------- push protecting ball for "hard" features */ - ball_data _ball ; - _ball._node[0] = _npos; - - _ball._pass = +0 ; - _ball._kind = _kind; - - _ball._ball[0] = _mesh. - _tria.node(_npos)->pval(0) ; - _ball._ball[1] = _mesh. - _tria.node(_npos)->pval(1) ; - - _ball._ball[2] = (real_type)0. ; - - _nbal += +1; - - _mesh.push_ball (_ball) ; - } - } - /* -------------------------------------------------------- * PUSH-EDGE: add new edge to restricted-tria. @@ -711,7 +666,6 @@ _tcpu.time_span(_ttic,_ttoc) ; # endif//__use_timers - iptr_type _nbal = +0 ; iptr_type _nedg = +0 ; iptr_type _ntri = +0 ; @@ -730,15 +684,12 @@ typename mesh_type::tria_pred(), +.8, _mesh._tset.get_alloc()) ; - iptr_list _tnew, _nnew ; + /*------------------------- DT cells to check for rDT */ + iptr_list _tnew ; _tnew.set_alloc ( _mesh._tria._tset.count()) ; - _nnew.set_alloc ( - _mesh._tria._nset.count()) ; - - /*------------------------- face in DT for rDT checks */ + iptr_type _tpos = +0 ; - iptr_type _npos = +0 ; for (auto _iter = _mesh._tria._tset.head() ; @@ -751,18 +702,6 @@ _tnew. push_tail( _tpos) ; } } - - for (auto _iter = - _mesh._tria._nset.head() ; - _iter != - _mesh._tria._nset.tend() ; - ++_iter, ++_npos) - { - if (_iter->mark() >= +0) - { - _nnew. push_tail( _npos) ; - } - } /*------------------------- calc. voronoi-dual points */ for( auto _iter = _tnew.head(); @@ -772,34 +711,6 @@ tria_circ(_mesh,*_iter) ; } - /*------------------------- test for restricted balls */ - if (_args.dims() >= 0 ) - { - - # ifdef __use_timers - _ttic = _time.now() ; - # endif//__use_timers - - for( auto _iter = _nnew.head(); - _iter != _nnew.tend(); - ++_iter ) - { - char_type _kind = feat_ball; - - init_ball(_mesh, _geom, - *_iter, - _kind, _nbal, - _args) ; - } - - # ifdef __use_timers - _ttoc = _time.now() ; - _tcpu._node_init += - _tcpu.time_span(_ttic,_ttoc) ; - # endif//__use_timers - - } - /*------------------------- test for restricted edges */ if (_args.dims() >= 1 && _geom.have_feat(1) ) @@ -896,10 +807,6 @@ std::to_string (_tcpu._mesh_seed)) ; _dump.push("\n") ; - _dump.push(" NODE-INIT = ") ; - _dump.push( - std::to_string (_tcpu._node_init)) ; - _dump.push("\n") ; _dump.push(" EDGE-INIT = ") ; _dump.push( std::to_string (_tcpu._edge_init)) ; @@ -909,10 +816,6 @@ std::to_string (_tcpu._tria_init)) ; _dump.push("\n") ; _dump.push("\n") ; - - _dump.push(" |rDEL-0| (node) = ") ; - _dump.push(std::to_string (_nbal)) ; - _dump.push("\n") ; _dump.push(" |rDEL-1| (edge) = ") ; _dump.push(std::to_string (_nedg)) ; diff --git a/src/libcpp/rdel_mesh/rdel_make_3.hpp b/src/libcpp/rdel_mesh/rdel_make_3.hpp index 8069f2b..d1329ec 100644 --- a/src/libcpp/rdel_mesh/rdel_make_3.hpp +++ b/src/libcpp/rdel_mesh/rdel_make_3.hpp @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 17 January, 2019 + * Last updated: 24 January, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -108,8 +108,6 @@ typedef typename mesh_type::node_data node_data ; - typedef typename - mesh_type::ball_data ball_data ; typedef typename mesh_type::edge_data edge_data ; typedef typename @@ -131,52 +129,7 @@ typedef mesh::rdel_timers < real_type , iptr_type > rdel_stat ; - - - /* - -------------------------------------------------------- - * INIT-BALL: add new ball to restricted-tria. - -------------------------------------------------------- - */ - - __static_call - __normal_call void_type init_ball ( - mesh_type &_mesh, - geom_type &_geom, - iptr_type _npos, - char_type _kind, - iptr_type &_nbal, - rdel_opts &_opts - ) - { - __unreferenced ( _geom ) ; - __unreferenced ( _opts ) ; - - if (_mesh._tria. - node(_npos)->feat()==hard_feat) - { - /*---------- push protecting ball for "hard" features */ - ball_data _ball ; - _ball._node[0] = _npos; - - _ball._pass = +0 ; - _ball._kind = _kind; - - _ball._ball[0] = _mesh. - _tria.node(_npos)->pval(0) ; - _ball._ball[1] = _mesh. - _tria.node(_npos)->pval(1) ; - _ball._ball[2] = _mesh. - _tria.node(_npos)->pval(2) ; - - _ball._ball[3] = (real_type)0. ; - - _nbal += +1; - - _mesh.push_ball (_ball) ; - } - } - + /* -------------------------------------------------------- * PUSH-EDGE: add new edge to restricted-tria. @@ -842,7 +795,6 @@ _tcpu.time_span(_ttic,_ttoc) ; # endif//__use_timers - iptr_type _nbal = +0 ; iptr_type _nedg = +0 ; iptr_type _nfac = +0 ; iptr_type _ntri = +0 ; @@ -868,15 +820,12 @@ typename mesh_type::tria_pred(), +.8, _mesh._tset.get_alloc()) ; - iptr_list _tnew, _nnew ; + /*------------------------- DT cells to check for rDT */ + iptr_list _tnew ; _tnew.set_alloc ( _mesh._tria._tset.count()) ; - _nnew.set_alloc ( - _mesh._tria._nset.count()) ; - - /*------------------------- face in DT for rDT checks */ + iptr_type _tpos = +0 ; - iptr_type _npos = +0 ; for (auto _iter = _mesh._tria._tset.head() ; @@ -890,18 +839,6 @@ } } - for (auto _iter = - _mesh._tria._nset.head() ; - _iter != - _mesh._tria._nset.tend() ; - ++_iter, ++_npos) - { - if (_iter->mark() >= +0) - { - _nnew. push_tail( _npos) ; - } - } - if (_geom.have_feat(1) || _geom.have_feat(2) ) { @@ -913,34 +850,6 @@ tria_circ(_mesh,*_iter) ; } } - - /*------------------------- test for restricted balls */ - if (_args.dims() >= 0 ) - { - - # ifdef __use_timers - _ttic = _time.now() ; - # endif//__use_timers - - for( auto _iter = _nnew.head(); - _iter != _nnew.tend(); - ++_iter ) - { - char_type _kind = feat_ball; - - init_ball(_mesh, _geom, - *_iter, - _kind, _nbal, - _args) ; - } - - # ifdef __use_timers - _ttoc = _time.now() ; - _tcpu._node_init += - _tcpu.time_span(_ttic,_ttoc) ; - # endif//__use_timers - - } /*------------------------- test for restricted edges */ if (_args.dims() >= 1 && @@ -1069,10 +978,6 @@ std::to_string (_tcpu._mesh_seed)) ; _dump.push("\n") ; - _dump.push(" NODE-INIT = ") ; - _dump.push( - std::to_string (_tcpu._node_init)) ; - _dump.push("\n") ; _dump.push(" EDGE-INIT = ") ; _dump.push( std::to_string (_tcpu._edge_init)) ; @@ -1087,10 +992,6 @@ _dump.push("\n") ; _dump.push("\n") ; - _dump.push(" |rDEL-0| (node) = ") ; - _dump.push(std::to_string (_nbal)) ; - _dump.push("\n") ; - _dump.push(" |rDEL-1| (edge) = ") ; _dump.push(std::to_string (_nedg)) ; _dump.push("\n") ; diff --git a/src/libcpp/rdel_mesh/rdel_mesh_2.hpp b/src/libcpp/rdel_mesh/rdel_mesh_2.hpp index ce85c70..8f96eee 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: 17 January, 2019 + * Last updated: 24 January, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -270,6 +270,7 @@ */ #include "rdel_update_face_2.inc" + #include "rdel_update_topo_2.inc" /* @@ -835,11 +836,13 @@ iptr_list _nnew, _nold ; iptr_list _tnew, _told ; + iptr_list _emrk ; + escr_list _escr ; tscr_list _tscr ; ball_list _bscr ; - edat_list _edat, _etmp ; + edat_list _edat, _eprv ; tdat_list _tdat ; ball_list _bdat ; @@ -853,19 +856,14 @@ _mesh._eset._lptr. set_count ( _mesh._tria._tset.count()*+3 , containers::loose_alloc, nullptr); + _mesh._tset._lptr. set_count ( _mesh._tria._tset.count()*+1 , containers::loose_alloc, nullptr); /*------------------------------ init. topo hash obj. */ typename - mesh_type::node_list _etin ( - typename mesh_type::node_hash(), - typename mesh_type::node_pred(), - +.8,_mesh._nset.get_alloc()) ; - - typename - mesh_type::edge_list _pedg ( + mesh_type::edge_list _epro ( typename mesh_type::edge_hash(), typename mesh_type::edge_pred(), +.8,_mesh._eset.get_alloc()) ; @@ -895,11 +893,11 @@ _tcpu.time_span(_ttic,_ttoc) ; # endif//__use_timers - /*-------------------- calc. size-func for seed nodes */ + /*------------------------------ calc. hfun. at seeds */ for (auto _node = - _mesh._tria._nset.head(); + _mesh._tria._nset.head() ; _node != - _mesh._tria._nset.tend(); + _mesh._tria._nset.tend() ; ++_node ) { if (_node->mark() >= +0) @@ -910,7 +908,7 @@ } /*-------------------- main: refine edges/faces/trias */ - iptr_type _pass = +0 ; + iptr_type _pass = +0 ; for(bool_type _done=false; !_done ; ) { @@ -924,22 +922,25 @@ if(++_pass>_args.iter()) break; /*------------------------- init. array workspace */ - _nnew.set_count( +0 ) ; + + _nnew.set_count( +0 ) ; // del-tri idx lists _nold.set_count( +0 ) ; _tnew.set_count( +0 ) ; _told.set_count( +0 ) ; - _escr.set_count( +0 ) ; + _escr.set_count( +0 ) ; // face "cost" lists _tscr.set_count( +0 ) ; _bscr.set_count( +0 ) ; - _etmp.set_count( +0 ) ; - _edat.set_count( +0 ) ; - _tdat.set_count( +0 ) ; - _bdat.set_count( +0 ) ; + _eprv.set_count( +0 ) ; // "old" edge in cav + _edat.set_count( +0 ) ; // "new" edge in cav + _tdat.set_count( +0 ) ; // "new" cell in cav + _bdat.set_count( +0 ) ; // "new" ball in cav /*--------- calc. "restricted-ness" incrementally */ + bool_type _irDT = false ; + if (_mode == null_mode ) { /*------------------------- init. protecting ball */ @@ -949,6 +950,8 @@ _mode = node_mode; + _irDT = true; // init. new face in rDT + init_rdel( _geom, _hfun, _mesh, false, _nnew, _tnew, @@ -976,6 +979,12 @@ _mode = edge_mode; + _irDT = true; // init. new face in rDT + + init_ball( _geom, _hfun, + _mesh, _epro, _pass, + _mode, _args) ; + init_rdel( _geom, _hfun, _mesh, true,// init. circum. for rDT _nnew, _tnew, @@ -996,6 +1005,7 @@ _edat. empty() ) { /*------------------------- init. restricted topo */ + _mode = etop_mode ; } @@ -1011,6 +1021,8 @@ _mode = tria_mode; + _irDT = true; // init. new face in rDT + init_rdel( _geom, _hfun, _mesh, false, _nnew, _tnew, @@ -1028,12 +1040,7 @@ /*------------- refine "bad" sub-faces until done */ - if ( _bscr.empty() && - _bdat.empty() && - _escr.empty() && - _edat.empty() && - _tscr.empty() && - _tdat.empty() ) + if (_irDT == false ) { char_type _tdim = -1; @@ -1049,9 +1056,9 @@ _kind =_bad_ball( _geom, _hfun, _mesh, _mode, - _pedg, _nnew, _nold, + _nnew, _nold, _tnew, _told, _nbpq, - _etmp, _edat, _escr, + _eprv, _edat, _escr, _tdat, _tscr, _bdat, _bscr, _tdim, _pass, _args) ; @@ -1072,9 +1079,9 @@ _kind =_bad_edge( _geom, _hfun, _mesh, _mode, - _pedg, _nnew, _nold, + _epro, _nnew, _nold, _tnew, _told, _eepq, - _etmp, _edat, _escr, + _eprv, _edat, _escr, _tdat, _tscr, _bdat, _bscr, _tdim, _pass, _args) ; @@ -1095,10 +1102,10 @@ _kind =_bad_etop( _geom, _hfun, _mesh, _mode, - _pedg, _nnew, _nold, + _epro, _nnew, _nold, _tnew, _told, - _etpq, _etin, - _etmp, _edat, _escr, + _etpq, _emrk, + _eprv, _edat, _escr, _tdat, _tscr, _bdat, _bscr, _tdim, _pass, _args) ; @@ -1119,9 +1126,9 @@ _kind =_bad_tria( _geom, _hfun, _mesh, _mode, - _pedg, _nnew, _nold, + _epro, _nnew, _nold, _tnew, _told, _ttpq, - _etmp, _edat, _escr, + _eprv, _edat, _escr, _tdat, _tscr, _bdat, _bscr, _tdim, _pass, _args) ; @@ -1135,27 +1142,31 @@ /*----------------------------- meshing converged */ else { _done = true ; } - if (_pass%_jlog_freq==+0 || _done ) + if (_pass%_jlog_freq==+0 || _done) { /*----------------------------- output to logfile */ std::stringstream _sstr; - _sstr << std::setw(11) << - _pass << std::setw(13) << - _mesh._eset.count() - << std::setw(13) << - _mesh._tset.count() - << "\n" ; + _sstr << std::setw(+11) << + _pass << std::setw(+13) << + _mesh._eset.count () + << std::setw(+13) << + _mesh._tset.count () + << "\n" ; _dump.push(_sstr.str()); } - if (_kind != rdel_opts::null_kind ) + if (_kind != rdel_opts::null_kind) { - /*--------------- update point-placement counters */ + /*----------------------------- inc. steiner info */ if (_tdim == +1) + { _enod[_kind] += +1 ; + } else if (_tdim == +2) + { _tnod[_kind] += +1 ; + } } } @@ -1163,64 +1174,34 @@ if (_pass%_trim_freq == +0 ) { /*--------------- trim null PQ items "on-the-fly" */ - //trim_nbpq(_mesh, _nbpq); - trim_eepq(_mesh, _eepq); - trim_ttpq(_mesh, _ttpq); - //trim_etpq(_mesh, _etpq); + //trim_nbpq( _mesh , + // _nbpq ) ; + trim_eepq( _mesh , + _eepq ) ; + trim_ttpq( _mesh , + _ttpq ) ; + //trim_etpq( _mesh , + // _etpq ) ; - trim_list(_nnew) ; - trim_list(_nold) ; - trim_list(_tnew) ; - trim_list(_told) ; + trim_list( _nnew ) ; + trim_list( _nold ) ; + trim_list( _tnew ) ; + trim_list( _told ) ; - trim_list(_etmp) ; - trim_list(_edat) ; - trim_list(_escr) ; - trim_list(_tdat) ; - trim_list(_tscr) ; - trim_list(_bscr) ; - trim_list(_bdat) ; + trim_list( _eprv ) ; + trim_list( _edat ) ; + trim_list( _escr ) ; + trim_list( _tdat ) ; + trim_list( _tscr ) ; + trim_list( _bscr ) ; + trim_list( _bdat ) ; } - - /*----------------------------- enqueue edge topo */ + + /*--------------- enqueue nodes for topol. checks */ - for (auto _edge = _edat.head(); - _edge != _edat.tend(); - ++_edge ) - { - /*------------------------- bypass "old" edge */ - typename mesh_type:: - edge_list:: - item_type *_eptr=nullptr; - if (_mesh. - find_edge(*_edge, _eptr)) - continue ; - - /*------------------------- enqueue edge node */ - iptr_type _npos; - for (_npos = +2; _npos-- != +0; ) - { - iptr_type _node = - _edge->_node[_npos] ; - - node_data _ndat; - _ndat._pass = _pass; - _ndat._node[0] = _node; - - typename mesh_type:: - node_list:: - item_type*_nptr=nullptr; - if(!_etin.find(_ndat, _nptr)) - { - if (_args.top1()) - { - /*--------------------- enqueue edge node */ - _etin.push(_ndat) ; - _etpq.push(_ndat) ; - } - } - } - } + fill_topo( _mesh, _pass, + _etpq, _emrk, + _edat, _eprv, _args) ; /*--------------- update restricted triangulation */ diff --git a/src/libcpp/rdel_mesh/rdel_mesh_3.hpp b/src/libcpp/rdel_mesh/rdel_mesh_3.hpp index 8765116..b4e7043 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: 17 January, 2019 + * Last updated: 24 January, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -303,6 +303,7 @@ */ #include "rdel_update_face_3.inc" + #include "rdel_update_topo_3.inc" /* @@ -958,13 +959,15 @@ iptr_list _nnew, _nold ; iptr_list _tnew, _told ; + iptr_list _emrk, _fmrk ; + escr_list _escr ; fscr_list _fscr ; tscr_list _tscr ; ball_list _bscr ; - edat_list _edat, _etmp ; - fdat_list _fdat, _ftmp ; + edat_list _edat, _eprv ; + fdat_list _fdat, _fprv ; tdat_list _tdat ; ball_list _bdat ; @@ -980,34 +983,24 @@ _mesh._eset._lptr. set_count ( _mesh._tria._tset.count()*+6 , containers::loose_alloc, nullptr) ; + _mesh._fset._lptr. set_count ( _mesh._tria._tset.count()*+4 , - containers::loose_alloc, nullptr) ; + containers::loose_alloc, nullptr) ; + _mesh._tset._lptr. set_count ( _mesh._tria._tset.count()*+1 , containers::loose_alloc, nullptr) ; /*------------------------------ init. topo hash obj. */ typename - mesh_type::node_list _etin ( - typename mesh_type::node_hash(), - typename mesh_type::node_pred(), - +.8,_mesh._nset.get_alloc()) ; - - typename - mesh_type::node_list _ftin ( - typename mesh_type::node_hash(), - typename mesh_type::node_pred(), - +.8,_mesh._nset.get_alloc()) ; - - typename - mesh_type::edge_list _pedg ( + mesh_type::edge_list _epro ( typename mesh_type::edge_hash(), typename mesh_type::edge_pred(), +.8,_mesh._eset.get_alloc()) ; typename - mesh_type::face_list _pfac ( + mesh_type::face_list _fpro ( typename mesh_type::face_hash(), typename mesh_type::face_pred(), +.8,_mesh._fset.get_alloc()) ; @@ -1050,7 +1043,7 @@ { if (_node->mark() >= +0) { - _node->idxh() = + _node->idxh() = hfun_type::null_hint () ; } } @@ -1070,25 +1063,28 @@ if(++_pass>_args.iter()) break; /*------------------------- init. array workspace */ - _nnew.set_count( +0 ) ; - _nold.set_count( +0 ) ; + + _nnew.set_count( +0 ) ; // del-tri idx lists + _nold.set_count( +0 ) ; _tnew.set_count( +0 ) ; _told.set_count( +0 ) ; - _escr.set_count( +0 ) ; + _escr.set_count( +0 ) ; // face "cost" lists _fscr.set_count( +0 ) ; _tscr.set_count( +0 ) ; _bscr.set_count( +0 ) ; - _etmp.set_count( +0 ) ; - _edat.set_count( +0 ) ; - _ftmp.set_count( +0 ) ; - _fdat.set_count( +0 ) ; - _tdat.set_count( +0 ) ; - _bdat.set_count( +0 ) ; + _eprv.set_count( +0 ) ; // "old" edge in cav + _edat.set_count( +0 ) ; // "new" edge in cav + _fprv.set_count( +0 ) ; // "old" face in cav + _fdat.set_count( +0 ) ; // "new" face in cav + _tdat.set_count( +0 ) ; // "new" cell in cav + _bdat.set_count( +0 ) ; // "new" ball in cav /*--------- calc. "restricted-ness" incrementally */ + bool_type _irDT = false ; + if (_mode == null_mode ) { /*------------------------- init. protecting ball */ @@ -1098,6 +1094,8 @@ _mode = node_mode; + _irDT = true; // init. new face in rDT + init_rdel( _geom, _hfun, _mesh, false, _nnew, _tnew, @@ -1126,6 +1124,13 @@ _mode = edge_mode; + _irDT = true; // init. new face in rDT + + init_ball( _geom, _hfun, + _mesh, + _epro, _fpro, _pass, + _mode, _args) ; + init_rdel( _geom, _hfun, _mesh, true,// init. circum. for rDT _nnew, _tnew, @@ -1162,6 +1167,8 @@ _mode = face_mode; + _irDT = true; // init. new face in rDT + init_rdel( _geom, _hfun, _mesh, false, _nnew, _tnew, @@ -1199,6 +1206,8 @@ _mode = tria_mode; + _irDT = true; // init. new face in rDT + init_rdel( _geom, _hfun, _mesh, false, _nnew, _tnew, @@ -1217,14 +1226,7 @@ /*------------- refine "bad" sub-faces until done */ - if ( _bscr.empty() && - _bdat.empty() && - _escr.empty() && - _edat.empty() && - _fscr.empty() && - _fdat.empty() && - _tscr.empty() && - _tdat.empty() ) + if (_irDT == false ) { char_type _tdim = -1; @@ -1240,11 +1242,10 @@ _kind =_bad_ball( _geom, _hfun, _mesh, _mode, - _pedg, _pfac, _nnew, _nold, _tnew, _told, _nbpq, - _etmp, _edat, _escr, - _ftmp, _fdat, _fscr, + _eprv, _edat, _escr, + _fprv, _fdat, _fscr, _tdat, _tscr, _bdat, _bscr, _tdim, _pass, _args) ; @@ -1265,11 +1266,11 @@ _kind =_bad_edge( _geom, _hfun, _mesh, _mode, - _pedg, _pfac, + _epro, _fpro, _nnew, _nold, _tnew, _told, _eepq, - _etmp, _edat, _escr, - _ftmp, _fdat, _fscr, + _eprv, _edat, _escr, + _fprv, _fdat, _fscr, _tdat, _tscr, _bdat, _bscr, _tdim, _pass, _args) ; @@ -1290,12 +1291,12 @@ _kind =_bad_etop( _geom, _hfun, _mesh, _mode, - _pedg, _pfac, + _epro, _fpro, _nnew, _nold, _tnew, _told, - _etpq, _etin, - _etmp, _edat, _escr, - _ftmp, _fdat, _fscr, + _etpq, _emrk, + _eprv, _edat, _escr, + _fprv, _fdat, _fscr, _tdat, _tscr, _bdat, _bscr, _tdim, _pass, _args) ; @@ -1316,11 +1317,11 @@ _kind =_bad_face( _geom, _hfun, _mesh, _mode, - _pedg, _pfac, + _epro, _fpro, _nnew, _nold, _tnew, _told, _ffpq, - _etmp, _edat, _escr, - _ftmp, _fdat, _fscr, + _eprv, _edat, _escr, + _fprv, _fdat, _fscr, _tdat, _tscr, _bdat, _bscr, _tdim, _pass, _args) ; @@ -1341,12 +1342,12 @@ _kind =_bad_ftop( _geom, _hfun, _mesh, _mode, - _pedg, _pfac, + _epro, _fpro, _nnew, _nold, _tnew, _told, - _ftpq, _ftin, - _etmp, _edat, _escr, - _ftmp, _fdat, _fscr, + _ftpq, _fmrk, + _eprv, _edat, _escr, + _fprv, _fdat, _fscr, _tdat, _tscr, _bdat, _bscr, _tdim, _pass, _args) ; @@ -1367,11 +1368,11 @@ _kind =_bad_tria( _geom, _hfun, _mesh, _mode, - _pedg, _pfac, + _epro, _fpro, _nnew, _nold, _tnew, _told, _ttpq, - _etmp, _edat, _escr, - _ftmp, _fdat, _fscr, + _eprv, _edat, _escr, + _fprv, _fdat, _fscr, _tdat, _tscr, _bdat, _bscr, _tdim, _pass, _args) ; @@ -1385,32 +1386,38 @@ /*----------------------------- meshing converged */ else { _done = true ; } - if (_pass%_jlog_freq==+0 || _done ) + if (_pass%_jlog_freq==+0 || _done) { /*----------------------------- output to logfile */ - std::stringstream _sstr ; - _sstr << std::setw(11) << - _pass << std::setw(13) << - _mesh._eset.count() - << std::setw(13) << - _mesh._fset.count() - << std::setw(13) << - _mesh._tset.count() - << "\n" ; - _dump.push(_sstr.str()) ; + std::stringstream _sstr; + _sstr << std::setw(+11) << + _pass << std::setw(+13) << + _mesh._eset.count () + << std::setw(+13) << + _mesh._fset.count () + << std::setw(+13) << + _mesh._tset.count () + << "\n" ; + _dump.push(_sstr.str()); } - if (_kind != rdel_opts::null_kind ) + if (_kind != rdel_opts::null_kind) { - /*--------------- update point-placement counters */ + /*----------------------------- inc. steiner info */ if (_tdim == +1) + { _enod[_kind] += +1 ; + } else if (_tdim == +2) + { _fnod[_kind] += +1 ; + } else if (_tdim == +3) + { _tnod[_kind] += +1 ; + } } } @@ -1418,109 +1425,43 @@ if (_pass%_trim_freq == +0 ) { /*--------------- trim null PQ items "on-the-fly" */ - //trim_nbpq(_mesh, _nbpq); - trim_eepq(_mesh, _eepq); - trim_ffpq(_mesh, _ffpq); - trim_ttpq(_mesh, _ttpq); - //trim_etpq(_mesh, _etpq); - //trim_ftpq(_mesh, _ftpq); + //trim_nbpq( _mesh , + // _nbpq ) ; + trim_eepq( _mesh , + _eepq ) ; + trim_ffpq( _mesh , + _ffpq ) ; + trim_ttpq( _mesh , + _ttpq ) ; + //trim_etpq( _mesh , + // _etpq ) ; + //trim_ftpq( _mesh , + // _ftpq ) ; - trim_list(_nnew) ; - trim_list(_nold) ; - trim_list(_tnew) ; - trim_list(_told) ; + trim_list( _nnew ) ; + trim_list( _nold ) ; + trim_list( _tnew ) ; + trim_list( _told ) ; - trim_list(_etmp) ; - trim_list(_edat) ; - trim_list(_escr) ; - trim_list(_ftmp) ; - trim_list(_fdat) ; - trim_list(_fscr) ; - trim_list(_tdat) ; - trim_list(_tscr) ; - trim_list(_bdat) ; - trim_list(_bscr) ; + trim_list( _eprv ) ; + trim_list( _edat ) ; + trim_list( _escr ) ; + trim_list( _fprv ) ; + trim_list( _fdat ) ; + trim_list( _fscr ) ; + trim_list( _tdat ) ; + trim_list( _tscr ) ; + trim_list( _bdat ) ; + trim_list( _bscr ) ; } - /*----------------------------- enqueue edge topo */ - - for (auto _edge = _edat.head(); - _edge != _edat.tend(); - ++_edge ) - { - /*------------------------- bypass "old" edge */ - typename mesh_type:: - edge_list:: - item_type *_eptr=nullptr; - if (_mesh. - find_edge(*_edge, _eptr)) - continue ; - - /*------------------------- enqueue edge node */ - iptr_type _npos; - for (_npos = +2; _npos-- != +0; ) - { - iptr_type _node = - _edge->_node[_npos]; - - node_data _ndat; - _ndat._pass = _pass ; - _ndat._node[0] = _node ; - - typename mesh_type:: - node_list:: - item_type*_nptr=nullptr; - if(!_etin.find(_ndat, _nptr)) - { - if (_args.top1()) - { - /*--------------------- enqueue edge node */ - _etin.push(_ndat) ; - _etpq.push(_ndat) ; - } - } - } - } - - /*----------------------------- enqueue face topo */ + /*--------------- enqueue nodes for topol. checks */ - for (auto _face = _fdat.head(); - _face != _fdat.tend(); - ++_face ) - { - /*------------------------- bypass "old" face */ - typename mesh_type:: - face_list:: - item_type *_fptr=nullptr; - if (_mesh. - find_face(*_face, _fptr)) - continue ; - - /*------------------------- enqueue face node */ - iptr_type _npos; - for (_npos = +3; _npos-- != +0; ) - { - iptr_type _node = - _face->_node[_npos]; - - node_data _ndat; - _ndat._pass = _pass ; - _ndat._node[0] = _node ; - - typename mesh_type:: - node_list:: - item_type*_nptr=nullptr; - if(!_ftin.find(_ndat, _nptr)) - { - if (_args.top2()) - { - /*--------------------- enqueue face node */ - _ftin.push(_ndat) ; - _ftpq.push(_ndat) ; - } - } - } - } + fill_topo( _mesh, _pass, + _etpq, _emrk, + _ftpq, _fmrk, + _edat, _eprv, + _fdat, _fprv, _args) ; /*--------------- update restricted triangulation */ diff --git a/src/libcpp/rdel_mesh/rdel_offh_delfront_2.inc b/src/libcpp/rdel_mesh/rdel_offh_delfront_2.inc index 3442e8c..0671dbe 100644 --- a/src/libcpp/rdel_mesh/rdel_offh_delfront_2.inc +++ b/src/libcpp/rdel_mesh/rdel_offh_delfront_2.inc @@ -31,9 +31,9 @@ * -------------------------------------------------------- * - * Last updated: 25 August, 2018 + * Last updated: 22 January, 2019 * - * Copyright 2013-2018 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -133,7 +133,7 @@ _pert += + 1 ; _pert *= _seps ; - for(iptr_type _iter = +0; _iter++ != +8 ; ) + for(iptr_type _iter = +0; _iter++ != +4 ; ) { _nsiz[ 2] = _hfun.eval (_ppos, _hint) ; @@ -160,14 +160,15 @@ _dist -= _pert * _dist ; + iptr_type _dual = -1 ; cosine_intersect _pred(_ipos,_dvec); _ball._rrad = _dist ; if (!_geom.intersect ( _ball, _pred) ) - return rdel_opts::null_kind; + return rdel_opts::null_kind ; if (!_pred. _find ) return rdel_opts::null_kind ; - + real_type _move = geometry::length_2d( _ppos,&_pred._proj.pval(0)); @@ -175,6 +176,14 @@ _ppos[ 0] = _pred._proj.pval(0); _ppos[ 1] = _pred._proj.pval(1); + _mesh._tria.find_node ( // test voro. limiter + &_ppos[0], _dual, _enod[ 0]) ; + + if (_dual != _enod[0] ) + { + return rdel_opts::null_kind ; + } + if (_move < _rtol * _hval) break ; } @@ -248,7 +257,7 @@ _pert += + 1 ; _pert *= _seps ; - for(iptr_type _iter = +0; _iter++ != +8 ; ) + for(iptr_type _iter = +0; _iter++ != +4 ; ) { _kind = rdel_opts::offh_kind ; diff --git a/src/libcpp/rdel_mesh/rdel_offh_delfront_3.inc b/src/libcpp/rdel_mesh/rdel_offh_delfront_3.inc index f186fb1..a79d96c 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: 09 January, 2019 + * Last updated: 22 January, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -141,7 +141,7 @@ _pert += + 1 ; _pert *= _seps ; - for(iptr_type _iter = +0; _iter++ != +8 ; ) + for(iptr_type _iter = +0; _iter++ != +4 ; ) { _nsiz[2] = _hfun.eval(_ppos, _hint) ; @@ -169,13 +169,14 @@ _dist -= _pert * _dist ; + iptr_type _dual = -1 ; cosine_intersect _pred(_ipos,_dvec); _ball._rrad = _dist ; if (!_geom.intersect ( _ball, _pred) ) - return rdel_opts::null_kind; + return rdel_opts::null_kind ; if (!_pred. _find ) - return rdel_opts::null_kind; + return rdel_opts::null_kind ; real_type _move = geometry::length_3d( @@ -184,8 +185,16 @@ _ppos[ 0] = _pred._proj.pval(0); _ppos[ 1] = _pred._proj.pval(1); _ppos[ 2] = _pred._proj.pval(2); + + _mesh._tria.find_node ( // test voro. limiter + &_ppos[0], _dual, _enod[ 0]) ; + + if (_dual != _enod[0] ) + { + return rdel_opts::null_kind ; + } - if (_move < _rtol * _hval) break; + if (_move < _rtol * _hval) break ; } return (_kind ) ; @@ -277,7 +286,7 @@ _pert += + 1 ; _pert *= _seps ; - for(iptr_type _iter = +0; _iter++ != +8 ; ) + for(iptr_type _iter = +0; _iter++ != +4 ; ) { _nsiz[2] = _hfun.eval(_ppos, _hint) ; @@ -338,14 +347,15 @@ _dist -= _pert * _dist ; + iptr_type _dual = -1 ; cosine_intersect _pred(_pmid,_dvec); _disc._rrad = _dist ; if (!_geom.intersect ( _disc, _sbal, _pred) ) - return rdel_opts::null_kind; + return rdel_opts::null_kind ; if (!_pred. _find ) - return rdel_opts::null_kind; + return rdel_opts::null_kind ; real_type _move = geometry::length_3d( @@ -354,8 +364,17 @@ _ppos[ 0] = _pred._proj.pval(0); _ppos[ 1] = _pred._proj.pval(1); _ppos[ 2] = _pred._proj.pval(2); + + _mesh._tria.find_node ( // test voro. limiter + &_ppos[0], _dual, _enod[ 0]) ; + + if (_dual != _enod[0] && + _dual != _enod[1] ) + { + return rdel_opts::null_kind ; + } - if (_move < _rtol * _hmin) break; + if (_move < _rtol * _hmin) break ; } return (_kind ) ; @@ -436,7 +455,7 @@ _pert += + 1 ; _pert *= _seps ; - for(iptr_type _iter = +0; _iter++ != +8 ; ) + for(iptr_type _iter = +0; _iter++ != +4 ; ) { _nsiz[3] = _hfun.eval(_ppos, _hint) ; diff --git a/src/libcpp/rdel_mesh/rdel_refine_ball_2.inc b/src/libcpp/rdel_mesh/rdel_refine_ball_2.inc index 260b8e3..255e137 100644 --- a/src/libcpp/rdel_mesh/rdel_refine_ball_2.inc +++ b/src/libcpp/rdel_mesh/rdel_refine_ball_2.inc @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 09 January, 2019 + * Last updated: 25 January, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -57,8 +57,6 @@ hfun_type &_hfun , mesh_type &_mesh , mode_type _mode , - typename - mesh_type::edge_list &_pedg , iptr_list &_nnew , iptr_list &_nold , iptr_list &_tnew , @@ -77,43 +75,51 @@ ) { class node_pred - { - /*----------------- find adj. set of tria-to-node */ - public : - typedef typename - mesh_type:: - tria_type tria_type ; - - public : - iptr_type _npos; + { + /*--------------------- find adj. set of tria-to-node */ + public : + + typedef typename + mesh_type::tria_type tria_type ; - public : - __inline_call node_pred ( - iptr_type _nsrc - ) : _npos(_nsrc) {} - /*----------------- find adj. set of tria-to-node */ - __inline_call - bool_type operator()( - tria_type&_tria, - iptr_type _tpos, - iptr_type - ) - { return ( - _tria.tria(_tpos)->node(+0) - == this->_npos || - _tria.tria(_tpos)->node(+1) - == this->_npos || - _tria.tria(_tpos)->node(+2) - == this->_npos ); - } - } ; + public : + iptr_type _npos; + + public : + __inline_call node_pred ( + iptr_type _nsrc + ) : _npos(_nsrc) {} + /*--------------------- find adj. set of tria-to-node */ + __inline_call bool_type operator() ( + tria_type&_tria, + iptr_type _tpos, + iptr_type _fpos + ) + { + iptr_type _tnod[3] = { + _tria.tria(_tpos)->node(0) , + _tria.tria(_tpos)->node(1) , + _tria.tria(_tpos)->node(2) + } ; + + __unreferenced(_fpos) ; + + if (_tnod[0] == this->_npos) + return true ; + if (_tnod[1] == this->_npos) + return true ; + if (_tnod[2] == this->_npos) + return true ; + + return false ; + } + } ; typename rdel_opts::node_kind _kind = rdel_opts::null_kind ; __unreferenced( _hfun) ; - __unreferenced( _mode) ; - __unreferenced( _eold) ; + __unreferenced( _mode) ; __unreferenced( _tdim) ; __unreferenced( _pass) ; __unreferenced( _args) ; @@ -128,20 +134,22 @@ _ecav.set_count( +0 ) ; _tscr.set_count( +0 ) ; _tcav.set_count( +0 ) ; + _bcav.set_count( +0 ) ; _bscr.set_count( +0 ) ; /*--------------------- pop() leading element from PQ */ - bool_type _find = false; + typename mesh_type:: + ball_list:: + item_type*_bptr = nullptr; ball_data _ball ; + + bool_type _find = false; for ( ; !_nbpq.empty() ; ) { /*--------------------- cycles until valid ball found */ _nbpq._pop_root(_ball) ; - typename mesh_type:: - ball_list:: - item_type *_bptr = nullptr ; if (_mesh.find_ball(_ball,_bptr)) { if (_bptr->_data._pass == @@ -149,16 +157,16 @@ { _find = true; break; } - } + } } - if (!_find) return (_kind) ; - + if ( !_find ) return _kind ; + /*--------------------- check dist. to existing nodes */ real_type static const _HTOL = - (real_type) +.33 * .33 ; + (real_type) +.333 * .333 ; real_type static const _STOL = - (real_type) +.50 * .50 ; + (real_type) +.500 * .500 ; iptr_list _tset; _tset.set_alloc( +32) ; @@ -217,12 +225,9 @@ _ball._ball[2] = std::sqrt ( _ball._ball [ 2]) ; - iptr_type _hpos = - _mesh._tria. - node(_nadj)->idxh() ; iptr_type _ndeg = _mesh._tria. - node(_nadj)->topo() ; + node(_nadj)->topo() ; /*--------------------- find intersections with GEOM. */ typename @@ -252,6 +257,10 @@ _hdat._rrad = (real_type)+1.5 * _ball._ball[2] ; + + _bptr->_data._ball[2] = // save to mesh "proper" + (real_type)+1.0 * _ball._ball[2] + * _ball._ball[2] ; _okay = _okay & _geom.intersect ( _pdat, _pred) ; // actual! @@ -286,41 +295,105 @@ (real_type) +.67 ; continue ; } - else { break ; } + else break ; } - + + return ( _kind ) ; + } + + /* + -------------------------------------------------------- + * INIT-BALL: setup ball "collar". + -------------------------------------------------------- + */ + + __static_call + __normal_call void_type init_ball ( + geom_type &_geom , + hfun_type &_hfun , + mesh_type &_mesh , + typename + mesh_type::edge_list &_epro , + iptr_type _mode , + iptr_type _pass , + rdel_opts &_args + ) + { /*--------------------- push finalised "collar" nodes */ - for (auto _iter = _pred._list.head() ; - _iter != _pred._list.tend() ; - ++_iter ) + __unreferenced(_hfun) ; + __unreferenced(_mode) ; + __unreferenced(_pass) ; + __unreferenced(_args) ; + + for (auto _bpos = + _mesh._bset._lptr.head () ; + _bpos != + _mesh._bset._lptr.tend () ; + ++_bpos ) { - iptr_type _hint = - _mesh._tria. - node(_nadj)->next(); - - iptr_type _node = -1; - if (_mesh._tria.push_node( - &_iter->pval(0), _node, - _hint ) ) + if ( *_bpos == nullptr) continue ; + + for (auto _ball = *_bpos ; + _ball != nullptr; + _ball = _ball->_next ) { - _mesh._tria.node( - _node)->idxh() = _hpos ; - _mesh._tria.node( - _node)->fdim() = 1 ; - _mesh._tria.node( - _node)->feat() = 0 ; - _mesh._tria.node( - _node)->topo() = 2 ; - edge_data _edat; - _edat._node[0] = _nadj ; - _edat._node[1] = _node ; + typename + geom_type::ball_type _hdat ; + _hdat._pmid[0] = + _ball->_data._ball[0] ; + _hdat._pmid[1] = + _ball->_data._ball[1] ; + + _hdat._rrad = + std::sqrt(_ball->_data._ball[2]) ; + + iptr_type _nadj = + _ball->_data._node[0] ; + + mesh::keep_all_2d < + real_type , + iptr_type > _halo ; + + _geom. + intersect(_hdat, _halo) ; + + for (auto _iter = _halo._list.head() ; + _iter != _halo._list.tend() ; + ++_iter ) + { + iptr_type _hint = + _mesh._tria. + node(_nadj)->next() ; + + 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() = 1 ; + + _mesh._tria.node( + _node)->feat() = 0 ; + + _mesh._tria.node( + _node)->topo() = 2 ; + + edge_data _edat; + _edat._node[0] = _nadj ; + _edat._node[1] = _node ; - _pedg.push(_edat) ; + _epro.push(_edat) ; + } + } + } - } - - return ( _kind ) ; + } } diff --git a/src/libcpp/rdel_mesh/rdel_refine_ball_3.inc b/src/libcpp/rdel_mesh/rdel_refine_ball_3.inc index d7c2243..42a6a5a 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: 09 January, 2019 + * Last updated: 25 January, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -57,10 +57,6 @@ hfun_type &_hfun , mesh_type &_mesh , mode_type _mode , - typename - mesh_type::edge_list &_pedg , - typename - mesh_type::face_list &_pfac , iptr_list &_nnew , iptr_list &_nold , iptr_list &_tnew , @@ -82,45 +78,54 @@ ) { class node_pred - { - /*----------------- find adj. set of tria-to-node */ - public : - typedef typename - mesh_type:: - tria_type tria_type ; - - public : - iptr_type _npos; + { + /*--------------------- find adj. set of tria-to-node */ + public : + + typedef typename + mesh_type::tria_type tria_type ; - public : - __inline_call node_pred ( - iptr_type _nsrc - ) : _npos(_nsrc) {} - /*----------------- find adj. set of tria-to-node */ - __inline_call - bool_type operator()( - tria_type&_tria, - iptr_type _tpos, - iptr_type - ) - { return ( - _tria.tria(_tpos)->node(+0) - == this->_npos || - _tria.tria(_tpos)->node(+1) - == this->_npos || - _tria.tria(_tpos)->node(+2) - == this->_npos ); - } - } ; + public : + iptr_type _npos; + + public : + __inline_call node_pred ( + iptr_type _nsrc + ) : _npos(_nsrc) {} + /*--------------------- find adj. set of tria-to-node */ + __inline_call bool_type operator() ( + tria_type&_tria, + iptr_type _tpos, + iptr_type _fpos + ) + { + iptr_type _tnod[4] = { + _tria.tria(_tpos)->node(0) , + _tria.tria(_tpos)->node(1) , + _tria.tria(_tpos)->node(2) , + _tria.tria(_tpos)->node(3) + } ; + + __unreferenced(_fpos) ; + + if (_tnod[0] == this->_npos) + return true ; + if (_tnod[1] == this->_npos) + return true ; + if (_tnod[2] == this->_npos) + return true ; + if (_tnod[3] == this->_npos) + return true ; + + return false ; + } + } ; typename rdel_opts::node_kind _kind = rdel_opts::null_kind ; __unreferenced( _hfun) ; __unreferenced( _mode) ; - __unreferenced( _pfac) ; - __unreferenced( _eold) ; - __unreferenced( _fold) ; __unreferenced( _tdim) ; __unreferenced( _pass) ; __unreferenced( _args) ; @@ -130,8 +135,10 @@ _tnew.set_count( +0 ) ; _told.set_count( +0 ) ; + _eold.set_count( +0 ) ; _escr.set_count( +0 ) ; _ecav.set_count( +0 ) ; + _fold.set_count( +0 ) ; _fscr.set_count( +0 ) ; _fcav.set_count( +0 ) ; _tscr.set_count( +0 ) ; @@ -141,16 +148,17 @@ _bcav.set_count( +0 ) ; /*--------------------- pop() leading element from PQ */ - bool_type _find = false; + typename mesh_type:: + ball_list:: + item_type*_bptr = nullptr; ball_data _ball ; + + bool_type _find = false; for ( ; !_nbpq.empty() ; ) { /*--------------------- cycles until valid ball found */ _nbpq._pop_root(_ball) ; - typename mesh_type:: - ball_list:: - item_type *_bptr = nullptr ; if (_mesh.find_ball(_ball,_bptr)) { if (_bptr->_data._pass == @@ -160,14 +168,14 @@ } } } - if (!_find) return (_kind) ; + if ( !_find ) return _kind ; /*--------------------- check dist. to existing nodes */ real_type static const _HTOL = - (real_type) +.33 * .33 ; + (real_type) +.333 * .333 ; real_type static const _STOL = - (real_type) +.50 * .50 ; + (real_type) +.500 * .500 ; iptr_list _tset; _tset.set_alloc( +32) ; @@ -219,16 +227,13 @@ } } - /*--------------------- otherwise push "collar" nodes */ + /*--------------------- otherwise form "collar" radii */ _ball._ball[3] = std:: min ( _ball._ball[3],_rmax) ; _ball._ball[3] = std::sqrt ( _ball._ball [ 3]) ; - iptr_type _hpos = - _mesh._tria. - node(_nadj)->idxh() ; iptr_type _ndeg = _mesh._tria. node(_nadj)->topo() ; @@ -263,6 +268,10 @@ _hdat._rrad = (real_type)+1.5 * _ball._ball[3] ; + + _bptr->_data._ball[3] = // save to mesh "proper" + (real_type)+1.0 * _ball._ball[3] + * _ball._ball[3] ; _okay = _okay & _geom.intersect ( _pdat, _pred) ; // actual! @@ -297,41 +306,110 @@ (real_type) +.67 ; continue ; } - else { break ; } + else break ; } - + + return ( _kind ) ; + } + + /* + -------------------------------------------------------- + * INIT-BALL: setup ball "collar". + -------------------------------------------------------- + */ + + __static_call + __normal_call void_type init_ball ( + 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 + ) + { /*--------------------- push finalised "collar" nodes */ - for (auto _iter = _pred._list.head() ; - _iter != _pred._list.tend() ; - ++_iter ) + __unreferenced(_hfun) ; + __unreferenced(_mode) ; + __unreferenced(_pass) ; + __unreferenced(_args) ; + __unreferenced(_fpro) ; + + for (auto _bpos = + _mesh._bset._lptr.head () ; + _bpos != + _mesh._bset._lptr.tend () ; + ++_bpos ) { - iptr_type _hint = - _mesh._tria. - node(_nadj)->next() ; - - iptr_type _node = -1; - if (_mesh._tria.push_node( - &_iter->pval(0), _node, - _hint ) ) + if ( *_bpos == nullptr) continue ; + + for (auto _ball = *_bpos ; + _ball != nullptr; + _ball = _ball->_next ) { - _mesh._tria.node( - _node)->idxh() = _hpos ; - _mesh._tria.node( - _node)->fdim() = 1 ; - _mesh._tria.node( - _node)->feat() = 0 ; - _mesh._tria.node( - _node)->topo() = 2 ; - edge_data _edat; - _edat._node[0] = _nadj ; - _edat._node[1] = _node ; + typename + geom_type::ball_type _hdat ; + _hdat._pmid[0] = + _ball->_data._ball[0] ; + _hdat._pmid[1] = + _ball->_data._ball[1] ; + _hdat._pmid[2] = + _ball->_data._ball[2] ; + + _hdat._rrad = + std::sqrt(_ball->_data._ball[3]) ; + + iptr_type _nadj = + _ball->_data._node[0] ; + + mesh::keep_all_3d < + real_type , + iptr_type > _halo ; + + _geom. + intersect(_hdat, _halo) ; + + for (auto _iter = _halo._list.head() ; + _iter != _halo._list.tend() ; + ++_iter ) + { + iptr_type _hint = + _mesh._tria. + node(_nadj)->next() ; + + 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() = 1 ; + + _mesh._tria.node( + _node)->feat() = 0 ; + + _mesh._tria.node( + _node)->topo() = 2 ; + + edge_data _edat; + _edat._node[0] = _nadj ; + _edat._node[1] = _node ; - _pedg.push(_edat) ; + _epro.push(_edat) ; + } + } + } } - - return ( _kind ) ; } diff --git a/src/libcpp/rdel_mesh/rdel_refine_base_2.inc b/src/libcpp/rdel_mesh/rdel_refine_base_2.inc index 4e646a2..857ba0c 100644 --- a/src/libcpp/rdel_mesh/rdel_refine_base_2.inc +++ b/src/libcpp/rdel_mesh/rdel_refine_base_2.inc @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 17 January, 2019 + * Last updated: 24 January, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -95,7 +95,7 @@ typename rdel_opts::node_kind _kind , typename - mesh_type::edge_list &_pedg , + mesh_type::edge_list &_epro , // "protected" 1-edge iptr_list &_nnew , iptr_list &_nold , iptr_list &_tnew , @@ -127,6 +127,9 @@ _hint, &_tnew, &_told)) { /*------------------------- bail if DT push fails */ + _kout = + rdel_opts ::null_kind; + _nnew.set_count(+0) ; _nold.set_count(+0) ; _tnew.set_count(+0) ; @@ -140,9 +143,6 @@ _bscr.set_count(+0) ; _bcav.set_count(+0) ; - _kout = - rdel_opts ::null_kind; - break ; } @@ -170,7 +170,6 @@ _told, _eold ) ; /*------------------------- form "new" rDT cavity */ - mode_type _dim0=null_mode; push_rdel( _geom, _hfun , _mesh, true, _nnew, _tnew, @@ -178,27 +177,27 @@ _tscr, _tcav, _bscr, _bcav, _sign, _pass, - _dim0, _mode, _args ) ; + null_mode, + _mode, _args ) ; - _mesh._tria. - node(_node)->fdim()=_tdim; + _mesh._tria.node( + _node)->fdim() = _tdim ; real_type _pnow[2] ; _pnow[0] = _ppos[0] ; _pnow[1] = _ppos[1] ; if (_cav_bnds( _geom, _mesh , - _told, _pedg, + _told, _epro, _ecav, _eold) ) { - /*-------- new node modifies restricted sub-faces */ - _mesh._tria.roll_back ( - _tnew, _told) ; + /*-------- new node modifies constraint sub-faces */ + _mesh. + _tria.roll_back(_tnew, _told); _tdim = +0 ; - _kout = - rdel_opts::null_kind ; + _kout = rdel_opts::null_kind ; _nnew.set_count(+0) ; _nold.set_count(+0) ; @@ -220,16 +219,15 @@ _ppos, _tdim) ) { /*-------- new node modifies restricted sub-faces */ - _mesh._tria.roll_back ( - _tnew, _told) ; + _mesh. + _tria.roll_back(_tnew, _told); _cav_bnds( _geom, _mesh , _told, _ecav, _eold, _pnow, _ppos, _tdim) ; - _kout = - rdel_opts::circ_kind ; + _kout = rdel_opts::circ_kind ; _nnew.set_count(+0) ; _nold.set_count(+0) ; diff --git a/src/libcpp/rdel_mesh/rdel_refine_base_3.inc b/src/libcpp/rdel_mesh/rdel_refine_base_3.inc index 89471d8..8c0b305 100644 --- a/src/libcpp/rdel_mesh/rdel_refine_base_3.inc +++ b/src/libcpp/rdel_mesh/rdel_refine_base_3.inc @@ -31,7 +31,7 @@ * -------------------------------------------------------- * - * Last updated: 17 January, 2019 + * Last updated: 24 January, 2019 * * Copyright 2013-2019 * Darren Engwirda @@ -97,9 +97,9 @@ typename rdel_opts::node_kind _kind , typename - mesh_type::edge_list &_pedg , + mesh_type::edge_list &_epro , typename - mesh_type::face_list &_pfac , + mesh_type::face_list &_fpro , iptr_list &_nnew , iptr_list &_nold , iptr_list &_tnew , @@ -182,7 +182,6 @@ _told, _eold, _fold ) ; /*------------------------- form "new" rDT cavity */ - mode_type _dim0=null_mode; push_rdel( _geom, _hfun , _mesh, true, _nnew, _tnew, @@ -191,10 +190,11 @@ _tscr, _tcav, _bscr, _bcav, _sign, _pass, - _dim0, _mode, _args ) ; + null_mode, + _mode, _args ) ; - _mesh._tria. - node(_node)->fdim()=_tdim; + _mesh._tria.node( + _node)->fdim() = _tdim; real_type _pnow[3] ; _pnow[0] = _ppos[0] ; @@ -202,19 +202,17 @@ _pnow[2] = _ppos[2] ; if (_cav_bnds( _geom, _mesh, - _told, _pedg, _pfac, + _told, _epro, _fpro, _ecav, _fcav, _eold, _fold) ) { - /*-------- new node modifies restricted sub-faces */ + /*-------- new node modifies constraint sub-faces */ _mesh. - _tria. - roll_back(_tnew, _told) ; + _tria.roll_back(_tnew, _told); _tdim = +0 ; - _kout = - rdel_opts::null_kind ; + _kout = rdel_opts::null_kind ; _nnew.set_count(+0) ; _nold.set_count(+0) ; @@ -232,26 +230,24 @@ _bscr.set_count(+0) ; _bcav.set_count(+0) ; - - break ; } else if (_cav_bnds( _geom, _mesh, - _told, _ecav, _fcav, - _pnow, _ppos, _tdim) ) + _told, + _ecav, _fcav, _pnow, + _ppos, _tdim) ) { /*-------- new node modifies restricted sub-faces */ _mesh. - _tria. - roll_back(_tnew, _told) ; + _tria.roll_back(_tnew, _told); _cav_bnds( _geom, _mesh, - _told, _ecav, _fcav, - _eold, _fold, - _pnow, _ppos, _tdim) ; + _told, _ecav, + _fcav, _eold, + _fold, _pnow, + _ppos, _tdim) ; - _kout = - rdel_opts::circ_kind ; + _kout = rdel_opts::circ_kind ; _nnew.set_count(+0) ; _nold.set_count(+0) ; diff --git a/src/libcpp/rdel_mesh/rdel_refine_face_2.inc b/src/libcpp/rdel_mesh/rdel_refine_face_2.inc index 3bf6c9a..60b6775 100644 --- a/src/libcpp/rdel_mesh/rdel_refine_face_2.inc +++ b/src/libcpp/rdel_mesh/rdel_refine_face_2.inc @@ -58,7 +58,7 @@ mesh_type &_mesh , mode_type _mode , typename - mesh_type::edge_list &_pedg , + mesh_type::edge_list &_epro , iptr_list &_nnew , iptr_list &_nold , iptr_list &_tnew , @@ -152,7 +152,7 @@ _kind = push_node( _geom , _hfun , _mesh, _mode , +1 , _ppos, _kind , - _pedg , _nnew, _nold , + _epro , _nnew, _nold , _tnew , _told, _eold , _ecav , _escr, _tcav , _tscr, @@ -185,7 +185,7 @@ mesh_type &_mesh , mode_type _mode , typename - mesh_type::edge_list &_pedg , + mesh_type::edge_list &_epro , iptr_list &_nnew , iptr_list &_nold , iptr_list &_tnew , @@ -278,7 +278,7 @@ _kind = push_node( _geom , _hfun , _mesh, _mode , +2 , _ppos, _kind , - _pedg , _nnew, _nold , + _epro , _nnew, _nold , _tnew , _told, _eold , _ecav , _escr, _tcav , _tscr, diff --git a/src/libcpp/rdel_mesh/rdel_refine_face_3.inc b/src/libcpp/rdel_mesh/rdel_refine_face_3.inc index 58c46e1..b3bb5a4 100644 --- a/src/libcpp/rdel_mesh/rdel_refine_face_3.inc +++ b/src/libcpp/rdel_mesh/rdel_refine_face_3.inc @@ -58,9 +58,9 @@ mesh_type &_mesh , mode_type _mode , typename - mesh_type::edge_list &_pedg , + mesh_type::edge_list &_epro , typename - mesh_type::face_list &_pfac , + mesh_type::face_list &_fpro , iptr_list &_nnew , iptr_list &_nold , iptr_list &_tnew , @@ -117,7 +117,7 @@ typename mesh_type:: edge_list:: item_type *_eptr = nullptr ; - if (_mesh.find_edge(_edat, _eptr)) + if (_mesh.find_edge(_edat,_eptr)) { if(_eptr->_data._pass == _qdat. _pass ) @@ -161,7 +161,7 @@ _kind = push_node( _geom , _hfun , _mesh, _mode , +1 , _ppos, _kind , - _pedg , _pfac, + _epro , _fpro, _nnew , _nold, _tnew , _told, _eold , _ecav, _escr , @@ -196,9 +196,9 @@ mesh_type &_mesh , mode_type _mode , typename - mesh_type::edge_list &_pedg , + mesh_type::edge_list &_epro , typename - mesh_type::face_list &_pfac , + mesh_type::face_list &_fpro , iptr_list &_nnew , iptr_list &_nold , iptr_list &_tnew , @@ -256,7 +256,7 @@ typename mesh_type:: face_list:: item_type *_fptr = nullptr ; - if (_mesh.find_face(_fdat, _fptr)) + if (_mesh.find_face(_fdat,_fptr)) { if(_fptr->_data._pass == _qdat. _pass ) @@ -299,7 +299,7 @@ _kind = push_node( _geom , _hfun , _mesh, _mode , +2 , _ppos, _kind , - _pedg , _pfac, + _epro , _fpro, _nnew , _nold, _tnew , _told, _eold , _ecav, _escr , @@ -334,9 +334,9 @@ mesh_type &_mesh , mode_type _mode , typename - mesh_type::edge_list &_pedg , + mesh_type::edge_list &_epro , typename - mesh_type::face_list &_pfac , + mesh_type::face_list &_fpro , iptr_list &_nnew , iptr_list &_nold , iptr_list &_tnew , @@ -395,7 +395,7 @@ typename mesh_type:: tria_list:: item_type *_tptr = nullptr ; - if (_mesh.find_tria(_tdat, _tptr)) + if (_mesh.find_tria(_tdat,_tptr)) { if(_tptr->_data._pass == _qdat. _pass ) @@ -437,7 +437,7 @@ _kind = push_node( _geom , _hfun , _mesh, _mode , +3 , _ppos, _kind , - _pedg , _pfac, + _epro , _fpro, _nnew , _nold, _tnew , _told, _eold , _ecav, _escr , diff --git a/src/libcpp/rdel_mesh/rdel_refine_topo_2.inc b/src/libcpp/rdel_mesh/rdel_refine_topo_2.inc index ca6a9e8..21f5887 100644 --- a/src/libcpp/rdel_mesh/rdel_refine_topo_2.inc +++ b/src/libcpp/rdel_mesh/rdel_refine_topo_2.inc @@ -31,9 +31,9 @@ * -------------------------------------------------------- * - * Last updated: 08 April, 2018 + * Last updated: 25 January, 2019 * - * Copyright 2013-2018 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -59,12 +59,12 @@ ) { typename - mesh_type::edge_list _hedg ( + mesh_type::edge_list _edup ( typename mesh_type::edge_hash(), typename mesh_type::edge_pred(), .8, _mesh._eset.get_alloc()) ; - _hedg._lptr.set_count ( + _edup._lptr.set_count ( _tset.count() * +3 , containers::loose_alloc, nullptr); @@ -115,10 +115,10 @@ if(!_mesh. find_edge(_edat, _mptr)) continue ; - if( _hedg.find(_edat, _eptr)) + if( _edup.find(_edat, _eptr)) continue ; else - { _hedg.push(_edat) ; + { _edup.push(_edat) ; } _eset.push_tail(_mptr->_data ) ; @@ -183,6 +183,7 @@ _pmax[1] = _sbal[1]; } } + /*--------------------- return steiner point kind */ return ( _kind ) ; } @@ -201,14 +202,13 @@ mesh_type &_mesh , mode_type _mode , typename - mesh_type::edge_list &_pedg , + mesh_type::edge_list &_epro , iptr_list &_nnew , iptr_list &_nold , iptr_list &_tnew , iptr_list &_told , node_heap &_etpq , - typename - mesh_type::node_list &_etin , + iptr_list &_emrk , edat_list &_etmp , edat_list &_ecav , escr_list &_escr , @@ -222,81 +222,95 @@ ) { class node_pred - { - /*----------------- find adj. set of tria-to-node */ - public : - typedef typename - mesh_type:: - tria_type tria_type ; - - public : - iptr_type _npos; + { + /*--------------------- find adj. set of tria-to-node */ + public : + + typedef typename + mesh_type::tria_type tria_type ; - public : - __inline_call node_pred ( - iptr_type _nsrc - ) : _npos(_nsrc) {} - /*----------------- find adj. set of tria-to-node */ - __inline_call - bool_type operator()( - tria_type&_tria, - iptr_type _tpos, - iptr_type - ) - { return ( - _tria.tria(_tpos)->node(+0) - == this->_npos || - _tria.tria(_tpos)->node(+1) - == this->_npos || - _tria.tria(_tpos)->node(+2) - == this->_npos ); - } - } ; + public : + iptr_type _npos; + + public : + __inline_call node_pred ( + iptr_type _nsrc + ) : _npos(_nsrc) {} + /*--------------------- find adj. set of tria-to-node */ + __inline_call bool_type operator() ( + tria_type&_tria, + iptr_type _tpos, + iptr_type _fpos + ) + { + iptr_type _tnod[3] = { + _tria.tria(_tpos)->node(0) , + _tria.tria(_tpos)->node(1) , + _tria.tria(_tpos)->node(2) + } ; + + __unreferenced(_fpos) ; + + if (_tnod[0] == this->_npos) + return true ; + if (_tnod[1] == this->_npos) + return true ; + if (_tnod[2] == this->_npos) + return true ; + + return false ; + } + } ; typename rdel_opts::node_kind _kind = rdel_opts::null_kind ; /*--------------------- pop() leading element from PQ */ - _nnew.set_count( +0) ; - _nold.set_count( +0) ; - _tnew.set_count( +0) ; - _told.set_count( +0) ; + _nnew.set_count( + 0) ; + _nold.set_count( + 0) ; + _tnew.set_count( + 0) ; + _told.set_count( + 0) ; - _etmp.set_count( +0) ; - _escr.set_count( +0) ; - _ecav.set_count( +0) ; - _tscr.set_count( +0) ; - _tcav.set_count( +0) ; - _bscr.set_count( +0) ; - _bcav.set_count( +0) ; - - iptr_list _tset; - _tset.set_alloc(+16) ; - edat_list _eset; - _eset.set_alloc(+16) ; + _etmp.set_count( + 0) ; + _escr.set_count( + 0) ; + _ecav.set_count( + 0) ; + _tscr.set_count( + 0) ; + _tcav.set_count( + 0) ; + + _bscr.set_count( + 0) ; + _bcav.set_count( + 0) ; + + containers:: + array _tset ; + containers:: + array _eset ; + + _tset.set_alloc( +32) ; + _eset.set_alloc( +32) ; /*--------------------- cycles until a bad disk found */ bool_type _find = false; iptr_type _npos = -1 ; - for ( ; !_etpq.empty() ; ) + for ( ; !_etpq.empty() ; ) { - node_data _ndat, _same ; - _etpq._pop_root( _ndat); - _etin._pop(_ndat,_same); + node_data _ndat; + _etpq._pop_root(_ndat) ; + + _npos = _ndat._node[0] ; - _npos = _ndat._node[ +0] ; + _emrk[_npos] = -1; _tset.set_count( +0) ; _eset.set_count( +0) ; _mesh._tria.walk_node( - _npos , - node_pred(_npos),_tset) ; + _npos , + node_pred(_npos),_tset) ; - if(!etop_disk(_mesh, _npos, - _tset, _eset) ) + if (!etop_disk(_mesh, + _npos, _tset, _eset) ) { - _find = true ; break ; + _find = true; break ; } } if (!_find) return _kind ; @@ -314,7 +328,7 @@ _kind = push_node( _geom, _hfun , _mesh, _mode, _tmax , _pmax, _kind, - _pedg , _nnew, _nold, + _epro , _nnew, _nold, _tnew , _told, _etmp, _ecav , _escr, _tcav , _tscr, diff --git a/src/libcpp/rdel_mesh/rdel_refine_topo_3.inc b/src/libcpp/rdel_mesh/rdel_refine_topo_3.inc index 8cac122..d6aa98d 100644 --- a/src/libcpp/rdel_mesh/rdel_refine_topo_3.inc +++ b/src/libcpp/rdel_mesh/rdel_refine_topo_3.inc @@ -31,9 +31,9 @@ * -------------------------------------------------------- * - * Last updated: 13 April, 2018 + * Last updated: 25 January, 2019 * - * Copyright 2013-2018 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -54,7 +54,8 @@ __normal_call iptr_type ftop_scan ( mesh_type &_mesh, fdat_list &_fset, - typename mesh_type::edge_list &_hedg + typename + mesh_type::edge_list &_edup ) { containers::array < @@ -73,7 +74,7 @@ _face != _fset.tend() ; ++_face, ++_fnum) { - if (_seen [_fnum] == +0 ) + if (_seen [_fnum] == +0) { _seen[_fnum] = +1 ; @@ -84,63 +85,66 @@ for ( ; !_fbfs.empty() ; ) { - _fbfs._pop_tail(_fpos) ; + + _fbfs._pop_tail(_fpos) ; + + for(auto _epos = +3; _epos-- != +0; ) + { + iptr_type _tadj = + _fset[_fpos]._tadj ; + iptr_type _fadj = + _fset[_fpos]._fadj ; + + iptr_type _fnod [ +4]; + mesh_type::tria_type::tria_type:: + face_node(_fnod, _fadj, 3, 2); + _fnod[0] = _mesh._tria. + tria( _tadj)->node(_fnod[0]); + _fnod[1] = _mesh._tria. + tria( _tadj)->node(_fnod[1]); + _fnod[2] = _mesh._tria. + tria( _tadj)->node(_fnod[2]); + + iptr_type _enod [ +3] ; + mesh_type::tria_type::tria_type:: + face_node(_enod, _epos, 2, 1); + _enod[0] = _fnod[_enod[ 0] ] ; + _enod[1] = _fnod[_enod[ 1] ] ; - for(auto _epos = +3; _epos-- != +0; ) - { - iptr_type _tadj = - _fset[_fpos]._tadj ; - iptr_type _fadj = - _fset[_fpos]._fadj ; + algorithms::isort ( + &_enod[0], &_enod[2], + std::less()) ; - iptr_type _fnod [ +4]; - mesh_type::tria_type::tria_type:: - face_node(_fnod, _fadj, 3, 2); - _fnod[0] = _mesh._tria. - tria( _tadj)->node(_fnod[0]); - _fnod[1] = _mesh._tria. - tria( _tadj)->node(_fnod[1]); - _fnod[2] = _mesh._tria. - tria( _tadj)->node(_fnod[2]); + edge_data _edat; + _edat._node[0] = _enod[0] ; + _edat._node[1] = _enod[1] ; - iptr_type _enod [ +3] ; - mesh_type::tria_type::tria_type:: - face_node(_enod, _epos, 2, 1); - _enod[0] = _fnod[_enod[ 0] ] ; - _enod[1] = _fnod[_enod[ 1] ] ; - - algorithms::isort ( - &_enod[0], &_enod[2], - std::less()) ; - - edge_data _edat; - _edat._node[0] = _enod[0] ; - _edat._node[1] = _enod[1] ; - - _aset.set_count(+0, - containers::loose_alloc); + _aset.set_count(+0, + containers::loose_alloc); + + if(!_edup.find(_edat, _aset)) + continue ; - if(!_hedg.find(_edat, _aset)) - continue ; + for (auto _iter = _aset.head() ; + _iter != _aset.tend() ; + ++_iter ) + { + iptr_type _iadj = + (*_iter)->_data._tadj; - for (auto _iter = _aset.head() ; - _iter != _aset.tend() ; - ++_iter ) + if (_seen[_iadj] == +0) { - iptr_type _iadj = - (*_iter)->_data._tadj; + _seen[_iadj] = +1 ; - if (_seen[_iadj] == +0) - { - _seen[_iadj] = +1 ; - _fbfs. - push_tail( _iadj ); - } + _fbfs. + push_tail( _iadj ); } } } } + + } } return ( _fcon ) ; @@ -161,12 +165,12 @@ ) { typename - mesh_type::edge_list _hedg ( + mesh_type::edge_list _edup ( typename mesh_type::edge_hash(), typename mesh_type::edge_pred(), +.8, _mesh._eset.get_alloc()) ; - _hedg._lptr.set_count ( + _edup._lptr.set_count ( _tset.count() * +6 , containers::loose_alloc, nullptr); @@ -217,10 +221,10 @@ if(!_mesh. find_edge(_edat, _mptr)) continue ; - if( _hedg.find(_edat, _eptr)) + if( _edup.find(_edat, _eptr)) continue ; else - { _hedg.push(_edat) ; + { _edup.push(_edat) ; } _eset.push_tail(_mptr->_data ) ; @@ -246,12 +250,12 @@ ) { typename - mesh_type::edge_list _hedg ( + mesh_type::edge_list _edup ( typename mesh_type::edge_hash(), typename mesh_type::edge_pred(), +.8, _mesh._eset.get_alloc()) ; typename - mesh_type::face_list _hfac ( + mesh_type::face_list _fdup ( typename mesh_type::face_hash(), typename mesh_type::face_pred(), +.8, _mesh._fset.get_alloc()) ; @@ -260,10 +264,11 @@ typename mesh_type:: edge_list::item_type*> _aset ; - _hedg._lptr.set_count ( + _edup._lptr.set_count ( _tset.count() * +6 , containers::loose_alloc, nullptr); - _hfac._lptr.set_count ( + + _fdup._lptr.set_count ( _tset.count() * +4 , containers::loose_alloc, nullptr); @@ -311,10 +316,10 @@ if(!_mesh. find_face(_fdat, _mptr)) continue ; - if( _hfac.find(_fdat, _fptr)) + if( _fdup.find(_fdat, _fptr)) continue ; else - { _hfac.push(_fdat) ; + { _fdup.push(_fdat) ; } _fset.push_tail(_mptr->_data ) ; @@ -371,15 +376,15 @@ typename mesh_type:: edge_list:: item_type*_eptr=nullptr; - if(!_hedg.find(_edat, _eptr)) + if(!_edup.find(_edat, _eptr)) { - _hedg.push(_edat) ; + _edup.push(_edat) ; _eset. push_tail(_edat) ; } else { - _hedg.push(_edat) ; + _edup.push(_edat) ; } } } @@ -388,10 +393,10 @@ _edge != _eset.tend() ; ++_edge ) { - _aset.set_count(+0, + _aset.set_count ( +0, containers::loose_alloc); - if(_hedg.find(*_edge, _aset)) + if(_edup.find(*_edge, _aset)) { /*------------------- "internal" topo. degree */ uint_type _topo = +2; @@ -413,7 +418,7 @@ } /*--------------------- count no. connected comp. */ - return ftop_scan(_mesh, _fset, _hedg) == 1; + return ftop_scan(_mesh, _fset, _edup) == 1; } /* @@ -527,16 +532,15 @@ mesh_type &_mesh , mode_type _mode , typename - mesh_type::edge_list &_pedg , + mesh_type::edge_list &_epro , typename - mesh_type::face_list &_pfac , + mesh_type::face_list &_fpro , iptr_list &_nnew , iptr_list &_nold , iptr_list &_tnew , iptr_list &_told , node_heap &_etpq , - typename - mesh_type::node_list &_etin , + iptr_list &_emrk , edat_list &_etmp , edat_list &_ecav , escr_list &_escr , @@ -553,59 +557,69 @@ ) { class node_pred - { - /*----------------- find adj. set of tria-to-node */ - public : - typedef typename - mesh_type:: - tria_type tria_type ; - - public : - iptr_type _npos; - - public : - __inline_call node_pred ( - iptr_type _nsrc - ) : _npos(_nsrc) {} - /*----------------- find adj. set of tria-to-node */ - __inline_call - bool_type operator()( - tria_type&_tria, - iptr_type _tpos, - iptr_type - ) - { return ( - _tria.tria(_tpos)->node(+0) - == this->_npos || - _tria.tria(_tpos)->node(+1) - == this->_npos || - _tria.tria(_tpos)->node(+2) - == this->_npos || - _tria.tria(_tpos)->node(+3) - == this->_npos ); - } - } ; + { + /*--------------------- find adj. set of tria-to-node */ + public : + + typedef typename + mesh_type::tria_type tria_type ; + + public : + iptr_type _npos; + + public : + __inline_call node_pred ( + iptr_type _nsrc + ) : _npos(_nsrc) {} + /*--------------------- find adj. set of tria-to-node */ + __inline_call bool_type operator() ( + tria_type&_tria, + iptr_type _tpos, + iptr_type _fpos + ) + { + iptr_type _tnod[4] = { + _tria.tria(_tpos)->node(0) , + _tria.tria(_tpos)->node(1) , + _tria.tria(_tpos)->node(2) , + _tria.tria(_tpos)->node(3) + } ; + + __unreferenced(_fpos) ; + + if (_tnod[0] == this->_npos) + return true ; + if (_tnod[1] == this->_npos) + return true ; + if (_tnod[2] == this->_npos) + return true ; + if (_tnod[3] == this->_npos) + return true ; + + return false ; + } + } ; typename rdel_opts::node_kind _kind = rdel_opts::null_kind ; /*--------------------- pop() leading element from PQ */ - _nnew.set_count( +0) ; - _nold.set_count( +0) ; - _tnew.set_count( +0) ; - _told.set_count( +0) ; - - _etmp.set_count( +0) ; - _escr.set_count( +0) ; - _ecav.set_count( +0) ; - _ftmp.set_count( +0) ; - _fscr.set_count( +0) ; - _fcav.set_count( +0) ; - _tscr.set_count( +0) ; - _tcav.set_count( +0) ; - - _bscr.set_count( +0) ; - _bcav.set_count( +0) ; + _nnew.set_count( + 0) ; + _nold.set_count( + 0) ; + _tnew.set_count( + 0) ; + _told.set_count( + 0) ; + + _etmp.set_count( + 0) ; + _escr.set_count( + 0) ; + _ecav.set_count( + 0) ; + _ftmp.set_count( + 0) ; + _fscr.set_count( + 0) ; + _fcav.set_count( + 0) ; + _tscr.set_count( + 0) ; + _tcav.set_count( + 0) ; + + _bscr.set_count( + 0) ; + _bcav.set_count( + 0) ; containers:: array _tset ; @@ -620,29 +634,29 @@ /*--------------------- cycles until a bad disk found */ bool_type _find = false; iptr_type _npos = -1 ; - for ( ; !_etpq.empty() ; ) + for ( ; !_etpq.empty() ; ) { - node_data _ndat, _same ; - _etpq._pop_root( _ndat); - _etin._pop(_ndat,_same); + node_data _ndat; + _etpq._pop_root(_ndat) ; + + _npos = _ndat._node[0] ; - _npos = - _ndat._node[ +0] ; + _emrk[_npos] = -1; _tset.set_count( +0) ; _eset.set_count( +0) ; _mesh._tria.walk_node( - _npos , - node_pred(_npos),_tset) ; + _npos , + node_pred(_npos),_tset) ; - if(!etop_disk(_mesh, _npos, - _tset, _eset) ) + if (!etop_disk(_mesh, + _npos, _tset, _eset) ) { _find = true; break ; } } - if (!_find) return ( _kind ); + if (!_find) return _kind ; /*------------------------- calc. best node in cavity */ real_type _pmax [ +4] ; @@ -657,7 +671,7 @@ _kind = push_node( _geom, _hfun , _mesh, _mode, _tmax , _pmax, _kind, - _pedg , _pfac, + _epro , _fpro, _nnew , _nold, _tnew , _told, _etmp , _ecav, _escr, @@ -683,16 +697,15 @@ mesh_type &_mesh , mode_type _mode , typename - mesh_type::edge_list &_pedg , + mesh_type::edge_list &_epro , typename - mesh_type::face_list &_pfac , + mesh_type::face_list &_fpro , iptr_list &_nnew , iptr_list &_nold , iptr_list &_tnew , iptr_list &_told , node_heap &_ftpq , - typename - mesh_type::node_list &_ftin , + iptr_list &_fmrk , edat_list &_etmp , edat_list &_ecav , escr_list &_escr , @@ -709,57 +722,67 @@ ) { class node_pred - { - /*----------------- find adj. set of tria-to-node */ - public : - typedef typename - mesh_type:: - tria_type tria_type ; - - public : - iptr_type _npos; - - public : - __inline_call node_pred ( - iptr_type _nsrc - ) : _npos(_nsrc) {} - /*----------------- find adj. set of tria-to-node */ - __inline_call - bool_type operator()( - tria_type&_tria, - iptr_type _tpos, - iptr_type - ) - { return ( - _tria.tria(_tpos)->node(+0) - == this->_npos || - _tria.tria(_tpos)->node(+1) - == this->_npos || - _tria.tria(_tpos)->node(+2) - == this->_npos || - _tria.tria(_tpos)->node(+3) - == this->_npos ); - } - } ; + { + /*--------------------- find adj. set of tria-to-node */ + public : + + typedef typename + mesh_type::tria_type tria_type ; + + public : + iptr_type _npos; + + public : + __inline_call node_pred ( + iptr_type _nsrc + ) : _npos(_nsrc) {} + /*--------------------- find adj. set of tria-to-node */ + __inline_call bool_type operator() ( + tria_type&_tria, + iptr_type _tpos, + iptr_type _fpos + ) + { + iptr_type _tnod[4] = { + _tria.tria(_tpos)->node(0) , + _tria.tria(_tpos)->node(1) , + _tria.tria(_tpos)->node(2) , + _tria.tria(_tpos)->node(3) + } ; + + __unreferenced(_fpos) ; + + if (_tnod[0] == this->_npos) + return true ; + if (_tnod[1] == this->_npos) + return true ; + if (_tnod[2] == this->_npos) + return true ; + if (_tnod[3] == this->_npos) + return true ; + + return false ; + } + } ; typename rdel_opts::node_kind _kind = rdel_opts::null_kind ; /*--------------------- pop() leading element from PQ */ - _tnew.set_count( +0) ; - _told.set_count( +0) ; - - _etmp.set_count( +0) ; - _escr.set_count( +0) ; - _ecav.set_count( +0) ; - _ftmp.set_count( +0) ; - _fscr.set_count( +0) ; - _fcav.set_count( +0) ; - _tscr.set_count( +0) ; - _tcav.set_count( +0) ; + _tnew.set_count( + 0) ; + _told.set_count( + 0) ; + + _etmp.set_count( + 0) ; + _escr.set_count( + 0) ; + _ecav.set_count( + 0) ; + _ftmp.set_count( + 0) ; + _fscr.set_count( + 0) ; + _fcav.set_count( + 0) ; + _tscr.set_count( + 0) ; + _tcav.set_count( + 0) ; - _bscr.set_count( +0) ; - _bcav.set_count( +0) ; + _bscr.set_count( + 0) ; + _bcav.set_count( + 0) ; containers:: //!! do we need these?? array _tset ; @@ -775,32 +798,32 @@ /*--------------------- cycles until a bad disk found */ bool_type _find = false; iptr_type _npos = -1 ; - for ( ; !_ftpq.empty() ; ) + for ( ; !_ftpq.empty() ; ) { - node_data _ndat, _same ; - _ftpq._pop_root( _ndat); - _ftin._pop(_ndat,_same); + node_data _ndat; + _ftpq._pop_root(_ndat) ; + + _npos = _ndat._node[0] ; - _npos = - _ndat._node[ +0] ; + _fmrk[_npos] = -1; _tset.set_count( +0) ; _eset.set_count( +0) ; _fset.set_count( +0) ; _mesh._tria.walk_node( - _npos , - node_pred(_npos),_tset) ; + _npos , + node_pred(_npos), _tset) ; - if(!ftop_disk(_mesh, _npos, - _tset, - _eset, _fset) ) + if (!ftop_disk(_mesh, + _npos, + _tset, _eset, _fset) ) { _find = true; break ; } } - if (!_find) return ( _kind ); - + if (!_find) return _kind ; + /*------------------------- calc. best node in cavity */ real_type _pmax [ +4] ; char_type _tmax = -1; @@ -816,7 +839,7 @@ _kind = push_node( _geom, _hfun , _mesh, _mode, _tmax , _pmax, _kind, - _pedg , _pfac, + _epro , _fpro, _nnew , _nold, _tnew , _told, _etmp , _ecav, _escr, diff --git a/src/libcpp/rdel_mesh/rdel_test_bounds_2.inc b/src/libcpp/rdel_mesh/rdel_test_bounds_2.inc index 5fe38df..3d7fd6c 100644 --- a/src/libcpp/rdel_mesh/rdel_test_bounds_2.inc +++ b/src/libcpp/rdel_mesh/rdel_test_bounds_2.inc @@ -31,9 +31,9 @@ * -------------------------------------------------------- * - * Last updated: 07 September, 2017 + * Last updated: 24 January, 2019 * - * Copyright 2013-2017 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -43,19 +43,9 @@ // from rdel_mesh_2.hpp - - /* - -------------------------------------------------------- - * GEOM-BNDS: full geometrical encroachment tests. - -------------------------------------------------------- - */ - // Generally, it seems that topology-only type encroa- - // chment is better ==> meshes are sparser, etc... - - // define __geomBNDS - - + # define __GEOMBNDS + /* -------------------------------------------------------- * FIND-RDEL: assemble restricted faces in cavity. @@ -89,8 +79,7 @@ ++_tpos ) { /*------------------------- extract restricted 1-face */ - iptr_type _epos ; - for (_epos = +3; _epos-- != +0; ) + for ( auto _epos = 3; _epos-- != 0; ) { iptr_type _tnod[ +3] ; mesh_type::tria_type::tria_type:: @@ -120,7 +109,7 @@ if (_eset.find(_fdat,_mptr)) continue ; if (_mesh. - find_edge(_fdat,_mptr)) + find_edge(_fdat,_mptr)) _ecav.push_tail(_mptr->_data); _eset.push( _fdat) ; @@ -193,16 +182,16 @@ { _pmax[2] = _sbal[2]; - _pmax[0] = _sbal[0]; - _pmax[1] = _sbal[1]; - - _mode = +1 ; + _mode = + 1 ; _bnds = true ; + + _pmax[0] = _sbal[0]; + _pmax[1] = _sbal[1]; } } - # ifdef __geomBNDS + # ifdef __GEOMBNDS if (_blen < _sbal[2]) { /*----------------- keep worst GEOM. encroachment */ @@ -210,15 +199,15 @@ { _pmax[2] = _sbal[2]; - _pmax[0] = _sbal[0]; - _pmax[1] = _sbal[1]; - - _mode = +1 ; + _mode = + 1 ; _bnds = true ; + + _pmax[0] = _sbal[0]; + _pmax[1] = _sbal[1]; } } - # endif + # endif//__GEOMBNDS } } @@ -306,16 +295,16 @@ { _pmax[2] = _sbal[2]; - _pmax[0] = _sbal[0]; - _pmax[1] = _sbal[1]; - - _mode = +1 ; + _mode = + 1 ; _bnds = true ; + + _pmax[0] = _sbal[0]; + _pmax[1] = _sbal[1]; } } - # ifdef __geomBNDS + # ifdef __GEOMBNDS if (_blen < _sbal[2]) { /*----------------- keep worst GEOM. encroachment */ @@ -323,15 +312,15 @@ { _pmax[2] = _sbal[2]; - _pmax[0] = _sbal[0]; - _pmax[1] = _sbal[1]; - - _mode = +1 ; + _mode = + 1 ; _bnds = true ; + + _pmax[0] = _sbal[0]; + _pmax[1] = _sbal[1]; } } - # endif + # endif//__GEOMBNDS } } @@ -350,7 +339,7 @@ mesh_type &_mesh , iptr_list &_told , typename - mesh_type::edge_list &_pedg , + mesh_type::edge_list &_epro , // "protected" edges edat_list &_enew , edat_list &_eold ) @@ -386,7 +375,7 @@ item_type *_mptr=nullptr; if(!_eset.find(*_epos, _mptr)) { - if( _pedg.find(*_epos, _mptr)) + if( _epro.find(*_epos, _mptr)) { /*----------------- bail on constraint violation! */ return true ; @@ -398,7 +387,7 @@ } - #undef __geomBNDS + #undef __GEOMBNDS diff --git a/src/libcpp/rdel_mesh/rdel_test_bounds_3.inc b/src/libcpp/rdel_mesh/rdel_test_bounds_3.inc index 02ed3af..1d166b0 100644 --- a/src/libcpp/rdel_mesh/rdel_test_bounds_3.inc +++ b/src/libcpp/rdel_mesh/rdel_test_bounds_3.inc @@ -31,9 +31,9 @@ * -------------------------------------------------------- * - * Last updated: 07 September, 2017 + * Last updated: 24 January, 2019 * - * Copyright 2013-2017 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -44,30 +44,7 @@ // from rdel_mesh_3.hpp - /* - -------------------------------------------------------- - * geom-BNDS: full geometrical encroachment tests. - -------------------------------------------------------- - */ - - // Generally, it seems that topology-only type encroa- - // chment is better ==> meshes are sparser, etc... - - // define __geomBNDS - - - /* - -------------------------------------------------------- - * prev-FACE: test against 2-face in prev. cavity. - -------------------------------------------------------- - */ - - // This doesn't seem to be a great idea ==> not a good - // way to split the high-quality surf.-tria. generated - // using DELFRONT... - - // define __prevFACE - + # define __GEOMBNDS /* -------------------------------------------------------- @@ -111,8 +88,7 @@ ++_tpos ) { /*------------------------- extract restricted 1-face */ - iptr_type _epos ; - for (_epos = +6; _epos-- != +0; ) + for ( auto _epos = 6; _epos-- != 0; ) { iptr_type _tnod[ +4] ; mesh_type::tria_type::tria_type:: @@ -148,8 +124,7 @@ _eset.push( _fdat) ; } /*------------------------- extract restricted 2-face */ - iptr_type _fpos ; - for (_fpos = +4; _fpos-- != +0; ) + for ( auto _fpos = 4; _fpos-- != 0; ) { iptr_type _tnod[ +4] ; mesh_type::tria_type::tria_type:: @@ -262,17 +237,17 @@ { _pmax[3] = _sbal[3]; + _mode = + 1 ; + + _bnds = true ; + _pmax[0] = _sbal[0]; _pmax[1] = _sbal[1]; _pmax[2] = _sbal[2]; - - _mode = +1; - - _bnds = true ; } } - # ifdef __geomBNDS + # ifdef __GEOMBNDS if (_blen < _sbal[3]) { /*----------------- keep worst GEOM. encroachment */ @@ -280,16 +255,16 @@ { _pmax[3] = _sbal[3]; + _mode = + 1 ; + + _bnds = true ; + _pmax[0] = _sbal[0]; _pmax[1] = _sbal[1]; _pmax[2] = _sbal[2]; - - _mode = +1; - - _bnds = true ; } } - # endif + # endif//__GEOMBNDS } } @@ -331,17 +306,17 @@ { _pmax[3] = _sbal[3]; + _mode = + 2 ; + + _bnds = true ; + _pmax[0] = _sbal[0]; _pmax[1] = _sbal[1]; _pmax[2] = _sbal[2]; - - _mode = +2; - - _bnds = true ; } } - # ifdef __geomBNDS + # ifdef __GEOMBNDS if (_blen < _sbal[3]) { /*----------------- keep worst GEOM. encroachment */ @@ -349,16 +324,16 @@ { _pmax[3] = _sbal[3]; + _mode = + 2 ; + + _bnds = true ; + _pmax[0] = _sbal[0]; _pmax[1] = _sbal[1]; _pmax[2] = _sbal[2]; - - _mode = +2; - - _bnds = true ; } } - # endif + # endif//__GEOMBNDS } } @@ -459,17 +434,17 @@ { _pmax[3] = _sbal[3]; + _mode = + 1 ; + + _bnds = true ; + _pmax[0] = _sbal[0]; _pmax[1] = _sbal[1]; _pmax[2] = _sbal[2]; - - _mode = +1; - - _bnds = true ; } } - # ifdef __geomBNDS + # ifdef __GEOMBNDS if (_blen < _sbal[3]) { /*----------------- keep worst GEOM. encroachment */ @@ -477,24 +452,22 @@ { _pmax[3] = _sbal[3]; + _mode = + 1 ; + + _bnds = true ; + _pmax[0] = _sbal[0]; _pmax[1] = _sbal[1]; _pmax[2] = _sbal[2]; - - _mode = +1; - - _bnds = true ; } } - # endif + # endif//__GEOMBNDS } } if (_mode > +1) { /*---------------------------------------- test faces */ - __unreferenced( _fold) ; - # ifdef __prevFACE for ( auto _fpos = _fnew.head(); _fpos != _fnew.tend(); ++_fpos ) @@ -535,17 +508,17 @@ { _pmax[3] = _sbal[3]; + _mode = + 2 ; + + _bnds = true ; + _pmax[0] = _sbal[0]; _pmax[1] = _sbal[1]; _pmax[2] = _sbal[2]; - - _mode = +2 ; - - _bnds = true ; } } - # ifdef __geomBNDS + # ifdef __GEOMBNDS if (_blen < _sbal[3]) { /*----------------- keep worst GEOM. encroachment */ @@ -553,18 +526,17 @@ { _pmax[3] = _sbal[3]; + _mode = + 2 ; + + _bnds = true ; + _pmax[0] = _sbal[0]; _pmax[1] = _sbal[1]; _pmax[2] = _sbal[2]; - - _mode = +2 ; - - _bnds = true ; } } - # endif + # endif//__GEOMBNDS } - # endif } return _bnds ; @@ -582,9 +554,9 @@ mesh_type &_mesh , iptr_list &_told , typename - mesh_type::edge_list &_pedg , + mesh_type::edge_list &_epro , // "protected" edges typename - mesh_type::face_list &_pfac , + mesh_type::face_list &_fpro , // "protected" faces edat_list &_enew , fdat_list &_fnew , edat_list &_eold , @@ -636,7 +608,7 @@ item_type *_mptr=nullptr; if(!_eset.find(*_epos, _mptr)) { - if( _pedg.find(*_epos, _mptr)) + if( _epro.find(*_epos, _mptr)) { /*----------------- bail on constraint violation! */ return true ; @@ -653,7 +625,7 @@ item_type *_mptr=nullptr; if(!_fset.find(*_fpos, _mptr)) { - if( _pfac.find(*_fpos, _mptr)) + if( _fpro.find(*_fpos, _mptr)) { /*----------------- bail on constraint violation! */ return true ; @@ -665,9 +637,7 @@ } - #undef __prevFACE - - #undef __geomBNDS + #undef __GEOMBNDS diff --git a/src/libcpp/rdel_mesh/rdel_update_topo_2.inc b/src/libcpp/rdel_mesh/rdel_update_topo_2.inc new file mode 100644 index 0000000..77ef20f --- /dev/null +++ b/src/libcpp/rdel_mesh/rdel_update_topo_2.inc @@ -0,0 +1,223 @@ + + /* + -------------------------------------------------------- + * RDEL-UPDATE-2: update nodes for topo. checks. + -------------------------------------------------------- + * + * 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: 24 January, 2019 + * + * Copyright 2013-2019 + * Darren Engwirda + * de2363@columbia.edu + * https://github.com/dengwirda/ + * + -------------------------------------------------------- + */ + + // from rdel_mesh_2.hpp + + + /* + -------------------------------------------------------- + * FILL-TOPO: enqueue node for topo. + -------------------------------------------------------- + */ + + __static_call + __normal_call void_type fill_topo ( + mesh_type &_mesh, + iptr_type _pass, + node_heap &_etpq, + iptr_list &_emrk, + edat_list &_enew, + edat_list &_eold, + rdel_opts &_opts + ) + { + if (_opts.top1()) + { + if (_eold.count() == +0) + { + _emrk.set_count ( + _mesh._tria._nset.count(), + containers::loose_alloc , -1) ; + + /*------------------------- add all: everything "new" */ + node_data _ndat ; + for (auto _iter = _enew.head(); + _iter != _enew.tend(); + ++_iter ) + { + _ndat._pass = _pass ; + + if (_emrk[_iter->_node[0]] == -1) + { + _emrk[_iter->_node[0]] = +1; + + _ndat. _node[0] = + _iter->_node[0] ; + + _etpq.push(_ndat) ; + } + + if (_emrk[_iter->_node[1]] == -1) + { + _emrk[_iter->_node[1]] = +1; + + _ndat. _node[0] = + _iter->_node[1] ; + + _etpq.push(_ndat) ; + } + } + } + else + { + _emrk.set_count ( + _mesh._tria._nset.count(), + containers::loose_alloc , -1) ; + + { + /*------------------------- init. for local hash obj. */ + typename + mesh_type::edge_list _eset ( + typename mesh_type::edge_hash(), + typename mesh_type::edge_pred(), + +.8, _mesh._eset.get_alloc()) ; + + /*------------------------- push alloc. for hash obj. */ + _eset._lptr.set_count ( + _enew.count() * 2 , + containers::loose_alloc, nullptr); + + /*------------------------- add "old" absent in "new" */ + node_data _ndat ; + for (auto _iter = _enew.head(); + _iter != _enew.tend(); + ++_iter ) + { + _eset.push(*_iter) ; + } + for (auto _iter = _eold.head(); + _iter != _eold.tend(); + ++_iter ) + { + _ndat._pass = _pass; + typename mesh_type:: + edge_list:: + item_type *_mptr = nullptr ; + if(!_eset.find(*_iter, _mptr)) + { + + if (_emrk[_iter->_node[0]] == -1) + { + _emrk[_iter->_node[0]] = +1; + + _ndat. _node[0] = + _iter->_node[0] ; + + _etpq.push(_ndat) ; + } + + if (_emrk[_iter->_node[1]] == -1) + { + _emrk[_iter->_node[1]] = +1; + + _ndat. _node[0] = + _iter->_node[1] ; + + _etpq.push(_ndat) ; + } + + } + } + } + + { + /*------------------------- init. for local hash obj. */ + typename + mesh_type::edge_list _eset ( + typename mesh_type::edge_hash(), + typename mesh_type::edge_pred(), + +.8, _mesh._eset.get_alloc()) ; + + /*------------------------- push alloc. for hash obj. */ + _eset._lptr.set_count ( + _eold.count() * 2 , + containers::loose_alloc, nullptr); + + /*------------------------- add "new" absent in "old" */ + node_data _ndat ; + for (auto _iter = _eold.head(); + _iter != _eold.tend(); + ++_iter ) + { + _eset.push(*_iter) ; + } + for (auto _iter = _enew.head(); + _iter != _enew.tend(); + ++_iter ) + { + _ndat._pass = _pass; + typename mesh_type:: + edge_list:: + item_type *_mptr = nullptr ; + if(!_eset.find(*_iter, _mptr)) + { + + if (_emrk[_iter->_node[0]] == -1) + { + _emrk[_iter->_node[0]] = +1; + + _ndat. _node[0] = + _iter->_node[0] ; + + _etpq.push(_ndat) ; + } + + if (_emrk[_iter->_node[1]] == -1) + { + _emrk[_iter->_node[1]] = +1; + + _ndat. _node[0] = + _iter->_node[1] ; + + _etpq.push(_ndat) ; + } + + } + } + } + } + } + } + + + diff --git a/src/libcpp/rdel_mesh/rdel_update_topo_3.inc b/src/libcpp/rdel_mesh/rdel_update_topo_3.inc new file mode 100644 index 0000000..7cdab0d --- /dev/null +++ b/src/libcpp/rdel_mesh/rdel_update_topo_3.inc @@ -0,0 +1,416 @@ + + /* + -------------------------------------------------------- + * RDEL-UPDATE-3: update nodes for topo. checks. + -------------------------------------------------------- + * + * 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: 24 January, 2019 + * + * Copyright 2013-2019 + * Darren Engwirda + * de2363@columbia.edu + * https://github.com/dengwirda/ + * + -------------------------------------------------------- + */ + + // from rdel_mesh_3.hpp + + + /* + -------------------------------------------------------- + * FILL-TOPO: enqueue node for topo. + -------------------------------------------------------- + */ + + __static_call + __normal_call void_type fill_topo ( + mesh_type &_mesh, + iptr_type _pass, + node_heap &_etpq, + iptr_list &_emrk, + node_heap &_ftpq, + iptr_list &_fmrk, + edat_list &_enew, + edat_list &_eold, + fdat_list &_fnew, + fdat_list &_fold, + rdel_opts &_opts + ) + { + if (_opts.top1()) + { + /*----------------------------- check nodes adj. to edges */ + if (_eold.count() == +0) + { + _emrk.set_count ( + _mesh._tria._nset.count(), + containers::loose_alloc , -1) ; + + /*------------------------- add all: everything "new" */ + node_data _ndat ; + for (auto _iter = _enew.head(); + _iter != _enew.tend(); + ++_iter ) + { + _ndat._pass = _pass ; + + if (_emrk[_iter->_node[0]] == -1) + { + _emrk[_iter->_node[0]] = +1; + + _ndat. _node[0] = + _iter->_node[0] ; + + _etpq.push(_ndat) ; + } + + if (_emrk[_iter->_node[1]] == -1) + { + _emrk[_iter->_node[1]] = +1; + + _ndat. _node[0] = + _iter->_node[1] ; + + _etpq.push(_ndat) ; + } + } + } + else + { + _emrk.set_count ( + _mesh._tria._nset.count(), + containers::loose_alloc , -1) ; + + { + /*------------------------- init. for local hash obj. */ + typename + mesh_type::edge_list _eset ( + typename mesh_type::edge_hash(), + typename mesh_type::edge_pred(), + +.8, _mesh._eset.get_alloc()) ; + + /*------------------------- push alloc. for hash obj. */ + _eset._lptr.set_count ( + _enew.count() * 2 , + containers::loose_alloc, nullptr); + + /*------------------------- add "old" absent in "new" */ + node_data _ndat ; + for (auto _iter = _enew.head(); + _iter != _enew.tend(); + ++_iter ) + { + _eset.push(*_iter) ; + } + for (auto _iter = _eold.head(); + _iter != _eold.tend(); + ++_iter ) + { + _ndat._pass = _pass; + typename mesh_type:: + edge_list:: + item_type *_mptr = nullptr ; + if(!_eset.find(*_iter, _mptr)) + { + + if (_emrk[_iter->_node[0]] == -1) + { + _emrk[_iter->_node[0]] = +1; + + _ndat. _node[0] = + _iter->_node[0] ; + + _etpq.push(_ndat) ; + } + + if (_emrk[_iter->_node[1]] == -1) + { + _emrk[_iter->_node[1]] = +1; + + _ndat. _node[0] = + _iter->_node[1] ; + + _etpq.push(_ndat) ; + } + + } + } + } + + { + /*------------------------- init. for local hash obj. */ + typename + mesh_type::edge_list _eset ( + typename mesh_type::edge_hash(), + typename mesh_type::edge_pred(), + +.8, _mesh._eset.get_alloc()) ; + + /*------------------------- push alloc. for hash obj. */ + _eset._lptr.set_count ( + _eold.count() * 2 , + containers::loose_alloc, nullptr); + + /*------------------------- add "new" absent in "old" */ + node_data _ndat ; + for (auto _iter = _eold.head(); + _iter != _eold.tend(); + ++_iter ) + { + _eset.push(*_iter) ; + } + for (auto _iter = _enew.head(); + _iter != _enew.tend(); + ++_iter ) + { + _ndat._pass = _pass; + typename mesh_type:: + edge_list:: + item_type *_mptr = nullptr ; + if(!_eset.find(*_iter, _mptr)) + { + + if (_emrk[_iter->_node[0]] == -1) + { + _emrk[_iter->_node[0]] = +1; + + _ndat. _node[0] = + _iter->_node[0] ; + + _etpq.push(_ndat) ; + } + + if (_emrk[_iter->_node[1]] == -1) + { + _emrk[_iter->_node[1]] = +1; + + _ndat. _node[0] = + _iter->_node[1] ; + + _etpq.push(_ndat) ; + } + + } + } + } + } + } + + if (_opts.top2()) + { + /*----------------------------- check nodes adj. to faces */ + if (_fold.count() == +0) + { + _fmrk.set_count ( + _mesh._tria._nset.count(), + containers::loose_alloc , -1) ; + + /*------------------------- add all: everything "new" */ + node_data _ndat ; + for (auto _iter = _fnew.head(); + _iter != _fnew.tend(); + ++_iter ) + { + _ndat._pass = _pass ; + + if (_fmrk[_iter->_node[0]] == -1) + { + _fmrk[_iter->_node[0]] = +1; + + _ndat. _node[0] = + _iter->_node[0] ; + + _ftpq.push(_ndat) ; + } + + if (_fmrk[_iter->_node[1]] == -1) + { + _fmrk[_iter->_node[1]] = +1; + + _ndat. _node[0] = + _iter->_node[1] ; + + _ftpq.push(_ndat) ; + } + + if (_fmrk[_iter->_node[2]] == -1) + { + _fmrk[_iter->_node[2]] = +1; + + _ndat. _node[0] = + _iter->_node[2] ; + + _ftpq.push(_ndat) ; + } + } + } + else + { + _fmrk.set_count ( + _mesh._tria._nset.count(), + containers::loose_alloc , -1) ; + + { + /*------------------------- init. for local hash obj. */ + typename + mesh_type::face_list _fset ( + typename mesh_type::face_hash(), + typename mesh_type::face_pred(), + +.8, _mesh._fset.get_alloc()) ; + + /*------------------------- push alloc. for hash obj. */ + _fset._lptr.set_count ( + _fnew.count() * 2 , + containers::loose_alloc, nullptr); + + /*------------------------- add "old" absent in "new" */ + node_data _ndat ; + for (auto _iter = _fnew.head(); + _iter != _fnew.tend(); + ++_iter ) + { + _fset.push(*_iter) ; + } + for (auto _iter = _fold.head(); + _iter != _fold.tend(); + ++_iter ) + { + _ndat._pass = _pass; + typename mesh_type:: + face_list:: + item_type *_mptr = nullptr ; + if(!_fset.find(*_iter, _mptr)) + { + + if (_fmrk[_iter->_node[0]] == -1) + { + _fmrk[_iter->_node[0]] = +1; + + _ndat. _node[0] = + _iter->_node[0] ; + + _ftpq.push(_ndat) ; + } + + if (_fmrk[_iter->_node[1]] == -1) + { + _fmrk[_iter->_node[1]] = +1; + + _ndat. _node[0] = + _iter->_node[1] ; + + _ftpq.push(_ndat) ; + } + + if (_fmrk[_iter->_node[2]] == -1) + { + _fmrk[_iter->_node[2]] = +1; + + _ndat. _node[0] = + _iter->_node[2] ; + + _ftpq.push(_ndat) ; + } + + } + } + } + + { + /*------------------------- init. for local hash obj. */ + typename + mesh_type::face_list _fset ( + typename mesh_type::face_hash(), + typename mesh_type::face_pred(), + +.8, _mesh._fset.get_alloc()) ; + + /*------------------------- push alloc. for hash obj. */ + _fset._lptr.set_count ( + _fold.count() * 2 , + containers::loose_alloc, nullptr); + + /*------------------------- add "new" absent in "old" */ + node_data _ndat ; + for (auto _iter = _fold.head(); + _iter != _fold.tend(); + ++_iter ) + { + _fset.push(*_iter) ; + } + for (auto _iter = _fnew.head(); + _iter != _fnew.tend(); + ++_iter ) + { + _ndat._pass = _pass; + typename mesh_type:: + face_list:: + item_type *_mptr = nullptr ; + if(!_fset.find(*_iter, _mptr)) + { + + if (_fmrk[_iter->_node[0]] == -1) + { + _fmrk[_iter->_node[0]] = +1; + + _ndat. _node[0] = + _iter->_node[0] ; + + _ftpq.push(_ndat) ; + } + + if (_fmrk[_iter->_node[1]] == -1) + { + _fmrk[_iter->_node[1]] = +1; + + _ndat. _node[0] = + _iter->_node[1] ; + + _ftpq.push(_ndat) ; + } + + if (_fmrk[_iter->_node[2]] == -1) + { + _fmrk[_iter->_node[2]] = +1; + + _ndat. _node[0] = + _iter->_node[2] ; + + _ftpq.push(_ndat) ; + } + + } + } + } + } + } + } + + + diff --git a/src/libcpp/tessellate/del_tri_euclidean_2.hpp b/src/libcpp/tessellate/del_tri_euclidean_2.hpp index 6d1aec0..0868a69 100644 --- a/src/libcpp/tessellate/del_tri_euclidean_2.hpp +++ b/src/libcpp/tessellate/del_tri_euclidean_2.hpp @@ -31,9 +31,9 @@ * ------------------------------------------------------------ * - * Last updated: 17 March, 2018 + * Last updated: 22 January, 2019 * - * Copyright 2013-2018 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -89,6 +89,56 @@ _jpos) ) ; } + template class near_pred + { +/*------------------------------------ walk--simplex test */ + public : + __inline_call bool_type operator() ( + mesh_type &_mesh, + __const_ptr ( real_type) _ppos, + real_type _dmin, + iptr_type &_nmin, + iptr_type _tpos, + iptr_type &_fadj + ) const + { + bool_type _done = true ; + iptr_type _fpos ; + + for(_fpos = 3; _fpos-- != 0; ) + { + /*--------------- check dist. wrt. k-th face apex */ + iptr_type _fnod[ 3] ; + mesh_type::tria_type:: + face_node(_fnod, _fpos, 2, 1) ; + /* + _fnod[0] = _mesh. + tria(_tpos)->node(_fnod[0]); + _fnod[1] = _mesh. + tria(_tpos)->node(_fnod[1]); + */ + _fnod[2] = _mesh. + tria(_tpos)->node(_fnod[2]); + + iptr_type _apex = _fnod[2] ; + + real_type _dsqr = lensqr_kd ( + _ppos, + &_mesh.node(_fnod[2])->pval(0)); + + if (_dsqr < _dmin) + { + _done = false; + _dmin = _dsqr; + _nmin = _apex; + _fadj = _fpos; + } + } + + return _done ; + } + } ; + template class walk_pred { /*------------------------------------ walk--simplex test */ @@ -107,7 +157,7 @@ for(_fpos = 3; _fpos-- != 0; ) { /*--------------- test orientation wrt. k-th face */ - iptr_type _fnod[3]; + iptr_type _fnod[ 3] ; mesh_type::tria_type:: face_node(_fnod, _fpos, 2, 1) ; _fnod[0] = _mesh. @@ -142,12 +192,12 @@ { /*------------------------------------ in-circumball test */ public : - __write_ptr ( real_type) _ppos ; + __const_ptr ( real_type) _ppos ; public : - __inline_call circ_pred ( - real_type *_psrc - ) : _ppos( _psrc){} + __inline_call circ_pred ( + __const_ptr ( real_type) _psrc + ) : _ppos( _psrc ) { ; } __inline_call bool_type operator()( mesh_type &_mesh, diff --git a/src/libcpp/tessellate/del_tri_euclidean_3.hpp b/src/libcpp/tessellate/del_tri_euclidean_3.hpp index 115fbb3..281aefb 100644 --- a/src/libcpp/tessellate/del_tri_euclidean_3.hpp +++ b/src/libcpp/tessellate/del_tri_euclidean_3.hpp @@ -31,9 +31,9 @@ * ------------------------------------------------------------ * - * Last updated: 17 March, 2018 + * Last updated: 22 January, 2019 * - * Copyright 2013-2018 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -89,6 +89,58 @@ _jpos) ) ; } + template class near_pred + { +/*------------------------------------ walk--simplex test */ + public : + __inline_call bool_type operator() ( + mesh_type &_mesh, + __const_ptr ( real_type) _ppos, + real_type _dmin, + iptr_type &_nmin, + iptr_type _tpos, + iptr_type &_fadj + ) const + { + bool_type _done = true ; + iptr_type _fpos ; + + for(_fpos = 4; _fpos-- != 0; ) + { + /*--------------- check dist. wrt. k-th face apex */ + iptr_type _fnod[ 4] ; + mesh_type::tria_type:: + face_node(_fnod, _fpos, 3, 2) ; + /* + _fnod[0] = _mesh. + tria(_tpos)->node(_fnod[0]); + _fnod[1] = _mesh. + tria(_tpos)->node(_fnod[1]); + _fnod[2] = _mesh. + tria(_tpos)->node(_fnod[2]); + */ + _fnod[3] = _mesh. + tria(_tpos)->node(_fnod[3]); + + iptr_type _apex = _fnod[3] ; + + real_type _dsqr = lensqr_kd ( + _ppos, + &_mesh.node(_fnod[3])->pval(0)); + + if (_dsqr < _dmin) + { + _done = false; + _dmin = _dsqr; + _nmin = _apex; + _fadj = _fpos; + } + } + + return _done ; + } + } ; + template class walk_pred { /*------------------------------------ walk--simplex test */ @@ -108,7 +160,7 @@ for(_fpos = 4; _fpos-- != 0; ) { /*--------------- test orientation wrt. k-th face */ - iptr_type _fnod[4]; + iptr_type _fnod[ 4] ; mesh_type::tria_type:: face_node(_fnod, _fpos, 3, 2) ; _fnod[0] = _mesh. @@ -153,12 +205,12 @@ { /*------------------------------------ in-circumball test */ public : - __write_ptr ( real_type) _ppos ; + __const_ptr ( real_type) _ppos ; public : - __inline_call circ_pred ( - real_type *_psrc - ) : _ppos( _psrc){} + __inline_call circ_pred ( + __const_ptr ( real_type) _psrc + ) : _ppos( _psrc ) { ; } __inline_call bool_type operator()( mesh_type &_mesh, diff --git a/src/libcpp/tessellate/delaunay_walk_tria.inc b/src/libcpp/tessellate/delaunay_scan_tria.inc similarity index 95% rename from src/libcpp/tessellate/delaunay_walk_tria.inc rename to src/libcpp/tessellate/delaunay_scan_tria.inc index b3644fa..3a2a140 100644 --- a/src/libcpp/tessellate/delaunay_walk_tria.inc +++ b/src/libcpp/tessellate/delaunay_scan_tria.inc @@ -1,7 +1,7 @@ /* -------------------------------------------------------- - * WALK-TRIA-LIST: breadth-first-search about tria. + * SCAN-TRIA-LIST: breadth-first-search about tria. -------------------------------------------------------- * * This program may be freely redistributed under the @@ -46,14 +46,14 @@ /* -------------------------------------------------------- - * WALK-TRIA-LIST: breadth-first-search about tria. + * SCAN-TRIA-LIST: breadth-first-search about tria. -------------------------------------------------------- */ template < typename func_type > - __normal_call void_type walk_tria_list ( + __normal_call void_type scan_tria_list ( iptr_type _elem, fptr_type _flag, func_type &_pred, diff --git a/src/libcpp/tessellate/delaunay_tri_k.hpp b/src/libcpp/tessellate/delaunay_tri_k.hpp index ebc8013..92815e1 100644 --- a/src/libcpp/tessellate/delaunay_tri_k.hpp +++ b/src/libcpp/tessellate/delaunay_tri_k.hpp @@ -31,9 +31,9 @@ * -------------------------------------------------------- * - * Last updated: 02 June, 2018 + * Last updated: 22 January, 2019 * - * Copyright 2013-2018 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -366,20 +366,21 @@ /* -------------------------------------------------------- - * WALK-TRIA-LIST: breadth-first-search about tria. + * SCAN-TRIA-LIST: breadth-first-search about tria. -------------------------------------------------------- */ -# include "delaunay_walk_tria.inc" +# include "delaunay_scan_tria.inc" /* -------------------------------------------------------- - * WALK-MESH-NODE: enclosing tria. via walking. + * WALK-TRIA-NEAR: "nearest" tria. via walking. + * WALK-NODE-NEAR: "nearest" vert. via walking. -------------------------------------------------------- */ -# include "delaunay_walk_node.inc" +# include "delaunay_walk_mesh.inc" /* @@ -423,8 +424,8 @@ */ __normal_call void_type push_root ( - real_type *_pmin, - real_type *_pmax + __const_ptr(real_type) _pmin, + __const_ptr(real_type) _pmax ) { /*---------------------------- de-alloc. existing */ @@ -487,19 +488,35 @@ /* -------------------------------------------------------- - * FIND-NODE: find enclosing triangle for a point. + * FIND-TRIA: find "nearest" triangle for a point. -------------------------------------------------------- */ - __inline_call bool_type find_node ( - real_type *_ppos, + __inline_call bool_type find_tria ( + __const_ptr(real_type) _ppos, iptr_type &_tpos, iptr_type _hint = null_flag() ) { - return walk_mesh_node ( + return walk_tria_near ( _ppos, _tpos, _hint) ; } + + /* + -------------------------------------------------------- + * FIND-NODE: find the "nearest" vert. to a point. + -------------------------------------------------------- + */ + + __inline_call bool_type find_node ( + __const_ptr(real_type) _ppos, + iptr_type &_npos, + iptr_type _hint = null_flag() + ) + { + return walk_node_near ( + _ppos, _npos, _hint) ; + } /* -------------------------------------------------------- @@ -511,7 +528,7 @@ typename list_type > __normal_call bool_type circ_list ( - real_type *_ppos, + __const_ptr(real_type) _ppos, iptr_type &_elem, list_type &_list, iptr_type _hint = null_flag() @@ -526,7 +543,7 @@ template circ_pred< self_type >_pred( _ppos) ; - walk_tria_list (_elem, +1, + scan_tria_list (_elem, +1, _pred, _work) ; /*------------------------------ push index lists */ @@ -546,7 +563,7 @@ */ __inline_call bool_type push_node ( - real_type *_ppos, + __const_ptr(real_type) _ppos, iptr_type &_node, iptr_type _hint = null_flag() ) @@ -564,7 +581,7 @@ typename list_type > __normal_call bool_type push_node ( - real_type *_ppos, + __const_ptr(real_type) _ppos, iptr_type &_node, iptr_type _hint = null_flag() , list_type *_tnew = nullptr , @@ -576,7 +593,7 @@ /*--------------------------- _find enclosing element */ iptr_type _elem = -1; - if (walk_mesh_node(_ppos, _elem, _hint)) + if (walk_tria_near(_ppos, _elem, _hint)) { /*--------------------------- push new node onto list */ @@ -621,7 +638,7 @@ self_type> _pred( _ppos) ; if (_circ == nullptr) - walk_tria_list(_elem, +1 , + scan_tria_list(_elem, +1 , _pred, _work) ; else _work.push_tail(_circ->head() , @@ -743,7 +760,7 @@ { this->_work.clear() ; /*----------------------------------- bfs about _node */ - walk_tria_list ( + scan_tria_list ( node(_node)->next(), +1 , _pred, _work) ; @@ -771,7 +788,7 @@ { this->_work.clear() ; /*----------------------------------- bfs about _tria */ - walk_tria_list (_tria, +1, + scan_tria_list (_tria, +1, _pred, _work) ; /*----------------------------------- push index list */ diff --git a/src/libcpp/tessellate/delaunay_walk_node.inc b/src/libcpp/tessellate/delaunay_walk_mesh.inc similarity index 55% rename from src/libcpp/tessellate/delaunay_walk_node.inc rename to src/libcpp/tessellate/delaunay_walk_mesh.inc index 5311675..5fb9bc2 100644 --- a/src/libcpp/tessellate/delaunay_walk_node.inc +++ b/src/libcpp/tessellate/delaunay_walk_mesh.inc @@ -1,7 +1,7 @@ /* -------------------------------------------------------- - * WALK-MESH-NODE: find enclosing tria. via traversal. + * DELAUNAY-WALK-MESH: find node/cell via traversal. -------------------------------------------------------- * * This program may be freely redistributed under the @@ -31,9 +31,9 @@ * -------------------------------------------------------- * - * Last updated: 12 May, 2017 + * Last updated: 22 January, 2019 * - * Copyright 2013-2017 + * Copyright 2013-2019 * Darren Engwirda * de2363@columbia.edu * https://github.com/dengwirda/ @@ -46,17 +46,17 @@ /* -------------------------------------------------------- - * WALK-MESH-NODE: find enclosing tria. via traversal. + * WALK-TRIA-NEAR: find "nearest" tria. via traversal. -------------------------------------------------------- */ - __normal_call bool_type walk_mesh_node ( - real_type *_ppos, + __normal_call bool_type walk_tria_near ( + __const_ptr(real_type) _ppos, iptr_type &_elem, iptr_type _hint = null_flag () ) { - /*------------- walk mesh until we find enclosing tri */ + /*------------- walk mesh until we find nearest tria. */ bool_type _done = false; bool_type _scan = false; iptr_type _ntri = @@ -133,7 +133,6 @@ } /*------------- walk elements toward _ppos from _elem */ - iptr_type _next = +0 ; iptr_type _face = +0 ; @@ -165,5 +164,135 @@ return ( _done ) ; } + /* + -------------------------------------------------------- + * WALK-NODE-NEAR: find "nearest" vert. via traversal. + -------------------------------------------------------- + */ + + __normal_call bool_type walk_node_near ( + __const_ptr(real_type) _ppos, + iptr_type &_node, + iptr_type _hint = null_flag () + ) + { + /*------------- walk mesh until we find nearest vert. */ + bool_type _done = false; + bool_type _scan = false; + + 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 ) + { + /*----------------------------- return when empty */ + return _done ; + } + if (_hint < +0 || _hint>=_nnod) + { + /*----------------------------- not a valid index */ + _node = +0 ; + _scan = true ; + } + else + if (node(_hint)->mark() <= -1 ) + { + /*----------------------------- not a valid _node */ + _node = +0 ; + _scan = true ; + } + else + { + _node = _hint ; + } + + /*------------- 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(); + + /*----------------------- reject "null" nodes */ + if (node(_next)->mark() < +0) + continue ; + + /*----------------------- calc. dist. to node */ + real_type _dist = + tria_pred::lensqr_kd(_ppos, + &node(_next)->pval(0)); + + if (_best > _dist) + { + /*----------------------- keep closest so far */ + _best = _dist; + _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)); + + iptr_type _next = +0 ; + iptr_type _face = +0 ; + + for (auto _iter = _ntri; + _iter-- != +0; _elem = _next) + { + /*--------------- test nearest-node condition */ + _done = typename tria_pred:: + template near_pred() ( + *this, _ppos, + _dmin, _node, + _elem, _face) ; + + if(!_done ) + { + /*--------------- push next element onto walk */ + _next = __unflip ( + tria(_elem)->next(_face)); + + if(_next == this->null_flag()) + { + /*--------------- stop at external boundaries */ + break ; + } + } + /*--------------- stop if dist. nondecreasing */ + else break ; + } + + return ( _done ) ; + } + diff --git a/src/marche.hpp b/src/marche.hpp new file mode 100644 index 0000000..c1d4e9e --- /dev/null +++ b/src/marche.hpp @@ -0,0 +1,279 @@ + + /* + -------------------------------------------------------- + * + * ,o, ,o, / + * ` ` e88~88e d88~\ /~~~8e Y88b e / + * 888 888 88 88 C888 88b Y88b d8b / + * 888 888 "8b_d8" Y88b e88~-888 Y888/Y88b/ + * 888 888 / 888D C88 888 Y8/ Y8/ + * 88P 888 Cb \_88P "8b_-888 Y Y + * \_8" Y8""8D + * + -------------------------------------------------------- + * MARCHE: "fast-marching" eikonal eqn. solver. + -------------------------------------------------------- + * + * Last updated: 22 January, 2019 + * + * Copyright 2013 -- 2019 + * Darren Engwirda + * darren.engwirda@columbia.edu + * https://github.com/dengwirda + * + -------------------------------------------------------- + * + * 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. + * + -------------------------------------------------------- + */ + + template < + typename jlog_data + > + __normal_call void_type marche_banner ( + jlog_data &_jlog + ) + { + /*-- NB: silliness re. escape sequences */ + _jlog.push ( + " \n" + "#------------------------------------------------------------\n" + "#\n" + "# ,o, ,o, / \n" + "# ` ` e88~88e d88~\\ /~~~8e Y88b e / \n" + "# 888 888 88 88 C888 88b Y88b d8b / \n" + "# 888 888 \"8b_d8\" Y88b e88~-888 Y888/Y88b/ \n" + "# 888 888 / 888D C88 888 Y8/ Y8/ \n" + "# 88P 888 Cb \\_88P \"8b_-888 Y Y \n" + "# \\_8\" Y8\"\"8D \n" + "#\n" + "#------------------------------------------------------------\n" + "# MARCHE: \"fast-marching\" eikonal eqn. solver. \n" + "#------------------------------------------------------------\n" + " \n" + " " __JGSWVSTR "\n\n" + ) ; + } + +# ifdef __lib_jigsaw + +# include "liblib/init_jig_t.hpp" +# include "liblib/init_msh_t.hpp" + +# include "liblib/load_jig_t.hpp" +# include "liblib/load_msh_t.hpp" + +# include "liblib/save_jig_t.hpp" +# include "liblib/save_msh_t.hpp" + + __normal_call iptr_type marche ( // lib-jigsaw + jigsaw_jig_t *_jjig , + jigsaw_msh_t *_fmsh + ) + { + iptr_type _retv = +0; + + hfun_data _ffun ; // FUNC data + jcfg_data _jcfg ; + +# ifdef __use_timers + typename std ::chrono:: + high_resolution_clock:: + time_point _ttic ; + typename std ::chrono:: + high_resolution_clock:: + time_point _ttoc ; + typename std ::chrono:: + high_resolution_clock _time; + + __unreferenced(_time) ; +# endif//__use_timers + + /*--------------------------------- setup *.JLOG data */ + if (_jjig != nullptr ) + { + _jcfg._verbosity = _jjig->_verbosity ; + } + jlog_null _jlog(_jcfg) ; + marche_banner (_jlog) ; + + /*--------------------------------- parse *.JCFG data */ + if (_jjig != nullptr ) + { + _jlog.push ( + " Reading CFG. data...\n\n" ) ; + +# ifdef __use_timers + _ttic = _time.now(); +# endif//__use_timers + + if ((_retv = copy_jcfg ( + _jcfg, + _jlog,*_jjig)) != __no_error) + { + return _retv ; + } + + if ((_retv = test_jcfg ( + _jcfg, _jlog)) != __no_error) + { + return _retv ; + } + + _jlog.push ( + " CFG. data summary...\n\n" ) ; + + if ((_retv = echo_jcfg ( + _jcfg, _jlog)) != __no_error) + { + return _retv ; + } + +# ifdef __use_timers + _ttoc = _time.now(); + _jlog.push(dump_time(_ttic, _ttoc)); +# endif//__use_timers + } + + /*-------------------------- success, if we got here! */ + + return ( _retv ) ; + } + +# else + +# if defined(__cmd_marche) + + __normal_call iptr_type main ( // cmd-marche + int _argc , + char **_argv + ) + { + hfun_data _ffun ; // FUNC data + +# ifdef __use_timers + typename std ::chrono:: + high_resolution_clock:: + time_point _ttic ; + typename std ::chrono:: + high_resolution_clock:: + time_point _ttoc ; + typename std ::chrono:: + high_resolution_clock _time; + + __unreferenced(_time) ; +# endif//__use_timers + + /*-------------------------- find *.JFCG file in args */ + iptr_type _retv = -1 ; + jcfg_data _jcfg ; + for (; _argc-- != +0; ) + { + std::string _ssrc(_argv[_argc]) ; + + std::string _path ; + std::string _name ; + std::string _fext ; + file_part ( _ssrc , + _path , _name , _fext) ; + + if (_ssrc.find("-whoami") == 0) + { + _retv = -2 ; + + std::cout << __JGSWVSTR ; + std::cout << std::endl ; + + break ; + } + + if (_fext.find("jig") == 0) + { + _retv = +0 ; + + _jcfg._jcfg_file =_ssrc ; + + _jcfg._file_path =_path ; + _jcfg._file_name =_name ; + + break ; + } + } + if (_retv != +0) return ( _retv ) ; + + /*--------------------------------- setup *.JLOG file */ + jlog_text _jlog(_jcfg) ; + marche_banner (_jlog) ; + + if(!_jcfg._jcfg_file.empty()) + { + /*--------------------------------- parse *.JCFG file */ + _jlog.push ( + " Reading CFG. file...\n\n" ) ; + +# ifdef __use_timers + _ttic = _time.now(); +# endif//__use_timers + + if ((_retv = read_jcfg ( + _jcfg, _jlog)) != __no_error) + { + return _retv ; + } + + if ((_retv = test_jcfg ( + _jcfg, _jlog)) != __no_error) + { + return _retv ; + } + + _jlog.push ( + " CFG. data summary...\n\n" ) ; + + if ((_retv = echo_jcfg ( + _jcfg, _jlog)) != __no_error) + { + return _retv ; + } + +# ifdef __use_timers + _ttoc = _time.now(); + _jlog.push(dump_time(_ttic, _ttoc)); +# endif//__use_timers + } + +/*-------------------------- success, if we got here! */ + + return ( _retv ) ; + } + +# endif //__cmd_marche + +# endif //__lib_jigsaw + + +