[Getfem-commits] (no subject)

[Yves Renard]

branch: mb_race_condition_convex_of_reference
commit 3bc266877f4c29d1cd9fe06a703f9b1a447814ce
Author: Yves Renard <address@hidden>
Date:   Wed Sep 20 10:01:56 2017 +0200

    use qhull to simplexify unrepertoried reference elements
---
 src/bgeot_convex_ref.cc                   | 121 +++++++++++++++++++++++++++---
 src/bgeot_convex_ref_simplexified.cc      |   3 +-
 src/getfem/bgeot_convex_ref.h             |   6 +-
 src/getfem/getfem_mesher.h                |   6 --
 src/getfem_contact_and_friction_common.cc |   2 +-
 src/getfem_mesh_level_set.cc              |   2 +-
 src/getfem_mesher.cc                      |  89 +---------------------
 7 files changed, 124 insertions(+), 105 deletions(-)

diff --git a/src/bgeot_convex_ref.cc b/src/bgeot_convex_ref.cc
index f6ad92f..700a526 100644
--- a/src/bgeot_convex_ref.cc
+++ b/src/bgeot_convex_ref.cc
@@ -25,10 +25,96 @@
 #include "getfem/bgeot_comma_init.h"
 
 namespace bgeot {
+  
+  // ******************************************************************
+  //    Interface with qhull
+  // ******************************************************************
+  
+# ifndef GETFEM_HAVE_LIBQHULL_QHULL_A_H
+  void qhull_delaunay(const std::vector<base_node> &,
+                gmm::dense_matrix<size_type>&) {
+    GMM_ASSERT1(false, "Qhull header files not installed. "
+                "Install qhull library and reinstall GetFEM++ library.");
+  }
+# else
+
+  extern "C" {
+    // #ifdef _MSC_VER
+# include <libqhull/qhull_a.h>
+    // #else
+    // # include <qhull/qhull.h>
+    // # include <qhull/mem.h>
+    // # include <qhull/qset.h>
+    // # include <qhull/geom.h>
+    // # include <qhull/merge.h>
+    // # include <qhull/poly.h>
+    // # include <qhull/io.h>
+    // # include <qhull/stat.h>
+    // #endif
+  }
+  
+  void qhull_delaunay(const std::vector<base_node> &pts,
+                     gmm::dense_matrix<size_type>& simplexes) {
+    // cout << "running delaunay with " << pts.size() << " points\n";
+    size_type dim = pts[0].size();   /* points dimension.           */
+    if (pts.size() <= dim) { gmm::resize(simplexes, dim+1, 0); return; }
+    if (pts.size() == dim+1) {
+      gmm::resize(simplexes, dim+1, 1);
+      for (size_type i=0; i <= dim; ++i) simplexes(i, 0) = i;
+      return;
+    }
+    std::vector<coordT> Pts(dim * pts.size());
+    for (size_type i=0; i < pts.size(); ++i)
+      gmm::copy(pts[i], gmm::sub_vector(Pts, gmm::sub_interval(i*dim, dim)));
+    boolT ismalloc=0;  /* True if qhull should free points in
+                        * qh_freeqhull() or reallocation      */
+    /* Be Aware: option QJ could destabilizate all, it can break everything. */
+    /* option Qbb -> QbB (????) */
+    /* option flags for qhull, see qh_opt.htm */
+    char flags[]= "qhull QJ d Qbb Pp T0"; //QJ s i TO";//"qhull Tv";
+    FILE *outfile= 0;    /* output from qh_produce_output()
+                          *  use NULL to skip qh_produce_output() */ 
+    FILE *errfile= stderr;    /* error messages from qhull code */ 
+    int exitcode;             /* 0 if no error from qhull */
+    facetT *facet;                  /* set by FORALLfacets */
+    int curlong, totlong;          /* memory remaining after qh_memfreeshort */
+    vertexT *vertex, **vertexp;
+    exitcode = qh_new_qhull (int(dim), int(pts.size()), &Pts[0], ismalloc,
+                             flags, outfile, errfile);
+    if (!exitcode) { /* if no error */ 
+      size_type nbf=0;
+      FORALLfacets { if (!facet->upperdelaunay) nbf++; }
+      gmm::resize(simplexes, dim+1, nbf);
+        /* 'qh facet_list' contains the convex hull */
+      nbf=0;
+      FORALLfacets {
+        if (!facet->upperdelaunay) {
+          size_type s=0;
+          FOREACHvertex_(facet->vertices) {
+            assert(s < (unsigned)(dim+1));
+            simplexes(s++,nbf) = qh_pointid(vertex->point);
+          }
+          nbf++;
+        }
+      }
+    }
+    qh_freeqhull(!qh_ALL);
+    qh_memfreeshort (&curlong, &totlong);
+    if (curlong || totlong)
+      cerr << "qhull internal warning (main): did not free " << totlong << 
+        " bytes of long memory (" << curlong << " pieces)\n";
+
+  }
+
+#endif
+
+  
 
   size_type simplexified_tab(pconvex_structure cvs, size_type **tab);
 
-  static void simplexify_convex(pconvex_structure cvs, mesh_structure &m) {
+  static void simplexify_convex(bgeot::convex_of_reference *cvr,
+                               mesh_structure &m) {
+    pconvex_structure cvs = cvr->structure();
     m.clear();
     auto basic_cvs = basic_structure(cvs);
     dim_type n = basic_cvs->dim();
@@ -40,9 +126,21 @@ namespace bgeot {
     else {
       size_type *tab;
       size_type nb = simplexified_tab(basic_cvs, &tab);
-      for (size_type nc = 0; nc < nb; ++nc) {
-        for (size_type i = 0; i <= n; ++i) ipts[i] = *tab++;
-        m.add_simplex(n, ipts.begin());
+      if (nb) {
+       for (size_type nc = 0; nc < nb; ++nc) {
+         for (size_type i = 0; i <= n; ++i) ipts[i] = *tab++;
+         m.add_simplex(n, ipts.begin());
+       }
+      }        else {
+#       ifdef GETFEM_HAVE_LIBQHULL_QHULL_A_H
+       cout << "COUCOU\n" << endl;
+       gmm::dense_matrix<size_type> t;
+       delaunay(cvr->points(), t);
+       for (size_type nc = 0; nc < gmm::mat_ncols(t); ++nc) {
+         for (size_type i = 0; i <= n; ++i) ipts[i] = t(i,nc);
+         m.add_simplex(n, ipts.begin());
+       }
+#       endif
       }
     }
   }
@@ -89,15 +187,14 @@ namespace bgeot {
     return p;
   }
 
-  convex_of_reference::convex_of_reference(
-    pconvex_structure cvs, bool auto_basic) :
-     convex<base_node>(move(cvs)), basic_convex_ref_(0), auto_basic(auto_basic)
-  {
+  convex_of_reference::convex_of_reference
+  (pconvex_structure cvs_, bool auto_basic_) :
+    convex<base_node>(move(cvs_)), basic_convex_ref_(0),
+    auto_basic(auto_basic_) {
     DAL_STORED_OBJECT_DEBUG_CREATED(this, "convex of refrence");
     psimplexified_convex = std::make_shared<mesh_structure>();
     // dal::singleton<cleanup_simplexified_convexes>::instance()
     //        .push_back(psimplexified_convex);
-    if (auto_basic) simplexify_convex(structure(), *psimplexified_convex);
   }
 
   bool convex_of_reference::is_basic() const {
@@ -199,6 +296,7 @@ namespace bgeot {
         }
       }
       ppoints = store_point_tab(convex<base_node>::points());
+      if (auto_basic) simplexify_convex(this, *psimplexified_convex);
     }
   };
 
@@ -363,6 +461,7 @@ namespace bgeot {
         convex<base_node>::points()[13] = base_node( 0.0,  0.0, 1.0);
       }
       ppoints = store_point_tab(convex<base_node>::points());
+      if (auto_basic) simplexify_convex(this, *psimplexified_convex);
     }
   };
 
@@ -416,6 +515,7 @@ namespace bgeot {
       convex<base_node>::points()[12] = base_node( 0.0,  0.0, 1.0);
 
       ppoints = store_point_tab(convex<base_node>::points());
+      if (auto_basic) simplexify_convex(this, *psimplexified_convex);
     }
   };
 
