fenics-buildbot team mailing list archive
-
fenics-buildbot team
-
Mailing list archive
-
Message #02037
buildbot failure in FEniCS Project on dolfin-master-full-raring-amd64
The Buildbot has detected a new failure on builder dolfin-master-full-raring-amd64 while building FEniCS Project.
Full details are available at:
http://fenicsproject.org:8010/builders/dolfin-master-full-raring-amd64/builds/286
Buildbot URL: http://fenicsproject.org:8010/
Buildslave for this Build: raring-amd64
Build Reason: 'try' job by user Johannes Ring
Build Source Stamp: HEAD (plus patch)
Blamelist: Johannes Ring
BUILD FAILED: failed make run_regressiontests
sincerely,
-The Buildbot
diff --git a/demo/undocumented/csg/2D/cpp/main.cpp b/demo/undocumented/csg/2D/cpp/main.cpp
index ed6ee3e..3f185e8 100644
--- a/demo/undocumented/csg/2D/cpp/main.cpp
+++ b/demo/undocumented/csg/2D/cpp/main.cpp
@@ -31,24 +31,43 @@ using namespace dolfin;
int main()
{
// Define 2D geometry
- Rectangle r(0.5, 0.5, 1.5, 1.5);
- Circle c(1, 1, 1);
- boost::shared_ptr<CSGGeometry> g2d = c - r;
+ // Rectangle r(0.5, 0.5, 1.5, 1.5);
+ // Circle c(1, 1, 1);
+ // boost::shared_ptr<CSGGeometry> g2d = c - r;
+
+ // Define 2D geometry
+ Rectangle r1(0., 0., 5., 5.);
+ Rectangle r2 (2., 1.25, 3., 1.75);
+ Circle c1(1, 4, .25);
+ Circle c2(4, 4, .25);
+ boost::shared_ptr<CSGGeometry> domain = r1 - r2 - c1 - c2;
+
+ Rectangle s1(1., 1., 4., 3.);
+ domain->set_subdomain(1, s1);
+
+ Rectangle s2(2., 2., 3., 4.);
+ domain->set_subdomain(2, s2);
+
// Test printing
info("\nCompact output of 2D geometry:");
- info(*g2d);
+ info(*domain);
info("");
info("\nVerbose output of 2D geometry:");
- info(*g2d, true);
+ info(*domain, true);
// Plot geometry
- plot(g2d, "2D Geometry (boundary)");
+ plot(domain, "2D Geometry (boundary)");
// Generate and plot mesh
- Mesh mesh2d(g2d, 100);
+ boost::shared_ptr<Mesh> mesh2d(new Mesh(domain, 45));
plot(mesh2d, "2D mesh");
+ // Convert mesh domains to mesh function for plotting
+ MeshFunction<std::size_t> mf(mesh2d, 2, mesh2d->domains());
+ plot(mf, "Subdomains");
+
+
interactive();
return 0;
}
diff --git a/demo/undocumented/csg/2D/python/demo_csg_2D.py b/demo/undocumented/csg/2D/python/demo_csg_2D.py
index 3e5ee3f..d89704d 100644
--- a/demo/undocumented/csg/2D/python/demo_csg_2D.py
+++ b/demo/undocumented/csg/2D/python/demo_csg_2D.py
@@ -26,22 +26,28 @@ if not has_cgal():
# Define 2D geometry
-r = Rectangle(0.5, 0.5, 1.5, 1.5);
-c = Circle (1, 1, 1);
-g2d = c - r;
+domain = Rectangle(0., 0., 5., 5.) - Rectangle(2., 1.25, 3., 1.75) - Circle(1, 4, .25) - Circle(4, 4, .25)
+domain.set_subdomain(1, Rectangle(1., 1., 4., 3.))
+domain.set_subdomain(2, Rectangle(2., 2., 3., 4.))
+
# Test printing
info("\nCompact output of 2D geometry:");
-info(g2d);
+info(domain);
info("");
info("\nVerbose output of 2D geometry:");
-info(g2d, True);
+info(domain, True);
# Plot geometry
-plot(g2d, "2D Geometry (boundary)");
+plot(domain, "2D Geometry (boundary)");
# Generate and plot mesh
-mesh2d = Mesh(g2d, 10);
-plot(mesh2d, "2D mesh");
+mesh2d = Mesh(domain, 45);
+plot(mesh2d, "2D mesh")
+
+# Convert subdomains to mesh function for plotting
+mf = MeshFunction("size_t", mesh2d, 2, mesh2d.domains())
+plot(mf, "Subdomains")
+
-interactive();
+interactive()
diff --git a/dolfin/generation/CSGCGALDomain2D.cpp b/dolfin/generation/CSGCGALDomain2D.cpp
new file mode 100644
index 0000000..8582947
--- /dev/null
+++ b/dolfin/generation/CSGCGALDomain2D.cpp
@@ -0,0 +1,345 @@
+// Copyright (C) 2013 Benjamin Kehlet
+//
+// This file is part of DOLFIN.
+//
+// DOLFIN is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// DOLFIN 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 Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
+//
+// First added: 2013-06-22
+// Last changed: 2013-08-06
+
+#include "CSGCGALDomain2D.h"
+
+#ifdef HAS_CGAL
+
+#include "CSGPrimitives2D.h"
+#include "CSGOperators.h"
+#include <dolfin/common/constants.h>
+#include <dolfin/log/LogStream.h>
+
+#include <CGAL/basic.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Boolean_set_operations_2.h>
+#include <CGAL/Polygon_set_2.h>
+
+#include <CGAL/Min_circle_2.h>
+#include <CGAL/Min_circle_2_traits_2.h>
+
+// Polygon typedefs
+typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_Kernel;
+typedef Exact_Kernel::Point_2 Point_2;
+typedef CGAL::Polygon_2<Exact_Kernel> Polygon_2;
+typedef Polygon_2::Vertex_const_iterator Vertex_const_iterator;
+typedef CGAL::Polygon_with_holes_2<Exact_Kernel> Polygon_with_holes_2;
+typedef Polygon_with_holes_2::Hole_const_iterator Hole_const_iterator;
+typedef CGAL::Polygon_set_2<Exact_Kernel> Polygon_set_2;
+
+// Min enclosing circle typedefs
+typedef CGAL::Min_circle_2_traits_2<Exact_Kernel> Min_Circle_Traits;
+typedef CGAL::Min_circle_2<Min_Circle_Traits> Min_circle;
+typedef CGAL::Circle_2<Exact_Kernel> CGAL_Circle;
+
+using namespace dolfin;
+
+struct CSGCGALDomain2DImpl
+{
+ Polygon_set_2 polygon_set;
+
+ CSGCGALDomain2DImpl(){}
+ CSGCGALDomain2DImpl(const Polygon_set_2& p)
+ : polygon_set(p) {}
+};
+//-----------------------------------------------------------------------------
+Polygon_2 make_circle(const Circle* c)
+{
+ std::vector<Point_2> pts;
+ pts.reserve(c->fragments());
+
+ for (std::size_t i = 0; i < c->fragments(); i++)
+ {
+ const double phi = (2*DOLFIN_PI*i) / c->fragments();
+ const double x = c->center().x() + c->radius()*cos(phi);
+ const double y = c->center().y() + c->radius()*sin(phi);
+ pts.push_back(Point_2(x, y));
+ }
+
+ return Polygon_2(pts.begin(), pts.end());
+}
+//-----------------------------------------------------------------------------
+Polygon_2 make_ellipse(const Ellipse* e)
+{
+ std::vector<Point_2> pts;
+
+ for (std::size_t i = 0; i < e->fragments(); i++)
+ {
+ const double phi = (2*DOLFIN_PI*i) / e->fragments();
+ const double x = e->center().x() + e->a()*cos(phi);
+ const double y = e->center().y() + e->b()*sin(phi);
+ pts.push_back(Point_2(x, y));
+ }
+
+ return Polygon_2(pts.begin(), pts.end());
+}
+//-----------------------------------------------------------------------------
+Polygon_2 make_rectangle(const Rectangle* r)
+{
+ const double x0 = std::min(r->first_corner().x(), r->second_corner().x());
+ const double y0 = std::min(r->first_corner().y(), r->second_corner().y());
+
+ const double x1 = std::max(r->first_corner().x(), r->second_corner().x());
+ const double y1 = std::max(r->first_corner().y(), r->second_corner().y());
+
+ std::vector<Point_2> pts;
+ pts.push_back(Point_2(x0, y0));
+ pts.push_back(Point_2(x1, y0));
+ pts.push_back(Point_2(x1, y1));
+ pts.push_back(Point_2(x0, y1));
+
+ Polygon_2 p(pts.begin(), pts.end());
+
+ return p;
+}
+//-----------------------------------------------------------------------------
+Polygon_2 make_polygon(const Polygon* p)
+{
+ std::vector<Point_2> pts;
+ std::vector<Point>::const_iterator v;
+ for (v = p->vertices().begin(); v != p->vertices().end(); ++v)
+ pts.push_back(Point_2(v->x(), v->y()));
+
+ return Polygon_2(pts.begin(), pts.end());
+}
+//-----------------------------------------------------------------------------
+CSGCGALDomain2D::CSGCGALDomain2D()
+ : impl(new CSGCGALDomain2DImpl)
+{
+
+}
+//-----------------------------------------------------------------------------
+CSGCGALDomain2D::~CSGCGALDomain2D()
+{
+}
+//-----------------------------------------------------------------------------
+CSGCGALDomain2D::CSGCGALDomain2D(const CSGGeometry *geometry)
+: impl(new CSGCGALDomain2DImpl)
+{
+ switch (geometry->getType())
+ {
+ case CSGGeometry::Union:
+ {
+ const CSGUnion *u = dynamic_cast<const CSGUnion*>(geometry);
+ dolfin_assert(u);
+
+ CSGCGALDomain2D a(u->_g0.get());
+ CSGCGALDomain2D b(u->_g1.get());
+
+ impl.swap(a.impl);
+ impl->polygon_set.join(b.impl->polygon_set);
+ break;
+ }
+ case CSGGeometry::Intersection:
+ {
+ const CSGIntersection* u = dynamic_cast<const CSGIntersection*>(geometry);
+ dolfin_assert(u);
+
+ CSGCGALDomain2D a(u->_g0.get());
+ CSGCGALDomain2D b(u->_g1.get());
+
+ impl.swap(a.impl);
+ impl->polygon_set.intersection(b.impl->polygon_set);
+ break;
+ }
+ case CSGGeometry::Difference:
+ {
+ const CSGDifference* u = dynamic_cast<const CSGDifference*>(geometry);
+ dolfin_assert(u);
+ CSGCGALDomain2D a(u->_g0.get());
+ CSGCGALDomain2D b(u->_g1.get());
+
+ impl.swap(a.impl);
+ impl->polygon_set.difference(b.impl->polygon_set);
+ break;
+ }
+ case CSGGeometry::Circle:
+ {
+ const Circle* c = dynamic_cast<const Circle*>(geometry);
+ dolfin_assert(c);
+ impl->polygon_set.insert(make_circle(c));
+ break;
+ }
+ case CSGGeometry::Ellipse:
+ {
+ const Ellipse* c = dynamic_cast<const Ellipse*>(geometry);
+ dolfin_assert(c);
+ impl->polygon_set.insert(make_ellipse(c));
+ break;
+ }
+ case CSGGeometry::Rectangle:
+ {
+ const Rectangle* r = dynamic_cast<const Rectangle*>(geometry);
+ dolfin_assert(r);
+ impl->polygon_set.insert(make_rectangle(r));
+ break;
+ }
+ case CSGGeometry::Polygon:
+ {
+ const Polygon* p = dynamic_cast<const Polygon*>(geometry);
+ dolfin_assert(p);
+ impl->polygon_set.insert(make_polygon(p));
+ break;
+ }
+ default:
+ dolfin_error("CSGCGALMeshGenerator2D.cpp",
+ "converting geometry to cgal polyhedron",
+ "Unhandled primitive type");
+ }
+}
+//-----------------------------------------------------------------------------
+CSGCGALDomain2D::CSGCGALDomain2D(const CSGCGALDomain2D &other)
+ : impl(new CSGCGALDomain2DImpl(other.impl->polygon_set))
+{
+}
+//-----------------------------------------------------------------------------
+CSGCGALDomain2D &CSGCGALDomain2D::operator=(const CSGCGALDomain2D &other)
+{
+ boost::scoped_ptr<CSGCGALDomain2DImpl> tmp(new CSGCGALDomain2DImpl(other.impl->polygon_set));
+
+ impl.swap(tmp);
+
+ return *this;
+}
+//-----------------------------------------------------------------------------
+double CSGCGALDomain2D::compute_boundingcircle_radius() const
+{
+ std::list<Polygon_with_holes_2> polygon_list;
+ impl->polygon_set.polygons_with_holes(std::back_inserter(polygon_list));
+
+ std::vector<Point_2> points;
+
+ for (std::list<Polygon_with_holes_2>::const_iterator pit = polygon_list.begin();
+ pit != polygon_list.end(); ++pit)
+ for (Polygon_2::Vertex_const_iterator vit = pit->outer_boundary().vertices_begin();
+ vit != pit->outer_boundary().vertices_end(); ++vit)
+ points.push_back(*vit);
+
+ Min_circle min_circle (points.begin(),
+ points.end(),
+ true); //randomize point order
+
+ return sqrt(CGAL::to_double(min_circle.circle().squared_radius()));
+}
+//-----------------------------------------------------------------------------
+void CSGCGALDomain2D::join_inplace(const CSGCGALDomain2D& other)
+{
+ impl->polygon_set.join(other.impl->polygon_set);
+}
+//-----------------------------------------------------------------------------
+void CSGCGALDomain2D::difference_inplace(const CSGCGALDomain2D& other)
+{
+ impl->polygon_set.difference(other.impl->polygon_set);
+}
+//-----------------------------------------------------------------------------
+void CSGCGALDomain2D::intersect_inplace(const CSGCGALDomain2D &other)
+{
+ impl->polygon_set.intersection(other.impl->polygon_set);
+}
+//-----------------------------------------------------------------------------
+bool CSGCGALDomain2D::point_in_domain(Point p) const
+{
+ const Point_2 p_(p.x(), p.y());
+ return impl->polygon_set.oriented_side(p_) == CGAL::ON_POSITIVE_SIDE;
+}
+//-----------------------------------------------------------------------------
+void CSGCGALDomain2D::get_vertices(std::list<std::vector<Point> >& l,
+ double truncate_threshold) const
+{
+ l.clear();
+
+ truncate_threshold *= truncate_threshold;
+
+ std::list<Polygon_with_holes_2> polygon_list;
+ impl->polygon_set.polygons_with_holes(std::back_inserter(polygon_list));
+
+ std::list<Polygon_with_holes_2>::const_iterator pit;
+ for (pit = polygon_list.begin(); pit != polygon_list.end(); ++pit)
+ {
+ const Polygon_2 &outer = pit->outer_boundary();
+
+ l.push_back(std::vector<Point>());
+ std::vector<Point> &v = l.back();
+ v.reserve(outer.size());
+
+ Polygon_2::Vertex_const_iterator prev = outer.vertices_begin();
+ Polygon_2::Vertex_const_iterator current = prev;
+ ++current;
+ for (; current != outer.vertices_end(); ++current)
+ {
+ if ( (*current - *prev).squared_length() < truncate_threshold)
+ continue;
+
+ const double x = CGAL::to_double(current->x());
+ const double y = CGAL::to_double(current->y());
+ v.push_back(Point(x, y));
+
+ prev = current;
+ }
+
+ current = outer.vertices_begin();
+ if ( (*current - *prev).squared_length() > truncate_threshold)
+ v.push_back(Point(CGAL::to_double(current->x()),
+ CGAL::to_double(current->y())));
+ }
+}
+//-----------------------------------------------------------------------------
+void CSGCGALDomain2D::get_holes(std::list<std::vector<Point> >& h,
+ double truncate_threshold) const
+{
+ h.clear();
+
+ std::list<Polygon_with_holes_2> polygon_list;
+ impl->polygon_set.polygons_with_holes(std::back_inserter(polygon_list));
+
+ std::list<Polygon_with_holes_2>::const_iterator pit;
+ for (pit = polygon_list.begin(); pit != polygon_list.end(); ++pit)
+ {
+ Hole_const_iterator hit;
+ for (hit = pit->holes_begin(); hit != pit->holes_end(); ++hit)
+ {
+ h.push_back(std::vector<Point>());
+ std::vector<Point> &v = h.back();
+ v.reserve(hit->size());
+
+ Polygon_2::Vertex_const_iterator prev = hit->vertices_begin();
+ Polygon_2::Vertex_const_iterator current = prev;
+ ++current;
+
+ for (; current != hit->vertices_end(); ++current)
+ {
+ if ( (*current - *prev).squared_length() < truncate_threshold)
+ continue;
+
+ const double x = CGAL::to_double(current->x());
+ const double y = CGAL::to_double(current->y());
+ v.push_back(Point(x, y));
+
+ prev = current;
+ }
+ current = hit->vertices_begin();
+ if ( (*current - *prev).squared_length() > truncate_threshold)
+ v.push_back(Point(CGAL::to_double(current->x()),
+ CGAL::to_double(current->y())));
+ }
+ }
+}
+
+#endif
diff --git a/dolfin/generation/CSGCGALDomain2D.h b/dolfin/generation/CSGCGALDomain2D.h
new file mode 100644
index 0000000..18990e8
--- /dev/null
+++ b/dolfin/generation/CSGCGALDomain2D.h
@@ -0,0 +1,71 @@
+// Copyright (C) 2013 Benjamin Kehlet
+//
+// This file is part of DOLFIN.
+//
+// DOLFIN is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// DOLFIN 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 Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
+//
+// First added: 2013-06-22
+// Last changed: 2013-08-06
+
+#ifdef HAS_CGAL
+
+#include "CSGGeometry.h"
+
+#include <dolfin/geometry/Point.h>
+#include <boost/scoped_ptr.hpp>
+
+struct CSGCGALDomain2DImpl;
+
+namespace dolfin
+{
+
+class CSGCGALDomain2D
+{
+ public:
+ // Create empty polygon
+ CSGCGALDomain2D();
+
+ // Construct polygon from Dolfin CSG geometry
+ CSGCGALDomain2D(const CSGGeometry *csg);
+
+ // Destructor
+ ~CSGCGALDomain2D();
+
+ // Copy constructor
+ CSGCGALDomain2D(const CSGCGALDomain2D &other);
+ CSGCGALDomain2D &operator=(const CSGCGALDomain2D &other);
+
+ // Boolean operators
+ void join_inplace(const CSGCGALDomain2D& other);
+ void intersect_inplace(const CSGCGALDomain2D& other);
+ void difference_inplace(const CSGCGALDomain2D& other);
+
+ bool point_in_domain(Point p) const;
+ double compute_boundingcircle_radius() const ;
+
+ // TODO: Replace this with a more C++-ish
+ // implementation, ie, take an outputiterator as arugment
+ // or define iterator
+ void get_vertices(std::list<std::vector<Point> >& v,
+ double truncate_threshold) const;
+
+ void get_holes(std::list<std::vector<Point> >& h,
+ double truncate_threshold) const;
+
+ boost::scoped_ptr<CSGCGALDomain2DImpl> impl;
+
+};
+}
+
+#endif
diff --git a/dolfin/generation/CSGCGALMeshGenerator2D.cpp b/dolfin/generation/CSGCGALMeshGenerator2D.cpp
index 8fb71d9..40c1824 100644
--- a/dolfin/generation/CSGCGALMeshGenerator2D.cpp
+++ b/dolfin/generation/CSGCGALMeshGenerator2D.cpp
@@ -19,34 +19,32 @@
// Modified by Benjamin Kehlet 2012-2013
//
// First added: 2012-05-10
-// Last changed: 2013-03-15
+// Last changed: 2013-08-06
-#include <vector>
-#include <cmath>
+#include "CSGCGALMeshGenerator2D.h"
+#include "CSGGeometry.h"
+#include "CSGOperators.h"
+#include "CSGPrimitives2D.h"
+#include "CSGCGALDomain2D.h"
#include <dolfin/common/constants.h>
#include <dolfin/math/basic.h>
#include <dolfin/mesh/Mesh.h>
#include <dolfin/mesh/MeshEditor.h>
+#include <dolfin/mesh/MeshFunction.h>
+#include <dolfin/mesh/MeshDomains.h>
+#include <dolfin/mesh/MeshValueCollection.h>
+#include <dolfin/log/log.h>
-#include "CSGCGALMeshGenerator2D.h"
-#include "CSGGeometry.h"
-#include "CSGOperators.h"
-#include "CSGPrimitives2D.h"
-#ifdef HAS_CGAL
-#include "CGALMeshBuilder.h"
+#include <vector>
+#include <cmath>
+#include <limits>
+
-#include <CGAL/basic.h>
-#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#ifdef HAS_CGAL
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Gmpq.h>
-#include <CGAL/Lazy_exact_nt.h>
-#include <CGAL/Simple_cartesian.h>
-#include <CGAL/Bounded_kernel.h>
-#include <CGAL/Nef_polyhedron_2.h>
-#include <CGAL/Polygon_2.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>
@@ -54,47 +52,21 @@
#include <CGAL/Delaunay_mesh_face_base_2.h>
#include <CGAL/Delaunay_mesh_size_criteria_2.h>
-#include <CGAL/Min_circle_2.h>
-#include <CGAL/Min_circle_2_traits_2.h>
-
-typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_Kernel;
typedef CGAL::Exact_predicates_inexact_constructions_kernel Inexact_Kernel;
-typedef CGAL::Lazy_exact_nt<CGAL::Gmpq> FT;
-typedef CGAL::Simple_cartesian<FT> EKernel;
-typedef CGAL::Bounded_kernel<EKernel> Extended_kernel;
-typedef CGAL::Nef_polyhedron_2<Extended_kernel> Nef_polyhedron_2;
-typedef Nef_polyhedron_2::Point Nef_point_2;
-
-typedef Nef_polyhedron_2::Explorer Explorer;
-typedef Explorer::Face_const_iterator Face_const_iterator;
-typedef Explorer::Hole_const_iterator Hole_const_iterator;
-typedef Explorer::Vertex_const_iterator Vertex_const_iterator;
-typedef Explorer::Halfedge_around_face_const_circulator
- Halfedge_around_face_const_circulator;
-typedef Explorer::Vertex_const_handle Vertex_const_handle;
-typedef Explorer::Halfedge_const_handle Halfedge_const_handle;
-
typedef CGAL::Triangulation_vertex_base_2<Inexact_Kernel> Vertex_base;
typedef CGAL::Constrained_triangulation_face_base_2<Inexact_Kernel> Face_base;
-// Min enclosing circle typedefs
-typedef CGAL::Min_circle_2_traits_2<Extended_kernel> Min_Circle_Traits;
-typedef CGAL::Min_circle_2<Min_Circle_Traits> Min_circle;
-typedef CGAL::Circle_2<Extended_kernel> CGAL_Circle;
-
-
-template <class Gt,
- class Fb >
-class Enriched_face_base_2 : public Fb
+template <class Gt, class Fb >
+class Enriched_face_base_2 : public Fb
{
-public:
+ public:
typedef Gt Geom_traits;
typedef typename Fb::Vertex_handle Vertex_handle;
typedef typename Fb::Face_handle Face_handle;
- template < typename TDS2 >
- struct Rebind_TDS
+ template <typename TDS2>
+ struct Rebind_TDS
{
typedef typename Fb::template Rebind_TDS<TDS2>::Other Fb2;
typedef Enriched_face_base_2<Gt,Fb2> Other;
@@ -102,24 +74,33 @@ public:
protected:
int status;
+ bool in_domain;
public:
Enriched_face_base_2(): Fb(), status(-1) {}
- Enriched_face_base_2(Vertex_handle v0, Vertex_handle v1, Vertex_handle v2)
- : Fb(v0,v1,v2), status(-1) {}
+ Enriched_face_base_2(Vertex_handle v0,
+ Vertex_handle v1,
+ Vertex_handle v2)
+ : Fb(v0,v1,v2), status(-1), in_domain(true) {}
- Enriched_face_base_2(Vertex_handle v0, Vertex_handle v1, Vertex_handle v2,
- Face_handle n0, Face_handle n1, Face_handle n2)
- : Fb(v0,v1,v2,n0,n1,n2), status(-1) {}
+ Enriched_face_base_2(Vertex_handle v0,
+ Vertex_handle v1,
+ Vertex_handle v2,
+ Face_handle n0,
+ Face_handle n1,
+ Face_handle n2)
+ : Fb(v0,v1,v2,n0,n1,n2), status(-1), in_domain(true) {}
inline
bool is_in_domain() const
- { return (status%2 == 1); }
+ //{ return (status%2 == 1); }
+ { return in_domain; }
inline
void set_in_domain(const bool b)
- { status = (b ? 1 : 0); }
+ //{ status = (b ? 1 : 0); }
+ { in_domain = b; }
inline
void set_counter(int i)
@@ -147,253 +128,274 @@ typedef CDT::Vertex_handle Vertex_handle;
typedef CDT::Face_handle Face_handle;
typedef CDT::All_faces_iterator All_faces_iterator;
-typedef CGAL::Polygon_2<Inexact_Kernel> Polygon_2;
typedef Inexact_Kernel::Point_2 Point_2;
using namespace dolfin;
//-----------------------------------------------------------------------------
CSGCGALMeshGenerator2D::CSGCGALMeshGenerator2D(const CSGGeometry& geometry)
- : _geometry(geometry)
+: geometry(geometry)
{
parameters = default_parameters();
+ //subdomains.push_back(reference_to_no_delete_pointer(geometry));
}
//-----------------------------------------------------------------------------
CSGCGALMeshGenerator2D::~CSGCGALMeshGenerator2D() {}
//-----------------------------------------------------------------------------
-Nef_polyhedron_2 make_circle(const Circle* c)
-{
- std::vector<Nef_point_2> pts;
- for (std::size_t i = 0; i < c->fragments(); i++)
+//-----------------------------------------------------------------------------
+static void print_edge(CDT::Edge edge)
+{
+ const int i = edge.second;
+ std::cout << "Edge: (" << edge.first->vertex( (i+1)%3 )->point()
+ << "), (" << edge.first->vertex( (i+2)%3 )->point() << ")" << std::endl;
+}
+//-----------------------------------------------------------------------------
+static void print_face(const CDT::Face_handle f)
+{
+ std::cout << "Face:" << std::endl;
+ std::cout << " " << f->vertex(0)->point() << std::endl;
+ std::cout << " " << f->vertex(1)->point() << std::endl;
+ std::cout << " " << f->vertex(2)->point() << std::endl;
+ std::cout << " " << (f->is_in_domain() ? "In domain" : "Not in domain") << std::endl;
+}
+//-----------------------------------------------------------------------------
+// print the faces of a triangulation
+static void print_triangulation(const CDT &cdt)
+{
+ for ( CDT::Finite_faces_iterator it = cdt.finite_faces_begin();
+ it != cdt.finite_faces_end(); ++it)
{
- const double phi = (2*DOLFIN_PI*i) / c->fragments();
- const double x = c->center().x() + c->radius()*cos(phi);
- const double y = c->center().y() + c->radius()*sin(phi);
- pts.push_back(Nef_point_2(x, y));
+ print_face(it);
}
-
- return Nef_polyhedron_2(pts.begin(), pts.end(), Nef_polyhedron_2::INCLUDED);
}
//-----------------------------------------------------------------------------
-Nef_polyhedron_2 make_ellipse(const Ellipse* e)
+static void print_constrained_edges(const CDT& cdt)
{
- std::vector<Nef_point_2> pts;
+ std::cout << "-- Constrained edges" << std::endl;
- for (std::size_t i = 0; i < e->fragments(); i++)
+ for (CDT::Finite_edges_iterator it = cdt.finite_edges_begin();
+ it != cdt.finite_edges_end();
+ it++)
{
- const double phi = (2*DOLFIN_PI*i) / e->fragments();
- const double x = e->center().x() + e->a()*cos(phi);
- const double y = e->center().y() + e->b()*sin(phi);
- pts.push_back(Nef_point_2(x, y));
+ if (cdt.is_constrained(*it))
+ print_edge(*it);
}
- return Nef_polyhedron_2(pts.begin(), pts.end(), Nef_polyhedron_2::INCLUDED);
+ std::cout << "--" << std::endl;
}
//-----------------------------------------------------------------------------
-Nef_polyhedron_2 make_rectangle(const Rectangle* r)
+void explore_subdomain(CDT &ct,
+ CDT::Face_handle start,
+ std::list<CDT::Face_handle>& other_domains)
{
- const double x0 = std::min(r->first_corner().x(), r->first_corner().y());
- const double y0 = std::max(r->first_corner().x(), r->first_corner().y());
+ std::list<Face_handle> queue;
+ queue.push_back(start);
- const double x1 = std::min(r->second_corner().x(), r->second_corner().y());
- const double y1 = std::max(r->second_corner().x(), r->second_corner().y());
+ while (!queue.empty())
+ {
+ CDT::Face_handle face = queue.front();
+ queue.pop_front();
- std::vector<Nef_point_2> pts;
- pts.push_back(Nef_point_2(x0, x1));
- pts.push_back(Nef_point_2(y0, x1));
- pts.push_back(Nef_point_2(y0, y1));
- pts.push_back(Nef_point_2(x0, y1));
+ for(int i = 0; i < 3; i++)
+ {
+ Face_handle n = face->neighbor(i);
+ if (ct.is_infinite(n))
+ continue;
- return Nef_polyhedron_2(pts.begin(), pts.end(), Nef_polyhedron_2::INCLUDED);
-}
-//-----------------------------------------------------------------------------
-Nef_polyhedron_2 make_polygon(const Polygon* p)
-{
- std::vector<Nef_point_2> pts;
- std::vector<Point>::const_iterator v;
- for (v = p->vertices().begin(); v != p->vertices().end(); ++v)
- pts.push_back(Nef_point_2(v->x(), v->y()));
+ const CDT::Edge e(face,i);
- return Nef_polyhedron_2(pts.begin(), pts.end(), Nef_polyhedron_2::INCLUDED);
+ if (n->counter() == -1)
+ {
+ if (ct.is_constrained(e))
+ {
+ // Reached a border
+ other_domains.push_back(n);
+ } else
+ {
+ // Set neighbor interface to the same and push to queue
+ n->set_counter(face->counter());
+ n->set_in_domain(face->is_in_domain());
+ queue.push_back(n);
+ }
+ }
+ }
+ }
}
//-----------------------------------------------------------------------------
-static Nef_polyhedron_2 convertSubTree(const CSGGeometry *geometry)
+// Set the member in_domain and counter for all faces in the cdt
+void explore_subdomains(CDT& cdt,
+ const CSGCGALDomain2D& total_domain,
+ const std::vector<std::pair<std::size_t, CSGCGALDomain2D> > &subdomain_geometries)
{
- switch (geometry->getType()) {
- case CSGGeometry::Union:
- {
- const CSGUnion* u = dynamic_cast<const CSGUnion*>(geometry);
- dolfin_assert(u);
- return convertSubTree(u->_g0.get()) + convertSubTree(u->_g1.get());
- break;
- }
- case CSGGeometry::Intersection:
- {
- const CSGIntersection* u = dynamic_cast<const CSGIntersection*>(geometry);
- dolfin_assert(u);
- return convertSubTree(u->_g0.get()) * convertSubTree(u->_g1.get());
- break;
- }
- case CSGGeometry::Difference:
- {
- const CSGDifference* u = dynamic_cast<const CSGDifference*>(geometry);
- dolfin_assert(u);
- return convertSubTree(u->_g0.get()) - convertSubTree(u->_g1.get());
- break;
- }
- case CSGGeometry::Circle:
- {
- const Circle* c = dynamic_cast<const Circle*>(geometry);
- dolfin_assert(c);
- return make_circle(c);
- break;
- }
- case CSGGeometry::Ellipse:
+ // Set counter to -1 for all faces
+ for (CDT::Finite_faces_iterator it = cdt.finite_faces_begin(); it != cdt.finite_faces_end(); ++it)
{
- const Ellipse* c = dynamic_cast<const Ellipse*>(geometry);
- dolfin_assert(c);
- return make_ellipse(c);
- break;
- }
- case CSGGeometry::Rectangle:
- {
- const Rectangle* r = dynamic_cast<const Rectangle*>(geometry);
- dolfin_assert(r);
- return make_rectangle(r);
- break;
+ it->set_counter(-1);
+ it->set_in_domain(false);
}
- case CSGGeometry::Polygon:
+
+ std::list<CDT::Face_handle> subdomains;
+ subdomains.push_back(cdt.finite_faces_begin());
+
+ //print_face(subdomains.front());
+ //dolfin_assert(face_in_domain(subdomains.front(), total_domain));
+
+ while (!subdomains.empty())
{
- const Polygon* p = dynamic_cast<const Polygon*>(geometry);
- dolfin_assert(p);
- return make_polygon(p);
- break;
- }
- default:
- dolfin_error("CSGCGALMeshGenerator2D.cpp",
- "converting geometry to cgal polyhedron",
- "Unhandled primitive type");
+ const CDT::Face_handle f = subdomains.front();
+ subdomains.pop_front();
+
+ if (f->counter() < 0)
+ {
+ // Set default marker (0)
+ f->set_counter(0);
+
+ const Point_2 p0 = f->vertex(0)->point();
+ const Point_2 p1 = f->vertex(1)->point();
+ const Point_2 p2 = f->vertex(2)->point();
+
+ Point p( (p0[0] + p1[0] + p2[0]) / 3.0,
+ (p0[1] + p1[1] + p2[1]) / 3.0 );
+
+ // Set the in_domain member (is face in the total domain)
+ f->set_in_domain(total_domain.point_in_domain(p));
+
+ for (int i = subdomain_geometries.size(); i > 0; --i)
+ {
+ if (subdomain_geometries[i-1].second.point_in_domain(p))
+ {
+ f->set_counter(subdomain_geometries[i-1].first);
+ break;
+ }
+ }
+
+ explore_subdomain(cdt, f, subdomains);
+ }
}
- return Nef_polyhedron_2();
}
//-----------------------------------------------------------------------------
-// Taken from examples/Triangulation_2/polygon_triangulation.cpp in the
-// CGAL source tree.
-//
-// Explore set of facets connected with non constrained edges, and
-// attribute to each such set a counter. We start from facets
-// incident to the infinite vertex, with a counter set to 0. Then we
-// recursively consider the non-explored facets incident to
-// constrained edges bounding the former set and increase thecounter
-// by 1. Facets in the domain are those with an odd nesting level.
-void mark_domains(CDT& ct,
- CDT::Face_handle start,
- int index,
- std::list<CDT::Edge>& border)
+ // Insert edges from polygon as constraints in the triangulation
+void add_subdomain(CDT& cdt, const CSGCGALDomain2D& cgal_geometry, double threshold)
{
- if (start->counter() != -1)
- return;
- std::list<CDT::Face_handle> queue;
- queue.push_back(start);
+ // Insert the outer boundaries of the domain
+ {
+ std::list<std::vector<Point> > v;
+ cgal_geometry.get_vertices(v, threshold);
- while (!queue.empty())
+ for (std::list<std::vector<Point> >::const_iterator pit = v.begin();
+ pit != v.end(); ++pit)
+ {
+ std::vector<Point>::const_iterator it = pit->begin();
+ Vertex_handle first = cdt.insert(Point_2(it->x(), it->y()));
+ Vertex_handle prev = first;
+ ++it;
+
+ for(; it != pit->end(); ++it)
+ {
+ Vertex_handle current = cdt.insert(Point_2(it->x(), it->y()));
+ cdt.insert_constraint(prev, current);
+ prev = current;
+ }
+
+ cdt.insert_constraint(first, prev);
+ }
+ }
+
+ // Insert holes
{
- CDT::Face_handle fh = queue.front();
- queue.pop_front();
- if (fh->counter() == -1)
+ std::list<std::vector<Point> > holes;
+ cgal_geometry.get_holes(holes, threshold);
+
+ for (std::list<std::vector<Point> >::const_iterator hit = holes.begin();
+ hit != holes.end(); ++hit)
{
- fh->counter() = index;
- fh->set_in_domain(index%2 == 1);
- for (int i = 0; i < 3; i++)
+
+ std::vector<Point>::const_iterator pit = hit->begin();
+ Vertex_handle first = cdt.insert(Point_2(pit->x(), pit->y()));
+ Vertex_handle prev = first;
+ ++pit;
+
+ for(; pit != hit->end(); ++pit)
{
- CDT::Edge e(fh, i);
- CDT::Face_handle n = fh->neighbor(i);
- if (n->counter() == -1)
- {
- if (ct.is_constrained(e))
- border.push_back(e);
- else
- queue.push_back(n);
- }
+ Vertex_handle current = cdt.insert(Point_2(pit->x(), pit->y()));
+ cdt.insert_constraint(prev, current);
+ prev = current;
}
+
+ cdt.insert_constraint(first, prev);
}
}
}
//-----------------------------------------------------------------------------
-void mark_domains(CDT& cdt)
+double shortest_constrained_edge(const CDT &cdt)
{
- for (CDT::All_faces_iterator it = cdt.all_faces_begin();
- it != cdt.all_faces_end(); ++it)
+ double min_length = std::numeric_limits<double>::max();
+ for (CDT::Finite_edges_iterator it = cdt.finite_edges_begin();
+ it != cdt.finite_edges_end();
+ it++)
{
- it->set_counter(-1);
- }
+ if (!cdt.is_constrained(*it))
+ continue;
- int index = 0;
- std::list<CDT::Edge> border;
- mark_domains(cdt, cdt.infinite_face(), index++, border);
- while (!border.empty())
- {
- CDT::Edge e = border.front();
- border.pop_front();
- CDT::Face_handle n = e.first->neighbor(e.second);
- if (n->counter() == -1)
- mark_domains(cdt, n, e.first->counter()+1, border);
+ // An edge is an std::pair<Face_handle, int>
+ // see CGAL/Triangulation_data_structure_2.h
+ CDT::Face_handle f = it->first;
+ const int i = it->second;
+
+ CDT::Point p1 = f->vertex( (i+1)%3 )->point();
+ CDT::Point p2 = f->vertex( (i+2)%3 )->point();
+
+ min_length = std::min(CGAL::to_double((p1-p2).squared_length()), min_length);
}
+
+ return min_length;
}
//-----------------------------------------------------------------------------
void CSGCGALMeshGenerator2D::generate(Mesh& mesh)
{
- Nef_polyhedron_2 cgal_geometry = convertSubTree(&_geometry);
-
// Create empty CGAL triangulation
CDT cdt;
- // Explore the Nef polyhedron and insert constraints in the
- // triangulation
- Explorer explorer = cgal_geometry.explorer();
- for (Face_const_iterator fit = explorer.faces_begin() ;
- fit != explorer.faces_end(); fit++)
+ // Convert the CSG geometry to a CGAL Polygon
+ log(TRACE, "Converting geometry to CGAL polygon");
+ CSGCGALDomain2D total_domain(&geometry);
+
+ add_subdomain(cdt, total_domain, parameters["edge_minimum"]);
+
+ // Empty polygon, will be populated when traversing the subdomains
+ CSGCGALDomain2D overlaying;
+
+ std::vector<std::pair<std::size_t, CSGCGALDomain2D> >
+ subdomain_geometries;
+
+ // Add the subdomains to the CDT. Traverse in reverse order to get the latest
+ // added subdomain on top
+ std::list<std::pair<std::size_t, boost::shared_ptr<const CSGGeometry> > >::const_reverse_iterator it;
+
+ if (!geometry.subdomains.empty())
+ log(TRACE, "Processing subdomains");
+
+ for (it = geometry.subdomains.rbegin(); it != geometry.subdomains.rend(); ++it)
{
- // Skip face if it is not part of polygon
- if (!explorer.mark(fit))
- continue;
+ const std::size_t current_index = it->first;
+ boost::shared_ptr<const CSGGeometry> current_subdomain = it->second;
- Halfedge_around_face_const_circulator hafc
- = explorer.face_cycle(fit), done(hafc);
- do
- {
- Vertex_handle va
- = cdt.insert(Point_2(to_double(hafc->vertex()->point().x()),
- to_double(hafc->vertex()->point().y())));
- Vertex_handle vb
- = cdt.insert(Point_2(to_double(hafc->next()->vertex()->point().x()),
- to_double(hafc->next()->vertex()->point().y())));
- cdt.insert_constraint(va, vb);
- hafc++;
- } while (hafc != done);
-
- Hole_const_iterator hit = explorer.holes_begin(fit);
- for (; hit != explorer.holes_end(fit); hit++)
- {
- Halfedge_around_face_const_circulator hafc1(hit), done1(hit);
- do
- {
- Vertex_handle va
- = cdt.insert(Point_2(to_double(hafc1->vertex()->point().x()),
- to_double(hafc1->vertex()->point().y())));
- Vertex_handle vb
- = cdt.insert(Point_2(to_double(hafc1->next()->vertex()->point().x()),
- to_double(hafc1->next()->vertex()->point().y())));
- cdt.insert_constraint(va, vb);
- hafc1++;
- } while (hafc1 != done1);
- }
+ CSGCGALDomain2D cgal_geometry(current_subdomain.get());
+ cgal_geometry.difference_inplace(overlaying);
+
+ subdomain_geometries.push_back(std::make_pair(current_index,
+ cgal_geometry));
+
+ add_subdomain(cdt, cgal_geometry, parameters["edge_minimum"]);
+
+ overlaying.join_inplace(cgal_geometry);
}
- // Mark parts that are inside and outside the domain
- mark_domains(cdt);
+ explore_subdomains(cdt, total_domain, subdomain_geometries);
+
+ log(TRACE, "Refining mesh");
// Create mesher
CGAL_Mesher_2 mesher(cdt);
@@ -411,34 +413,25 @@ void CSGCGALMeshGenerator2D::generate(Mesh& mesh)
Point_2 p2 = fit->vertex(2)->point();
const double x = (p0[0] + p1[0] + p2[0]) / 3;
const double y = (p0[1] + p1[1] + p2[1]) / 3;
+
list_of_seeds.push_back(Point_2(x, y));
}
}
+
mesher.set_seeds(list_of_seeds.begin(), list_of_seeds.end(), true);
// Set shape and size criteria
const int mesh_resolution = parameters["mesh_resolution"];
if (mesh_resolution > 0)
{
- cout << "Mesh resolution set" << endl;
-
- std::vector<Nef_point_2> points;
- for (Vertex_const_iterator it = explorer.vertices_begin();
- it != explorer.vertices_end();
- it++)
- points.push_back(it->point());
+ const double min_radius = total_domain.compute_boundingcircle_radius();
+ const double cell_size = 2.0*min_radius/mesh_resolution;
- Min_circle min_circle (points.begin(),
- points.end(),
- true); //randomize point order
-
- const double cell_size
- = 2.0*sqrt(CGAL::to_double(min_circle.circle().squared_radius()))/mesh_resolution;
Mesh_criteria_2 criteria(parameters["triangle_shape_bound"],
cell_size);
mesher.set_criteria(criteria);
- }
+ }
else
{
// Set shape and size criteria
@@ -453,10 +446,12 @@ void CSGCGALMeshGenerator2D::generate(Mesh& mesh)
// Make sure triangulation is valid
dolfin_assert(cdt.is_valid());
+ // Mark the subdomains
+ explore_subdomains(cdt, total_domain, subdomain_geometries);
+
// Clear mesh
mesh.clear();
- // Get various dimensions
const std::size_t gdim = cdt.finite_vertices_begin()->point().dimension();
const std::size_t tdim = cdt.dimension();
const std::size_t num_vertices = cdt.number_of_vertices();
@@ -464,7 +459,7 @@ void CSGCGALMeshGenerator2D::generate(Mesh& mesh)
// Count valid cells
std::size_t num_cells = 0;
CDT::Finite_faces_iterator cgal_cell;
- for (cgal_cell = cdt.finite_faces_begin();
+ for (cgal_cell = cdt.finite_faces_begin();
cgal_cell != cdt.finite_faces_end(); ++cgal_cell)
{
// Add cell if it is in the domain
@@ -490,8 +485,6 @@ void CSGCGALMeshGenerator2D::generate(Mesh& mesh)
Point p;
p[0] = cgal_vertex->point()[0];
p[1] = cgal_vertex->point()[1];
- if (gdim == 3)
- p[2] = cgal_vertex->point()[2];
// Add mesh vertex
mesh_editor.add_vertex(vertex_index, p);
@@ -499,36 +492,37 @@ void CSGCGALMeshGenerator2D::generate(Mesh& mesh)
// Attach index to vertex and increment
cgal_vertex->info() = vertex_index++;
}
+
dolfin_assert(vertex_index == num_vertices);
- // Add cells to mesh
+ // Add cells to mesh and build domain marker mesh function
+ MeshDomains &domain_markers = mesh.domains();
std::size_t cell_index = 0;
- for (cgal_cell = cdt.finite_faces_begin();
- cgal_cell != cdt.finite_faces_end(); ++cgal_cell)
+ for (cgal_cell = cdt.finite_faces_begin(); cgal_cell != cdt.finite_faces_end(); ++cgal_cell)
{
// Add cell if it is in the domain
if (cgal_cell->is_in_domain())
{
- mesh_editor.add_cell(cell_index++,
+ mesh_editor.add_cell(cell_index,
cgal_cell->vertex(0)->info(),
cgal_cell->vertex(1)->info(),
cgal_cell->vertex(2)->info());
+
+ domain_markers.set_marker(std::make_pair(cell_index, cgal_cell->counter()), 2);
+ ++cell_index;
}
}
dolfin_assert(cell_index == num_cells);
// Close mesh editor
mesh_editor.close();
-
- // Build DOLFIN mesh from CGAL triangulation
- // FIXME: Why does this not work correctly?
- //CGALMeshBuilder::build(mesh, cdt);
}
#else
namespace dolfin
{
CSGCGALMeshGenerator2D::CSGCGALMeshGenerator2D(const CSGGeometry& geometry)
+ : geometry(geometry)
{
dolfin_error("CSGCGALMeshGenerator2D.cpp",
"Create mesh generator",
diff --git a/dolfin/generation/CSGCGALMeshGenerator2D.h b/dolfin/generation/CSGCGALMeshGenerator2D.h
index 7845f2a..1f644d2 100644
--- a/dolfin/generation/CSGCGALMeshGenerator2D.h
+++ b/dolfin/generation/CSGCGALMeshGenerator2D.h
@@ -18,7 +18,7 @@
// Modified by Johannes Ring, 2012
//
// First added: 2012-05-10
-// Last changed: 2013-03-15
+// Last changed: 2013-04-05
#ifndef __CSG_CGAL_MESH_GENERATOR2D_H
#define __CSG_CGAL_MESH_GENERATOR2D_H
@@ -40,6 +40,7 @@ namespace dolfin
public :
CSGCGALMeshGenerator2D(const CSGGeometry& geometry);
+ //CSGCGALMeshGenerator2D(const std::vector<boost::shared_ptr<const CSGGeometry> >& subdomains);
~CSGCGALMeshGenerator2D();
@@ -52,16 +53,15 @@ namespace dolfin
p.add("mesh_resolution", 64);
p.add("triangle_shape_bound", 0.125);
p.add("cell_size", 0.25);
+
+ // shorter edges in the domain will be collapsed before meshing
+ p.add("edge_minimum", 10e-5);
return p;
}
private:
-
- #ifdef HAS_CGAL
- const CSGGeometry& _geometry;
- #endif
-
+ const CSGGeometry &geometry;
};
}
diff --git a/dolfin/generation/CSGGeometry.cpp b/dolfin/generation/CSGGeometry.cpp
index 9b54812..3fb4230 100644
--- a/dolfin/generation/CSGGeometry.cpp
+++ b/dolfin/generation/CSGGeometry.cpp
@@ -15,10 +15,13 @@
// You should have received a copy of the GNU Lesser General Public License
// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
//
+// Modified by Benjamin Kehlet, 2013
+//
// First added: 2012-04-11
-// Last changed: 2012-04-13
+// Last changed: 2013-04-18
#include "CSGGeometry.h"
+#include <dolfin/common/NoDeleter.h>
using namespace dolfin;
@@ -33,3 +36,34 @@ CSGGeometry::~CSGGeometry()
// Do nothing
}
//-----------------------------------------------------------------------------
+void CSGGeometry::set_subdomain(std::size_t i, boost::shared_ptr<CSGGeometry> s)
+{
+ dolfin_assert(dim() == s->dim());
+
+ if (i == 0)
+ {
+ error("Setting reserved CSG subdomain (0)");
+ }
+
+ // Check if i already used
+ std::list<std::pair<std::size_t, boost::shared_ptr<const CSGGeometry> > >::iterator it = subdomains.begin();
+ while (it != subdomains.end())
+ {
+ if (it->first == i)
+ {
+ warning("Double declaration of CSG subdomain with index %u.", i);
+
+ // Remove existing declaration
+ it = subdomains.erase(it);
+ }
+ else
+ ++it;
+ }
+
+ subdomains.push_back(std::make_pair(i, s));
+}
+//-----------------------------------------------------------------------------
+void CSGGeometry::set_subdomain(std::size_t i, CSGGeometry& s)
+{
+ set_subdomain(i, reference_to_no_delete_pointer(s));
+}
diff --git a/dolfin/generation/CSGGeometry.h b/dolfin/generation/CSGGeometry.h
index d4d6981..30a976a 100644
--- a/dolfin/generation/CSGGeometry.h
+++ b/dolfin/generation/CSGGeometry.h
@@ -25,6 +25,9 @@
#define __CSG_GEOMETRY_H
#include <cstddef>
+#include <vector>
+#include <list>
+#include <boost/shared_ptr.hpp>
#include <dolfin/common/Variable.h>
namespace dolfin
@@ -49,12 +52,19 @@ namespace dolfin
/// Informal string representation
virtual std::string str(bool verbose) const = 0;
+ void set_subdomain(std::size_t i, boost::shared_ptr<CSGGeometry> s);
+ void set_subdomain(std::size_t i, CSGGeometry& s);
+
enum Type { Box, Sphere, Cone, Tetrahedron, Surface3D, Circle, Ellipse, Rectangle, Polygon, Union, Intersection, Difference };
virtual Type getType() const = 0;
virtual bool is_operator() const = 0;
+
+ std::list<std::pair<std::size_t, boost::shared_ptr<const CSGGeometry> > > subdomains;
};
+
+
}
#endif
diff --git a/dolfin/generation/CSGOperators.cpp b/dolfin/generation/CSGOperators.cpp
index a23097d..896ba91 100644
--- a/dolfin/generation/CSGOperators.cpp
+++ b/dolfin/generation/CSGOperators.cpp
@@ -16,9 +16,10 @@
// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
//
// Modified by Johannes Ring, 2012
+// Modified by Benjamin Kehlet, 2013
//
// First added: 2012-04-13
-// Last changed: 2012-05-03
+// Last changed: 2013-06-24
#include <sstream>
@@ -46,12 +47,8 @@ CSGUnion::CSGUnion(boost::shared_ptr<CSGGeometry> g0,
"Dimensions of geomestries don't match (%d vs %d)",
g0->dim(), g1->dim());
}
-}
-//-----------------------------------------------------------------------------
-std::size_t CSGUnion::dim() const
-{
- assert(_g0->dim() == _g1->dim());
- return _g0->dim();
+
+ dim_ = g0->dim();
}
//-----------------------------------------------------------------------------
std::string CSGUnion::str(bool verbose) const
@@ -95,12 +92,8 @@ CSGDifference::CSGDifference(boost::shared_ptr<CSGGeometry> g0,
"Dimensions of geomestries don't match (%d vs %d)",
g0->dim(), g1->dim());
}
-}
-//-----------------------------------------------------------------------------
-std::size_t CSGDifference::dim() const
-{
- assert(_g0->dim() == _g1->dim());
- return _g0->dim();
+
+ dim_ = g0->dim();
}
//-----------------------------------------------------------------------------
std::string CSGDifference::str(bool verbose) const
@@ -144,12 +137,8 @@ CSGIntersection::CSGIntersection(boost::shared_ptr<CSGGeometry> g0,
"Dimensions of geomestries don't match (%d vs %d)",
g0->dim(), g1->dim());
}
-}
-//-----------------------------------------------------------------------------
-std::size_t CSGIntersection::dim() const
-{
- assert(_g0->dim() == _g1->dim());
- return _g0->dim();
+
+ dim_ = g0->dim();
}
//-----------------------------------------------------------------------------
std::string CSGIntersection::str(bool verbose) const
diff --git a/dolfin/generation/CSGOperators.h b/dolfin/generation/CSGOperators.h
index 2ada53f..8645d30 100644
--- a/dolfin/generation/CSGOperators.h
+++ b/dolfin/generation/CSGOperators.h
@@ -15,10 +15,10 @@
// You should have received a copy of the GNU Lesser General Public License
// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
//
-// Modified by Benjamin Kehlet, 2012
+// Modified by Benjamin Kehlet, 2012-2013
//
// First added: 2012-04-13
-// Last changed: 2012-11-05
+// Last changed: 2013-06-24
#ifndef __CSG_OPERATORS_H
#define __CSG_OPERATORS_H
@@ -36,6 +36,10 @@ namespace dolfin
{
public:
virtual bool is_operator() const { return true; }
+ std::size_t dim() const { return dim_;}
+
+ protected :
+ std::size_t dim_;
};
/// Union of CSG geometries
@@ -47,9 +51,6 @@ namespace dolfin
CSGUnion(boost::shared_ptr<CSGGeometry> g0,
boost::shared_ptr<CSGGeometry> g1);
- /// Return dimension of geometry
- std::size_t dim() const;
-
/// Informal string representation
std::string str(bool verbose) const;
@@ -69,9 +70,6 @@ namespace dolfin
CSGDifference(boost::shared_ptr<CSGGeometry> g0,
boost::shared_ptr<CSGGeometry> g1);
- /// Return dimension of geometry
- std::size_t dim() const;
-
/// Informal string representation
std::string str(bool verbose) const;
@@ -92,9 +90,6 @@ namespace dolfin
CSGIntersection(boost::shared_ptr<CSGGeometry> g0,
boost::shared_ptr<CSGGeometry> g1);
- /// Return dimension of geometry
- std::size_t dim() const;
-
/// Informal string representation
std::string str(bool verbose) const;
diff --git a/dolfin/generation/CSGPrimitive.h b/dolfin/generation/CSGPrimitive.h
index 0dc81d2..3b02024 100644
--- a/dolfin/generation/CSGPrimitive.h
+++ b/dolfin/generation/CSGPrimitive.h
@@ -15,10 +15,10 @@
// You should have received a copy of the GNU Lesser General Public License
// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
//
-// Modified by Benjamin Kehlet, 2012
+// Modified by Benjamin Kehlet, 2012-2013
//
// First added: 2012-04-11
-// Last changed: 2012-11-05
+// Last changed: 2013-06-24
#ifndef __CSG_PRIMITIVE_H
#define __CSG_PRIMITIVE_H
@@ -34,6 +34,7 @@ namespace dolfin
{
public:
virtual bool is_operator() const { return false; }
+ virtual std::size_t dim() const = 0;
};
}
diff --git a/dolfin/generation/CSGPrimitives2D.h b/dolfin/generation/CSGPrimitives2D.h
index e5af011..1e79e49 100644
--- a/dolfin/generation/CSGPrimitives2D.h
+++ b/dolfin/generation/CSGPrimitives2D.h
@@ -151,10 +151,10 @@ namespace dolfin
Type getType() const { return CSGGeometry::Rectangle; }
/// Return first corner
- Point first_corner() const { return Point(_x0, _y0); }
+ Point first_corner() const { return Point(_x0, _x1); }
/// Return second corner
- Point second_corner() const { return Point(_x1, _y1); }
+ Point second_corner() const { return Point(_y0, _y1); }
private:
