[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
From: |
Yves Renard |
Subject: |
[Getfem-commits] (no subject) |
Date: |
Wed, 20 Sep 2017 04:02:53 -0400 (EDT) |
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
-
}