@@ -474,6 +574,7 @@ namespace bgeot {
       convex<base_node>::points()[14] = base_node(0.0, 1.0, 1.0);
 
       ppoints = store_point_tab(convex<base_node>::points());
+      if (auto_basic) simplexify_convex(this, *psimplexified_convex);
     }
   };
 
@@ -554,6 +655,7 @@ namespace bgeot {
       if (basic_convex_ref(a) != a || basic_convex_ref(b) != b)
         basic_convex_ref_ = convex_ref_product(basic_convex_ref(a),
                                                basic_convex_ref(b));
+      if (auto_basic) simplexify_convex(this, *psimplexified_convex);
     }
   };
 
@@ -633,6 +735,7 @@ namespace bgeot {
         gmm::scale(normals_[f], 1/gmm::vect_norm2(normals_[f]));
       }
       ppoints = store_point_tab(convex<base_node>::points());
+      if (auto_basic) simplexify_convex(this, *psimplexified_convex);
     }
   };
 
diff --git a/src/bgeot_convex_ref_simplexified.cc 
b/src/bgeot_convex_ref_simplexified.cc
index 4e1a230..4157730 100644
--- a/src/bgeot_convex_ref_simplexified.cc
+++ b/src/bgeot_convex_ref_simplexified.cc
@@ -310,7 +310,8 @@ namespace bgeot {
       return simplexified_pyramid_nb;
     }
 
