[mshr] 01/03: New upstream version 2017.1.0
Johannes Ring
johannr-guest at moszumanska.debian.org
Thu May 11 11:55:06 UTC 2017
This is an automated email from the git hooks/post-receive script.
johannr-guest pushed a commit to branch experimental
in repository mshr.
commit 8240543801c270929637ebd8a5cb700bda41b469
Author: Johannes Ring <johannr at simula.no>
Date: Thu May 11 13:00:32 2017 +0200
New upstream version 2017.1.0
---
.gitignore | 1 +
CMakeLists.txt | 4 +-
ChangeLog | 13 -----
ChangeLog.rst | 32 ++++++++++++
include/mshr/CSGCGALDomain2D.h | 4 +-
include/mshr/CSGGeometry.h | 19 +++++--
include/mshr/CSGOperators.h | 44 +++++++++++-----
include/mshr/CSGPrimitives2D.h | 15 +++++-
include/mshr/CSGPrimitives3D.h | 9 +++-
release.conf | 2 +-
src/CSGCGALDomain2D.cpp | 72 ++++++++++++++------------
src/CSGCGALDomain3D.cpp | 8 ++-
src/CSGCGALMeshGenerator2D.cpp | 14 +++--
src/CSGGeometry.cpp | 62 ++++++++++++++++++++--
src/CSGOperators.cpp | 115 ++++++++++++++++++++++++++++++++++++++++-
src/CSGPrimitives2D.cpp | 82 +++++++++++++++++++++++++++--
src/CSGPrimitives3D.cpp | 14 ++++-
src/MeshGenerator.cpp | 30 +++--------
swig/CMakeLists.txt | 9 ++--
test/CMakeLists.txt | 10 ++++
test/test-csg-predicates.py | 34 ++++++++++++
test/test-num-segments-2d.py | 13 +++++
22 files changed, 491 insertions(+), 115 deletions(-)
diff --git a/.gitignore b/.gitignore
index 1207296..c63635d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,3 @@
build/
build.*/
+release/
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 12c290b..e24e473 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,7 +1,7 @@
project( MSHR )
set(MSHR_VERSION_RELEASE 1)
-set(MSHR_VERSION_MAJOR "2016")
-set(MSHR_VERSION_MINOR "2")
+set(MSHR_VERSION_MAJOR "2017")
+set(MSHR_VERSION_MINOR "1")
set(MSHR_VERSION_MICRO "0")
set(MSHR_VERSION "${MSHR_VERSION_MAJOR}.${MSHR_VERSION_MINOR}.${MSHR_VERSION_MICRO}")
if (NOT MSHR_VERSION_RELEASE)
diff --git a/ChangeLog b/ChangeLog
deleted file mode 100644
index 6c31a5b..0000000
--- a/ChangeLog
+++ /dev/null
@@ -1,13 +0,0 @@
-2016.2.0 [2016-11-30]
- - Improvements and bugfixes
- - Upgrade to CGAL 4.9 and use CGAL as header only library
-2016.1.0 [2016-06-23]
- - Major improvements, in particular boundary conditions
-1.6.0 [2015-07-28]
- - Upgrade to CGAL 4.6
- - Add more demos
- - Faster and more memory efficient implementation of csg operations
- with -DENABLE_EXPERIMENTAL=True
- - Bugfixes and cleanups
-1.5.0 [2015-02-04]
- - Initial release of mshr.
diff --git a/ChangeLog.rst b/ChangeLog.rst
new file mode 100644
index 0000000..601f374
--- /dev/null
+++ b/ChangeLog.rst
@@ -0,0 +1,32 @@
+Changelog
+=========
+
+2017.1.0 (2017-05-11)
+---------------------
+
+- Bugfixes
+
+2016.2.0 (2016-11-30)
+---------------------
+
+- Improvements and bugfixes
+- Upgrade to CGAL 4.9 and use CGAL as header only library
+
+2016.1.0 (2016-06-23)
+---------------------
+
+- Major improvements, in particular boundary conditions
+
+1.6.0 (2015-07-28)
+------------------
+
+- Upgrade to CGAL 4.6
+- Add more demos
+- Faster and more memory efficient implementation of csg operations
+ with -DENABLE_EXPERIMENTAL=True
+- Bugfixes and cleanups
+
+1.5.0 (2015-02-04)
+------------------
+
+- Initial release of mshr.
diff --git a/include/mshr/CSGCGALDomain2D.h b/include/mshr/CSGCGALDomain2D.h
index ba3f2a8..35bca27 100644
--- a/include/mshr/CSGCGALDomain2D.h
+++ b/include/mshr/CSGCGALDomain2D.h
@@ -1,4 +1,4 @@
-// Copyright (C) 2013-2014 Benjamin Kehlet
+// Copyright (C) 2013-2017 Benjamin Kehlet
//
// This file is part of mshr.
//
@@ -43,7 +43,7 @@ class CSGCGALDomain2D : public dolfin::Variable
/// @brief Construct polygon from Dolfin CSG geometry
CSGCGALDomain2D(
std::shared_ptr<const CSGGeometry> geometry,
- double segment_granularity=0.1
+ double segment_granularity
);
/// @brief Destructor
diff --git a/include/mshr/CSGGeometry.h b/include/mshr/CSGGeometry.h
index 1c59df0..cdbb5fd 100644
--- a/include/mshr/CSGGeometry.h
+++ b/include/mshr/CSGGeometry.h
@@ -1,4 +1,4 @@
-// Copyright (C) 2012 Anders Logg and 2013-2014 Benjamin Kehlet
+// Copyright (C) 2012 Anders Logg and 2013-2017 Benjamin Kehlet
//
// This file is part of mshr.
//
@@ -6,12 +6,12 @@
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
-//
+//
// mshr is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
-//
+//
// You should have received a copy of the GNU General Public License
// along with mshr. If not, see <http://www.gnu.org/licenses/>.
//
@@ -27,6 +27,7 @@
#include <list>
#include <dolfin/common/Variable.h>
+#include <dolfin/geometry/Point.h>
namespace mshr
{
@@ -66,6 +67,18 @@ namespace mshr
/// @brief Return const list of subdomain geometries
const std::list<std::pair<std::size_t, std::shared_ptr<const CSGGeometry>>>& get_subdomains() const { return subdomains; }
+ // These functions are (for now) implemented for 2D only.
+ std::pair<dolfin::Point, double> estimate_bounding_sphere(std::size_t numSamples=100) const;
+
+ // Compute axis aligned box that is guaranteed to bound the geometry. May overshoot.
+ virtual std::pair<dolfin::Point, dolfin::Point> bounding_box() const = 0;
+
+ // Test if given point is inside geometry. 2D only (for now)
+ virtual bool inside(dolfin::Point p) const = 0;
+
+ // Test if line segment is entirely inside geometry. 2D only (for now)
+ virtual bool inside(dolfin::Point p1, dolfin::Point p2) const;
+
enum Type { Box,
Sphere,
Cylinder,
diff --git a/include/mshr/CSGOperators.h b/include/mshr/CSGOperators.h
index 8592037..35a58ac 100644
--- a/include/mshr/CSGOperators.h
+++ b/include/mshr/CSGOperators.h
@@ -1,4 +1,4 @@
-// Copyright (C) 2012 Anders Logg, 2012-2014 Benjamin Kehlet
+// Copyright (C) 2012 Anders Logg, 2012-2017 Benjamin Kehlet
//
// This file is part of mshr.
//
@@ -6,7 +6,7 @@
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
-//
+//
// mshr is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
@@ -57,6 +57,9 @@ namespace mshr
/// @param verbose vervosity level
std::string str(bool verbose) const;
+ std::pair<dolfin::Point, dolfin::Point> bounding_box() const;
+ bool inside(dolfin::Point p) const;
+
Type getType() const { return CSGGeometry::Union; }
const std::shared_ptr<CSGGeometry> _g0;
@@ -77,6 +80,9 @@ namespace mshr
/// @brief get informal string representation
std::string str(bool verbose) const;
+ std::pair<dolfin::Point, dolfin::Point> bounding_box() const;
+ bool inside(dolfin::Point p) const;
+
Type getType() const { return CSGGeometry::Difference; }
const std::shared_ptr<CSGGeometry> _g0;
@@ -98,6 +104,9 @@ namespace mshr
/// @brief get informal string representation
std::string str(bool verbose) const;
+ std::pair<dolfin::Point, dolfin::Point> bounding_box() const;
+ bool inside(dolfin::Point p) const;
+
Type getType() const { return CSGGeometry::Intersection; }
const std::shared_ptr<CSGGeometry> _g0;
@@ -112,7 +121,7 @@ namespace mshr
{
public:
- /// @brief create translation
+ /// @brief create translation
/// @param g a CSG geometry
/// @param t the translation vector
CSGTranslation(std::shared_ptr<CSGGeometry> g,
@@ -121,6 +130,9 @@ namespace mshr
/// @brief get informal string representation
std::string str(bool verbose) const;
+ std::pair<dolfin::Point, dolfin::Point> bounding_box() const;
+ bool inside(dolfin::Point p) const;
+
Type getType() const { return CSGGeometry::Translation; }
const std::shared_ptr<CSGGeometry> g;
@@ -136,22 +148,23 @@ namespace mshr
/// @brief Create scaling
/// @param g a CSG geometry
- /// @param scale_factor
+ /// @param scale_factor
CSGScaling(std::shared_ptr<CSGGeometry> g,
double scale_factor);
/// @brief Scale (translated) geometry. Geometry will be translated, scaled and translated back
/// @param g a CSG geometry
- /// @param c translation
+ /// @param c translation
/// @param scale_factor the scale factor
CSGScaling(std::shared_ptr<CSGGeometry> g,
dolfin::Point c,
double scale_factor);
-
-
std::string str(bool verbose) const;
+ std::pair<dolfin::Point, dolfin::Point> bounding_box() const;
+ bool inside(dolfin::Point p) const;
+
Type getType() const { return CSGGeometry::Scaling; }
const std::shared_ptr<CSGGeometry> g;
@@ -180,7 +193,7 @@ namespace mshr
dolfin::Point v,
double theta);
- /// @brief create 3D rotation
+ /// @brief create 3D rotation
/// @param g A CSG geometry
/// @param rot_axis The rotation axis
/// @param rot_center The rotation center
@@ -192,6 +205,9 @@ namespace mshr
Type getType() const { return CSGGeometry::Rotation; }
std::string str(bool verbose) const;
+ std::pair<dolfin::Point, dolfin::Point> bounding_box() const;
+ bool inside(dolfin::Point p) const;
+
const std::shared_ptr<CSGGeometry> g;
const dolfin::Point rot_axis;
const dolfin::Point c;
@@ -264,28 +280,28 @@ namespace mshr
/// Create intersection of two geometries
inline std::shared_ptr<CSGIntersection> operator*(std::shared_ptr<CSGGeometry> g0,
- std::shared_ptr<CSGGeometry> g1)
+ std::shared_ptr<CSGGeometry> g1)
{
return std::shared_ptr<CSGIntersection>(new CSGIntersection(g0, g1));
}
/// Create intersection of two geometries
inline std::shared_ptr<CSGIntersection> operator*(CSGGeometry& g0,
- std::shared_ptr<CSGGeometry> g1)
+ std::shared_ptr<CSGGeometry> g1)
{
return reference_to_no_delete_pointer(g0) * g1;
}
/// Create intersection of two geometries
inline std::shared_ptr<CSGIntersection> operator*(std::shared_ptr<CSGGeometry> g0,
- CSGGeometry& g1)
+ CSGGeometry& g1)
{
return g0 * reference_to_no_delete_pointer(g1);
}
/// Create intersection of two geometries
inline std::shared_ptr<CSGIntersection> operator*(CSGGeometry& g0,
- CSGGeometry& g1)
+ CSGGeometry& g1)
{
return reference_to_no_delete_pointer(g0) * reference_to_no_delete_pointer(g1);
}
@@ -299,9 +315,9 @@ namespace mshr
//--- Scaling operators ---
inline std::shared_ptr<CSGScaling> operator*(std::shared_ptr<CSGGeometry> g,
- double s)
+ double s)
{
- return std::shared_ptr<CSGScaling>(new CSGScaling(g, s, false));
+ return std::shared_ptr<CSGScaling>(new CSGScaling(g, s));
}
}
diff --git a/include/mshr/CSGPrimitives2D.h b/include/mshr/CSGPrimitives2D.h
index 6496bc9..8789bef 100644
--- a/include/mshr/CSGPrimitives2D.h
+++ b/include/mshr/CSGPrimitives2D.h
@@ -1,4 +1,4 @@
-// Copyright (C) 2012 Anders Logg, 2012-2015 Benjamin Kehlet
+// Copyright (C) 2012 Anders Logg, 2012-2017 Benjamin Kehlet
//
// This file is part of mshr.
//
@@ -71,6 +71,9 @@ namespace mshr
/// @brief get number of segments used when computing polygonal approximation
std::size_t segments() const { return _segments; }
+ std::pair<dolfin::Point, dolfin::Point> bounding_box() const;
+ bool inside(dolfin::Point p) const;
+
private:
dolfin::Point c;
double _r;
@@ -111,6 +114,10 @@ namespace mshr
/// @brief get resolution when computing polygonal approximation
std::size_t segments() const { return _segments; }
+ std::pair<dolfin::Point, dolfin::Point> bounding_box() const;
+ bool inside(dolfin::Point p) const;
+
+
private:
dolfin::Point c;
double _a, _b;
@@ -142,6 +149,9 @@ namespace mshr
/// @brief get second corner
dolfin::Point second_corner() const { return b; }
+ std::pair<dolfin::Point, dolfin::Point> bounding_box() const;
+ bool inside(dolfin::Point p) const;
+
private:
const dolfin::Point a, b;
};
@@ -170,6 +180,9 @@ namespace mshr
/// @brief get vertices in polygon. Can be modified inplace.
const std::vector<dolfin::Point>& vertices() const { return _vertices; }
+ std::pair<dolfin::Point, dolfin::Point> bounding_box() const;
+ bool inside(dolfin::Point p) const;
+
private:
const std::vector<dolfin::Point> _vertices;
};
diff --git a/include/mshr/CSGPrimitives3D.h b/include/mshr/CSGPrimitives3D.h
index 66da580..433957a 100644
--- a/include/mshr/CSGPrimitives3D.h
+++ b/include/mshr/CSGPrimitives3D.h
@@ -1,4 +1,4 @@
-// Copyright (C) 2012 Anders Logg, 2012-2015 Benjamin Kehlet
+// Copyright (C) 2012 Anders Logg, 2012-2017 Benjamin Kehlet
//
// This file is part of mshr.
//
@@ -39,6 +39,13 @@ namespace mshr
/// @return get dimension of geometry
std::size_t dim() const { return 3; }
+ // Compute axis aligned box that is guaranteed to bound the geometry. May overshoot.
+ std::pair<dolfin::Point, dolfin::Point> bounding_box() const;
+
+ // Test if given point is inside geometry. 2D only (for now)
+ bool inside(dolfin::Point p) const;
+
+
};
// TODO: Make sphere a special case of ellipsoid.
diff --git a/release.conf b/release.conf
index 39f7747..309f558 100644
--- a/release.conf
+++ b/release.conf
@@ -2,4 +2,4 @@
PACKAGE="mshr"
BRANCH="master"
-FILES="ChangeLog CMakeLists.txt"
+FILES="ChangeLog.rst CMakeLists.txt"
diff --git a/src/CSGCGALDomain2D.cpp b/src/CSGCGALDomain2D.cpp
index 638efc2..061a145 100644
--- a/src/CSGCGALDomain2D.cpp
+++ b/src/CSGCGALDomain2D.cpp
@@ -224,15 +224,20 @@ Polygon_2 make_circle(
)
{
unsigned int num_segments = 0;
- if (c->segments() > 0) {
+ if (c->segments() > 0)
+ {
num_segments = c->segments();
} else {
dolfin_assert(segment_granularity > 0.0);
+
num_segments = std::round(
(2 * DOLFIN_PI * c->radius()) / segment_granularity
);
}
+ // A polygon with less than 3 segments is degenerate
+ num_segments = std::max(num_segments, 3u);
+
std::vector<Point_2> pts;
pts.reserve(num_segments);
@@ -247,15 +252,15 @@ Polygon_2 make_circle(
return Polygon_2(pts.begin(), pts.end());
}
//-----------------------------------------------------------------------------
-Polygon_2 make_ellipse(
- const Ellipse* e,
- const double segment_granularity
- )
+Polygon_2 make_ellipse(const Ellipse* e,
+ const double segment_granularity)
{
unsigned int num_segments = 0;
if (e->segments() > 0) {
num_segments = e->segments();
- } else {
+ }
+ else
+ {
dolfin_assert(segment_granularity > 0.0);
// https://en.wikipedia.org/wiki/Ellipse#Circumference
const double a_min_b = e->a() - e->b();
@@ -267,6 +272,10 @@ Polygon_2 make_ellipse(
num_segments = std::round(arc_length_ellipse_approx / segment_granularity);
}
+ // A polygon with less segments than 3 is degenerate
+ // FIXME: Should this result in an empty polygon?
+ num_segments = std::max(num_segments, 3u);
+
std::vector<Point_2> pts;
pts.reserve(num_segments);
@@ -347,9 +356,8 @@ CSGCGALDomain2D::~CSGCGALDomain2D()
{
}
//-----------------------------------------------------------------------------
-CSGCGALDomain2D::CSGCGALDomain2D(
- std::shared_ptr<const CSGGeometry> geometry,
- double segment_granularity)
+CSGCGALDomain2D::CSGCGALDomain2D(std::shared_ptr<const CSGGeometry> geometry,
+ double segment_granularity)
: impl(new CSGCGALDomain2DImpl)
{
if (geometry->dim() != 2)
@@ -362,11 +370,11 @@ CSGCGALDomain2D::CSGCGALDomain2D(
{
case CSGGeometry::Union:
{
- const CSGUnion *u = dynamic_cast<const CSGUnion*>(geometry.get());
+ std::shared_ptr<const CSGUnion> u = std::dynamic_pointer_cast<const CSGUnion, const CSGGeometry>(geometry);
dolfin_assert(u);
- CSGCGALDomain2D a(u->_g0);
- CSGCGALDomain2D b(u->_g1);
+ CSGCGALDomain2D a(u->_g0, segment_granularity);
+ CSGCGALDomain2D b(u->_g1, segment_granularity);
impl.swap(a.impl);
impl->polygon_set.join(b.impl->polygon_set);
@@ -374,11 +382,11 @@ CSGCGALDomain2D::CSGCGALDomain2D(
}
case CSGGeometry::Intersection:
{
- const CSGIntersection* u = dynamic_cast<const CSGIntersection*>(geometry.get());
+ auto u = std::dynamic_pointer_cast<const CSGIntersection, const CSGGeometry>(geometry);
dolfin_assert(u);
- CSGCGALDomain2D a(u->_g0);
- CSGCGALDomain2D b(u->_g1);
+ CSGCGALDomain2D a(u->_g0, segment_granularity);
+ CSGCGALDomain2D b(u->_g1, segment_granularity);
impl.swap(a.impl);
impl->polygon_set.intersection(b.impl->polygon_set);
@@ -386,10 +394,10 @@ CSGCGALDomain2D::CSGCGALDomain2D(
}
case CSGGeometry::Difference:
{
- const CSGDifference* u = dynamic_cast<const CSGDifference*>(geometry.get());
+ auto u = std::dynamic_pointer_cast<const CSGDifference, const CSGGeometry>(geometry);
dolfin_assert(u);
- CSGCGALDomain2D a(u->_g0);
- CSGCGALDomain2D b(u->_g1);
+ CSGCGALDomain2D a(u->_g0, segment_granularity);
+ CSGCGALDomain2D b(u->_g1, segment_granularity);
impl.swap(a.impl);
impl->polygon_set.difference(b.impl->polygon_set);
@@ -397,9 +405,9 @@ CSGCGALDomain2D::CSGCGALDomain2D(
}
case CSGGeometry::Translation :
{
- const CSGTranslation* t = dynamic_cast<const CSGTranslation*>(geometry.get());
+ auto t = std::dynamic_pointer_cast<const CSGTranslation, const CSGGeometry>(geometry);
dolfin_assert(t);
- CSGCGALDomain2D a(t->g);
+ CSGCGALDomain2D a(t->g, segment_granularity);
Exact_Kernel::Aff_transformation_2 translation(CGAL::TRANSLATION, Vector_2(t->t.x(), t->t.y()));
std::unique_ptr<CSGCGALDomain2DImpl> transformed = do_transformation(a.impl->polygon_set, translation);
impl.swap(transformed);
@@ -407,9 +415,9 @@ CSGCGALDomain2D::CSGCGALDomain2D(
}
case CSGGeometry::Scaling :
{
- const CSGScaling* t = dynamic_cast<const CSGScaling*>(geometry.get());
+ auto t = std::dynamic_pointer_cast<const CSGScaling, const CSGGeometry>(geometry);
dolfin_assert(t);
- CSGCGALDomain2D a(t->g);
+ CSGCGALDomain2D a(t->g, segment_granularity);
Exact_Kernel::Aff_transformation_2 tr(CGAL::IDENTITY);
// Translate if requested
@@ -431,9 +439,9 @@ CSGCGALDomain2D::CSGCGALDomain2D(
}
case CSGGeometry::Rotation :
{
- const CSGRotation* t = dynamic_cast<const CSGRotation*>(geometry.get());
+ auto t = std::dynamic_pointer_cast<const CSGRotation, const CSGGeometry>(geometry);
dolfin_assert(t);
- CSGCGALDomain2D a(t->g);
+ CSGCGALDomain2D a(t->g, segment_granularity);
Exact_Kernel::Aff_transformation_2 tr(CGAL::IDENTITY);
// Translate if requested
@@ -455,30 +463,30 @@ CSGCGALDomain2D::CSGCGALDomain2D(
}
case CSGGeometry::Circle:
{
- const Circle* c = dynamic_cast<const Circle*>(geometry.get());
+ auto c = std::dynamic_pointer_cast<const Circle, const CSGGeometry>(geometry);
dolfin_assert(c);
- impl->polygon_set.insert(make_circle(c, segment_granularity));
+ impl->polygon_set.insert(make_circle(c.get(), segment_granularity));
break;
}
case CSGGeometry::Ellipse:
{
- const Ellipse* c = dynamic_cast<const Ellipse*>(geometry.get());
+ auto c = std::dynamic_pointer_cast<const Ellipse, const CSGGeometry>(geometry);
dolfin_assert(c);
- impl->polygon_set.insert(make_ellipse(c, segment_granularity));
+ impl->polygon_set.insert(make_ellipse(c.get(), segment_granularity));
break;
}
case CSGGeometry::Rectangle:
{
- const Rectangle* r = dynamic_cast<const Rectangle*>(geometry.get());
+ auto r = std::dynamic_pointer_cast<const Rectangle, const CSGGeometry>(geometry);
dolfin_assert(r);
- impl->polygon_set.insert(make_rectangle(r));
+ impl->polygon_set.insert(make_rectangle(r.get()));
break;
}
case CSGGeometry::Polygon:
{
- const Polygon* p = dynamic_cast<const Polygon*>(geometry.get());
+ auto p = std::dynamic_pointer_cast<const Polygon, const CSGGeometry>(geometry);
dolfin_assert(p);
- impl->polygon_set.insert(make_polygon(p));
+ impl->polygon_set.insert(make_polygon(p.get()));
break;
}
default:
diff --git a/src/CSGCGALDomain3D.cpp b/src/CSGCGALDomain3D.cpp
index 3ef7742..9ac4f26 100644
--- a/src/CSGCGALDomain3D.cpp
+++ b/src/CSGCGALDomain3D.cpp
@@ -843,9 +843,13 @@ namespace
generator.parameters["mesh_resolution"] = 2.0;
generator.parameters["partition"] = false;
+ const std::pair<dolfin::Point, double> bounding_circle = polygon.estimate_bounding_sphere();
+
// Generate 2D mesh (on all nodes if in parallel)
std::shared_ptr<mshr::CSGCGALDomain2D>
- d2(new mshr::CSGCGALDomain2D(dolfin::reference_to_no_delete_pointer(polygon)));
+ d2(new mshr::CSGCGALDomain2D(dolfin::reference_to_no_delete_pointer(polygon),
+ bounding_circle.second/6.));
+
mesh2d = generator.generate(d2);
}
@@ -1506,7 +1510,7 @@ namespace
parameters = default_parameters();
// Create the polyhedron
- BuildFromFacetList<Exact_HalfedgeDS> builder(vertices, facets, {});
+ BuildFromFacetList<Exact_HalfedgeDS> builder(vertices, facets, std::set<std::size_t>{});
impl->p.delegate(builder);
}
//-----------------------------------------------------------------------------
diff --git a/src/CSGCGALMeshGenerator2D.cpp b/src/CSGCGALMeshGenerator2D.cpp
index 6716966..c30f569 100644
--- a/src/CSGCGALMeshGenerator2D.cpp
+++ b/src/CSGCGALMeshGenerator2D.cpp
@@ -1,4 +1,4 @@
-// Copyright (C) 2012 Johannes Ring, 2012-2016 Benjamin Kehlet
+// Copyright (C) 2012 Johannes Ring, 2012-2017 Benjamin Kehlet
//
// This file is part of mshr.
//
@@ -217,17 +217,15 @@ std::shared_ptr<dolfin::Mesh> CSGCGALMeshGenerator2D::generate(const std::shared
// Empty polygon, will be populated when traversing the subdomains
CSGCGALDomain2D overlaying;
- // Add the subdomains to the PSLG. Traverse in reverse order to get the latest
- // added subdomain on top
- std::list<std::pair<std::size_t, std::shared_ptr<const CSGGeometry> > >::const_reverse_iterator it;
-
if (!subdomains.empty())
log(dolfin::TRACE, "Processing subdomains");
- for (const std::pair<std::size_t, std::shared_ptr<const CSGCGALDomain2D>>& current_subdomain : subdomains)
+ // Add the subdomains to the PSLG. Traverse in reverse order to get the latest
+ // added subdomain on top
+ for (auto rit = subdomains.rbegin(); rit != subdomains.rend(); rit++)
{
- const std::size_t current_index = current_subdomain.first;
- CSGCGALDomain2D cgal_geometry = *current_subdomain.second;
+ const std::size_t current_index = rit->first;
+ CSGCGALDomain2D cgal_geometry(*(rit->second));
// Only the part inside the total domain
cgal_geometry.intersect_inplace(*total_domain, 1e-15);
diff --git a/src/CSGGeometry.cpp b/src/CSGGeometry.cpp
index 0f521cf..17d522b 100644
--- a/src/CSGGeometry.cpp
+++ b/src/CSGGeometry.cpp
@@ -1,4 +1,4 @@
-// Copyright (C) 2012 Anders Logg
+// Copyright (C) 2012 Anders Logg, Benjamin Kehlet 2013-2017
//
// This file is part of mshr.
//
@@ -6,21 +6,25 @@
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
-//
+//
// mshr is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
-//
+//
// You should have received a copy of the GNU General Public License
// along with mshr. If not, see <http://www.gnu.org/licenses/>.
//
-// Modified by Benjamin Kehlet, 2013
-
#include <dolfin/common/NoDeleter.h>
+#include <dolfin/log/LogStream.h>
#include <mshr/CSGGeometry.h>
+// Bounding sphere computation
+#include <CGAL/Cartesian.h>
+#include <CGAL/Min_sphere_of_spheres_d.h>
+#include <CGAL/Min_sphere_of_spheres_d_traits_3.h>
+
namespace mshr
{
@@ -81,5 +85,53 @@ bool CSGGeometry::has_subdomains() const
{
return subdomains.size() > 0;
}
+//-----------------------------------------------------------------------------
+bool CSGGeometry::inside(dolfin::Point p1, dolfin::Point p2) const
+{
+ dolfin_not_implemented();
+ return false;
+}
+//-----------------------------------------------------------------------------
+std::pair<dolfin::Point, double> CSGGeometry::estimate_bounding_sphere(std::size_t numSamples) const
+{
+ std::pair<dolfin::Point, dolfin::Point> aabb = bounding_box();
+ const dolfin::Point min = aabb.first;
+ const dolfin::Point max = aabb.second;
+
+ //typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+ typedef CGAL::Cartesian<double> K;
+ typedef K::Point_2 Point_2;
+ typedef CGAL::Min_sphere_of_spheres_d_traits_2<K, K::FT> MinSphereTraits;
+ typedef CGAL::Min_sphere_of_spheres_d<MinSphereTraits> Min_sphere;
+ typedef MinSphereTraits::Sphere Sphere;
+
+ std::vector<Sphere> S;
+ S.reserve(numSamples);
+ const std::size_t N = static_cast<std::size_t>(std::sqrt(static_cast<double>(numSamples)+.5));
+ const double dX = (max.x()-min.x())/N;
+ const double dY = (max.y()-min.y())/N;
+ // const double dTheta = 2*DOLFIN_PI/N;
+ for (std::size_t i = 0; i < N; i++)
+ {
+ for (std::size_t j = 0; j < N; j++)
+ {
+ const dolfin::Point p(min.x() + i*dX, min.y() + j*dY);
+ if (inside(p))
+ S.push_back(Sphere(Point_2(p.x(), p.y()), 0.0));
+ }
+ }
+
+ Min_sphere ms(S.begin(), S.end());
+ CGAL_assertion(ms.is_valid());
+
+ auto it = ms.center_cartesian_begin();
+ const double center_x = CGAL::to_double(*it);
+ it++;
+ const double center_y = CGAL::to_double(*it);
+ const double radius = CGAL::to_double(ms.radius());
+
+ return std::make_pair(dolfin::Point(center_x, center_y),
+ radius);
+}
}
diff --git a/src/CSGOperators.cpp b/src/CSGOperators.cpp
index e9aa69f..a01e48e 100644
--- a/src/CSGOperators.cpp
+++ b/src/CSGOperators.cpp
@@ -1,4 +1,4 @@
-// Copyright (C) 2012 Anders Logg
+// Copyright (C) 2012 Anders Logg, Benjamin Kehlet 2013-2017
//
// This file is part of mshr.
//
@@ -16,7 +16,6 @@
// along with mshr. If not, see <http://www.gnu.org/licenses/>.
//
// Modified by Johannes Ring, 2012
-// Modified by Benjamin Kehlet, 2013
#include <mshr/CSGOperators.h>
@@ -81,6 +80,25 @@ std::string CSGUnion::str(bool verbose) const
return s.str();
}
//-----------------------------------------------------------------------------
+std::pair<dolfin::Point, dolfin::Point> CSGUnion::bounding_box() const
+{
+ std::pair<dolfin::Point, dolfin::Point> a = _g0->bounding_box();
+ std::pair<dolfin::Point, dolfin::Point> b = _g1->bounding_box();
+
+ return std::make_pair(dolfin::Point(std::min(a.first.x(), b.first.x()),
+ std::min(a.first.y(), b.first.y()),
+ std::min(a.first.z(), b.first.z())),
+ dolfin::Point(std::max(a.second.x(), b.second.x()),
+ std::max(a.second.y(), b.second.y()),
+ std::max(a.second.z(), b.second.z())));
+}
+//-----------------------------------------------------------------------------
+bool CSGUnion::inside(dolfin::Point p) const
+{
+ return _g0->inside(p) || _g1->inside(p);
+}
+
+//-----------------------------------------------------------------------------
// CSGDifference
//-----------------------------------------------------------------------------
CSGDifference::CSGDifference(std::shared_ptr<CSGGeometry> g0,
@@ -126,6 +144,18 @@ std::string CSGDifference::str(bool verbose) const
return s.str();
}
//-----------------------------------------------------------------------------
+std::pair<dolfin::Point, dolfin::Point> CSGDifference::bounding_box() const
+{
+ // FIXME: This is where the bounding_box implementation may overshoot
+ return _g0->bounding_box();
+}
+//-----------------------------------------------------------------------------
+bool CSGDifference::inside(dolfin::Point p) const
+{
+ return _g0->inside(p) && !_g1->inside(p);
+}
+
+//-----------------------------------------------------------------------------
// CSGIntersection
//-----------------------------------------------------------------------------
CSGIntersection::CSGIntersection(std::shared_ptr<CSGGeometry> g0,
@@ -171,6 +201,25 @@ std::string CSGIntersection::str(bool verbose) const
return s.str();
}
//-----------------------------------------------------------------------------
+std::pair<dolfin::Point, dolfin::Point> CSGIntersection::bounding_box() const
+{
+ std::pair<dolfin::Point, dolfin::Point> a = _g0->bounding_box();
+ std::pair<dolfin::Point, dolfin::Point> b = _g1->bounding_box();
+
+ return std::make_pair(dolfin::Point(std::min(a.first.x(), b.first.x()),
+ std::min(a.first.y(), b.first.y()),
+ std::min(a.first.z(), b.first.z())),
+ dolfin::Point(std::max(a.second.x(), b.second.x()),
+ std::max(a.second.y(), b.second.y()),
+ std::max(a.second.z(), b.second.z())));
+}
+//-----------------------------------------------------------------------------
+bool CSGIntersection::inside(dolfin::Point p) const
+{
+ return _g0->inside(p) && _g1->inside(p);
+}
+
+//-----------------------------------------------------------------------------
// CSGTranslation
//-----------------------------------------------------------------------------
CSGTranslation::CSGTranslation(std::shared_ptr<CSGGeometry> g,
@@ -201,6 +250,19 @@ std::string CSGTranslation::str(bool verbose) const
return s.str();
}
//-----------------------------------------------------------------------------
+std::pair<dolfin::Point, dolfin::Point> CSGTranslation::bounding_box() const
+{
+ std::pair<dolfin::Point, dolfin::Point> aabb = g->bounding_box();
+
+ return std::make_pair(aabb.first+t, aabb.second+t);
+}
+//-----------------------------------------------------------------------------
+bool CSGTranslation::inside(dolfin::Point p) const
+{
+ return g->inside(p-t);
+}
+
+//-----------------------------------------------------------------------------
// CSGScaling
//-----------------------------------------------------------------------------
CSGScaling::CSGScaling(std::shared_ptr<CSGGeometry> g,
@@ -248,6 +310,26 @@ std::string CSGScaling::str(bool verbose) const
return ss.str();
}
//-----------------------------------------------------------------------------
+std::pair<dolfin::Point, dolfin::Point> CSGScaling::bounding_box() const
+{
+ std::pair<dolfin::Point, dolfin::Point> aabb = g->bounding_box();
+
+ const dolfin::Point scaled_c = (translate ? (1-s)*c : dolfin::Point(0,0,0));
+
+ if (translate)
+ return std::make_pair((aabb.first-c)*s + c, (aabb.second-c)*s + c);
+ else
+ return std::make_pair(aabb.first*s, aabb.second*s);
+
+}
+//-----------------------------------------------------------------------------
+bool CSGScaling::inside(dolfin::Point p) const
+{
+ const dolfin::Point p_scaled = (translate ? (p-c)*s + c : p*s);
+ return g->inside(p_scaled);
+}
+
+//-----------------------------------------------------------------------------
// CSGRotation
//-----------------------------------------------------------------------------
CSGRotation::CSGRotation(std::shared_ptr<CSGGeometry> g,
@@ -330,4 +412,33 @@ std::string CSGRotation::str(bool verbose) const
return ss.str();
}
//-----------------------------------------------------------------------------
+std::pair<dolfin::Point, dolfin::Point> CSGRotation::bounding_box() const
+{
+ if (dim() != 2)
+ {
+ dolfin_not_implemented();
+ }
+
+ const std::pair<dolfin::Point, dolfin::Point> aabb = g->bounding_box();
+ const dolfin::Point axis(0,0,1);
+
+ // Rotate all corners
+ const dolfin::Point a = aabb.first.rotate(axis, theta);
+ const dolfin::Point b = dolfin::Point(aabb.first.x(), aabb.second.y()).rotate(axis, theta);
+ const dolfin::Point c = dolfin::Point(aabb.second.x(), aabb.first.y()).rotate(axis, theta);
+ const dolfin::Point d = aabb.second.rotate(axis, theta);
+
+ return std::make_pair(dolfin::Point(std::min(a.x(), std::min(b.x(), std::min(c.x(), d.x()))),
+ std::min(a.y(), std::min(b.y(), std::min(c.y(), d.y())))),
+ dolfin::Point(std::max(a.x(), std::max(b.x(), std::max(c.x(), d.x()))),
+ std::max(a.y(), std::max(b.y(), std::max(c.y(), d.y())))));
+
+}
+//-----------------------------------------------------------------------------
+bool CSGRotation::inside(dolfin::Point p) const
+{
+ const dolfin::Point p_rotated = translate ? (p-c).rotate(dolfin::Point(0,0,1), -theta) + c : p.rotate(dolfin::Point(0,0,1), -theta);
+ return g->inside(p_rotated);
+}
+
}
diff --git a/src/CSGPrimitives2D.cpp b/src/CSGPrimitives2D.cpp
index 110d3b8..05fa47e 100644
--- a/src/CSGPrimitives2D.cpp
+++ b/src/CSGPrimitives2D.cpp
@@ -1,4 +1,4 @@
-// Copyright (C) 2012 Anders Logg, 2014-2015 Benjamin Kehlet
+// Copyright (C) 2012 Anders Logg, 2014-2017 Benjamin Kehlet
//
// This file is part of mshr.
//
@@ -16,8 +16,6 @@
// along with mshr. If not, see <http://www.gnu.org/licenses/>.
//
// Modified by Johannes Ring, 2012
-// Modified by Benjamin Kehlet, 2012-2013
-
#include <mshr/CSGPrimitives2D.h>
#include <dolfin/math/basic.h>
@@ -27,6 +25,8 @@
#include <limits>
#include <algorithm>
+#include <CGAL/Cartesian.h>
+#include <CGAL/Polygon_2.h>
namespace mshr
{
@@ -75,6 +75,20 @@ std::string Circle::str(bool verbose) const
return s.str();
}
//-----------------------------------------------------------------------------
+std::pair<dolfin::Point, dolfin::Point> Circle::bounding_box() const
+{
+ return std::make_pair(dolfin::Point(c.x()-radius(), c.y()-radius()),
+ dolfin::Point(c.x()+radius(), c.y()+radius()));
+}
+//-----------------------------------------------------------------------------
+bool Circle::inside(dolfin::Point p) const
+{
+ // dolfin::cout << "Inside: " << c << " r=" << _r << " ::: " << p << "(" << ((p-c).squared_norm() <= _r*_r ? "Inside" : "Outside") << dolfin::endl;
+ return (p-c).squared_norm() <= _r*_r;
+}
+
+
+//-----------------------------------------------------------------------------
// Ellipse
//-----------------------------------------------------------------------------
Ellipse::Ellipse(dolfin::Point c, double a, double b,
@@ -115,6 +129,20 @@ std::string Ellipse::str(bool verbose) const
return s.str();
}
//-----------------------------------------------------------------------------
+std::pair<dolfin::Point, dolfin::Point> Ellipse::bounding_box() const
+{
+ return std::make_pair(dolfin::Point(c.x()-a(), c.y()-b()),
+ dolfin::Point(c.x()+a(), c.y()+b()));
+}
+//-----------------------------------------------------------------------------
+bool Ellipse::inside(dolfin::Point p) const
+{
+ const double x = (c.x()-p.x())/_a;
+ const double y = (c.y()-p.y())/_b;
+ return x*x + y*y <= 1;
+}
+
+//-----------------------------------------------------------------------------
// Rectangle
//-----------------------------------------------------------------------------
Rectangle::Rectangle(dolfin::Point a, dolfin::Point b)
@@ -147,6 +175,18 @@ std::string Rectangle::str(bool verbose) const
return s.str();
}
//-----------------------------------------------------------------------------
+std::pair<dolfin::Point, dolfin::Point> Rectangle::bounding_box() const
+{
+ return std::make_pair(a, b);
+}
+//-----------------------------------------------------------------------------
+bool Rectangle::inside(dolfin::Point p) const
+{
+ return p.x() >= a.x() && p.x() <= b.x() && p.y() >= a.y() && p.y() <= b.y();
+}
+
+
+//-----------------------------------------------------------------------------
// Polygon
//-----------------------------------------------------------------------------
Polygon::Polygon(const std::vector<dolfin::Point>& vertices)
@@ -213,4 +253,40 @@ bool Polygon::ccw() const
return signed_area > 0;
}
+//-----------------------------------------------------------------------------
+std::pair<dolfin::Point, dolfin::Point> Polygon::bounding_box() const
+{
+ std::vector<dolfin::Point>::const_iterator it = _vertices.begin();
+ dolfin::Point min = *it;
+ dolfin::Point max = *it;
+ it++;
+
+ for (; it != _vertices.end(); it++)
+ {
+ const dolfin::Point& p = *it;
+ min[0] = std::min(min[0], p[0]);
+ min[1] = std::min(min[1], p[1]);
+ min[2] = std::min(min[2], p[2]);
+
+ max[0] = std::max(max[0], p[0]);
+ max[1] = std::max(max[1], p[1]);
+ max[2] = std::max(max[2], p[2]);
+ }
+
+ return std::make_pair(min, max);
+}
+//-----------------------------------------------------------------------------
+bool Polygon::inside(dolfin::Point p) const
+{
+ typedef CGAL::Cartesian<double> K;
+ typedef CGAL::Point_2<K> Point_2;
+ typedef CGAL::Polygon_2<K> Polygon_2;
+
+ Polygon_2 polygon;
+ for (const dolfin::Point& vertex : _vertices)
+ polygon.push_back(Point_2(vertex.x(), vertex.y()));
+
+ return polygon.has_on_bounded_side(Point_2(p.x(), p.y()));
+}
+
}
diff --git a/src/CSGPrimitives3D.cpp b/src/CSGPrimitives3D.cpp
index 02fd67c..a83fbf4 100644
--- a/src/CSGPrimitives3D.cpp
+++ b/src/CSGPrimitives3D.cpp
@@ -1,4 +1,4 @@
-// Copyright (C) 2012 Anders Logg and 2012, 2014-2015 Benjamin Kehlet
+// Copyright (C) 2012 Anders Logg and 2012, 2014-2017 Benjamin Kehlet
//
// This file is part of mshr.
//
@@ -28,6 +28,18 @@ namespace mshr
CSGPrimitive3D::CSGPrimitive3D()
{}
//-----------------------------------------------------------------------------
+std::pair<dolfin::Point, dolfin::Point> CSGPrimitive3D::bounding_box() const
+{
+ dolfin_not_implemented();
+ return std::make_pair(dolfin::Point(0., 0., 0.), dolfin::Point(0., 0., 0.));
+}
+//-----------------------------------------------------------------------------
+bool CSGPrimitive3D::inside(dolfin::Point p) const
+{
+ dolfin_not_implemented();
+ return false;
+}
+
//-----------------------------------------------------------------------------
diff --git a/src/MeshGenerator.cpp b/src/MeshGenerator.cpp
index 556fa59..8a6d896 100644
--- a/src/MeshGenerator.cpp
+++ b/src/MeshGenerator.cpp
@@ -55,36 +55,22 @@ std::shared_ptr<dolfin::Mesh> generate_mesh(std::shared_ptr<const CSGGeometry> g
CSGCGALMeshGenerator2D generator;
- // Compute cell size
- // Chicken-and-egg problem:
- // The cell size is needed to transform the circles into polygons, but
- // the cell size can only be determined via the bounding circle radius
- // which can only be computed if we already have a polygon.
- //
- // Proper solution:
- // Compute the bounding circle radius from the geometry information
- // alone.
- //
- // Workaround:
- // Polygonize with some arbitrary cell size, compute the bounding circle
- // of the resulting polygonal object, and take the cell size from that.
- double cell_size;
- {
- const double tmp_cell_size = 0.1;
- CSGCGALDomain2D total_domain_coarse(geometry, tmp_cell_size);
- const double min_radius = total_domain_coarse.compute_boundingcircle_radius();
- cell_size = 2.0*min_radius/resolution;
- }
+ // Estimate the bounding circle radius from the geometry
+ const std::pair<dolfin::Point, double> bounding_circle = geometry->estimate_bounding_sphere();
+ const double segment_granularity = 2*bounding_circle.second/resolution;
+ const double cell_size = 2.0*bounding_circle.second/resolution;
+
+ // Calculate surface representation
log(dolfin::TRACE, "Request cell size: %f", cell_size);
generator.parameters["mesh_resolution"] = -1.0;
generator.parameters["cell_size"] = cell_size;
- std::shared_ptr<CSGCGALDomain2D> total_domain(new CSGCGALDomain2D(geometry));
+ std::shared_ptr<CSGCGALDomain2D> total_domain(new CSGCGALDomain2D(geometry, segment_granularity));
std::vector<std::pair<std::size_t, std::shared_ptr<const CSGCGALDomain2D>>> subdomain_geometries;
for (const std::pair<std::size_t, std::shared_ptr<const CSGGeometry>>& subdomain : geometry->get_subdomains())
{
subdomain_geometries.push_back(std::make_pair(subdomain.first,
- std::shared_ptr<CSGCGALDomain2D>(new CSGCGALDomain2D(subdomain.second))));
+ std::shared_ptr<CSGCGALDomain2D>(new CSGCGALDomain2D(subdomain.second, segment_granularity))));
}
return generator.generate(total_domain, subdomain_geometries);
}
diff --git a/swig/CMakeLists.txt b/swig/CMakeLists.txt
index 832318a..2e4753e 100644
--- a/swig/CMakeLists.txt
+++ b/swig/CMakeLists.txt
@@ -10,9 +10,12 @@ endif()
find_package(SWIG REQUIRED)
include(${SWIG_USE_FILE})
-# FIXME: Instead use DOLFIN_PYTHON_EXECUTABLE when
-# this has been added to DOLFINConfig.cmake
-find_package(PythonInterp 2)
+if (NOT DEFINED DOLFIN_PYTHON_EXECUTABLE )
+ message(WARNING "DOLFIN_PYTHON_EXECUTABLE not defined. This should be in DOLFINConfig.cmake")
+ find_package(PythonInterp)
+else()
+ set(PYTHON_EXECUTABLE ${DOLFIN_PYTHON_EXECUTABLE})
+endif()
# Don't use these. Instad get the values from Dolfin
# find_package(PythonLibs)
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 633a021..20834d0 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -66,7 +66,17 @@ set_property(TEST "Python-ASCFileReader" PROPERTY
ENVIRONMENT "PYTHONPATH=${PYTHON_ENVIR}"
)
+############# Test the ASCFileReader ##########################################################
+add_test("Python-NumSegments2D" "${PYTHON_EXECUTABLE}" "${CMAKE_SOURCE_DIR}/test/test-num-segments-2d.py")
+set_property(TEST "Python-NumSegments2D" PROPERTY
+ ENVIRONMENT "PYTHONPATH=${PYTHON_ENVIR}"
+)
+############# Test the ASCFileReader ##########################################################
+add_test("Python-CSGPredicates2D" "${PYTHON_EXECUTABLE}" "${CMAKE_SOURCE_DIR}/test/test-csg-predicates.py")
+set_property(TEST "Python-CSGPredicates2D" PROPERTY
+ ENVIRONMENT "PYTHONPATH=${PYTHON_ENVIR}"
+)
############# Try meshing some surface files found in the cgal source ########################
# TODO: Print mesh quality measure from program and check the value
set(CGAL_DATA_DIR "${CMAKE_SOURCE_DIR}/3rdparty/CGAL/demo/Polyhedron/data")
diff --git a/test/test-csg-predicates.py b/test/test-csg-predicates.py
new file mode 100644
index 0000000..e1cee1a
--- /dev/null
+++ b/test/test-csg-predicates.py
@@ -0,0 +1,34 @@
+from dolfin import *
+import mshr
+import math
+
+
+r = mshr.Rectangle(Point(0,-.5), Point(4,.5))
+
+# Trivial inside outside
+assert r.inside(Point( .5, .0))
+assert not r.inside(Point(1.5, 1.5))
+
+rotated = mshr.CSGRotation(r, math.pi/2)
+assert not rotated.inside(Point(1.5, 1.5))
+assert rotated.inside(Point(0., 1.5))
+assert not rotated.inside(Point(0., 4.5))
+assert not rotated.inside(Point(0., -.5))
+
+# Translation
+translated = mshr.CSGTranslation(rotated, Point(2, 1))
+assert translated.inside(Point(2, 1.5))
+assert not translated.inside(Point(0, 1.5))
+
+# Circle
+c = mshr.Circle(Point(2, 1), .5)
+assert c.inside(Point(2, 1))
+assert not c.inside(Point(3, 1))
+
+c_rotated = mshr.CSGRotation(c, Point(2,1), math.pi)
+assert c_rotated.inside(Point(2,1))
+assert not c_rotated.inside(Point(3, 1))
+
+c_rotated_2 = mshr.CSGRotation(c_rotated, pi/2)
+assert not c_rotated_2.inside(Point(2,1))
+assert c_rotated_2.inside(Point(-1,2))
diff --git a/test/test-num-segments-2d.py b/test/test-num-segments-2d.py
new file mode 100644
index 0000000..e1c2273
--- /dev/null
+++ b/test/test-num-segments-2d.py
@@ -0,0 +1,13 @@
+from dolfin import *
+import mshr
+
+# The ellipse is completely contained in the circle and will not be
+# visible in the resulting mesh, but it is so small that with
+# mesh_resolution = 10 that when converting it to a polygon the number
+# of segments in the polygon will be 0 (ie. the polygon is degenerate)
+# if not handled correctly.
+c = mshr.Circle(Point(0, 0), 1.5)
+e = mshr.Ellipse(Point(.05, 0), .01, .05)
+
+mesh = mshr.generate_mesh(c+e, 10)
+
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/fenics/mshr.git
More information about the debian-science-commits
mailing list