Merge pull request #13 from dengwirda/dev

0.9.8 updates: "Off-centre" fix -- dual cell checks

Author: dengwirda

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):
 #   ----------------------

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
+        ) ;
 
     /* 
     --------------------------------------------------------

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"
 

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!
+    
+    }
+     
+    } ;
+    
+    
+    

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__)

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
+     * [email protected]
+     * 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
                 <real_type>::epsilon(),+.67);
 
@@ -417,17 +397,20 @@
                   _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 = 
                 geometry::lensqr_2d(
                     _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 ;