-    GMM_ASSERT1(false, "No simplexification for this element");
+    *tab = nullptr; // Not taken into account
+    return 0;
   }
 
 
diff --git a/src/getfem/bgeot_convex_ref.h b/src/getfem/bgeot_convex_ref.h
index a6795cd..b69cebe 100644
--- a/src/getfem/bgeot_convex_ref.h
+++ b/src/getfem/bgeot_convex_ref.h
@@ -131,7 +131,11 @@ namespace bgeot {
   inline pconvex_ref basic_convex_ref(pconvex_ref cvr)
   { return cvr->auto_basic ? cvr : cvr->basic_convex_ref_; }
 
-
+  // interface with qhull
+  void qhull_delaunay(const std::vector<base_node> &pts,
+                     gmm::dense_matrix<size_type>& simplexes);
+  
+  
   //@name public functions for obtaining a convex of reference
   //@{
 
diff --git a/src/getfem/getfem_mesher.h b/src/getfem/getfem_mesher.h
index c702f4c..ffd5c38 100644
--- a/src/getfem/getfem_mesher.h
+++ b/src/getfem/getfem_mesher.h
@@ -1067,12 +1067,6 @@ namespace getfem {
 
   inline pmesher_signed_distance new_mesher_torus(scalar_type R, scalar_type r)
   { return std::make_shared<mesher_torus>(R,r); }
-
-
-  // interface with qhull
-  void delaunay(const std::vector<base_node> &pts,
-               gmm::dense_matrix<size_type>& simplexes);
-
   
   // mesher
   void build_mesh(mesh &m, const pmesher_signed_distance& dist_,
diff --git a/src/getfem_contact_and_friction_common.cc 
b/src/getfem_contact_and_friction_common.cc
index 385ce8c..ed7fc19 100644
--- a/src/getfem_contact_and_friction_common.cc
+++ b/src/getfem_contact_and_friction_common.cc
@@ -460,7 +460,7 @@ namespace getfem {
     //   gmm::fill_random(rr);
     //   boundary_points[i] += 1E-9*rr;
     // }
-    getfem::delaunay(boundary_points, simplexes);
+    bgeot::qhull_delaunay(boundary_points, simplexes);
 
     // connectivity analysis
     for (size_type is = 0; is < gmm::mat_ncols(simplexes); ++is) {
diff --git a/src/getfem_mesh_level_set.cc b/src/getfem_mesh_level_set.cc
index 0f2917f..98cc1eb 100644
--- a/src/getfem_mesh_level_set.cc
+++ b/src/getfem_mesh_level_set.cc
@@ -306,7 +306,7 @@ struct Chrono {
     double t0=gmm::uclock_sec();
     if (noisy) cout << "running delaunay with " << fixed_points.size()
                    << " points.." << std::flush;
-    delaunay(fixed_points, simplexes);
+    bgeot::qhull_delaunay(fixed_points, simplexes);
     if (noisy) cout << " -> " << gmm::mat_ncols(simplexes)
                    << " simplexes [" << gmm::uclock_sec()-t0 << "sec]\n";
   }
diff --git a/src/getfem_mesher.cc b/src/getfem_mesher.cc
index 0127b2c..b3bca0f 100644
--- a/src/getfem_mesher.cc
+++ b/src/getfem_mesher.cc
@@ -1171,7 +1171,7 @@ namespace getfem {
         cout << "NEW DELAUNAY, running on " << pts.size() << " points\n";
       size_type nbpt = pts.size();
       add_point_hull();
-      delaunay(pts, t);
+      bgeot::qhull_delaunay(pts, t);
       pts.resize(nbpt);
       if (noisy > 1) cout << "number of elements before selection = "
                           << gmm::mat_ncols(t) << "\n";
@@ -1309,7 +1309,7 @@ namespace getfem {
           control_mesh_surface();
           size_type nbpt = pts.size();
           add_point_hull();          
-          delaunay(pts, t);
+         bgeot::qhull_delaunay(pts, t);
           pts.resize(nbpt);
           select_elements((prefind == 3) ? 0 : 1);
           suppress_flat_boundary_elements();
@@ -1343,7 +1343,7 @@ namespace getfem {
             control_mesh_surface();
             nbpt = pts.size();
             add_point_hull();
-            delaunay(pts, t);
+           bgeot::qhull_delaunay(pts, t);
             pts.resize(nbpt);
             select_elements((prefind == 3) ? 0 : 1);
             suppress_flat_boundary_elements();
@@ -1426,88 +1426,5 @@ namespace getfem {
               iter_max, prefind, dist_point_hull, boundary_threshold_flatness);
   }
   
-
-  // ******************************************************************
-  //    Interface with qhull
-  // ******************************************************************
-
-# ifndef GETFEM_HAVE_LIBQHULL_QHULL_A_H
-  void delaunay(const std::vector<base_node> &,
-                gmm::dense_matrix<size_type>&) {
-    GMM_ASSERT1(false, "Qhull header files not installed. "
-                "Install qhull library and reinstall GetFEM++ library.");
-  }
-# else
-
-  extern "C" {
-// #ifdef _MSC_VER
-# include <libqhull/qhull_a.h>
-// #else
-// # include <qhull/qhull.h>
-// # include <qhull/mem.h>
-// # include <qhull/qset.h>
-// # include <qhull/geom.h>
-// # include <qhull/merge.h>
-// # include <qhull/poly.h>
-// # include <qhull/io.h>
-// # include <qhull/stat.h>
-// #endif
-}
-
-  void delaunay(const std::vector<base_node> &pts,
-                gmm::dense_matrix<size_type>& simplexes) {
-    // cout << "running delaunay with " << pts.size() << " points\n";
-    size_type dim = pts[0].size();   /* points dimension.           */
-    if (pts.size() <= dim) { gmm::resize(simplexes, dim+1, 0); return; }
-    if (pts.size() == dim+1) {
-      gmm::resize(simplexes, dim+1, 1);
-      for (size_type i=0; i <= dim; ++i) simplexes(i, 0) = i;
-      return;
-    }
-    std::vector<coordT> Pts(dim * pts.size());
-    for (size_type i=0; i < pts.size(); ++i)
-      gmm::copy(pts[i], gmm::sub_vector(Pts, gmm::sub_interval(i*dim, dim)));
-    boolT ismalloc=0;  /* True if qhull should free points in
-                        * qh_freeqhull() or reallocation      */
-    /* Be Aware: option QJ could destabilizate all, it can break everything. */
-    /* option Qbb -> QbB (????) */
-    /* option flags for qhull, see qh_opt.htm */
-    char flags[]= "qhull QJ d Qbb Pp T0"; //QJ s i TO";//"qhull Tv";
-    FILE *outfile= 0;    /* output from qh_produce_output()
-                          *  use NULL to skip qh_produce_output() */ 
-    FILE *errfile= stderr;    /* error messages from qhull code */ 
-    int exitcode;             /* 0 if no error from qhull */
-    facetT *facet;                  /* set by FORALLfacets */
-    int curlong, totlong;          /* memory remaining after qh_memfreeshort */
-    vertexT *vertex, **vertexp;
-    exitcode = qh_new_qhull (int(dim), int(pts.size()), &Pts[0], ismalloc,
-                             flags, outfile, errfile);
-    if (!exitcode) { /* if no error */ 
-      size_type nbf=0;
-      FORALLfacets { if (!facet->upperdelaunay) nbf++; }
-      gmm::resize(simplexes, dim+1, nbf);
-        /* 'qh facet_list' contains the convex hull */
-      nbf=0;
-      FORALLfacets {
-        if (!facet->upperdelaunay) {
-          size_type s=0;
-          FOREACHvertex_(facet->vertices) {
-            assert(s < (unsigned)(dim+1));
-            simplexes(s++,nbf) = qh_pointid(vertex->point);
-          }
-          nbf++;
-        }
-      }
-    }
-    qh_freeqhull(!qh_ALL);
-    qh_memfreeshort (&curlong, &totlong);
-    if (curlong || totlong)
-      cerr << "qhull internal warning (main): did not free " << totlong << 
-        " bytes of long memory (" << curlong << " pieces)\n";
-
-  }
-
-#endif
-
 }