dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #16776
[Fwd: [Merge] lp:~massing/dolfin/cgal_branch into lp:dolfin]
I suggest that anyone who would like to see what the 'merge proposal'
interface looks like on Launchpad should follow the link
https://code.launchpad.net/dolfin/+activereviews
before the merge is approved.
Garth
-------- Original Message --------
Subject: [Merge] lp:~massing/dolfin/cgal_branch into lp:dolfin
Date: Mon, 30 Nov 2009 01:20:28 -0000
From: Andre Massing <massing@xxxxxxxxx>
Reply-To: mp+15390@xxxxxxxxxxxxxxxxxx
To: mp+15390@xxxxxxxxxxxxxxxxxx
Andre Massing has proposed merging lp:~massing/dolfin/cgal_branch into
lp:dolfin.
Requested reviews:
DOLFIN Core Team (dolfin-core)
Related bugs:
#489867 CGAL exception bug reappeared after migration to bzr and
cleaning up templates.
https://bugs.launchpad.net/bugs/489867
This branch contains a CGAL based interface for some mesh intersection
functionality. It extends and replaces the old GTS based implementation.
Please consider it :) I haven't yet completely removed the old GTS
interface code and the intersect member function, spread in some
CellType classes in order to run some tests and benchmarks. If the this
codes made it into the trunk I would suggest to keep this left-overs
maybe for one release cycle to facilitate bug fixing (who knows :),
benchmarks etc.
Summary of changes:
- Intersection detection for mesh-point (1d-3d) and mesh-mesh
intersection (2d-1d, 2d-2d, 3d-2d,3d-3d)
(computing all cell ids for one mesh, now *uniquely* collected in a
std:set)
- Computing first intersected cell for mesh-point intersection,
accelerates significantly function
eval
ADDED:
- IntersectionOperator class: Provides the interface for the new
intersection
functionality based on the CGAL library (computational geometry).
Implementation via the pImpl idiom.
- IntersectionOperatorImplementation class: Abstract base class for the
pImpl
implementation. Only inherited classes contain the search tree, which type
depends is based on the geometric primitive it includes.
- IntersectionOperatorImplementation_d class: Contains search tree for
the mesh
and executes search operations.
- MeshPrimitive class: Implementation of the concept "AABBPrimitive"
required by the
CGAL AABB tree. Relies on the Primitive_Traits class to perform data
conversion,
- Primitive_Traits template class: Traits class which provides conversion
functionality from MeshEntity (s.a. IntervalCell, TetrahedronCell
etc.) -->
CGAL geometric primitive (Segment_3, Tetrahedron_3).
- demo/mesh/intersection/python/demo_3d.py demonstrates intersection of
3D cube
with 3D sphere running through the cube.
- added intersection functionality to complement CGAL interface
(Point_3_*_intersection.h) and add new intersection detection
(Tetrahedron -
Tetrahedron) and a template specialization for the Triangle
Tetrahedron do_intersect function to overwrite its CGAL default
implementation
(fix CGAL roundoffs exceptions due inexact geometry kernel).
CHANGED:
- MeshEntity class: Introduced default constructor (requires to change
private
member const Mesh & _mesh to const Mesh * _mesh), changed _mesh. to
_mesh-> in
(derived) classes when necessary.
- added comparison operator
- MeshEntityInterator class: Introduced default constructor, added
comparison
operator, added set_end function to create iterators which points past
the end
(more STL like, necessary for the AABB tree).
- Mesh class: Provide appropriate delegate functions.
- Function class: Replace GTS based intersection detector.
- swig: Added wrapper for std::set<uint> and a has_cgal function.
--
https://code.launchpad.net/~massing/dolfin/cgal_branch/+merge/15390
You are subscribed to branch lp:dolfin.
=== modified file '.bzrignore'
--- .bzrignore 2009-11-24 16:28:00 +0000
+++ .bzrignore 2009-11-30 01:20:26 +0000
@@ -38,3 +38,4 @@
RE:demo/.*/(cpp|python)/.*\.(pvd|vtu|pvtu|xml)
RE:demo/.*/cpp/demo
RE:test/.*/cpp/test
+*.os
=== added directory 'bench/mesh/basic'
=== renamed file 'bench/mesh/SConstruct' => 'bench/mesh/basic/SConstruct'
=== renamed file 'bench/mesh/bench.log' => 'bench/mesh/basic/bench.log'
=== renamed file 'bench/mesh/bench.py' => 'bench/mesh/basic/bench.py'
=== renamed file 'bench/mesh/main.cpp' => 'bench/mesh/basic/main.cpp'
=== modified file 'demo/mesh/intersection/python/demo.py'
--- demo/mesh/intersection/python/demo.py 2009-09-03 09:05:39 +0000
+++ demo/mesh/intersection/python/demo.py 2009-11-30 01:20:26 +0000
@@ -10,8 +10,8 @@
from dolfin import *
from numpy import *
-if not has_gts():
- print "DOLFIN must be compiled with GTS to run this demo."
+if not has_cgal():
+ print "DOLFIN must be compiled with CGAL to run this demo."
exit(0)
# Create meshes (omega0 overlapped by omega1)
@@ -34,7 +34,7 @@
# Compute intersection with boundary of square
boundary = BoundaryMesh(omega1)
- cells = omega0.intersection(boundary, False)
+ cells = omega0.all_intersected_entities(boundary)
# Mark intersected values
intersection.values()[:] = 0
@@ -49,13 +49,15 @@
else:
plot(intersection)
+# interactive()
# Rotate circle around (0.5, 0.5)
xr = x[:, 0] - 0.5
yr = x[:, 1] - 0.5
x[:,0] = 0.5 + (cos(dtheta)*xr - sin(dtheta)*yr)
x[:,1] = 0.5 + (sin(dtheta)*xr + cos(dtheta)*yr)
+ omega0.intersection_operator().clear()
theta += dtheta
# Hold plot
-interactive()
+#interactive()
=== added file 'demo/mesh/intersection/python/demo_3d.py'
--- demo/mesh/intersection/python/demo_3d.py 1970-01-01 00:00:00 +0000
+++ demo/mesh/intersection/python/demo_3d.py 2009-11-30 01:20:26 +0000
@@ -0,0 +1,72 @@
+"""This demo shows the intersection of the boundary of a unit square
+(omega1) with a unit circle (omega2) rotating around the center of the
+square.
+@todo Change camera perspective/ viewpoint to improve intersection visibility.
+"""
+
+__author__ = "Andre Massing (massing@xxxxxxxxx)"
+__date__ = "2008-11-17"
+__copyright__ = "Copyright (C) 2009 Andre Massing"
+__license__ = "GNU LGPL Version 2.1"
+
+from dolfin import *
+from numpy import *
+
+if not has_cgal():
+ print "DOLFIN must be compiled with CGAL to run this demo."
+ exit(0)
+
+sphere = UnitSphere(20)
+cube = UnitCube(20, 20, 20)
+#cube = UnitCube(10, 10, 10)
+
+# Access mesh geometry
+x = sphere.coordinates()
+
+# Start center and propagtion speed for the sphere.
+dt = 0.1
+t = -0.71
+
+# Scale and move the circle.
+x *= 0.4
+x[:] += t
+
+intersection = MeshFunction("uint", cube, cube.topology().dim())
+_first = True
+
+counter = 0
+while t < 1.4 :
+
+ boundary = BoundaryMesh(sphere)
+ # Compute intersection with boundary of square
+ cells = cube.all_intersected_entities(boundary)
+ print cells
+
+ # Mark intersected values
+ intersection.values()[:] = 0
+ intersection.values()[cells] = 1
+
+ counter +=1
+
+ # Plot intersection
+ if _first:
+ p = plot(intersection, rescale=False)
+ p.ren.SetViewPoint(-1,-1,-1)
+ p.ren.ResetCamera()
+ _first = False
+ if counter == 2:
+ interactive()
+ p.movie("movie")
+
+ else:
+ plot(intersection)
+
+ #Propagate sphere along the line t(1,1,1).
+ x[:,0] += dt
+ x[:,1] += dt
+ x[:,2] += dt
+
+ t += dt
+
+# Hold plot
+interactive()
=== modified file 'dolfin/adaptivity/AdaptiveObjects.cpp'
--- dolfin/adaptivity/AdaptiveObjects.cpp 2009-11-25 22:07:01 +0000
+++ dolfin/adaptivity/AdaptiveObjects.cpp 2009-11-30 01:20:26 +0000
@@ -6,7 +6,7 @@
#include <dolfin/common/NoDeleter.h>
#include <dolfin/la/GenericVector.h>
-#include <dolfin/mesh/IntersectionDetector.h>
+#include <dolfin/mesh/IntersectionOperator.h>
#include <dolfin/mesh/UniformMeshRefinement.h>
#include <dolfin/mesh/LocalMeshRefinement.h>
#include <dolfin/function/FunctionSpace.h>
@@ -154,6 +154,7 @@
// Copy data from new mesh to old
*mesh = new_mesh;
+
}
//-----------------------------------------------------------------------------
void AdaptiveObjects::refine(FunctionSpace* function_space,
@@ -198,9 +199,6 @@
// Copy vector from new function to old
*function->_vector = *new_function._vector;
- // FIXME: This won't be needed later with new CGAL interface
- function->intersection_detector.reset();
-
// FIXME: Might need to touch/update other data for function space here!
}
//-----------------------------------------------------------------------------
=== modified file 'dolfin/common/types.h'
--- dolfin/common/types.h 2008-10-02 19:11:34 +0000
+++ dolfin/common/types.h 2009-11-30 01:20:26 +0000
@@ -2,7 +2,7 @@
// Licensed under the GNU LGPL Version 2.1.
//
// First added: 2008-04-22
-// Last changed: 2008-10-02
+// Last changed: 2009-11-27
//
// This file provides DOLFIN typedefs for basic types.
@@ -10,6 +10,7 @@
#define __DOLFIN_TYPES_H
#include <complex>
+#include <set>
namespace dolfin
{
@@ -20,6 +21,9 @@
// Complex numbers
typedef std::complex<double> complex;
+ // (Ordered) set of uints
+ typedef std::set<uint> uint_set;
+
}
#endif
=== modified file 'dolfin/function/Function.cpp'
--- dolfin/function/Function.cpp 2009-11-29 00:43:25 +0000
+++ dolfin/function/Function.cpp 2009-11-30 01:20:26 +0000
@@ -3,9 +3,10 @@
//
// Modified by Garth N. Wells, 2005-2009.
// Modified by Martin Sandve Alnes, 2008.
+// Modified by Andre Massing, 2009.
//
// First added: 2003-11-28
-// Last changed: 2009-11-10
+// Last changed: 2009-11-29
#include <algorithm>
#include <boost/assign/list_of.hpp>
@@ -19,7 +20,7 @@
#include <dolfin/fem/FiniteElement.h>
#include <dolfin/fem/DofMap.h>
#include <dolfin/fem/UFC.h>
-#include <dolfin/mesh/IntersectionDetector.h>
+#include <dolfin/mesh/IntersectionOperator.h>
#include <dolfin/mesh/Vertex.h>
#include <dolfin/adaptivity/AdaptiveObjects.h>
#include "Data.h"
@@ -281,26 +282,14 @@
not_working_in_parallel("Function::eval at arbitray points.");
- // Initialize intersection detector if not done before
- if (!intersection_detector)
- {
- dolfin_debug("Initializing intersection detector");
- intersection_detector.reset(new IntersectionDetector(_function_space->mesh()));
- }
-
// Find the cell that contains x
const double* _x = &x[0];
Point point(_function_space->mesh().geometry().dim(), _x);
- std::vector<uint> cells;
- intersection_detector->intersection(point, cells);
- if (cells.size() < 1)
- {
- error("Unable to evaluate function at x = %s, not inside domain.",
- to_string(_x, geometric_dimension()).c_str());
- }
+ int id = _function_space->mesh().any_intersected_entity(point);
+ if (id == -1)
+ error("Unable to evaluate function at given point (not inside domain).");
- // Create cell
- Cell cell(_function_space->mesh(), cells[0]);
+ Cell cell(_function_space->mesh(), id);
UFCCell ufc_cell(cell);
// Evaluate function
=== modified file 'dolfin/function/Function.h'
--- dolfin/function/Function.h 2009-11-29 21:40:18 +0000
+++ dolfin/function/Function.h 2009-11-30 01:20:26 +0000
@@ -4,9 +4,10 @@
// Modified by Garth N. Wells, 2005-2009.
// Modified by Kristian B. Oelgaard, 2007.
// Modified by Martin Sandve Alnes, 2008.
+// Modified by Andre Massing, 2009.
//
// First added: 2003-11-28
-// Last changed: 2009-11-09
+// Last changed: 2009-11-29
#ifndef __FUNCTION_H
#define __FUNCTION_H
@@ -32,7 +33,6 @@
class FunctionSpace;
class GenericVector;
class Data;
- class IntersectionDetector;
/// This class represents a function u_h in a finite element
/// function space V_h, given by
@@ -177,9 +177,6 @@
mutable std::map<uint, uint> global_to_local;
mutable std::vector<uint> _off_process_dofs;
- // Intersection detector, used for evaluation at arbitrary points
- mutable boost::scoped_ptr<IntersectionDetector> intersection_detector;
-
// Scratch space, used for storing temporary local data
class LocalScratch
{
=== modified file 'dolfin/function/FunctionSpace.cpp'
--- dolfin/function/FunctionSpace.cpp 2009-11-17 13:03:46 +0000
+++ dolfin/function/FunctionSpace.cpp 2009-11-30 01:20:26 +0000
@@ -18,8 +18,6 @@
#include <dolfin/fem/UFC.h>
#include <dolfin/log/log.h>
#include <dolfin/mesh/Vertex.h>
-#include <dolfin/mesh/IntersectionDetector.h>
-#include <dolfin/mesh/Cell.h>
#include <dolfin/mesh/Mesh.h>
#include <dolfin/mesh/MeshFunction.h>
#include <dolfin/mesh/MeshPartitioning.h>
=== modified file 'dolfin/function/FunctionSpace.h'
--- dolfin/function/FunctionSpace.h 2009-11-09 21:12:36 +0000
+++ dolfin/function/FunctionSpace.h 2009-11-30 01:20:26 +0000
@@ -28,7 +28,6 @@
class DofMap;
class Function;
class GenericFunction;
- class IntersectionDetector;
class GenericVector;
template <class T> class MeshFunction;
=== modified file 'dolfin/mesh/Cell.h'
--- dolfin/mesh/Cell.h 2009-10-02 15:00:52 +0000
+++ dolfin/mesh/Cell.h 2009-11-30 01:20:26 +0000
@@ -2,9 +2,10 @@
// Licensed under the GNU LGPL Version 2.1.
//
// Modified by Johan Hoffman 2006.
+// Modified by Andre Massing 2009.
//
// First added: 2006-06-01
-// Last changed: 2009-08-06
+// Last changed: 2009-11-22
#ifndef __CELL_H
#define __CELL_H
@@ -24,74 +25,81 @@
public:
/// Constructor
+ Cell() : MeshEntity() {}
Cell(const Mesh& mesh, uint index) : MeshEntity(mesh, mesh.topology().dim(), index) {}
/// Destructor
~Cell() {}
/// Return type of cell
- inline CellType::Type type() const { return _mesh.type().cell_type(); }
+ inline CellType::Type type() const { return _mesh->type().cell_type(); }
/// Compute orientation of cell (0 is right, 1 is left)
inline double orientation() const
- { return _mesh.type().orientation(*this); }
+ { return _mesh->type().orientation(*this); }
/// Compute (generalized) volume of cell
inline double volume() const
- { return _mesh.type().volume(*this); }
+ { return _mesh->type().volume(*this); }
/// Compute diameter of cell
inline double diameter() const
- { return _mesh.type().diameter(*this); }
+ { return _mesh->type().diameter(*this); }
/// Compute midpoint of cell
Point midpoint() const;
/// Compute component i of normal of given facet with respect to the cell
inline double normal(uint facet, uint i) const
- { return _mesh.type().normal(*this, facet, i); }
+ { return _mesh->type().normal(*this, facet, i); }
/// Compute normal of given facet with respect to the cell
inline Point normal(uint facet) const
- { return _mesh.type().normal(*this, facet); }
+ { return _mesh->type().normal(*this, facet); }
/// Compute the area/length of given facet with respect to the cell
inline double facet_area(uint facet) const
- { return _mesh.type().facet_area(*this, facet); }
+ { return _mesh->type().facet_area(*this, facet); }
/// Order entities locally
inline void order(MeshFunction<uint>* global_vertex_indices)
- { _mesh.type().order(*this, global_vertex_indices); }
+ { _mesh->type().order(*this, global_vertex_indices); }
/// Check if entities are ordered
inline bool ordered(MeshFunction<uint>* global_vertex_indices)
- { return _mesh.type().ordered(*this, global_vertex_indices); }
+ { return _mesh->type().ordered(*this, global_vertex_indices); }
/// Check for intersection with point
inline bool intersects(const Point& p) const
- { return _mesh.type().intersects(*this, p); }
+ { return _mesh->type().intersects(*this, p); }
/// Check for intersection with line defined by points
inline bool intersects(const Point& p0, const Point& p1) const
- { return _mesh.type().intersects(*this, p0, p1); }
+ { return _mesh->type().intersects(*this, p0, p1); }
/// Check for intersection with cell
inline bool intersects(const Cell& cell)
- { return _mesh.type().intersects(*this, cell); }
+ { return _mesh->type().intersects(*this, cell); }
};
/// A CellIterator is a MeshEntityIterator of topological codimension 0.
-
class CellIterator : public MeshEntityIterator
{
public:
+ CellIterator() : MeshEntityIterator() {}
CellIterator(const Mesh& mesh) : MeshEntityIterator(mesh, mesh.topology().dim()) {}
CellIterator(const MeshEntity& entity) : MeshEntityIterator(entity, entity.mesh().topology().dim()) {}
+// CellIterator(const MeshEntityIterator& it):
+
inline Cell& operator*() { return *operator->(); }
inline Cell* operator->() { return static_cast<Cell*>(MeshEntityIterator::operator->()); }
+
+ //Can I do this. That looks quite dirty :)
+// CellIterator end_iterator() {return *static_cast<CellIterator*>(&MeshEntityIterator::end_iterator());}
+
};
}
=== modified file 'dolfin/mesh/CellType.h'
--- dolfin/mesh/CellType.h 2009-09-08 13:03:05 +0000
+++ dolfin/mesh/CellType.h 2009-11-30 01:20:26 +0000
@@ -2,9 +2,10 @@
// Licensed under the GNU LGPL Version 2.1.
//
// Modified by Kristoffer Selim, 2008.
+// Modified by Andre Massing, 2009.
//
// First added: 2006-06-05
-// Last changed: 2009-09-08
+// Last changed: 2009-11-10
#ifndef __CELL_TYPE_H
#define __CELL_TYPE_H
=== added file 'dolfin/mesh/Changesets_description.txt'
--- dolfin/mesh/Changesets_description.txt 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/Changesets_description.txt 2009-11-30 01:20:26 +0000
@@ -0,0 +1,31 @@
+Added:
+
+Mesh_d class: Template class to wrap the dimension independent mesh class into a
+dimension dependent type.
+
+CellBoostIterator class: More STL compliant iterator over of cells of mesh.
+
+IntersectionOperator class: Provides the interface for the new intersection functionality based
+on the CGAL library (computational geometry). Implementation via the pImpl
+idiom.
+
+IntersectionOperatorImplementation class: Abstract base class for the pImpl
+implementation. Only inherited classes contain the search tree, which type depends
+is based on the geometric primitive it includes.
+
+IntersectionOperatorImplementation_d class: Contains search tree for the mesh
+and executes search operations.
+
+Changed:
+
+MeshEntity class:
+- Introduced default constructor (requires to change private member const
+Mesh & _mesh to const Mesh * _mesh)
+- added comparison operator
+
+MeshEntityInterator class:
+- Introduced default constructor
+- added comparison operator
+- added set_end function to create iterators which points past the end.
+- synchronize immediately internal setting of position with entity index
+ (otherwise ugly bugs could occur by using the comparison operator)
=== modified file 'dolfin/mesh/Edge.cpp'
--- dolfin/mesh/Edge.cpp 2009-07-10 12:16:02 +0000
+++ dolfin/mesh/Edge.cpp 2009-11-30 01:20:26 +0000
@@ -19,8 +19,8 @@
const uint* vertices = entities(0);
assert(vertices);
- const Vertex v0(_mesh, vertices[0]);
- const Vertex v1(_mesh, vertices[1]);
+ const Vertex v0(*_mesh, vertices[0]);
+ const Vertex v1(*_mesh, vertices[1]);
const Point p0 = v0.point();
const Point p1 = v1.point();
@@ -37,8 +37,8 @@
const uint* vertices = entities(0);
assert(vertices);
- const Vertex v0(_mesh, vertices[0]);
- const Vertex v1(_mesh, vertices[1]);
+ const Vertex v0(*_mesh, vertices[0]);
+ const Vertex v1(*_mesh, vertices[1]);
const Point p0 = v0.point();
const Point p1 = v1.point();
=== added file 'dolfin/mesh/IntersectionOperator.cpp'
--- dolfin/mesh/IntersectionOperator.cpp 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/IntersectionOperator.cpp 2009-11-30 01:20:26 +0000
@@ -0,0 +1,142 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-01
+// Last changed: 2009-11-28
+
+#include <algorithm>
+#include <map>
+#include <string>
+
+#include <dolfin/log/dolfin_log.h>
+#include <dolfin/common/types.h>
+#include <dolfin/common/NoDeleter.h>
+
+#include "Mesh.h"
+#include "Edge.h"
+#include "Facet.h"
+#include "Vertex.h"
+#include "Cell.h"
+#include "IntersectionOperator.h"
+#include "IntersectionOperatorImplementation.h"
+#include "MeshPrimitive.h"
+#include "Primitive_Traits.h"
+
+using namespace dolfin;
+
+#ifdef HAS_CGAL
+
+IntersectionOperator::IntersectionOperator(const Mesh & mesh, const std::string & kernel_type)
+ : _mesh(reference_to_no_delete_pointer(mesh)),_kernel_type(kernel_type) {}
+
+IntersectionOperator::IntersectionOperator(boost::shared_ptr<const Mesh> mesh, const std::string & kernel_type)
+ : _mesh(mesh), _kernel_type(kernel_type) {}
+
+IntersectionOperator::~IntersectionOperator() {}
+
+void IntersectionOperator::all_intersected_entities(const Point & point, uint_set & ids_result) const
+{ rImpl().all_intersected_entities(point,ids_result);}
+
+void IntersectionOperator::all_intersected_entities(const std::vector<Point> & points, uint_set & ids_result) const
+{ rImpl().all_intersected_entities(points,ids_result);}
+
+void IntersectionOperator::all_intersected_entities(const Mesh & another_mesh, uint_set & ids_result) const
+{ rImpl().all_intersected_entities(another_mesh, ids_result);}
+
+int IntersectionOperator::any_intersected_entity(const Point & point) const
+{ return rImpl().any_intersected_entity(point);}
+
+void IntersectionOperator::reset_kernel(const std::string & kernel_type)
+{ _pImpl.reset(create_intersection_operator(_mesh,kernel_type)); }
+
+void IntersectionOperator::clear()
+{ _pImpl.reset(); }
+
+const Mesh& IntersectionOperator::mesh() const
+{
+ assert(_mesh);
+ return *_mesh;
+}
+
+boost::shared_ptr<const Mesh> IntersectionOperator::mesh_ptr()
+{
+ assert(_mesh);
+ return _mesh;
+}
+
+IntersectionOperatorImplementation * IntersectionOperator::create_intersection_operator(boost::shared_ptr<const Mesh> mesh, const std::string & kernel_type = "SimpleCartesian")
+{
+ if (kernel_type == "ExactPredicates")
+ {
+ switch( mesh->type().cell_type())
+ {
+ case CellType::point : return new IntersectionOperatorImplementation_d< Primitive_Traits<PointCell, EPICK> >(mesh);
+ case CellType::interval : return new IntersectionOperatorImplementation_d< Primitive_Traits<IntervalCell, EPICK> >(mesh);
+ case CellType::triangle : return new IntersectionOperatorImplementation_d< Primitive_Traits<TriangleCell, EPICK> >(mesh);
+ case CellType::tetrahedron : return new IntersectionOperatorImplementation_d< Primitive_Traits<TetrahedronCell, EPICK> >(mesh);
+ default: error("DOLFIN IntersectionOperator::create_intersection_operator: \n Mesh CellType is not known."); return 0;
+ }
+ }
+ //default is SimpleCartesion
+ else
+ {
+// if (kernel_type != "SimpleCartesian")
+// warning("Type %s of geometry kernel is not known. Falling back to SimpleCartesian.",kernel_type);
+ switch( mesh->type().cell_type())
+ {
+ case CellType::point : return new IntersectionOperatorImplementation_d< Primitive_Traits<PointCell, SCK > >(mesh);
+ case CellType::interval : return new IntersectionOperatorImplementation_d< Primitive_Traits<IntervalCell, SCK > >(mesh);
+ case CellType::triangle : return new IntersectionOperatorImplementation_d< Primitive_Traits<TriangleCell, SCK> >(mesh);
+ case CellType::tetrahedron : return new IntersectionOperatorImplementation_d< Primitive_Traits<TetrahedronCell, SCK > >(mesh);
+ default: error("DOLFIN IntersectionOperator::create_intersection_operator: \n Mesh CellType is not known."); return 0;
+ }
+ }
+}
+
+const IntersectionOperatorImplementation& IntersectionOperator::rImpl() const
+{
+ if (!_pImpl)
+ _pImpl.reset(const_cast<IntersectionOperator *>(this)->create_intersection_operator(_mesh,_kernel_type));
+ return *_pImpl;
+}
+
+#else
+
+IntersectionOperator::IntersectionOperator(const Mesh & _mesh)
+{
+ error("DOLFIN has been compiled without CGAL, IntersectionOperator is not available.");
+}
+
+IntersectionOperator::IntersectionOperator(boost::shared_ptr<const Mesh> _mesh, std::string &)
+{
+ error("DOLFIN has been compiled without CGAL, IntersectionOperator is not available.");
+}
+
+IntersectionOperator::~IntersectionOperator() {}
+
+void IntersectionOperator::all_intersected_entities(const Point & point, uint_set & ids_result) {}
+
+void IntersectionOperator::all_intersected_entities(const std::vector<Point> & points, uint_set & ids_result) {}
+
+void IntersectionOperator::all_intersected_entities(const Mesh & another_mesh, uint_set & ids_result) {}
+
+int IntersectionOperator::any_intersected_entity(const Point & point) {return -1;}
+
+void IntersectionOperator::clear() {}
+
+const Mesh& IntersectionOperator::mesh() const
+{
+ assert(_mesh);
+ return *_mesh;
+}
+
+boost::shared_ptr<Mesh> IntersectionOperator::mesh_ptr()
+{
+ assert(_mesh);
+ return _mesh;
+}
+
+void IntersectionOperator::reset_kernel(const std::string & kernel_type) {}
+void IntersectionOperator::reset() {}
+
+#endif
=== added file 'dolfin/mesh/IntersectionOperator.h'
--- dolfin/mesh/IntersectionOperator.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/IntersectionOperator.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,98 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-01
+// Last changed: 2009-11-27
+
+#ifndef __INTERSECTIONOPERATOR_H
+#define __INTERSECTIONOPERATOR_H
+
+#include <vector>
+#include <string>
+
+#include <boost/scoped_ptr.hpp>
+#include <boost/shared_ptr.hpp>
+
+#include <dolfin/common/types.h>
+
+namespace dolfin
+{
+
+// typedef std::set<uint> uint_set;
+ // Forward declarations
+ class MeshEntity;
+ class Mesh;
+ class Point;
+ class IntersectionOperatorImplementation;
+
+ class IntersectionOperator
+ {
+ public:
+
+ ///Create intersection detector for the mesh \em mesh.
+ ///@param kernel_type The CGAL geometric kernel is used to compute predicates,
+ ///intersections and such. Depending on this choice the kernel
+ ///(kernel_type = "ExcactPredicates") can compute predicates excactly
+ ///(without roundoff error) or only approximately (default, kernel_type =
+ ///"SimpleCartesian").
+ IntersectionOperator(const Mesh & _mesh,const std::string & kernel_type = "SimpleCartesian");
+ IntersectionOperator(boost::shared_ptr<const Mesh> _mesh, const std::string & kernel_type = "SimpleCartesian");
+
+ ///Destructor. Needed be explicit written, otherwise default inline here, with prohibits
+ ///pImpl with scoped_ptr.
+ ~IntersectionOperator();
+
+ ///Compute all id of all cells which are intersects by a \em point.
+ ///\param[out] ids_result The ids of the intersected entities are saved in a set for efficienty
+ ///reasons, to avoid to sort out duplicates later on.
+ void all_intersected_entities(const Point & point, uint_set & ids_result) const ;
+
+ ///Compute all id of all cells which are intersects any point in \em points.
+ ///\param[out] ids_result The ids of the intersected entities are saved in a set for efficienty
+ ///reasons, to avoid to sort out duplicates later on.
+ void all_intersected_entities(const std::vector<Point> & points, uint_set & ids_result) const;
+
+ ///Compute all id of all cells which are intersects by the given mesh \em another_mesh;
+ ///\param[out] ids_result The ids of the intersected entities are saved in a set for efficienty
+ ///reasons, to avoid to sort out duplicates later on.
+ void all_intersected_entities(const Mesh & another_mesh, uint_set & ids_result) const;
+
+ ///Computes only the first id of the entity, which contains the point. Returns -1 if no cell is intersected.
+ ///@internal @remark This makes the function evaluation significantly faster.
+ int any_intersected_entity(const Point & point) const;
+
+ ///Rebuilds the underlying search structure from scratch and uses the kernel kernel_type
+ ///underlying CGAL Geometry kernel.
+ void reset_kernel(const std::string & kernel_type = "SimpleCartesian");
+
+ ///clears search structure. Should be used if the mesh has changed
+ void clear();
+
+ const Mesh& mesh() const;
+ boost::shared_ptr<const Mesh> mesh_ptr();
+
+ private:
+
+// IntersectionOperatorImplementation& rImpl();
+
+ ///Helper function to introduce lazy initialization.
+ const IntersectionOperatorImplementation& rImpl() const;
+
+ ///Factory function to create the dimension dependent intersectionoperator
+ ///implementation.
+ IntersectionOperatorImplementation * create_intersection_operator(boost::shared_ptr<const Mesh> mesh,const std::string & kernel_type);
+
+ ///Pointer to implementation. Mutable to enable lazy initialization.
+ mutable boost::scoped_ptr<IntersectionOperatorImplementation> _pImpl;
+
+ ///Pointer to mesh.
+ boost::shared_ptr<const Mesh> _mesh;
+
+ ///String description of the used geometry kernel.
+ std::string _kernel_type;
+
+
+ };
+}
+
+#endif
=== added file 'dolfin/mesh/IntersectionOperatorImplementation.h'
--- dolfin/mesh/IntersectionOperatorImplementation.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/IntersectionOperatorImplementation.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,176 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-11
+// Last changed: 2009-11-29
+
+#ifndef __INTERSECTIONOPERATORIMPLEMENTATION_H
+#define __INTERSECTIONOPERATORIMPLEMENTATION_H
+
+#ifdef HAS_CGAL
+
+#include <vector>
+
+#include <boost/scoped_ptr.hpp>
+#include <boost/shared_ptr.hpp>
+#include <boost/bimap.hpp>
+#include <boost/optional.hpp>
+
+#include <dolfin/common/types.h>
+
+#include "added_intersection_3.h" //additional intersection functionality, *Must* include before the AABB_tree!
+
+#include <CGAL/AABB_tree.h> // *Must* be inserted before kernel!
+#include <CGAL/AABB_traits.h>
+
+#include <CGAL/Simple_cartesian.h>
+#include "Triangle_3_Tetrahedron_3_do_intersect_SCK.h" //template specialization for Simple_cartesian kernel
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Bbox_3.h>
+#include <CGAL/Point_3.h>
+
+
+#include "Primitive_Traits.h"
+#include "MeshPrimitive.h"
+
+
+typedef CGAL::Simple_cartesian<double> SCK;
+typedef CGAL::Exact_predicates_inexact_constructions_kernel EPICK;
+
+namespace dolfin
+{
+
+ class Point;
+ class Mesh;
+
+ ///@brief Interface class for the actual implementation of the intersection operations.
+ ///
+ ///@internal This is neccessary since the search tree has a dimension dependent type, hence encapsulates in the
+ ///inheriting IntersectionOperatorImplementation_d! It provides the glue level between the dimension independent implementation
+ ///of the mesh class and the dimension dependent search structures in CGAL.
+ class IntersectionOperatorImplementation
+ {
+ public:
+ //Only default constructor, since the search tree has a dimension dependent type, hence encapsulates in the
+ //inheriting IntersectionOperatorImplementation_d!
+
+ virtual void all_intersected_entities(const Point & point, uint_set & ids_result) const = 0;
+ virtual void all_intersected_entities(const std::vector<Point> & points, uint_set & ids_result) const = 0;
+ virtual void all_intersected_entities(const Mesh & another_mesh, uint_set & ids_result) const = 0;
+ virtual int any_intersected_entity(const Point & point) const = 0;
+
+ };
+
+ template <class PrimitiveTrait>
+ class IntersectionOperatorImplementation_d : public IntersectionOperatorImplementation
+ {
+ typedef PrimitiveTrait PT;
+ typedef typename PT::K K;
+ typedef MeshPrimitive<PrimitiveTrait> CellPrimitive;
+
+// typedef Tree_Traits<K,CellPrimitive> AABB_PrimitiveTraits;
+ typedef CGAL::AABB_traits<K,CellPrimitive> AABB_PrimitiveTraits;
+ typedef CGAL::AABB_tree<AABB_PrimitiveTraits> Tree;
+
+ public:
+ ///Constructor.
+ IntersectionOperatorImplementation_d(boost::shared_ptr<const Mesh> _mesh) : _mesh(_mesh)
+ {
+ build_tree();
+ }
+
+ virtual void all_intersected_entities(const Point & point, uint_set & ids_result) const;
+ virtual void all_intersected_entities(const std::vector<Point> & points, uint_set & ids_result) const;
+ virtual void all_intersected_entities(const Mesh & another_mesh, uint_set & ids_result) const;
+
+ virtual int any_intersected_entity(const Point & point) const;
+
+ ///Topological dimension of the mesh.
+ static const uint dim = PrimitiveTrait::dim;
+
+ private:
+
+ void build_tree();
+ boost::shared_ptr<const Mesh> _mesh;
+ boost::scoped_ptr<Tree> tree;
+
+ };
+
+ template <class PT>
+ void IntersectionOperatorImplementation_d<PT>::all_intersected_entities(const Point & point, uint_set & ids_result) const
+ {
+ //@remark For a set the start iterator required by the insert_iterator constructor does not really matter.
+ std::insert_iterator< uint_set > output_it(ids_result,ids_result.end());
+ tree->all_intersected_primitives(Primitive_Traits<PointPrimitive,K>::datum(point), output_it);
+ }
+
+ template <class PT>
+ void IntersectionOperatorImplementation_d<PT>::all_intersected_entities(const std::vector<Point> & points, uint_set & ids_result) const
+ {
+ //@remark For a set the start iterator required by the insert_iterator constructor does not really matter.
+ std::insert_iterator< uint_set > output_it(ids_result,ids_result.end());
+ for (std::vector<Point>::const_iterator p = points.begin(); p != points.end(); ++p)
+ {
+ tree->all_intersected_primitives(Primitive_Traits<PointPrimitive,K>::datum(*p), output_it);
+ }
+ }
+
+// template< template <class PT>
+// void IntersectionOperatorImplementation_d<PT>::all_intersected_entities(const Point & point, uint_set & ids_result) const
+
+ template <class PT>
+ void IntersectionOperatorImplementation_d<PT>::all_intersected_entities(const Mesh & another_mesh, uint_set & ids_result) const
+ {
+ typedef typename PT::K K;
+ //Avoid instantiation of an insert_iterator for each cell.
+ std::insert_iterator< uint_set > output_it(ids_result,ids_result.end());
+ switch( another_mesh.type().cell_type())
+ {
+ case CellType::point :
+ for (CellIterator cell(another_mesh); !cell.end(); ++cell)
+ tree->all_intersected_primitives(Primitive_Traits<PointCell,K>::datum(*cell),output_it); break;
+ case CellType::interval :
+ if (dim == 1 || dim == 3)
+ dolfin_not_implemented();
+ else
+ for (CellIterator cell(another_mesh); !cell.end(); ++cell)
+ tree->all_intersected_primitives(Primitive_Traits<IntervalCell,K>::datum(*cell),output_it); break;
+ case CellType::triangle :
+ for (CellIterator cell(another_mesh); !cell.end(); ++cell)
+ tree->all_intersected_primitives(Primitive_Traits<TriangleCell,K>::datum(*cell),output_it); break;
+ case CellType::tetrahedron :
+ for (CellIterator cell(another_mesh); !cell.end(); ++cell)
+ {
+ tree->all_intersected_primitives(Primitive_Traits<TetrahedronCell,K>::datum(*cell),output_it);
+ }
+ break;
+ default: error("DOLFIN IntersectionOperatorImplementation::all_intersected_entities: \n Mesh CellType is not known.");
+ }
+ }
+
+ template <class PT>
+ int IntersectionOperatorImplementation_d<PT>::any_intersected_entity(const Point & point) const
+ {
+ boost::optional<uint> id = tree->any_intersected_primitive(Primitive_Traits<PointPrimitive,K>::datum(point));
+ if (id)
+ return *id;
+ else
+ return -1;
+ }
+
+ template <class PT>
+ void IntersectionOperatorImplementation_d<PT>::build_tree()
+ {
+ if (_mesh)
+ {
+ MeshEntityIterator cell_iter(*_mesh,dim);
+ tree.reset(new Tree(cell_iter,cell_iter.end_iterator()));
+ }
+ }
+
+}
+
+#endif
+#endif
=== modified file 'dolfin/mesh/Mesh.cpp'
--- dolfin/mesh/Mesh.cpp 2009-11-11 19:45:40 +0000
+++ dolfin/mesh/Mesh.cpp 2009-11-30 01:20:26 +0000
@@ -5,9 +5,10 @@
// Modified by Garth N. Wells 2007.
// Modified by Niclas Jansson 2008.
// Modified by Kristoffer Selim 2008.
+// Modified by Andre Massing, 2009.
//
// First added: 2006-05-09
-// Last changed: 2009-11-11
+// Last changed: 2009-11-27
#include <sstream>
@@ -22,7 +23,7 @@
#include "UniformMeshRefinement.h"
#include "LocalMeshRefinement.h"
#include "LocalMeshCoarsening.h"
-#include "IntersectionDetector.h"
+#include "IntersectionOperator.h"
#include "TopologyComputation.h"
#include "MeshSmoothing.h"
#include "MeshOrdering.h"
@@ -39,21 +40,21 @@
//-----------------------------------------------------------------------------
Mesh::Mesh()
: Variable("mesh", "DOLFIN mesh"),
- _data(*this), _cell_type(0), detector(0), _ordered(false)
+ _data(*this), _cell_type(0), _intersection_operator(*this), _ordered(false)
{
// Do nothing
}
//-----------------------------------------------------------------------------
Mesh::Mesh(const Mesh& mesh)
: Variable("mesh", "DOLFIN mesh"),
- _data(*this), _cell_type(0), detector(0), _ordered(false)
+ _data(*this), _cell_type(0), _intersection_operator(*this), _ordered(false)
{
*this = mesh;
}
//-----------------------------------------------------------------------------
Mesh::Mesh(std::string filename)
: Variable("mesh", "DOLFIN mesh"),
- _data(*this), _cell_type(0), detector(0), _ordered(false)
+ _data(*this), _cell_type(0), _intersection_operator(*this), _ordered(false)
{
if (MPI::num_processes() > 1)
{
@@ -94,7 +95,7 @@
_ordered = mesh._ordered;
- return *this;
+ return *this;
}
//-----------------------------------------------------------------------------
dolfin::uint Mesh::init(uint dim) const
@@ -166,8 +167,7 @@
_data.clear();
delete _cell_type;
_cell_type = 0;
- delete detector;
- detector = 0;
+ _intersection_operator.clear();
_ordered = false;
}
//-----------------------------------------------------------------------------
@@ -247,84 +247,34 @@
MeshSmoothing::smooth(*this);
}
//-----------------------------------------------------------------------------
-void Mesh::intersection(const Point& p, std::vector<uint>& cells, bool fixed_mesh)
-{
- // Don't reuse detector if mesh has moved
- if (!fixed_mesh)
- {
- delete detector;
- detector = 0;
- }
-
- // Create detector if necessary
- if (!detector)
- detector = new IntersectionDetector(*this);
-
- detector->intersection(p, cells);
-}
-//-----------------------------------------------------------------------------
-void Mesh::intersection(const Point& p1, const Point& p2, std::vector<uint>& cells, bool fixed_mesh)
-{
- // Don't reuse detector if the mesh has moved
- if (!fixed_mesh)
- {
- delete detector;
- detector = 0;
- }
-
- // Create detector if necessary
- if (!detector)
- detector = new IntersectionDetector(*this);
-
- detector->intersection(p1, p2, cells);
-}
-//-----------------------------------------------------------------------------
-void Mesh::intersection(Cell& cell, std::vector<uint>& cells, bool fixed_mesh)
-{
- // Don't reuse detector if the has has moved
- if (!fixed_mesh)
- {
- delete detector;
- detector = 0;
- }
-
- // Create detector if necessary
- if (!detector)
- detector = new IntersectionDetector(*this);
-
- detector->intersection(cell, cells);
-}
-//-----------------------------------------------------------------------------
-void Mesh::intersection(std::vector<Point>& points, std::vector<uint>& intersection, bool fixed_mesh)
-{
- // Don't reuse detector if the mesh has moved
- if (!fixed_mesh)
- {
- delete detector;
- detector = 0;
- }
-
- // Create detector if necessary
- if (!detector)
- detector = new IntersectionDetector(*this);
-
- detector->intersection(points, intersection);
-}
-//-----------------------------------------------------------------------------
-void Mesh::intersection(Mesh& mesh, std::vector<uint>& intersection, bool fixed_mesh)
-{
- // Don't reuse detector if the mesh has moved
- if (!fixed_mesh)
- {
- delete detector;
- detector = 0;
- }
-
- // Create detector if necessary
- if (!detector)
- detector = new IntersectionDetector(*this);
-
- detector->intersection(mesh, intersection);
+void Mesh::all_intersected_entities(const Point & point, uint_set & ids_result) const
+{
+ _intersection_operator.all_intersected_entities(point, ids_result);
+}
+//-----------------------------------------------------------------------------
+void Mesh::all_intersected_entities(const std::vector<Point> & points, uint_set & ids_result) const
+{
+ _intersection_operator.all_intersected_entities(points, ids_result);
+}
+//-----------------------------------------------------------------------------
+void Mesh::all_intersected_entities(const Mesh & another_mesh, uint_set & ids_result) const
+{
+ _intersection_operator.all_intersected_entities(another_mesh, ids_result);
+}
+//-----------------------------------------------------------------------------
+int Mesh::any_intersected_entity(const Point & point) const
+{
+ return _intersection_operator.any_intersected_entity(point);
+}
+//-----------------------------------------------------------------------------
+IntersectionOperator& Mesh::intersection_operator()
+{
+ return _intersection_operator;
+}
+//-----------------------------------------------------------------------------
+const IntersectionOperator& Mesh::intersection_operator() const
+{
+ return _intersection_operator;
}
//-----------------------------------------------------------------------------
double Mesh::hmin() const
=== modified file 'dolfin/mesh/Mesh.h'
--- dolfin/mesh/Mesh.h 2009-11-11 19:45:40 +0000
+++ dolfin/mesh/Mesh.h 2009-11-30 01:20:26 +0000
@@ -6,27 +6,30 @@
// Modified by Garth N. Wells, 2007.
// Modified by Niclas Jansson, 2008.
// Modified by Kristoffer Selim, 2008.
+// Modified by Andre Massing, 2009.
//
// First added: 2006-05-08
-// Last changed: 2009-11-11
+// Last changed: 2009-11-27
#ifndef __MESH_H
#define __MESH_H
#include <string>
+
#include <dolfin/common/types.h>
#include <dolfin/common/Variable.h>
#include <dolfin/ale/ALEType.h>
#include "MeshTopology.h"
#include "MeshGeometry.h"
#include "MeshData.h"
+#include "IntersectionOperator.h"
#include "CellType.h"
namespace dolfin
{
-
+
template <class T> class MeshFunction;
- class IntersectionDetector;
+ class IntersectionOperator;
class Function;
class BoundaryMesh;
class XMLMesh;
@@ -119,6 +122,12 @@
/// Return mesh geometry (const version)
inline const MeshGeometry& geometry() const { return _geometry; }
+ ///Return intersectionoperator (const version);
+ const IntersectionOperator& intersection_operator() const;
+
+ ///Return intersectionoperator (non-const version);
+ IntersectionOperator& intersection_operator();
+
/// Return mesh data (non-const version)
MeshData& data() { return _data; }
@@ -173,20 +182,24 @@
/// Smooth mesh using Lagrangian mesh smoothing
void smooth(uint num_smoothings=1);
- /// Compute cells intersecting point
- void intersection(const Point& p, std::vector<uint>& cells, bool fixed_mesh=true);
-
- /// Compute cells overlapping line defined by points
- void intersection(const Point& p1, const Point& p2, std::vector<uint>& cells, bool fixed_mesh=true);
-
- /// Compute cells overlapping cell
- void intersection(Cell& cell, std::vector<uint>& cells, bool fixed_mesh=true);
-
- /// Compute intersection with curve defined by points
- void intersection(std::vector<Point>& points, std::vector<uint>& intersection, bool fixed_mesh=true);
-
- /// Compute intersection with mesh
- void intersection(Mesh& mesh, std::vector<uint>& cells, bool fixed_mesh=true);
+ ///Compute all id of all cells which are intersects by a \em point.
+ ///\param[out] ids_result The ids of the intersected entities are saved in a set for efficienty
+ ///reasons, to avoid to sort out duplicates later on.
+ void all_intersected_entities(const Point & point, uint_set & ids_result) const;
+
+ ///Compute all id of all cells which are intersects any point in \em points.
+ ///\param[out] ids_result The ids of the intersected entities are saved in a set for efficienty
+ ///reasons, to avoid to sort out duplicates later on.
+ void all_intersected_entities(const std::vector<Point> & points, uint_set & ids_result) const;
+
+ ///Compute all id of all cells which are intersects by the given mesh \em another_mesh;
+ ///\param[out] ids_result The ids of the intersected entities are saved in a set for efficienty
+ ///reasons, to avoid to sort out duplicates later on.
+ void all_intersected_entities(const Mesh & another_mesh, uint_set & ids_result) const;
+
+ ///Computes only the first id of the entity, which contains the point. Returns -1 if no cell is intersected.
+ ///@internal @remark This makes the function evaluation significantly faster.
+ int any_intersected_entity(const Point & point) const;
/// Compute minimum cell diameter
double hmin() const;
@@ -221,13 +234,12 @@
CellType* _cell_type;
// Intersection detector
- IntersectionDetector* detector;
+ IntersectionOperator _intersection_operator;
// True if mesh has been ordered
mutable bool _ordered;
};
-
}
#endif
=== modified file 'dolfin/mesh/MeshEntity.cpp'
--- dolfin/mesh/MeshEntity.cpp 2009-08-11 12:14:45 +0000
+++ dolfin/mesh/MeshEntity.cpp 2009-11-30 01:20:26 +0000
@@ -1,32 +1,35 @@
// Copyright (C) 2006-2009 Anders Logg.
// Licensed under the GNU LGPL Version 2.1.
//
+// Modified by Andre Massing, 2009.
+//
// First added: 2006-05-11
-// Last changed: 2009-08-11
+// Last changed: 2009-11-16
#include <dolfin/log/dolfin_log.h>
#include "MeshEntity.h"
+#include "Vertex.h"
using namespace dolfin;
//-----------------------------------------------------------------------------
MeshEntity::MeshEntity(const Mesh& mesh, uint dim, uint index)
- : _mesh(mesh), _dim(dim), _index(index)
+ : _mesh(&mesh), _dim(dim), _index(index)
{
// Check index range
- if (index < mesh.num_entities(dim))
+ if (index < _mesh->num_entities(dim))
return;
// Initialize mesh entities
mesh.init(dim);
// Check index range again
- if (index < mesh.num_entities(dim))
+ if (index < _mesh->num_entities(dim))
return;
// Illegal index range
error("Mesh entity index %d out of range [0, %d] for entity of dimension %d.",
- index, mesh.num_entities(dim), dim);
+ index,_mesh->num_entities(dim), dim);
}
//-----------------------------------------------------------------------------
MeshEntity::~MeshEntity()
@@ -37,12 +40,12 @@
bool MeshEntity::incident(const MeshEntity& entity) const
{
// Must be in the same mesh to be incident
- if ( &_mesh != &entity._mesh )
+ if ( _mesh != entity._mesh )
return false;
// Get list of entities for given topological dimension
- const uint* entities = _mesh.topology()(_dim, entity._dim)(_index);
- const uint num_entities = _mesh.topology()(_dim, entity._dim).size(_index);
+ const uint* entities = _mesh->topology()(_dim, entity._dim)(_index);
+ const uint num_entities = _mesh->topology()(_dim, entity._dim).size(_index);
// Check if any entity matches
for (uint i = 0; i < num_entities; ++i)
@@ -56,12 +59,12 @@
dolfin::uint MeshEntity::index(const MeshEntity& entity) const
{
// Must be in the same mesh to be incident
- if ( &_mesh != &entity._mesh )
+ if ( _mesh != entity._mesh )
error("Unable to compute index of given entity defined on a different mesh.");
// Get list of entities for given topological dimension
- const uint* entities = _mesh.topology()(_dim, entity._dim)(_index);
- const uint num_entities = _mesh.topology()(_dim, entity._dim).size(_index);
+ const uint* entities = _mesh->topology()(_dim, entity._dim)(_index);
+ const uint num_entities = _mesh->topology()(_dim, entity._dim).size(_index);
// Check if any entity matches
for (uint i = 0; i < num_entities; ++i)
@@ -74,6 +77,19 @@
return 0;
}
//-----------------------------------------------------------------------------
+
+#ifdef HAS_CGAL
+template <typename K>
+CGAL::Bbox_3 MeshEntity::bbox () const
+{
+ VertexIterator v(*this);
+ CGAL::Bbox_3 box(v->point().bbox<K>());
+ for (++v; !v.end(); ++v)
+ box = box + v->point().bbox<K>();
+ return box;
+}
+#endif
+
std::string MeshEntity::str(bool verbose) const
{
if (verbose)
=== modified file 'dolfin/mesh/MeshEntity.h'
--- dolfin/mesh/MeshEntity.h 2009-09-08 13:03:05 +0000
+++ dolfin/mesh/MeshEntity.h 2009-11-30 01:20:26 +0000
@@ -1,34 +1,50 @@
// Copyright (C) 2006-2009 Anders Logg.
// Licensed under the GNU LGPL Version 2.1.
//
+// Modified by Andre Massing, 2009.
+//
// First added: 2006-05-11
-// Last changed: 2009-09-08
+// Last changed: 2009-11-16
#ifndef __MESH_ENTITY_H
#define __MESH_ENTITY_H
+#include <iostream>
+
+#ifdef HAS_CGAL
+#include <CGAL/Bbox_3.h>
+#endif
+
#include <dolfin/common/types.h>
#include <dolfin/log/dolfin_log.h>
#include "Mesh.h"
namespace dolfin
{
+
/// A MeshEntity represents a mesh entity associated with
/// a specific topological dimension of some mesh.
-
class MeshEntity
{
public:
+ /// Default Constructor
+ MeshEntity() :_mesh(0), _dim(0), _index(0) {}
+
/// Constructor
MeshEntity(const Mesh& mesh, uint dim, uint index);
/// Destructor
virtual ~MeshEntity();
+ ///Comparision Operator
+ bool operator==(const MeshEntity& another) const { return (_mesh == another._mesh && _dim == another._dim && _index == another._index); }
+
+ bool operator!=(const MeshEntity& another) const { return !operator==(another); }
+
/// Return mesh associated with mesh entity
- inline const Mesh& mesh() const { return _mesh; }
+ inline const Mesh& mesh() const { return *_mesh; }
/// Return topological dimension
inline uint dim() const { return _dim; }
@@ -37,10 +53,10 @@
inline uint index() const { return _index; }
/// Return number of incident mesh entities of given topological dimension
- inline uint num_entities(uint dim) const { return _mesh.topology()(_dim, dim).size(_index); }
+ inline uint num_entities(uint dim) const { return _mesh->topology()(_dim, dim).size(_index); }
/// Return array of indices for incident mesh entitites of given topological dimension
- inline const uint* entities(uint dim) const { return _mesh.topology()(_dim, dim)(_index); }
+ inline const uint* entities(uint dim) const { return _mesh->topology()(_dim, dim)(_index); }
/// Check if given entity is indicent
bool incident(const MeshEntity& entity) const;
@@ -48,8 +64,12 @@
/// Compute local index of given incident entity (error if not found)
uint index(const MeshEntity& entity) const;
+#ifdef HAS_CGAL
+ ///Returns a 3D bounding box of the mesh entity. For lower dimension it may be a degenerated box.
+ template <typename K> CGAL::Bbox_3 bbox() const;
+#endif
+
// Note: Not a subclass of Variable for efficiency!
-
/// Return informal string representation (pretty-print)
std::string str(bool verbose) const;
@@ -58,8 +78,8 @@
// Friends
friend class MeshEntityIterator;
- // The mesh
- const Mesh& _mesh;
+ //The mesh
+ Mesh const * _mesh;
// Topological dimension
uint _dim;
=== modified file 'dolfin/mesh/MeshEntityIterator.h'
--- dolfin/mesh/MeshEntityIterator.h 2009-11-27 06:43:01 +0000
+++ dolfin/mesh/MeshEntityIterator.h 2009-11-30 01:20:26 +0000
@@ -1,8 +1,10 @@
// Copyright (C) 2006-2008 Anders Logg.
// Licensed under the GNU LGPL Version 2.1.
//
+// Modified by Andre Massing, 2009.
+//
// First added: 2006-05-09
-// Last changed: 2009-11-25
+// Last changed: 2009-11-27
#ifndef __MESH_ENTITY_ITERATOR_H
#define __MESH_ENTITY_ITERATOR_H
@@ -42,6 +44,10 @@
class MeshEntityIterator
{
public:
+ ///Default constructor
+ MeshEntityIterator()
+ : _pos(0), pos_end(0), index(0)
+ {}
/// Create iterator for mesh entities over given topological dimension
MeshEntityIterator(const Mesh& mesh, uint dim)
@@ -79,7 +85,10 @@
/// Destructor
virtual ~MeshEntityIterator() {}
- /// Step to next mesh entity (prefix increment)
+ /// Copy Constructor
+ MeshEntityIterator(const MeshEntityIterator& it) : entity(it.entity), _pos(it._pos), pos_end(it.pos_end), index(it.index) {};
+
+ ///Step to next mesh entity (prefix increment)
MeshEntityIterator& operator++() { ++_pos; return *this; }
/// Step to the previous mesh entity (prefix decrease)
@@ -88,8 +97,18 @@
/// Return current position
inline uint pos() const { return _pos; }
- /// Check if iterator has reached the end
- inline bool end() const { return _pos >= pos_end; }
+ ///Comparison operator.
+ ///@internal
+ ///Uncommenting following results into the warning message:
+ //dolfin/mesh/MeshEntityIterator.h:94: Warning|508| Declaration of 'operator ==' shadows declaration accessible via operator->(),
+ //Use const_cast to use operator* inside comparison, which automatically
+ //updates the entity index corresponding to pos *before* comparison (since
+ //update of entity delays until request for entity)
+ bool operator==(const MeshEntityIterator & it) const {return (
+ (const_cast<MeshEntityIterator *>(this))->operator*() == (const_cast<MeshEntityIterator *>(&it))->operator*()
+ && _pos == it._pos && index == it.index);}
+
+ bool operator!=(const MeshEntityIterator & it) const { return !operator==(it); }
/// Dereference operator
inline MeshEntity& operator*() { return *operator->(); }
@@ -97,6 +116,17 @@
/// Member access operator
inline MeshEntity* operator->() { entity._index = (index ? index[_pos] : _pos); return &entity; }
+
+ /// Check if iterator has reached the end
+ inline bool end() const { return _pos >= pos_end; }
+
+ ///Provide a safeguard iterator pointing beyond the end of an iteration
+ ///process, either iterating over the mesh /or incident entities. Added to
+ ///be bit more like STL iteratoren, since many algorithms rely on a kind of
+ ///beyond iterator.
+ MeshEntityIterator end_iterator() {MeshEntityIterator sg(*this); sg.set_end(); return sg;}
+
+
// Note: Not a subclass of Variable for efficiency!
// Commented out to avoid warning about shadowing str() for MeshEntity
/// Return informal string representation (pretty-print)
@@ -114,8 +144,11 @@
/// c1 looks to be an iterator over the entities around c0 when it is in
/// fact a copy of c0.
- MeshEntityIterator(const MeshEntityIterator& entity) : entity(entity.entity.mesh(), 0, 0), _pos(0)
- { error("Illegal use of mesh entity iterator."); }
+// MeshEntityIterator(const MeshEntityIterator& entity) : entity(entity.entity.mesh(), 0, 0), _pos(0)
+// { error("Illegal use of mesh entity iterator."); }
+
+ ///Set pos to end position. To create a kind of mesh.end() iterator.
+ void set_end() { _pos = pos_end; }
// Mesh entity
MeshEntity entity;
=== added file 'dolfin/mesh/MeshPrimitive.h'
--- dolfin/mesh/MeshPrimitive.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/MeshPrimitive.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,81 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-11
+// Last changed: 2009-11-11
+
+#ifndef meshprimitive_INC
+#define meshprimitive_INC
+
+#ifdef HAS_CGAL
+
+#include <CGAL/Bbox_3.h>
+
+#include "Cell.h"
+#include "Mesh.h"
+#include "MeshEntityIterator.h"
+#include "Vertex.h"
+
+#include "Primitive_Traits.h"
+
+namespace dolfin {
+
+template <typename PrimitiveTrait>
+class MeshPrimitive
+{
+public:
+ typedef uint Id;
+ typedef typename PrimitiveTrait::Datum Datum;
+ typedef typename PrimitiveTrait::K K;
+ typedef typename PrimitiveTrait::Primitive Primitive;
+ typedef typename K::Point_3 Point_3;
+
+ //topological dimension of the Cell;
+ static const uint dim = PrimitiveTrait::dim;
+
+ //static function, so that only reference to a mesh and cell index have to
+ //be saved.
+ static Cell getEntity(const MeshPrimitive & p) {return Cell(p.mesh(),p.index);}
+
+ const Mesh& mesh() const { return *_mesh; }
+
+private:
+ Id index; // this is what the AABB tree stores internally
+ const Mesh * _mesh;
+
+public:
+ MeshPrimitive(): index(0), _mesh(0) {} // default constructor needed
+
+ // the following constructor is the one that receives the iterators from the
+ // iterator range given as input to the AABB_tree
+ MeshPrimitive(MeshEntityIterator cell)
+ : index(cell->index()), _mesh(&(cell->mesh())) {}
+
+ Id id() const { return index;}
+
+ ///*Not* required by the CGAL primitive concept, but added for efficieny
+ //reasons. Explanation: We use a modified AABB_tree, in which the local BBox
+ //functor class has been redefined to use the bbox function of dolfin mesh entities.
+ ///Otherwise the bbox function of the Datum object (see below) would have
+ //been used, which means that we would have had to convert dolfin cells into
+ //CGAL primitives only to initialize the tree, which is probably very costly
+ //for 1 million of triangles.
+// CGAL::Bbox_3 bbox () const {return MeshPrimitive<PrimitiveTrait>::getEntity(*this).bbox<K>();}
+
+ ///Provides a reference point required by the Primitive concept of CGAL
+ //AABB_tree. Uses conversion operator in dolfin::Point to create a certain
+ //CGAL Point_3 type.
+ Point_3 reference_point() const
+ {
+ return VertexIterator(MeshPrimitive<PrimitiveTrait>::getEntity(*this))->point();
+ }
+
+ //First line compiles but not the second..?
+ Datum datum() const { return Primitive_Traits<Primitive,K>::datum(MeshPrimitive<PrimitiveTrait>::getEntity(*this));}
+// Datum datum() const { return PrimitiveTrait::datum(MeshPrimitive<PrimitiveTrait>::getEntity(*this));}
+};
+
+}
+
+#endif
+#endif /* ----- #ifndef meshprimitive_INC ----- */
=== modified file 'dolfin/mesh/Point.h'
--- dolfin/mesh/Point.h 2009-10-09 21:02:24 +0000
+++ dolfin/mesh/Point.h 2009-11-30 01:20:26 +0000
@@ -2,21 +2,27 @@
// Licensed under the GNU LGPL Version 2.1.
//
// Modified by Garth N. Wells, 2006.
+// Modified by Andre Massing, 2009.
//
// First added: 2006-06-12
-// Last changed: 2009-09-08
+// Last changed: 2009-11-11
#ifndef __POINT_H
#define __POINT_H
#include <iostream>
#include <dolfin/log/dolfin_log.h>
-#include <dolfin/log/dolfin_log.h>
#include <dolfin/common/types.h>
+#ifdef HAS_CGAL
+// #include "CGAL_includes.h"
+ #include <CGAL/Bbox_3.h>
+ #include <CGAL/Point_3.h>
+// #include <CGAL/Simple_cartesian.h>
+#endif
+
namespace dolfin
{
-
/// A Point represents a point in R^3 with coordinates x, y, z, or,
/// alternatively, a vector in R^3, supporting standard operations
/// like the norm, distances, scalar and vector products etc.
@@ -41,12 +47,10 @@
~Point() {};
/// Return address of coordinate in direction i
- inline double& operator[] (uint i)
- { assert(i < 3); return _x[i]; }
+ inline double& operator[] (uint i) { assert(i < 3); return _x[i]; }
/// Return coordinate in direction i
- inline double operator[] (uint i) const
- { assert(i < 3); return _x[i]; }
+ inline double operator[] (uint i) const { assert(i < 3); return _x[i]; }
/// Return x-coordinate
inline double x() const { return _x[0]; }
@@ -58,44 +62,44 @@
inline double z() const { return _x[2]; }
/// Return coordinate array
- inline const double* coordinates() const
- { return _x; }
+ inline const double* coordinates() const { return _x; }
/// Compute sum of two points
- Point operator+ (const Point& p) const
- { Point q(_x[0] + p._x[0], _x[1] + p._x[1], _x[2] + p._x[2]); return q; }
+ Point operator+ (const Point& p) const { Point q(_x[0] + p._x[0], _x[1] + p._x[1], _x[2] + p._x[2]); return q; }
/// Compute difference of two points
- Point operator- (const Point& p) const
- { Point q(_x[0] - p._x[0], _x[1] - p._x[1], _x[2] - p._x[2]); return q; }
+ Point operator- (const Point& p) const { Point q(_x[0] - p._x[0], _x[1] - p._x[1], _x[2] - p._x[2]); return q; }
/// Add given point
- const Point& operator+= (const Point& p)
- { _x[0] += p._x[0]; _x[1] += p._x[1]; _x[2] += p._x[2]; return *this; }
+ const Point& operator+= (const Point& p) { _x[0] += p._x[0]; _x[1] += p._x[1]; _x[2] += p._x[2]; return *this; }
/// Subtract given point
- const Point& operator-= (const Point& p)
- { _x[0] -= p._x[0]; _x[1] -= p._x[1]; _x[2] -= p._x[2]; return *this; }
+ const Point& operator-= (const Point& p) { _x[0] -= p._x[0]; _x[1] -= p._x[1]; _x[2] -= p._x[2]; return *this; }
/// Multiplication with scalar
- Point operator* (double a) const
- { Point p(a*_x[0], a*_x[1], a*_x[2]); return p; }
+ Point operator* (double a) const { Point p(a*_x[0], a*_x[1], a*_x[2]); return p; }
/// Incremental multiplication with scalar
- const Point& operator*= (double a)
- { _x[0] *= a; _x[1] *= a; _x[2] *= a; return *this; }
+ const Point& operator*= (double a) { _x[0] *= a; _x[1] *= a; _x[2] *= a; return *this; }
/// Division by scalar
- Point operator/ (double a) const
- { Point p(_x[0]/a, _x[1]/a, _x[2]/a); return p; }
+ Point operator/ (double a) const { Point p(_x[0]/a, _x[1]/a, _x[2]/a); return p; }
/// Incremental division by scalar
- const Point& operator/= (double a)
- { _x[0] /= a; _x[1] /= a; _x[2] /= a; return *this; }
+ const Point& operator/= (double a) { _x[0] /= a; _x[1] /= a; _x[2] /= a; return *this; }
/// Assignment operator
- const Point& operator= (const Point& p)
- { _x[0] = p._x[0]; _x[1] = p._x[1]; _x[2] = p._x[2]; return *this; }
+ const Point& operator= (const Point& p) { _x[0] = p._x[0]; _x[1] = p._x[1]; _x[2] = p._x[2]; return *this; }
+
+#ifdef HAS_CGAL
+ ///Conversion operator to appropriate CGAL Point_3 class.
+ template <typename Kernel>
+ operator CGAL::Point_3<Kernel>() const { return CGAL::Point_3<Kernel>(_x[0],_x[1],_x[2]); }
+
+ ///Provides a CGAL bounding box, using conversion operator.
+ template <typename Kernel>
+ CGAL::Bbox_3 bbox() {return CGAL::Point_3<Kernel>(*this).bbox();}
+#endif
/// Compute distance to given point
double distance(const Point& p) const;
=== added file 'dolfin/mesh/Point_3_Bbox_3_intersection.h'
--- dolfin/mesh/Point_3_Bbox_3_intersection.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/Point_3_Bbox_3_intersection.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,107 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-11
+// Last changed: 2009-11-10
+
+#ifndef POINT_3_BBOX_3_INTERSECTION_H
+#define POINT_3_BBOX_3_INTERSECTION_H
+
+#include <CGAL/Bbox_3.h>
+#include <CGAL/Point_3.h>
+#include <CGAL/Object.h>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+ template <class K>
+ inline
+ bool do_intersect(const typename K::Point_3& pt,
+ const CGAL::Bbox_3& bbox,
+ const K &)
+ {
+ for(int i = 0; i < 3; ++i)
+ if(bbox.max(i) < pt[i] || bbox.min(i) > pt[i])
+ return false;
+ return true;
+ }
+
+ template <class K>
+ inline
+ bool do_intersect(const CGAL::Bbox_3& bbox,
+ const typename K::Point_3& pt,
+ const K &)
+ {
+ for(int i = 0; i < 3; ++i)
+ if(bbox.max(i) < pt[i] || bbox.min(i) > pt[i])
+ return false;
+ return true;
+ }
+
+template <class K>
+inline
+Object
+intersection(const typename K::Point_3 &pt,
+ const typename K::Bbox_3 &bbox,
+ const K&)
+{
+ if (do_intersect(pt,bbox)) {
+ return make_object(pt);
+ }
+ return Object();
+}
+
+template <class K>
+inline
+Object
+intersection( const typename K::Bbox_3 &bbox,
+ const typename K::Point_3 &pt,
+ const K&)
+{
+ if (do_intersect(pt,bbox)) {
+ return make_object(pt);
+ }
+ return Object();
+}
+
+
+} //namespace CGALi
+
+template <class K>
+inline
+bool do_intersect(const CGAL::Point_3<K>& point,
+ const CGAL::Bbox_3& bbox)
+{
+ return typename K::Do_intersect_3()(point, bbox);
+// return CGALi::do_intersect(bbox, point, K());
+}
+
+template <class K>
+inline
+bool do_intersect(const CGAL::Bbox_3& bbox,
+ const CGAL::Point_3<K>& point)
+{
+ return typename K::Do_intersect_3()(point, bbox);
+// return CGALi::do_intersect(bbox, point, K());
+}
+
+template <class K>
+inline Object
+intersection(const Bbox_3 & bbox, const Point_3<K> & pt)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(pt, bbox);
+}
+
+template <class K>
+inline Object
+intersection(const Point_3<K> & pt, const Bbox_3 & bbox)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(pt, bbox);
+}
+
+CGAL_END_NAMESPACE
+
+#endif /* ----- #ifndef POINT_3_BBOX_3_INTERSECTION_H ----- */
=== added file 'dolfin/mesh/Point_3_Iso_Cuboid_3_intersection.h'
--- dolfin/mesh/Point_3_Iso_Cuboid_3_intersection.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/Point_3_Iso_Cuboid_3_intersection.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,108 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-11
+// Last changed: 2009-11-10
+
+#ifndef CGAL_POINT_3_ISO_CUBOID_3_INTERSECTION_H
+#define CGAL_POINT_3_ISO_CUBOID_3_INTERSECTION_H
+
+#include <CGAL/Iso_cuboid_3.h>
+#include <CGAL/Point_3.h>
+#include <CGAL/Object.h>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+template <class K>
+inline
+bool
+do_intersect(const typename K::Point_3 &pt,
+ const typename K::Iso_cuboid_3 &iso,
+ const K&)
+{
+ return !iso.has_on_unbounded_side(pt);
+}
+
+template <class K>
+inline
+bool
+do_intersect(const typename K::Iso_cuboid_3 &iso,
+ const typename K::Point_3 &pt,
+ const K&)
+{
+ return !iso.has_on_unbounded_side(pt);
+}
+
+template <class K>
+Object
+intersection(const typename K::Point_3 &pt,
+ const typename K::Iso_cuboid_3 &iso,
+ const K& k)
+{
+ if (CGALi::do_intersect(pt,iso,k)) {
+ return make_object(pt);
+ }
+ return Object();
+}
+
+
+template <class K>
+Object
+intersection(const typename K::Iso_cuboid_3 &iso,
+ const typename K::Point_3 &pt,
+ const K& k)
+{
+ if (CGALi::do_intersect(pt,iso,k)) {
+ return make_object(pt);
+ }
+ return Object();
+}
+
+} // namespace CGALi
+
+
+template <class K>
+inline
+bool
+do_intersect(const Iso_cuboid_3<K> &iso,
+ const Point_3<K> &pt)
+{
+ typedef typename K::Do_intersect_3 Do_intersect;
+ return Do_intersect()(pt, iso);
+}
+
+template <class K>
+inline
+bool
+do_intersect(const Point_3<K> &pt,
+ const Iso_cuboid_3<K> &iso)
+{
+ typedef typename K::Do_intersect_3 Do_intersect;
+ return Do_intersect()(pt, iso);
+}
+
+template <class K>
+inline
+Object
+intersection(const Iso_cuboid_3<K> &iso,
+ const Point_3<K> &pt)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(pt, iso);
+}
+
+template <class K>
+inline
+Object
+intersection(const Point_3<K> &pt,
+ const Iso_cuboid_3<K> &iso)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(pt, iso);
+}
+
+CGAL_END_NAMESPACE
+
+#endif
=== added file 'dolfin/mesh/Point_3_Line_3_intersection.h'
--- dolfin/mesh/Point_3_Line_3_intersection.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/Point_3_Line_3_intersection.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,104 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-11
+// Last changed: 2009-11-10
+
+#ifndef CGAL_POINT_3_LINE_3_INTERSECTION_H
+#define CGAL_POINT_3_LINE_3_INTERSECTION_H
+
+#include <CGAL/Line_3.h>
+#include <CGAL/Point_3.h>
+#include <CGAL/Object.h>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+template <class K>
+inline bool
+do_intersect(const typename K::Point_3 &pt,
+ const typename K::Line_3 &line,
+ const K&)
+{
+ return line.has_on(pt);
+}
+
+template <class K>
+inline bool
+do_intersect(const typename K::Line_3 &line,
+ const typename K::Point_3 &pt,
+ const K&)
+{
+ return line.has_on(pt);
+}
+
+template <class K>
+Object
+intersection(const typename K::Point_3 &pt,
+ const typename K::Line_3 &line,
+ const K& k)
+{
+ if (do_intersect(pt,line, k)) {
+ return make_object(pt);
+ }
+ return Object();
+}
+
+template <class K>
+Object
+intersection(const typename K::Line_3 &line,
+ const typename K::Point_3 &pt,
+ const K& k)
+{
+ if (do_intersect(pt,line, k)) {
+ return make_object(pt);
+ }
+ return Object();
+}
+
+} // namespace CGALi
+
+template <class K>
+inline
+bool
+do_intersect(const Line_3<K> &line,
+ const Point_3<K> &pt)
+{
+ typedef typename K::Do_intersect_3 Do_intersect;
+ return Do_intersect()(pt, line);
+}
+
+template <class K>
+inline
+bool
+do_intersect(const Point_3<K> &pt,
+ const Line_3<K> &line)
+{
+ typedef typename K::Do_intersect_3 Do_intersect;
+ return Do_intersect()(pt, line);
+}
+
+template <class K>
+inline
+Object
+intersection(const Line_3<K> &line,
+ const Point_3<K> &pt)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(pt, line);
+}
+
+template <class K>
+inline
+Object
+intersection(const Point_3<K> &pt,
+ const Line_3<K> &line)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(pt, line);
+}
+
+CGAL_END_NAMESPACE
+
+#endif
=== added file 'dolfin/mesh/Point_3_Point_3_intersection.h'
--- dolfin/mesh/Point_3_Point_3_intersection.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/Point_3_Point_3_intersection.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,63 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-11
+// Last changed: 2009-11-10
+
+#ifndef CGAL_POINT_3_POINT_3_INTERSECTION_H
+#define CGAL_POINT_3_POINT_3_INTERSECTION_H
+
+#include <CGAL/Point_3.h>
+#include <CGAL/Object.h>
+
+///@file This file contains some small extension to the CGAL library, for
+//instance unifying their do_intersect functions to also deal with Point_3
+//and Primitive intersections or some additional intersection collision test.
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+ template <class K >
+ inline bool do_intersect(const typename K::Point_3 & pt1,
+ const typename K::Point_3 & pt2,
+ const K & k)
+ {
+ return pt1 == pt1;
+ }
+
+ template <class K>
+ Object
+ intersection(const typename K::Point_3 &pt1,
+ const typename K::Point_3 &pt2)
+ {
+ if (pt1 == pt2) {
+ return make_object(pt1);
+ }
+ return Object();
+ }
+
+
+}// namespace CGALi
+
+template <class K>
+inline
+bool
+do_intersect(const Point_3<K> &pt1, const Point_3<K> &pt2)
+{
+ typedef typename K::Do_intersect_3 Do_intersect;
+ return Do_intersect()(pt1, pt2);
+}
+
+template <class K>
+inline
+Object
+intersection(const Point_3<K> &pt1, const Point_3<K> &pt2)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(pt1, pt2);
+}
+
+CGAL_END_NAMESPACE
+
+#endif // CGAL_POINT_3_POINT_3_INTERSECTION_H
=== added file 'dolfin/mesh/Point_3_Ray_3_intersection.h'
--- dolfin/mesh/Point_3_Ray_3_intersection.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/Point_3_Ray_3_intersection.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,104 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-11
+// Last changed: 2009-11-10
+
+#ifndef CGAL_POINT_3_RAY_3_INTERSECTION_H
+#define CGAL_POINT_3_RAY_3_INTERSECTION_H
+
+#include <CGAL/Ray_3.h>
+#include <CGAL/Point_3.h>
+#include <CGAL/Object.h>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+template <class K>
+inline
+bool
+do_intersect(const typename K::Point_3 &pt,
+ const typename K::Ray_3 &ray,
+ const K&)
+{
+ return ray.has_on(pt);
+}
+
+
+template <class K>
+inline
+bool
+do_intersect(const typename K::Ray_3 &ray,
+ const typename K::Point_3 &pt,
+ const K&)
+{
+ return ray.has_on(pt);
+}
+
+
+template <class K>
+Object
+intersection(const typename K::Point_3 &pt,
+ const typename K::Ray_3 &ray,
+ const K& k)
+{
+ if (do_intersect(pt,ray, k)) {
+ return make_object(pt);
+ }
+ return Object();
+}
+
+template <class K>
+Object
+intersection(const typename K::Ray_3 &ray,
+ const typename K::Point_3 &pt,
+ const K& k)
+{
+ if (do_intersect(pt,ray, k)) {
+ return make_object(pt);
+ }
+ return Object();
+}
+
+} // namespace CGALi
+
+
+template <class K>
+inline
+bool
+do_intersect(const Ray_3<K> &ray, const Point_3<K> &pt)
+{
+ typedef typename K::Do_intersect_3 Do_intersect;
+ return Do_intersect()(pt, ray);
+}
+
+template <class K>
+inline
+bool
+do_intersect(const Point_3<K> &pt, const Ray_3<K> &ray)
+{
+ typedef typename K::Do_intersect_3 Do_intersect;
+ return Do_intersect()(pt, ray);
+}
+
+
+template <class K>
+inline Object
+intersection(const Ray_3<K> &ray, const Point_3<K> &pt)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(pt, ray);
+}
+
+template <class K>
+inline Object
+intersection(const Point_3<K> &pt, const Ray_3<K> &ray)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(pt, ray);
+}
+
+CGAL_END_NAMESPACE
+
+#endif
=== added file 'dolfin/mesh/Point_3_Segment_3_intersection.h'
--- dolfin/mesh/Point_3_Segment_3_intersection.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/Point_3_Segment_3_intersection.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,103 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-11
+// Last changed: 2009-11-10
+
+#ifndef CGAL_POINT_3_SEGMENT_3_INTERSECTION_H
+#define CGAL_POINT_3_SEGMENT_2_INTERSECTION_H
+
+#include <CGAL/Segment_3.h>
+#include <CGAL/Point_3.h>
+#include <CGAL/Object.h>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+template <class K>
+inline
+bool
+do_intersect(const typename K::Point_3 &pt,
+ const typename K::Segment_3 &seg,
+ const K&)
+{
+ return seg.has_on(pt);
+}
+
+template <class K>
+inline
+bool
+do_intersect(const typename K::Segment_3 &seg,
+ const typename K::Point_3 &pt,
+ const K&)
+{
+ return seg.has_on(pt);
+}
+
+
+template <class K>
+inline
+Object
+intersection(const typename K::Point_3 &pt,
+ const typename K::Segment_3 &seg,
+ const K&)
+{
+ if (do_intersect(pt,seg)) {
+ return make_object(pt);
+ }
+ return Object();
+}
+
+template <class K>
+inline
+Object
+intersection( const typename K::Segment_3 &seg,
+ const typename K::Point_3 &pt,
+ const K&)
+{
+ if (do_intersect(pt,seg)) {
+ return make_object(pt);
+ }
+ return Object();
+}
+
+} // namespace CGALi
+
+
+template <class K>
+inline bool
+do_intersect(const Segment_3<K> &seg, const Point_3<K> &pt)
+{
+ typedef typename K::Do_intersect_3 Do_intersect;
+ return Do_intersect()(pt, seg);
+}
+
+template <class K>
+inline bool
+do_intersect(const Point_3<K> &pt, const Segment_3<K> &seg)
+{
+ typedef typename K::Do_intersect_3 Do_intersect;
+ return Do_intersect()(pt, seg);
+}
+
+
+template <class K>
+inline Object
+intersection(const Segment_3<K> &seg, const Point_3<K> &pt)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(pt, seg);
+}
+
+template <class K>
+inline Object
+intersection(const Point_3<K> &pt, const Segment_3<K> &seg)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(pt, seg);
+}
+
+CGAL_END_NAMESPACE
+
+#endif
=== added file 'dolfin/mesh/Point_3_Tetrahedron_3_intersection.h'
--- dolfin/mesh/Point_3_Tetrahedron_3_intersection.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/Point_3_Tetrahedron_3_intersection.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,104 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-11
+// Last changed: 2009-11-10
+
+#ifndef CGAL_POINT_3_TETRAHEDRON_3_INTERSECTION_H
+#define CGAL_POINT_3_TETRAHEDRON_3_INTERSECTION_H
+
+#include <CGAL/Point_3.h>
+#include <CGAL/Object.h>
+#include <CGAL/Tetrahedron_3.h>
+
+CGAL_BEGIN_NAMESPACE
+
+
+namespace CGALi {
+
+template <class K>
+inline
+bool
+do_intersect(const typename K::Point_3 &pt,
+ const typename K::Tetrahedron_3 &tet,
+ const K&)
+{
+ return !tet.has_on_unbounded_side(pt);
+}
+
+template <class K>
+inline
+bool
+do_intersect(const typename K::Tetrahedron_3 &tet,
+ const typename K::Point_3 &pt,
+ const K&)
+{
+ return !tet.has_on_unbounded_side(pt);
+}
+
+
+template <class K>
+inline
+Object
+intersection(const typename K::Point_3 &pt,
+ const typename K::Tetrahedron_3 &tet,
+ const K&)
+{
+ if (do_intersect(pt,tet)) {
+ return make_object(pt);
+ }
+ return Object();
+}
+
+template <class K>
+inline
+Object
+intersection( const typename K::Tetrahedron_3 &tet,
+ const typename K::Point_3 &pt,
+ const K&)
+{
+ if (do_intersect(pt,tet)) {
+ return make_object(pt);
+ }
+ return Object();
+}
+
+} // namespace CGALi
+
+
+template <class K>
+inline bool
+do_intersect(const Tetrahedron_3<K> &tet, const Point_3<K> &pt)
+{
+ typedef typename K::Do_intersect_3 Do_intersect;
+ return Do_intersect()(pt, tet);
+}
+
+template <class K>
+inline bool
+do_intersect(const Point_3<K> &pt, const Tetrahedron_3<K> &tet)
+{
+ typedef typename K::Do_intersect_3 Do_intersect;
+ return Do_intersect()(pt, tet);
+}
+
+
+template <class K>
+inline Object
+intersection(const Tetrahedron_3<K> &tet, const Point_3<K> &pt)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(pt, tet);
+}
+
+template <class K>
+inline Object
+intersection(const Point_3<K> &pt, const Tetrahedron_3<K> &tet)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(pt, tet);
+}
+
+CGAL_END_NAMESPACE
+
+#endif
=== added file 'dolfin/mesh/Point_3_Triangle_3_intersection.h'
--- dolfin/mesh/Point_3_Triangle_3_intersection.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/Point_3_Triangle_3_intersection.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,109 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-11
+// Last changed: 2009-11-10
+
+#ifndef CGAL_POINT_3_TRIANGLE_3_INTERSECTION_H
+#define CGAL_POINT_3_TRIANGLE_3_INTERSECTION_H
+
+#include <CGAL/Point_3.h>
+#include <CGAL/Triangle_3.h>
+#include <CGAL/Object.h>
+
+//already defined in original
+//CGAL, but not with p.has_on(pt) function??
+
+#ifdef SUBSTITUTE_ORIGINAL_CGAL_FUNCTION
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+template <class K>
+inline
+bool
+do_intersect(const typename K::Point_3 &pt,
+ const typename K::Triangle_3 &tri,
+ const K&)
+{
+ return tri.has_on(pt);
+}
+
+template <class K>
+inline
+bool
+do_intersect(const typename K::Triangle_3 &tri,
+ const typename K::Point_3 &pt,
+ const K&)
+{
+ return tri.has_on(pt);
+}
+
+
+template <class K>
+inline
+Object
+intersection(const typename K::Point_3 &pt,
+ const typename K::Triangle_3 &tri,
+ const K&)
+{
+ if (do_intersect(pt,tri)) {
+ return make_object(pt);
+ }
+ return Object();
+}
+
+template <class K>
+inline
+Object
+intersection( const typename K::Triangle_3 &tri,
+ const typename K::Point_3 &pt,
+ const K&)
+{
+ if (do_intersect(pt,tri)) {
+ return make_object(pt);
+ }
+ return Object();
+}
+
+} // namespace CGALi
+
+
+template <class K>
+inline bool
+do_intersect(const Triangle_3<K> &tri, const Point_3<K> &pt)
+{
+ typedef typename K::Do_intersect_3 Do_intersect;
+ return Do_intersect()(pt, tri);
+}
+
+template <class K>
+inline bool
+do_intersect(const Point_3<K> &pt, const Triangle_3<K> &tri)
+{
+ typedef typename K::Do_intersect_3 Do_intersect;
+ return Do_intersect()(pt, tri);
+}
+
+
+template <class K>
+inline Object
+intersection(const Triangle_3<K> &tri, const Point_3<K> &pt)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(pt, tri);
+}
+
+template <class K>
+inline Object
+intersection(const Point_3<K> &pt, const Triangle_3<K> &tri)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(pt, tri);
+}
+
+CGAL_END_NAMESPACE
+
+#endif
+#endif
=== added file 'dolfin/mesh/Primitive_Traits.h'
--- dolfin/mesh/Primitive_Traits.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/Primitive_Traits.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,112 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-16
+// Last changed: 2009-11-28
+
+#ifndef primitives_traits_INC
+#define primitives_traits_INC
+
+#ifdef HAS_CGAL
+
+#include "Vertex.h"
+#include "MeshEntityIterator.h"
+#include "PointCell.h"
+#include "IntervalCell.h"
+#include "TriangleCell.h"
+#include "TetrahedronCell.h"
+
+namespace dolfin {
+
+struct PointPrimitive {};
+
+//struct PointCellPrimitive {};
+
+//struct IntervalCellPrimitive {};
+
+//struct TriangleCellPrimitive {};
+
+//struct TetrahedronCellPrimitive {};
+
+///Forward declaration for a general Traits class. This traits class is
+///supposed to provide a datum function, which returns a geometric primitive
+///object, which type corresponds to the primitive type (Point, PointCell,
+///Tetrahedron(Cell) etc.) and the passed geometric CGAL kernel.
+template <typename Primitive, typename Kernel> struct Primitive_Traits;
+
+template <typename Kernel> struct Primitive_Traits<PointPrimitive,Kernel> {
+ typedef Kernel K;
+ typedef PointPrimitive Primitive;
+ typedef typename K::Point_3 Datum;
+ static const int dim = 0;
+ static Datum datum(const Point & point) {
+ return Datum(point);
+ }
+};
+
+template <typename Kernel> struct Primitive_Traits<PointCell,Kernel> {
+ typedef Kernel K;
+ typedef PointCell Primitive;
+ typedef typename K::Point_3 Datum;
+ static const int dim = 0;
+ static Datum datum(const Cell & cell) {
+ VertexIterator v(cell);
+ return Datum(v->point());
+ }
+};
+
+template <typename Kernel> struct Primitive_Traits<IntervalCell,Kernel> {
+ typedef Kernel K;
+ typedef IntervalCell Primitive;
+ typedef typename K::Point_3 Point_3;
+ typedef typename K::Segment_3 Datum;
+ static const int dim = 1;
+ static Datum datum(const Cell & cell) {
+ VertexIterator v(cell);
+ Point_3 p1(v->point());
+ ++v;
+ Point_3 p2(v->point());
+ return Datum(p1,p2);
+ }
+};
+
+template <typename Kernel> struct Primitive_Traits<TriangleCell,Kernel> {
+ typedef Kernel K;
+ typedef TriangleCell Primitive;
+ typedef typename K::Point_3 Point_3;
+ typedef typename K::Triangle_3 Datum;
+ static const int dim = 2;
+ static Datum datum(const Cell & cell) {
+ VertexIterator v(cell);
+ Point_3 p1(v->point());
+ ++v;
+ Point_3 p2(v->point());
+ ++v;
+ Point_3 p3(v->point());
+ return Datum(p1,p2,p3);
+ }
+};
+
+template <typename Kernel> struct Primitive_Traits<TetrahedronCell,Kernel> {
+ typedef Kernel K;
+ typedef TetrahedronCell Primitive;
+ typedef typename K::Point_3 Point_3;
+ typedef typename K::Tetrahedron_3 Datum;
+ static const int dim = 3;
+ static Datum datum(const Cell & cell) {
+ VertexIterator v(cell);
+ Point_3 p1(v->point());
+ ++v;
+ Point_3 p2(v->point());
+ ++v;
+ Point_3 p3(v->point());
+ ++v;
+ Point_3 p4(v->point());
+ return Datum(p1,p2,p3,p4);
+ }
+};
+
+} //end namespace dolfin
+
+#endif
+#endif /* ----- #ifndef primitives_traits_INC ----- */
=== added file 'dolfin/mesh/Segment_3_Segment_3_intersection.h'
--- dolfin/mesh/Segment_3_Segment_3_intersection.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/Segment_3_Segment_3_intersection.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,69 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-11
+// Last changed: 2009-11-10
+
+#ifndef segment_3_segment_3_intersection_INC
+#define segment_3_segment_3_intersection_INC
+
+
+#include <CGAL/Segment_3.h>
+#include <CGAL/Object.h>
+
+///@file This file contains some small extension to the CGAL library, for
+//instance unifying their do_intersect functions to also deal with Segment_3
+//and Primitive intersections or some additional intersection collision test.
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+ template <class K >
+ inline bool do_intersect(const typename K::Segment_3 & s1,
+ const typename K::Segment_3 & s2,
+ const K & k)
+ {
+ //throw exception
+ return false;
+ }
+
+template <class K>
+inline
+Object
+intersection(const typename K::Segment_3 &s1,
+ const typename K::Segment_3 &s2,
+ const K&)
+{
+ //throw exception
+ if (do_intersect(s1,s2)) {
+ return Object();
+ }
+ return Object();
+}
+
+
+}// namespace CGALi
+
+template <class K>
+inline
+bool
+do_intersect(const Segment_3<K> &s1, const Segment_3<K> &s2)
+{
+ typedef typename K::Do_intersect_3 Do_intersect;
+ return Do_intersect()(s1, s2);
+}
+
+template <class K>
+inline
+Object
+intersection(const Segment_3<K> &s1, const Segment_3<K> &s2)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(s1, s2);
+}
+
+CGAL_END_NAMESPACE
+
+
+#endif /* ----- #ifndef segment_3_segment_3_intersection_INC ----- */
=== added file 'dolfin/mesh/Segment_3_Tetrahedron_3_intersection.h'
--- dolfin/mesh/Segment_3_Tetrahedron_3_intersection.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/Segment_3_Tetrahedron_3_intersection.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,108 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-11
+// Last changed: 2009-11-10
+
+#ifndef segment_3_tetrahedron_3_intersection_INC
+#define segment_3_tetrahedron_3_intersection_INC
+
+#include <CGAL/Segment_3.h>
+#include <CGAL/Tetrahedron_3.h>
+#include <CGAL/Object.h>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+template <class K>
+inline
+bool
+do_intersect(const typename K::Tetrahedron_3 &tet,
+ const typename K::Segment_3 &seg,
+ const K&)
+{
+ //throw exception!
+ return false;
+}
+
+template <class K>
+inline
+bool
+do_intersect(const typename K::Segment_3 &seg,
+ const typename K::Tetrahedron_3 &tet,
+ const K&)
+{
+ //throw exception!
+ return false;
+}
+
+
+template <class K>
+inline
+Object
+intersection(const typename K::Tetrahedron_3 &tet,
+ const typename K::Segment_3 &seg,
+ const K&)
+{
+ //throw exception! If intersect either point or segment.
+ if (do_intersect(tet,seg)) {
+ return Object();
+ }
+ return Object();
+}
+
+template <class K>
+inline
+Object
+intersection( const typename K::Segment_3 &seg,
+ const typename K::Tetrahedron_3 &tet,
+ const K&)
+{
+ //throw exception!
+ if (do_intersect(tet,seg)) {
+ return Object();
+ }
+ return Object();
+}
+
+} // namespace CGALi
+
+
+template <class K>
+inline bool
+do_intersect(const Segment_3<K> &seg, const Tetrahedron_3<K> &tet)
+{
+ typedef typename K::Do_intersect_3 Do_intersect;
+ return Do_intersect()(tet, seg);
+}
+
+template <class K>
+inline bool
+do_intersect(const Tetrahedron_3<K> &tet, const Segment_3<K> &seg)
+{
+ typedef typename K::Do_intersect_3 Do_intersect;
+ return Do_intersect()(tet, seg);
+}
+
+
+template <class K>
+inline Object
+intersection(const Segment_3<K> &seg, const Tetrahedron_3<K> &tet)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(tet, seg);
+}
+
+template <class K>
+inline Object
+intersection(const Tetrahedron_3<K> &tet, const Segment_3<K> &seg)
+{
+ typedef typename K::Intersect_3 Intersect;
+ return Intersect()(tet, seg);
+}
+
+CGAL_END_NAMESPACE
+
+
+#endif /* ----- #ifndef segment_3_tetrahedron_3_intersection_INC ----- */
=== added file 'dolfin/mesh/Tetrahedron_3_Bbox_3_intersection.h'
--- dolfin/mesh/Tetrahedron_3_Bbox_3_intersection.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/Tetrahedron_3_Bbox_3_intersection.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,115 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-11
+// Last changed: 2009-11-11
+
+#ifndef TETRAHEDRON_3_BBOX_3_INTERSECTION_INC
+#define TETRAHEDRON_3_BBOX_3_INTERSECTION_INC
+
+#include <CGAL/Bbox_3.h>
+#include <CGAL/Tetrahedron_3.h>
+#include <CGAL/AABB_intersections/Bbox_3_triangle_3_do_intersect.h>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+ //This code is not optimized!!
+ template <class K>
+ inline
+ bool do_intersect(const typename K::Tetrahedron_3& tet,
+ const CGAL::Bbox_3& bbox,
+ const K & k)
+ {
+ typedef typename K::Point_3 Point;
+ typedef typename K::Triangle_3 Triangle;
+
+ //Check first whether on point of one primitive intersect the other
+ if (do_intersect(tet[0],bbox,k)) return true;
+ if (!k.has_on_unbounded_side_3_object()(tet,Point(bbox.xmin(),bbox.ymin(),bbox.zmin()))) return true;
+// if (!tet.has_on_unbounded_side(Point(bbox.xmin(),bbox.ymin(),bbox.zmin()))) return true;
+
+ //Otherwise one tetrahedron face must intersect the bbox in order to intersect.
+ if (do_intersect(bbox, Triangle(tet[0], tet[1], tet[2]))) return true;
+ if (do_intersect(bbox, Triangle(tet[0], tet[1], tet[3]))) return true;
+ if (do_intersect(bbox, Triangle(tet[0], tet[2], tet[3]))) return true;
+ if (do_intersect(bbox, Triangle(tet[1], tet[2], tet[3]))) return true;
+
+ return false;
+ }
+
+ template <class K>
+ inline
+ bool do_intersect(const CGAL::Bbox_3& bbox,
+ const typename K::Tetrahedron_3& tet,
+ const K & k)
+ {
+ return do_intersect(tet, bbox, k);
+ }
+
+//template <class K>
+//inline
+//Object
+//intersection(const typename K::Tetrahedron_3 &tet,
+// const typename K::Bbox_3 &bbox,
+// const K&)
+//{
+// if (do_intersect(tet,bbox)) {
+// return Object();
+// }
+// return Object();
+//}
+
+//template <class K>
+//inline
+//Object
+//intersection( const typename K::Bbox_3 &bbox,
+// const typename K::Tetrahedron_3 &tet,
+// const K&)
+//{
+// if (do_intersect(tet,bbox)) {
+// return Object();
+// }
+// return Object();
+//}
+
+
+} //namespace CGALi
+
+template <class K>
+inline
+bool do_intersect(const CGAL::Tetrahedron_3<K>& point,
+ const CGAL::Bbox_3& bbox)
+{
+ return typename K::Do_intersect_3()(point, bbox);
+}
+
+template <class K>
+inline
+bool do_intersect(const CGAL::Bbox_3& bbox,
+ const CGAL::Tetrahedron_3<K>& point)
+{
+ return typename K::Do_intersect_3()(point, bbox);
+}
+
+//template <class K>
+//inline Object
+//intersection(const Bbox_3 & bbox, const Tetrahedron_3<K> & tet)
+//{
+// typedef typename K::Intersect_3 Intersect;
+// return Intersect()(tet, bbox);
+//}
+
+//template <class K>
+//inline Object
+//intersection(const Tetrahedron_3<K> & tet, const Bbox_3 & bbox)
+//{
+// typedef typename K::Intersect_3 Intersect;
+// return Intersect()(tet, bbox);
+//}
+
+CGAL_END_NAMESPACE
+
+
+#endif /* ----- #ifndef TETRAHEDRON_3_BBOX_3_INTERSECTION_INC ----- */
=== added file 'dolfin/mesh/Tetrahedron_3_Tetrahedron_3_intersection.h'
--- dolfin/mesh/Tetrahedron_3_Tetrahedron_3_intersection.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/Tetrahedron_3_Tetrahedron_3_intersection.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,63 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-11
+// Last changed: 2009-11-10
+
+#ifndef TETRAHEDRON_3_TETRAHEDRON_3_INTERSECTION_INC
+#define TETRAHEDRON_3_TETRAHEDRON_3_INTERSECTION_INC
+
+#include <CGAL/Tetrahedron_3.h>
+//#include <CGAL/Triangle_3_Tetrahedron_3_do_intersect.h>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+ //This code is not optimized!!
+ template <class K >
+ bool
+ do_intersect(const typename K::Tetrahedron_3 & tet1,
+ const typename K::Tetrahedron_3 & tet2,
+ const K & k)
+ {
+ typedef typename K::Triangle_3 Triangle;
+
+ //Check first whether on point of one primitive intersect the other
+ //Check for all points, might be more efficient...?
+ if (!tet1.has_on_unbounded_side(tet2[0])) return true;
+ if (!tet1.has_on_unbounded_side(tet2[1])) return true;
+ if (!tet1.has_on_unbounded_side(tet2[2])) return true;
+ if (!tet1.has_on_unbounded_side(tet2[3])) return true;
+
+ if (!tet2.has_on_unbounded_side(tet1[0])) return true;
+ if (!tet2.has_on_unbounded_side(tet1[1])) return true;
+ if (!tet2.has_on_unbounded_side(tet1[2])) return true;
+ if (!tet2.has_on_unbounded_side(tet1[3])) return true;
+
+// if (!k.has_on_unbounded_side_3_object()(tet1,tet2[0])) return true;
+// if (!k.has_on_unbounded_side_3_object()(tet2,tet1[0])) return true;
+
+ //Otherwise one tetrahedron face must intersect the bbox in order to intersect.
+ if (do_intersect(tet1, Triangle(tet2[0], tet2[1], tet2[2]))) return true;
+ if (do_intersect(tet1, Triangle(tet2[0], tet2[1], tet2[3]))) return true;
+ if (do_intersect(tet1, Triangle(tet2[0], tet2[2], tet2[3]))) return true;
+ if (do_intersect(tet1, Triangle(tet2[1], tet2[2], tet2[3]))) return true;
+
+ return false;
+ }
+
+}// namespace CGALi
+
+template <class K>
+inline bool
+do_intersect(const Tetrahedron_3<K> &tet1, const Tetrahedron_3<K> &tet2)
+{
+ typedef typename K::Do_intersect_3 Do_intersect;
+ return Do_intersect()(tet1, tet2);
+}
+
+CGAL_END_NAMESPACE
+
+
+#endif /* ----- #ifndef TETRAHEDRON_3_TETRAHEDRON_3_INTERSECTION_INC ----- */
=== added file 'dolfin/mesh/Tree_Traits.h'
--- dolfin/mesh/Tree_Traits.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/Tree_Traits.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,59 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-10-12
+// Last changed: 2009-11-10
+
+#ifndef __TREE_TRAITS_H
+#define __TREE_TRAITS_H
+
+#include <CGAL/AABB_traits.h>
+
+namespace dolfin
+{
+
+ template<typename GeomTraits, typename AABB_primitive>
+ class Tree_Traits: public CGAL::AABB_traits<GeomTraits, AABB_primitive>
+ {
+
+ public:
+ typedef CGAL::AABB_traits<GeomTraits, AABB_primitive> AT;
+ typedef typename CGAL::Bbox_3 Bounding_box;
+ typedef AABB_primitive Primitive;
+
+
+ //Redefine this class in order to overwrite static compute_bbox methods
+ //and to use our own (calling directly the bbox of primitive, which is not
+ //required by th CGAL primitive concept.
+ class Compute_bbox {
+ public:
+ template<typename ConstPrimitiveIterator>
+ typename AT::Bounding_box operator()(ConstPrimitiveIterator first,
+ ConstPrimitiveIterator beyond) const
+ {
+ typename AT::Bounding_box bbox = compute_bbox(*first);
+ for(++first; first != beyond; ++first)
+ {
+ bbox = bbox + compute_bbox(*first);
+ }
+ return bbox;
+ }
+ };
+
+ Compute_bbox compute_bbox_object() {return Compute_bbox();}
+
+ private:
+ /**
+ * @brief Computes bounding box of one primitive
+ * @param pr the primitive
+ * @return the bounding box of the primitive \c pr
+ */
+ static Bounding_box compute_bbox(const Primitive& pr)
+ {
+ return pr.bbox();
+ }
+ };
+
+}
+
+#endif // ----- #ifndef __TREE_TRAITS_H -----
=== added file 'dolfin/mesh/Triangle_3_Tetrahedron_3_do_intersect_SCK.h'
--- dolfin/mesh/Triangle_3_Tetrahedron_3_do_intersect_SCK.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/Triangle_3_Tetrahedron_3_do_intersect_SCK.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,95 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-11
+// Last changed: 2009-11-29
+
+#ifndef CGAL_TRIANGLE_3_TETRAHEDRON_3_DO_INTERSECT_SCK_H
+#define CGAL_TRIANGLE_3_TETRAHEDRON_3_DO_INTERSECT_SCK_H
+
+#include <CGAL/Triangle_3_Triangle_3_do_intersect.h>
+#include <CGAL/enum.h>
+#include <CGAL/Simple_cartesian.h>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+template <>
+Simple_cartesian<double>::Boolean
+do_intersect<Simple_cartesian<double> >(const Simple_cartesian<double>::Triangle_3 &tr,
+ const Simple_cartesian<double>::Tetrahedron_3 &tet,
+ const Simple_cartesian<double> & k)
+{
+ typedef Simple_cartesian<double>::Triangle_3 Triangle;
+ typedef Simple_cartesian<double>::Point_3 Point;
+
+ CGAL_kernel_precondition( ! k.is_degenerate_3_object() (tr) );
+ CGAL_kernel_precondition( ! k.is_degenerate_3_object() (tet) );
+
+ if (!tet.has_on_unbounded_side(tr[0])) return true;
+ if (!tet.has_on_unbounded_side(tr[1])) return true;
+ if (!tet.has_on_unbounded_side(tr[2])) return true;
+
+ if (do_intersect(tr, Triangle(tet[0], tet[1], tet[2]), k)) return true;
+ if (do_intersect(tr, Triangle(tet[0], tet[1], tet[3]), k)) return true;
+ if (do_intersect(tr, Triangle(tet[0], tet[2], tet[3]), k)) return true;
+ if (do_intersect(tr, Triangle(tet[1], tet[2], tet[3]), k)) return true;
+
+ return false;
+}
+
+template <>
+Simple_cartesian<double>::Boolean
+do_intersect<Simple_cartesian<double> >(const Simple_cartesian<double>::Tetrahedron_3 &tet,
+ const Simple_cartesian<double>::Triangle_3 &tr,
+ const Simple_cartesian<double> & k)
+{
+ return do_intersect(tr, tet, k);
+}
+
+} //namespace CGALi
+
+
+template <>
+inline bool do_intersect<Simple_cartesian<double> >(const Tetrahedron_3<Simple_cartesian<double> > &tet,
+ const Triangle_3<Simple_cartesian<double> > &tr)
+{
+ return Simple_cartesian<double> ::Do_intersect_3()(tr,tet);
+}
+
+template <>
+inline bool do_intersect<Simple_cartesian<double> >(const Triangle_3<Simple_cartesian<double> > &tr,
+ const Tetrahedron_3<Simple_cartesian<double> > &tet)
+{
+ return Simple_cartesian<double> ::Do_intersect_3()(tr,tet);
+}
+
+//template <>
+//inline bool do_intersect<CGAL::Simple_cartesian<double> > (const Tetrahedron_3<CGAL::Simple_cartesian<double> > &tet,
+// const Triangle_3<CGAL::Simple_cartesian<double> > &tr)
+//{
+// typedef CGAL::Simple_cartesian<double>::Triangle_3 Triangle;
+
+// if (!tet.has_on_unbounded_side(tr[0])) return true;
+// if (!tet.has_on_unbounded_side(tr[1])) return true;
+// if (!tet.has_on_unbounded_side(tr[2])) return true;
+
+// if (do_intersect(tr, Triangle(tet[0], tet[1], tet[2]))) return true;
+// if (do_intersect(tr, Triangle(tet[0], tet[1], tet[3]))) return true;
+// if (do_intersect(tr, Triangle(tet[0], tet[2], tet[3]))) return true;
+// if (do_intersect(tr, Triangle(tet[1], tet[2], tet[3]))) return true;
+//
+// return false;
+//}
+
+//template <>
+//inline bool do_intersect<CGAL::Simple_cartesian<double> > (const Triangle_3<CGAL::Simple_cartesian<double> > &tr,
+// const Tetrahedron_3<CGAL::Simple_cartesian<double> > &tet)
+//{
+// return do_intersect<CGAL::Simple_cartesian<double> > (tet,tr);
+//}
+
+CGAL_END_NAMESPACE
+
+#endif // CGAL_TRIANGLE_3_TETRAHEDRON_3_DO_INTERSECT_SCK_H
=== modified file 'dolfin/mesh/Vertex.h'
--- dolfin/mesh/Vertex.h 2009-04-27 13:01:54 +0000
+++ dolfin/mesh/Vertex.h 2009-11-30 01:20:26 +0000
@@ -30,16 +30,16 @@
~Vertex() {}
/// Return value of vertex coordinate i
- inline double x(uint i) const { return _mesh.geometry().x(_index, i); }
+ inline double x(uint i) const { return _mesh->geometry().x(_index, i); }
/// Return vertex coordinates as a 3D point value
- inline Point point() const { return _mesh.geometry().point(_index); }
+ inline Point point() const { return _mesh->geometry().point(_index); }
/// Return array of vertex coordinates
- //double* x() { return _mesh.geometry().x(_index); }
+ //double* x() { return _mesh->geometry().x(_index); }
/// Return array of vertex coordinates (const version)
- const double* x() const { return _mesh.geometry().x(_index); }
+ const double* x() const { return _mesh->geometry().x(_index); }
};
=== added file 'dolfin/mesh/added_intersection_3.h'
--- dolfin/mesh/added_intersection_3.h 1970-01-01 00:00:00 +0000
+++ dolfin/mesh/added_intersection_3.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,25 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-09-11
+// Last changed: 2009-11-10
+
+#ifndef ADDED_INTERSECTION_3_INC
+#define ADDED_INTERSECTION_3_INC
+
+#include "Point_3_Point_3_intersection.h"
+#include "Point_3_Segment_3_intersection.h"
+//#include "Point_3_Triangle_3_intersection.h" //already defined in original
+//CGAL, but not with p.has_on(pt) function??
+#include "Point_3_Tetrahedron_3_intersection.h"
+#include "Point_3_Iso_Cuboid_3_intersection.h"
+#include "Point_3_Line_3_intersection.h"
+#include "Point_3_Ray_3_intersection.h"
+#include "Point_3_Bbox_3_intersection.h"
+#include "Segment_3_Segment_3_intersection.h"
+#include "Segment_3_Tetrahedron_3_intersection.h"
+#include "Tetrahedron_3_Tetrahedron_3_intersection.h"
+//#include "Triangle_3_Tetrahedron_3_do_intersect_SCK.h"
+#include "Tetrahedron_3_Bbox_3_intersection.h"
+
+#endif /* ----- #ifndef ADDED_INTERSECTION_3_INC ----- */
=== modified file 'dolfin/mesh/dolfin_mesh.h'
--- dolfin/mesh/dolfin_mesh.h 2009-08-10 06:41:26 +0000
+++ dolfin/mesh/dolfin_mesh.h 2009-11-30 01:20:26 +0000
@@ -22,6 +22,8 @@
#include <dolfin/mesh/MeshFunction.h>
#include <dolfin/mesh/Mesh.h>
#include <dolfin/mesh/MeshPartitioning.h>
+#include <dolfin/mesh/MeshPrimitive.h>
+#include <dolfin/mesh/Primitive_Traits.h>
#include <dolfin/mesh/LocalMeshData.h>
#include <dolfin/mesh/SubDomain.h>
#include <dolfin/mesh/SubMesh.h>
@@ -36,5 +38,7 @@
#include <dolfin/mesh/Rectangle.h>
#include <dolfin/mesh/UnitSphere.h>
#include <dolfin/mesh/IntersectionDetector.h>
+#include <dolfin/mesh/IntersectionOperator.h>
+#include <dolfin/mesh/added_intersection_3.h>
#endif
=== modified file 'dolfin/swig/common_post.i'
--- dolfin/swig/common_post.i 2009-09-09 13:22:16 +0000
+++ dolfin/swig/common_post.i 2009-11-30 01:20:26 +0000
@@ -7,6 +7,8 @@
//%template(STLVectorString) std::vector<std::string>;
//%template(STLPairUInt) std::pair<dolfin::uint,dolfin::uint>;
+/*%template(BOOSTUnorderSetUInt) boost::unordered_set<dolfin::uint>;*/
+
%extend dolfin::Variable
{
std::string __str__() const
=== modified file 'dolfin/swig/defines.i'
--- dolfin/swig/defines.i 2009-09-04 16:38:32 +0000
+++ dolfin/swig/defines.i 2009-11-30 01:20:26 +0000
@@ -46,6 +46,15 @@
#endif
}
+bool has_cgal()
+{
+#ifdef HAS_CGAL
+ return true;
+#else
+ return false;
+#endif
+}
+
bool has_umfpack()
{
#ifdef HAS_UMFPACK
=== modified file 'dolfin/swig/dolfin.i'
--- dolfin/swig/dolfin.i 2009-10-07 18:45:48 +0000
+++ dolfin/swig/dolfin.i 2009-11-30 01:20:26 +0000
@@ -9,7 +9,7 @@
// Modified by Garth N. Wells, 2009.
//
// First added: 2005-10-24
-// Last changed: 2009-10-07
+// Last changed: 2009-11-27
// The PyDOLFIN extension module
%module(package="dolfin", directors="1") cpp
@@ -31,6 +31,7 @@
%include "dolfin/swig/typemaps.i"
%include "dolfin/swig/numpy_typemaps.i"
%include "dolfin/swig/std_vector_typemaps.i"
+%include "dolfin/swig/std_set_typemaps.i"
// Global exceptions
%include <exception.i>
=== modified file 'dolfin/swig/generate.py'
--- dolfin/swig/generate.py 2009-10-07 18:45:48 +0000
+++ dolfin/swig/generate.py 2009-11-30 01:20:26 +0000
@@ -16,7 +16,8 @@
# List of headers to exclude (add more here)
excludes = ["plot.h", "ParameterSystem.h", "ParameterList.h", \
- "ConvectionMatrix.h", "MassMatrix.h", "StiffnessMatrix.h", "LoadVector.h"]
+ "ConvectionMatrix.h", "MassMatrix.h", "StiffnessMatrix.h", "LoadVector.h", \
+ "IntersectionOperatorImplementation.h" ]
# Name of SWIG interface file to be generated
interface_file = "kernel_modules.i"
=== modified file 'dolfin/swig/kernel_modules.i'
--- dolfin/swig/kernel_modules.i 2009-11-11 19:59:27 +0000
+++ dolfin/swig/kernel_modules.i 2009-11-30 01:20:26 +0000
@@ -117,6 +117,7 @@
%include "dolfin/mesh/Rectangle.h"
%include "dolfin/mesh/UnitSphere.h"
%include "dolfin/mesh/IntersectionDetector.h"
+%include "dolfin/mesh/IntersectionOperator.h"
%include "dolfin/swig/mesh_post.i"
// DOLFIN headers included from function
=== modified file 'dolfin/swig/mesh_pre.i'
--- dolfin/swig/mesh_pre.i 2009-11-27 06:43:01 +0000
+++ dolfin/swig/mesh_pre.i 2009-11-30 01:20:26 +0000
@@ -8,7 +8,7 @@
// Modified by Johan Hake 2008-2009
//
// First added: 2006-09-20
-// Last changed: 2009-11-25
+// Last changed: 2009-11-30
//=============================================================================
// SWIG directives for the DOLFIN Mesh kernel module (pre)
@@ -107,6 +107,8 @@
%ignore dolfin::MeshGeometry::operator=;
%ignore dolfin::MeshTopology::operator=;
%ignore dolfin::MeshConnectivity::operator=;
+%ignore dolfin::MeshEntityIterator::operator->;
+
//-----------------------------------------------------------------------------
// Map increment, decrease and dereference operators for iterators
=== added file 'dolfin/swig/std_set_typemaps.i'
--- dolfin/swig/std_set_typemaps.i 1970-01-01 00:00:00 +0000
+++ dolfin/swig/std_set_typemaps.i 2009-11-30 01:20:26 +0000
@@ -0,0 +1,63 @@
+namespace std
+{
+ template <class T> class set
+ {
+ };
+}
+
+
+%define ARGOUT_TYPEMAP_BOOST_UNORDERED_SET_OF_PRIMITIVES(TYPE, TYPE_UPPER, ARG_NAME, NUMPY_TYPE)
+//-----------------------------------------------------------------------------
+// In typemap removing the argument from the expected in list
+//-----------------------------------------------------------------------------
+%typemap (in,numinputs=0) std::set<TYPE>& ARG_NAME (std::set<TYPE> set_temp)
+{
+ $1 = &set_temp;
+}
+
+//-----------------------------------------------------------------------------
+// Argout typemap, returning a NumPy array for the boost::unordered_set<TYPE>
+//-----------------------------------------------------------------------------
+%typemap(argout) std::set<TYPE> & ARG_NAME
+{
+ PyObject* o0 = 0;
+ PyObject* o1 = 0;
+ PyObject* o2 = 0;
+ npy_intp size = $1->size();
+ PyArrayObject *ret = reinterpret_cast<PyArrayObject*>(PyArray_SimpleNew(1, &size, NUMPY_TYPE));
+ TYPE* data = static_cast<TYPE*>(PyArray_DATA(ret));
+
+ int i = 0;
+ for (std::set<TYPE>::const_iterator it = (*$1).begin(); it != (*$1).end(); ++it)
+ {
+ data[i] = *it;
+ ++i;
+ }
+ o0 = PyArray_Return(ret);
+ // If the $result is not already set
+ if ((!$result) || ($result == Py_None))
+ {
+ $result = o0;
+ }
+ // If the result is set by another out typemap build a tuple of arguments
+ else
+ {
+ // If the the argument is set but is not a tuple make one and put the result in it
+ if (!PyTuple_Check($result))
+ {
+ o1 = $result;
+ $result = PyTuple_New(1);
+ PyTuple_SetItem($result, 0, o1);
+ }
+ o2 = PyTuple_New(1);
+ PyTuple_SetItem(o2, 0, o0);
+ o1 = $result;
+ $result = PySequence_Concat(o1, o2);
+ Py_DECREF(o1);
+ Py_DECREF(o2);
+ }
+}
+
+%enddef
+
+ARGOUT_TYPEMAP_BOOST_UNORDERED_SET_OF_PRIMITIVES(dolfin::uint, INT32, ids_result, NPY_INT)
=== added directory 'sandbox/cgal'
=== added file 'sandbox/cgal/README'
--- sandbox/cgal/README 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/README 2009-11-30 01:20:26 +0000
@@ -0,0 +1,5 @@
+Most of these programms are only tiny code-snippet to try out some C++ features
+and restrictions (I <irony>love</irony> templates) and in particular to get so
+experience with the CGAL library. It was also used in the early phase of
+testing new intersection class, so there might some old snippets lying around,
+which does not compile. Have a look at the demos to see how the interface is working now.
=== added directory 'sandbox/cgal/benchmark'
=== added file 'sandbox/cgal/benchmark/SConstruct'
--- sandbox/cgal/benchmark/SConstruct 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/benchmark/SConstruct 2009-11-30 01:20:26 +0000
@@ -0,0 +1,28 @@
+import os, commands
+
+# Get compiler from pkg-config
+compiler = commands.getoutput('pkg-config --variable=compiler dolfin')
+
+# Create a SCons Environment based on the main os environment
+env = Environment(ENV=os.environ, CXX=compiler)
+
+# Get compiler flags from pkg-config
+env.ParseConfig('pkg-config --cflags --libs dolfin')
+
+#Use user environment as base.
+# env = Environment(ENV = os.environ);
+#local_path = os.environ['LOCALPATH']
+#
+#env.Append(CCFLAGS = '-frounding-math') #necessary for CGAL
+#lib_path = local_path + os.sep + 'lib'
+#include_path = local_path + os.sep + 'include'
+#libs = Split('CGAL CGAL_Core')
+#
+#env.Append(LIBS = libs)
+#env.Append(LIBPATH = lib_path)
+#env.Append(CPPPATH = include_path)
+
+#Program
+env.Program(['check_eval_results.cpp'])
+env.Program(['point_detection.cpp'])
+env.Program(['mesh_intersection_detection.cpp'])
=== added file 'sandbox/cgal/benchmark/check_eval_results.cpp'
--- sandbox/cgal/benchmark/check_eval_results.cpp 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/benchmark/check_eval_results.cpp 2009-11-30 01:20:26 +0000
@@ -0,0 +1,172 @@
+
+#include <cstdlib>
+#include <iostream>
+#include <climits>
+
+//#include <boost/foreach.hpp>
+#include <boost/progress.hpp>
+
+#include <dolfin.h>
+#include "Projection.h"
+
+using namespace dolfin;
+
+using boost::timer;
+using boost::progress_timer;
+using boost::progress_display;
+//using std::cout;
+//using std::endl;
+
+
+
+#ifdef HAS_GTS
+
+class F_1 : public Function
+{
+public:
+ void eval(double* values, const double* x) const
+ {
+ values[0] = x[0];
+ }
+};
+
+class F_2 : public Function
+{
+public:
+ void eval(double* values, const double* x) const
+ {
+ values[0] = x[1];
+ }
+};
+
+class F_3 : public Function
+{
+public:
+ void eval(double* values, const double* x) const
+ {
+ values[0] = x[2];
+ }
+};
+
+class F_4 : public Function
+{
+public:
+ void eval(double* values, const double* x) const
+ {
+ values[0] = sin(3.0*x[0])*sin(3.0*x[1]);
+ }
+};
+
+class F_5 : public Function
+{
+public:
+ void eval(double* values, const double* x) const
+ {
+ values[0] = sin(3.0*x[0])*sin(3.0*x[1])*sin(3.0*x[2]);
+ }
+};
+
+int main()
+{
+
+// bool write = false;
+
+// std::cout << "Should function eval results be written to stdout ? " << std::endl;
+// std::cin >> write;
+
+ //intialize the random number generator
+ srand(1);
+
+ unsigned int N_max = 100;
+ for (unsigned int N =10; N <= N_max; )
+ {
+ UnitCube mesh(N,N,N);
+
+ Projection::FunctionSpace V(mesh);
+ Projection::BilinearForm a(V, V);
+ Projection::LinearForm L(V);
+
+ // F_1 f_1;
+ // F_2 f_2;
+ // F_3 f_3;
+ // F_4 f_4;
+
+ // L.f = f_1;
+ // VariationalProblem pde_1(a, L);
+ // Function Pf_1;
+ // pde_1.solve(Pf_1);
+
+ // L.f = f_2;
+ // VariationalProblem pde_2(a, L);
+ // Function Pf_2;
+ // pde_2.solve(Pf_2);
+
+ // L.f = f_3;
+ // VariationalProblem pde_3(a, L);
+ // Function Pf_3;
+ // pde_3.solve(Pf_3);
+
+ // L.f = f_4;
+ // VariationalProblem pde_4(a, L);
+ // Function Pf_4;
+ // pde_4.solve(Pf_4);
+
+ F_5 f_5;
+ L.f = f_5;
+ VariationalProblem pde_5(a, L);
+ Function Pf_5;
+ pde_5.solve(Pf_5);
+
+ int max = 100000;
+
+ double X[3] = {0.0, 0.0, 0.0};
+ double value = 0.0;
+
+ std::cout <<"#Starting function evaluation of random points" <<std::endl;
+ timer t0;
+
+ for (int i = 1; i <= max; ++i)
+ {
+
+ X[0] = std::rand()/static_cast<double>(RAND_MAX);
+ X[1] = std::rand()/static_cast<double>(RAND_MAX);
+ X[2] = std::rand()/static_cast<double>(RAND_MAX);
+
+ // cout << "X[0] : " << X[0];
+ // cout << "X[1] : " << X[1];
+ // cout << "X[2] : " << X[2];
+ // f_1.eval(&value, X);
+ // f_2.eval(&value, X);
+ // f_3.eval(&value, X);
+ // f_4.eval(&value, X);
+ // f_5.eval(&value, X);
+ // info("f_5(x) = %g", value);
+
+ Pf_5.eval(&value, X);
+ // info("Pf_5(x) = %g", value);
+
+ //increase number of points by 10
+ // i = 10;
+
+ }
+
+ std::cout <<"Time elapsed: " << t0.elapsed() <<std::endl;
+
+ //Increase mesh size
+ N += 10;
+ }
+
+ return 0;
+}
+
+
+#else
+
+int main()
+{
+ info("DOLFIN must be compiled with GTS to run this demo.");
+ return 0;
+}
+
+#endif
+
=== added file 'sandbox/cgal/benchmark/mesh_intersection_detection.cpp'
--- sandbox/cgal/benchmark/mesh_intersection_detection.cpp 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/benchmark/mesh_intersection_detection.cpp 2009-11-30 01:20:26 +0000
@@ -0,0 +1,81 @@
+#include <cstdlib>
+#include <iostream>
+#include <climits>
+
+//#include <boost/foreach.hpp>
+#include <boost/progress.hpp>
+
+#include <dolfin.h>
+#include "Projection.h"
+
+using namespace dolfin;
+
+using boost::timer;
+using boost::progress_timer;
+using boost::progress_display;
+
+
+int main()
+{
+ unsigned int N_max = 1000;
+ for (unsigned int N =1000; N <= N_max; )
+ {
+// UnitCube mesh0(N,N,N);
+// UnitCube mesh1(N,N,N);
+ UnitSquare mesh0(N,N);
+// UnitSquare mesh1(100,100);
+ UnitCircle mesh1(50);
+
+ MeshGeometry & geo1 = mesh1.geometry();
+ // Move and scale circle
+ for (VertexIterator vertex(mesh1); !vertex.end(); ++vertex)
+ {
+ double* x = geo1.x(vertex->index());
+ x[0] = 0.2*x[0] + 0.5;
+ x[1] = 0.2*x[1] + 0.5;
+ }
+
+ std::cout <<"#Creating CGAL intersection operator ...";
+ timer t0;
+ IntersectionOperator new_ic(mesh0);
+ std::cout <<"Time elapsed: " << t0.elapsed() <<std::endl;
+
+ unordered_set s_cells;
+
+ std::cout <<"#Starting intersection with mesh1 using CGAL Interface\t" ;
+ t0.restart();
+ new_ic.all_intersected_entities(mesh1,s_cells);
+
+ std::cout <<"Time elapsed: " << t0.elapsed() <<std::endl;
+ std::cout <<"Lenght of cell container: " <<s_cells.size() <<std::endl;
+
+// for (unordered_set::const_iterator ic = s_cells.begin(); ic != s_cells.end(); ++ic)
+// {
+// std::cout <<"Id: " <<*ic <<std::endl;
+// }
+
+
+ std::cout <<"#Creating GTS intersection operator ..." ;
+ t0.restart();
+ IntersectionDetector old_ic(mesh0);
+ std::cout <<"Time elapsed: " << t0.elapsed() <<std::endl;
+
+// std::vector<unsigned int> v_cells;
+// std::cout <<"#Starting intersection with mesh1 using GTS Interface\t" ;
+
+// t0.restart();
+// old_ic.intersection(mesh1,v_cells);
+// std::cout <<"Time elapsed: " << t0.elapsed() <<std::endl;
+
+// std::cout <<"Lenght of cell container: " <<v_cells.size() <<std::endl;
+
+// for (std::vector<unsigned int>::const_iterator ic = v_cells.begin(); ic != v_cells.end(); ++ic)
+// {
+// std::cout <<"Id: " <<*ic <<std::endl;
+// }
+
+ N += 10;
+ }
+
+ return 0;
+}
=== added file 'sandbox/cgal/benchmark/point_detection.cpp'
--- sandbox/cgal/benchmark/point_detection.cpp 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/benchmark/point_detection.cpp 2009-11-30 01:20:26 +0000
@@ -0,0 +1,77 @@
+#include <cstdlib>
+#include <iostream>
+#include <climits>
+
+//#include <boost/foreach.hpp>
+#include <boost/progress.hpp>
+
+#include <dolfin.h>
+#include "Projection.h"
+
+using namespace dolfin;
+
+using boost::timer;
+using boost::progress_timer;
+using boost::progress_display;
+
+
+int main()
+{
+ unsigned int N_max = 150;
+ for (unsigned int N =10; N <= N_max; )
+ {
+// UnitCube mesh(N,N,N);
+ UnitSquare mesh(N,N);
+ IntersectionDetector old_ic(mesh);
+ IntersectionOperator new_ic(mesh);
+
+ int max = 3000000;
+
+ double X[3] = {0.0, 0.0, 0.0};
+
+ srand(1);
+
+ timer t0;
+
+ std::cout <<"#Starting intersection with random points using CGAL Interface\t" ;
+ unordered_set s_cells;
+ for (int i = 1; i <= max; ++i)
+ {
+
+ X[0] = std::rand()/static_cast<double>(RAND_MAX);
+ X[1] = std::rand()/static_cast<double>(RAND_MAX);
+// X[2] = std::rand()/static_cast<double>(RAND_MAX);
+
+ Point p(X[0],X[1],X[2]);
+
+ new_ic.all_intersected_entities(p,s_cells);
+ }
+
+ std::cout <<"Time elapsed: " << t0.elapsed() <<std::endl;
+// std::cout <<"Lenght of cell container: " <<s_cells.size() <<std::endl;
+
+ srand(1);
+ t0.restart();
+
+ std::cout <<"#Starting intersection with random points using GTS Interface\t";
+ std::vector<unsigned int> v_cells;
+ for (int i = 1; i <= max; ++i)
+ {
+
+ X[0] = std::rand()/static_cast<double>(RAND_MAX);
+ X[1] = std::rand()/static_cast<double>(RAND_MAX);
+// X[2] = std::rand()/static_cast<double>(RAND_MAX);
+
+ Point p(X[0],X[1],X[2]);
+
+ old_ic.intersection(p,v_cells);
+ }
+
+ std::cout <<"Time elapsed: " << t0.elapsed() <<std::endl;
+// std::cout <<"Lenght of cell container: " <<v_cells.size() <<std::endl;
+
+ N += 10;
+ }
+
+ return 0;
+}
=== added directory 'sandbox/cgal/mesh_intersection'
=== added directory 'sandbox/cgal/mesh_intersection/cpp'
=== added file 'sandbox/cgal/mesh_intersection/cpp/SConstruct'
--- sandbox/cgal/mesh_intersection/cpp/SConstruct 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/mesh_intersection/cpp/SConstruct 2009-11-30 01:20:26 +0000
@@ -0,0 +1,18 @@
+# This is an example of a simple SConstruct file for building
+# programs against DOLFIN. To build this demo, just type 'scons'.
+
+import os, commands
+
+# Get compiler from pkg-config
+compiler = commands.getoutput('pkg-config --variable=compiler dolfin')
+
+# Create a SCons Environment based on the main os environment
+env = Environment(ENV=os.environ, CXX=compiler)
+
+# Get compiler flags from pkg-config
+env.ParseConfig('pkg-config --cflags --libs dolfin')
+
+#env.Program(['tetrahedron_triangle_intersection_test.cpp'])
+env.Program(['intersection_test.cpp'])
+env.Program(['intersection_test_2.cpp'])
+env.Program(['intersection_operator_test.cpp'])
=== added file 'sandbox/cgal/mesh_intersection/cpp/intersection_operator_test.cpp'
--- sandbox/cgal/mesh_intersection/cpp/intersection_operator_test.cpp 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/mesh_intersection/cpp/intersection_operator_test.cpp 2009-11-30 01:20:26 +0000
@@ -0,0 +1,122 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU GPL Version 2.1.
+//
+// First added: 2008-10-08
+// Last changed: 2009-11-22
+
+#include <set>
+
+#include <dolfin/mesh/dolfin_mesh.h> //Should be included before the CGAL Kernel.
+#include <math.h>
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+typedef CGAL::Simple_cartesian<double> K;
+//typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef K::Point_3 Point_3;
+typedef K::Segment_3 Segment_3;
+typedef K::Triangle_3 Triangle_3;
+typedef K::Tetrahedron_3 Tetrahedron_3;
+
+
+using namespace dolfin;
+
+//CGAL points
+typedef std::vector<Point_3>::iterator P3_iterator;
+//dolfin points
+typedef std::vector<Point>::iterator P_iterator;
+
+typedef std::vector<Segment_3>::iterator Seg_iterator;
+typedef std::vector<Triangle_3>::iterator Tri_iterator;
+typedef std::vector<Tetrahedron_3>::iterator Tet_iterator;
+
+int main()
+{
+ Point e0(0.0,0.0,0.0);
+ Point e1(1.0,0.0,0.0);
+ Point e2(0.0,1.0,0.0);
+ Point e3(0.0,0.0,1.0);
+
+ Point a(0.1, 0.0, 0.0);
+ Point b(0.0, 0.1, 0.0);
+ Point c(0.0, 0.0, 0.1);
+ Point d(0.0, 1.0/3.0, 0.0);
+
+ double height = 1.2;
+ Point e(1.0/3.0, height, 0.0);
+ Point f(2.0/3.0, height, 0.0);
+ Point g(0.5, 0.0, 0.0);
+
+ //points shifted by (2,0,0)
+ Point A(3.0, 0.0, 0.0);
+ Point B(2.0, 1.0, 0.0);
+ Point C(2.0, 0.0, 1.0);
+ Point D(2.0, 0.0, 0.0);
+
+ //points shifted by (-1,0,0)
+ Point x(-1.0, 0.0, 0.0);
+ Point y(-1.0, 1.0, 0.0);
+ Point z(-1.0, 0.0, 1.0);
+
+ std::vector<Point> points;
+ points.push_back(e1);
+ points.push_back(e2);
+ points.push_back(e3);
+ points.push_back(a);
+ points.push_back(b);
+ points.push_back(c);
+ points.push_back(d);
+ points.push_back(e0);
+ points.push_back(A);
+ points.push_back(B);
+ points.push_back(C);
+ points.push_back(D);
+
+ UnitCube cube(3,3,2);
+ cout <<"Total number of cells in Cube:" << cube.num_cells() <<endl;
+
+ UnitSphere sphere(3);
+ cout <<"Total number of cells in Sphere:" << sphere.num_cells() <<endl;
+
+ IntersectionOperator io(cube);
+// IntersectionOperator io(cube,"SimpleCartesian"); //the same as before
+// IntersectionOperator io(cube,"ExactPredicates");
+
+ int counter = 0;
+
+ //Point intersection, point by point
+ for (P_iterator i = points.begin(); i != points.end(); ++i)
+ {
+ ++counter;
+ std::set<unsigned int> cells;
+ io.all_intersected_entities(*i,cells);
+ std::cout <<"Intersection with point :" << counter <<std::endl;
+ for (std::set<unsigned int>::const_iterator ic = cells.begin(); ic != cells.end(); ++ic)
+ {
+ std::cout <<"Id: " <<*ic <<std::endl;
+ }
+ }
+
+ //Entire point set at all
+ std::set<unsigned int> cells;
+ io.all_intersected_entities(points,cells);
+ cout <<"Intersection with pointlist:" <<endl;
+ for (std::set<unsigned int>::const_iterator ic = cells.begin(); ic != cells.end(); ++ic)
+ {
+ std::cout <<"Id: " <<*ic <<std::endl;
+ }
+
+ //Mesh intersection
+ cells.clear();
+ io.all_intersected_entities(sphere,cells);
+ cout <<"Intersection with sphere:" <<endl;
+ for (std::set<unsigned int>::const_iterator ic = cells.begin(); ic != cells.end(); ++ic)
+ {
+ std::cout <<"Id: " <<*ic <<std::endl;
+ }
+
+ return 0;
+}
+
+
=== added file 'sandbox/cgal/mesh_intersection/cpp/intersection_test.cpp'
--- sandbox/cgal/mesh_intersection/cpp/intersection_test.cpp 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/mesh_intersection/cpp/intersection_test.cpp 2009-11-30 01:20:26 +0000
@@ -0,0 +1,184 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU GPL Version 2.1.
+//
+// First added: 2008-10-08
+// Last changed: 2009-11-22
+//
+// The IntersectionOperator class relies on a CGAL-based search tree structure
+// This file illustrates how to use this search structure in combination with
+// CGAL primitives.
+
+
+#include <dolfin/mesh/dolfin_mesh.h> //Should be included before the CGAL Kernel.
+
+#include <CGAL/AABB_tree.h> // *Must* be inserted before kernel!
+#include <CGAL/AABB_traits.h>
+
+#include <CGAL/Simple_cartesian.h>
+//#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Bbox_3.h>
+
+#include <math.h>
+
+using namespace dolfin;
+
+typedef CGAL::Simple_cartesian<double> K;
+//Alternative kernel
+//typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+typedef K::Point_3 Point_3;
+typedef K::Segment_3 Segment_3;
+typedef K::Triangle_3 Triangle_3;
+typedef K::Tetrahedron_3 Tetrahedron_3;
+
+typedef TetrahedronCellPrimitive Primitive;
+typedef Primitive_Traits<Primitive,K> PT;
+typedef MeshPrimitive<PT> CellPrimitive;
+
+typedef CGAL::AABB_traits<K, CellPrimitive> AABB_triangle_traits;
+typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+
+typedef Tree::Object_and_primitive_id Object_and_primitive_id;
+typedef Tree::Primitive_id Primitive_id;
+
+typedef std::vector<Point_3>::iterator P_iterator;
+typedef std::vector<Segment_3>::iterator Seg_iterator;
+typedef std::vector<Triangle_3>::iterator Tri_iterator;
+typedef std::vector<Tetrahedron_3>::iterator Tet_iterator;
+
+//testing intersection with dolfin points.
+
+typedef std::vector<Point> Points;
+typedef Points::iterator Poi_Iterator;
+
+void print_point(const Point & p)
+{
+ std::cout.precision(16);
+ std::cout << "<Point x = " << p.x() << " y = " << p.y() << " z = " << p.z() << ">" <<std::endl;
+}
+
+
+int main()
+{
+ //Building points
+ //Unit points
+ Point_3 e0(0.0,0.0,0.0);
+ Point_3 e1(1.0,0.0,0.0);
+ Point_3 e2(0.0,1.0,0.0);
+ Point_3 e3(0.0,0.0,1.0);
+
+ Point_3 a(0.1, 0.0, 0.0);
+ Point_3 b(0.0, 0.1, 0.0);
+ Point_3 c(0.0, 0.0, 0.1);
+ Point_3 d(0.0, 1.0/3.0, 0.0);
+
+ double height = 1.2;
+ Point_3 e(1.0/3.0, height, 0.0);
+ Point_3 f(2.0/3.0, height, 0.0);
+ Point_3 g(0.5, 0.0, 0.0);
+
+ //points shifted by (2,0,0)
+ Point_3 A(3.0, 0.0, 0.0);
+ Point_3 B(2.0, 1.0, 0.0);
+ Point_3 C(2.0, 0.0, 1.0);
+ Point_3 D(2.0, 0.0, 0.0);
+
+ //points shifted by (-1,0,0)
+ Point_3 x(-1.0, 0.0, 0.0);
+ Point_3 y(-1.0, 1.0, 0.0);
+ Point_3 z(-1.0, 0.0, 1.0);
+
+ std::vector<Point_3> points;
+ points.push_back(a);
+ points.push_back(b);
+ points.push_back(c);
+ points.push_back(d);
+ points.push_back(A);
+ points.push_back(B);
+ points.push_back(C);
+ points.push_back(D);
+
+ std::vector<Triangle_3> triangles;
+ triangles.push_back(Triangle_3(e0,a,b));
+ triangles.push_back(Triangle_3(e0,e1,b));
+ triangles.push_back(Triangle_3(e0,e1,c));
+ triangles.push_back(Triangle_3(x,y,z));
+ triangles.push_back(Triangle_3(e0,x,y));
+ triangles.push_back(Triangle_3(d,x,y));
+ triangles.push_back(Triangle_3(e,f,g));
+ triangles.push_back(Triangle_3(A,B,C));
+ triangles.push_back(Triangle_3(x,y,d));
+ triangles.push_back(Triangle_3(x,y,D));
+
+ std::vector<Tetrahedron_3> tetrahedrons;
+ tetrahedrons.push_back(Tetrahedron_3(a,b,c,d));
+ tetrahedrons.push_back(Tetrahedron_3(A,B,C,D));
+
+ UnitCube mesh(3,3,3);
+
+ for (MeshEntityIterator check_iter(mesh,3); !check_iter.end(); ++check_iter)
+ {
+ cout <<"Cell Index: " <<check_iter->index() << endl;
+ for (VertexIterator v(*check_iter); !v.end(); ++v)
+ {
+ cout <<"vertex index = " << v->index() <<"\t"
+ << "coordinates:\t" << v->point() << endl;
+ std::cout <<"\t\t\t\t\t";
+ print_point(v->point());
+ }
+ }
+
+ MeshEntityIterator cell_iter(mesh,3);
+
+ //Build the search tree.
+ Tree tree(cell_iter,cell_iter.end_iterator());
+
+ //Point intersection
+// int pj = 0;
+// for (P_iterator i = points.begin(); i != points.end(); ++i)
+// {
+// ++pj;
+// if(tree.do_intersect(*i))
+// std::cout << std::endl << "intersection(s) with point " << pj;
+
+// std::cout << std::endl << tree.number_of_intersected_primitives(*i)
+// << " intersection(s) with point number " << pj << std::endl;
+
+ //Compute the cell ids of the corresponding intersection operations.
+// std::vector<Primitive_id> primitives;
+// tree.all_intersected_primitives(*i,std::back_inserter(primitives));
+// std::cout <<"Cell ids for point " << pj << ":" <<std::endl;
+// for (std::vector<Primitive_id>::iterator id_iter = primitives.begin();
+// id_iter != primitives.end(); ++id_iter)
+// {
+// std::cout <<"Id : " <<*id_iter <<std::endl;
+// }
+// }
+
+ //counts intersections for each triangle.
+// cout << endl <<"Compute all intersection ids for each triangle: ";
+// int trij = 0;
+// for (Tri_iterator i = triangles.begin(); i != triangles.end(); ++i)
+// {
+// ++trij;
+// if(tree.do_intersect(*i))
+// std::cout << std::endl <<"intersection(s) with triangle " << trij;
+// std::cout << std::endl << tree.number_of_intersected_primitives(*i)
+// << " intersection(s) with triangle number " << trij << std::endl;
+
+ //Compute the cell ids of the corresponding intersection operation.
+// std::vector<Primitive_id> primitives;
+// tree.all_intersected_primitives(*i,std::back_inserter(primitives));
+// std::cout <<"Cell ids for triangle " << trij << ":" <<std::endl;
+// for (std::vector<Primitive_id>::iterator id_iter = primitives.begin();
+// id_iter != primitives.end(); ++id_iter)
+// {
+// std::cout <<"Id : " <<*id_iter <<std::endl;
+// }
+// }
+
+ return 0;
+}
+
+
=== added file 'sandbox/cgal/mesh_intersection/cpp/intersection_test_2.cpp'
--- sandbox/cgal/mesh_intersection/cpp/intersection_test_2.cpp 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/mesh_intersection/cpp/intersection_test_2.cpp 2009-11-30 01:20:26 +0000
@@ -0,0 +1,96 @@
+#include <dolfin.h>
+#include <math.h>
+
+#include <dolfin/mesh/Point.h>
+
+#include <CGAL/Bbox_3.h>
+
+#include <CGAL/Simple_cartesian.h> //used kernel
+typedef CGAL::Simple_cartesian<double> K;
+
+typedef K::Point_3 Point_3;
+typedef K::Segment_3 Segment_3;
+typedef K::Triangle_3 Triangle_3;
+typedef K::Tetrahedron_3 Tetrahedron_3;
+typedef K::Iso_cuboid_3 Iso_cuboid_3;
+typedef K::Ray_3 Ray_3;
+typedef K::Line_3 Line_3;
+
+typedef CGAL::Bbox_3 Bbox_3;
+
+using namespace dolfin;
+
+void print_point(const Point & p)
+{
+ std::cout.precision(16);
+ std::cout << "<Point x = " << p.x() << " y = " << p.y() << " z = " << p.z() << ">" <<std::endl;
+}
+
+
+int main()
+{
+
+ UnitCube cube(3,3,2);
+ UnitSphere sphere(3);
+
+ Cell cell_1(cube,102);
+ Cell cell_2(sphere, 79);
+
+ cout <<"Cube Cell:\n";
+ for (VertexIterator v(cell_1); !v.end(); ++v)
+ {
+ cout <<"vertex index = " << v->index() <<"\t"
+ << "coordinates:\t" << v->point() << endl;
+ std::cout <<"\t\t\t\t\t";
+ print_point(v->point());
+ }
+
+ cout <<"Sphere Cell:\n";
+ for (VertexIterator v(cell_2); !v.end(); ++v)
+ {
+ cout <<"vertex index = " << v->index() <<"\t"
+ << "coordinates:\t" << v->point() << endl;
+ std::cout <<"\t\t\t\t\t";
+ print_point(v->point());
+ }
+
+ if (CGAL::do_intersect(Primitive_Traits<TetrahedronCellPrimitive,K>::datum(cell_1),
+ (Primitive_Traits<TetrahedronCellPrimitive,K>::datum(cell_2))))
+ cout <<"Cells do intersect." << endl;
+
+ Point_3 a(0.6666666666666666,0.6666666666666666,0.5);
+ Point_3 b(1,0.6666666666666666,0.5);
+ Point_3 c(1.0,1.0,0.5);
+ Point_3 d(1.0,1.0,1.0);
+
+ Point_3 A(0.9037749551350623, 0.9037749551350623, 0.9037749551350623);
+ Point_3 B(1.096225044864938, 0.9037749551350623, 0.9037749551350623);
+ Point_3 C(1.096225044864938, 0.9037749551350623,1.096225044864938);
+ Point_3 D(1.096225044864938, 1.096225044864938, 1.096225044864938);
+
+ Tetrahedron_3 tet_1(a,b,c,d);
+ Triangle_3 tri_1(a,b,d);
+ Tetrahedron_3 tet_2(A,B,C,D);
+
+ if (CGAL::do_intersect(tet_2, Triangle_3(tet_1[0],tet_1[1],tet_1[2])))
+ cout <<"Intersection of Tetraeder 1 with face 1" << endl;
+ else
+ cout <<"Intersection of Tetraeder 1 with face 1" << endl;
+
+ if (CGAL::do_intersect(tet_2, Triangle_3(tet_1[0],tet_1[1],tet_1[3])))
+ cout <<"Intersection of Tetraeder 1 with face 2" << endl;
+ else
+ cout <<"Intersection of Tetraeder 1 with face 2" << endl;
+
+ if (CGAL::do_intersect(tet_2, Triangle_3(tet_1[0],tet_1[2],tet_1[3])))
+ cout <<"Intersection of Tetraeder 1 with face 3" << endl;
+ else
+ cout <<"Intersection of Tetraeder 1 with face 3" << endl;
+
+ if (CGAL::do_intersect(tet_2, Triangle_3(tet_1[1],tet_1[2],tet_1[3])))
+ cout <<"Intersection of Tetraeder 1 with face 4" << endl;
+ else
+ cout <<"Intersection of Tetraeder 1 with face 4" << endl;
+
+ return 0;
+}
=== added directory 'sandbox/cgal/mesh_intersection/python'
=== added file 'sandbox/cgal/mesh_intersection/python/demo.py'
--- sandbox/cgal/mesh_intersection/python/demo.py 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/mesh_intersection/python/demo.py 2009-11-30 01:20:26 +0000
@@ -0,0 +1,91 @@
+"""This demo shows the intersection of the boundary of a unit square
+(omega1) with a unit circle (omega2) rotating around the center of the
+square."""
+
+__author__ = "Andre Massing (masssing@xxxxxxxxx)"
+__date__ = ""
+__copyright__ = "Copyright (C) 2009 Andre Massing"
+__license__ = "GNU LGPL Version 2.1"
+
+from dolfin import *
+from numpy import *
+
+#modify to check for CGAL
+#if not has_gts():
+# print "DOLFIN must be compiled with GTS to run this demo."
+# exit(0)
+
+# Create meshes (omega0 overlapped by omega1)
+#mesh0 = UnitSquare(3, 3)
+#mesh1 = UnitCircle(5)
+mesh0 = UnitCube(3, 3, 2)
+mesh1 = UnitSphere(3)
+
+print "Cube"
+print "Cells:"
+#print mesh0.cells()
+#print mesh0.coordinates()
+print mesh0.coordinates()[mesh0.cells()[102]]
+
+print
+
+print "Sphere"
+print "Cells:"
+print mesh1.coordinates()[mesh1.cells()[79]]
+#print mesh1.cells()
+#print mesh1.coordinates()
+
+#Create intersection operator for mesh0
+iop = IntersectionOperator(mesh0)
+
+# Access mesh geometry
+x = mesh1.coordinates()
+
+# Move and scale circle
+x *= 0.5
+x += 1.0
+
+# Iterate over angle
+theta = 0.0
+dtheta = 0.01*DOLFIN_PI
+intersection = MeshFunction("uint", mesh0, mesh0.topology().dim())
+_first = False
+
+cells = iop.all_intersected_entities(mesh1)
+p = plot(intersection, rescale=False)
+p.add_polygon([[0, 0, -0.01], [1, 0, -0.01], [1, 1, -0.01], [0, 1, -0.01], [0, 0, -0.01]])
+p.ren.ResetCamera()
+
+while theta > 2*DOLFIN_PI:
+
+ # Compute intersection with boundary of square
+ cells = iop.all_intersected_entities(mesh1)
+# boundary =BoundaryMesh(mesh1)
+# cells = iop.all_intersected_entities(boundary)
+
+ # Mark intersected values
+ intersection.values()[:] = 0
+ intersection.values()[cells] = 1
+
+ # Plot intersection
+ if _first:
+ p = plot(intersection, rescale=False)
+ p.add_polygon([[0, 0, -0.01], [1, 0, -0.01], [1, 1, -0.01], [0, 1, -0.01], [0, 0, -0.01]])
+ p.ren.ResetCamera()
+ _first = False
+ else:
+ plot(intersection)
+
+# interactive()
+
+
+ # Rotate circle around (0.5, 0.5)
+ xr = x[:, 0] - 0.5
+ yr = x[:, 1] - 0.5
+ x[:,0] = 0.5 + (cos(dtheta)*xr - sin(dtheta)*yr)
+ x[:,1] = 0.5 + (sin(dtheta)*xr + cos(dtheta)*yr)
+
+ theta += dtheta
+
+# Hold plot
+interactive()
=== added file 'sandbox/cgal/mesh_intersection/python/intersection_operator_test.py'
--- sandbox/cgal/mesh_intersection/python/intersection_operator_test.py 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/mesh_intersection/python/intersection_operator_test.py 2009-11-30 01:20:26 +0000
@@ -0,0 +1,48 @@
+"""This demo shows the intersection of the boundary of a unit square
+(omega1) with a unit circle (omega2) rotating around the center of the
+square."""
+
+__author__ = "Andre Massing (masssing@xxxxxxxxx)"
+__date__ = ""
+__copyright__ = "Copyright (C) 2009 Andre Massing"
+__license__ = "GNU LGPL Version 2.1"
+
+from dolfin import *
+from numpy import *
+
+#modify to check for CGAL
+#if not has_gts():
+# print "DOLFIN must be compiled with GTS to run this demo."
+# exit(0)
+
+# Create meshes (omega0 overlapped by omega1)
+#mesh0 = UnitSquare(3, 3)
+#mesh1 = UnitCircle(5)
+mesh0 = UnitCube(3, 3, 2)
+mesh1 = UnitSphere(3)
+
+print "Cube"
+print "Cells:"
+#print mesh0.cells()
+#print mesh0.coordinates()
+print mesh0.coordinates()[mesh0.cells()[102]]
+
+print
+
+#Create intersection operator for mesh0
+iop = IntersectionOperator(mesh0)
+
+# Access mesh geometry
+x = mesh1.coordinates()
+
+# Move and scale circle
+x *= 0.5
+x += 1.0
+
+print "Sphere"
+print "Cells:"
+print mesh1.coordinates()[mesh1.cells()[79]]
+print mesh1.cells()
+print mesh1.coordinates()
+
+cells = iop.all_intersected_entities(mesh1)
=== added directory 'sandbox/cgal/primitive_intersection'
=== added file 'sandbox/cgal/primitive_intersection/AABB_custom_indexed_triangle_set_example.cpp'
--- sandbox/cgal/primitive_intersection/AABB_custom_indexed_triangle_set_example.cpp 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/primitive_intersection/AABB_custom_indexed_triangle_set_example.cpp 2009-11-30 01:20:26 +0000
@@ -0,0 +1,170 @@
+// Copyright (c) 2009 INRIA Sophia-Antipolis (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org); you may redistribute it under
+// the terms of the Q Public License version 1.0.
+// See the file LICENSE.QPL distributed with CGAL.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL: $
+// $Id: $
+//
+//
+// Author(s) : Camille Wormser, Pierre Alliez
+//
+//******************************************************************************
+// File Description : Example of AABB tree used with indexed triangle set
+//
+//******************************************************************************
+
+#include <iostream>
+#include <list>
+#include <boost/iterator.hpp>
+
+#include <CGAL/AABB_tree.h> // must be inserted before kernel
+#include <CGAL/AABB_traits.h>
+
+#include <CGAL/Simple_cartesian.h>
+
+typedef CGAL::Simple_cartesian<double> K;
+
+
+// My own point type:
+struct My_point {
+ double x;
+ double y;
+ double z;
+
+ My_point (double _x, double _y, double _z)
+ : x(_x), y(_y), z(_z) {}
+};
+
+// The triangles are stored in a flat array of indices
+// referring to an array of points: three consecutive
+// indices represent a triangle.
+typedef std::vector<size_t>::const_iterator Point_index_iterator;
+
+// Let us now define the iterator on triangles that the tree needs:
+class Triangle_iterator
+ : public boost::iterator_adaptor<
+ Triangle_iterator // Derived
+ , Point_index_iterator // Base
+ , boost::use_default // Value
+ , boost::forward_traversal_tag // CategoryOrTraversal
+ >
+{
+public:
+ Triangle_iterator()
+ : Triangle_iterator::iterator_adaptor_() {}
+
+ explicit Triangle_iterator(Point_index_iterator p)
+ : Triangle_iterator::iterator_adaptor_(p) {}
+
+private:
+ friend class boost::iterator_core_access;
+ void increment() { this->base_reference() += 3; }
+};
+
+// The following primitive provides the conversion facilities between
+// my own triangle and point types and the CGAL ones
+struct My_triangle_primitive {
+public:
+ typedef int Id;
+ typedef Triangle_iterator Iter;
+
+ // the CGAL types returned
+ typedef K::Point_3 Point;
+ typedef K::Triangle_3 Datum;
+
+ // a static pointer to the vector containing the points
+ // is needed to build the triangles on the fly:
+ static const std::vector<My_point>* point_container;
+
+private:
+ static int counter;
+ Iter m_it; // this is what the AABB tree stores internally
+ Id index;
+
+public:
+ My_triangle_primitive() {} // default constructor needed
+
+ // the following constructor is the one that receives the iterators from the
+ // iterator range given as input to the AABB_tree
+ My_triangle_primitive(Triangle_iterator a)
+ : m_it(a), index(++counter) { std::cout <<"Copy number " << index; }
+
+ Id id() const { return index; }
+
+ // on the fly conversion from the internal data to the CGAL types
+ Datum datum() const
+ {
+ Point_index_iterator p_it = m_it.base();
+ const My_point& mp = (*point_container)[*p_it];
+ Point p(mp.x, mp.y, mp.z);
+ ++p_it;
+ const My_point& mq = (*point_container)[*p_it];
+ Point q(mq.x, mq.y, mq.z);
+ ++p_it;
+ const My_point& mr = (*point_container)[*p_it];
+ Point r(mr.x, mr.y, mr.z);
+
+ return Datum(p, q, r); // assembles triangle from three points
+ }
+
+ // one point which must be on the primitive
+ Point reference_point() const
+ {
+ const My_point& mp = (*point_container)[*m_it];
+ return Point(mp.x, mp.y, mp.z);
+ }
+};
+
+int My_triangle_primitive::counter = 0;
+
+
+// types
+typedef CGAL::AABB_traits<K, My_triangle_primitive> My_AABB_traits;
+typedef CGAL::AABB_tree<My_AABB_traits> Tree;
+const std::vector<My_point>* My_triangle_primitive::point_container = 0;
+
+int main()
+{
+ // generates point set
+ My_point a(1.0, 0.0, 0.0);
+ My_point b(0.0, 1.0, 0.0);
+ My_point c(0.0, 0.0, 1.0);
+ My_point d(0.0, 0.0, 0.0);
+
+ std::vector<My_point> points;
+ My_triangle_primitive::point_container = &points;
+ points.push_back(a);
+ points.push_back(b);
+ points.push_back(c);
+ points.push_back(d);
+
+ // generates indexed triangle set
+ std::vector<size_t> triangles;
+ triangles.push_back(0); triangles.push_back(1); triangles.push_back(2);
+ triangles.push_back(0); triangles.push_back(1); triangles.push_back(3);
+ triangles.push_back(0); triangles.push_back(3); triangles.push_back(2);
+
+ // constructs AABB tree
+ Tree tree(Triangle_iterator(triangles.begin()),
+ Triangle_iterator(triangles.end()));
+
+ // counts #intersections
+ K::Ray_3 ray_query(K::Point_3(0.2, 0.2, 0.2), K::Point_3(0.0, 1.0, 0.0));
+ std::cout << tree.number_of_intersected_primitives(ray_query)
+ << " intersections(s) with ray query" << std::endl;
+
+ // computes closest point
+ K::Point_3 point_query(2.0, 2.0, 2.0);
+ K::Point_3 closest_point = tree.closest_point(point_query);
+
+ return EXIT_SUCCESS;
+}
=== added file 'sandbox/cgal/primitive_intersection/SConstruct'
--- sandbox/cgal/primitive_intersection/SConstruct 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/primitive_intersection/SConstruct 2009-11-30 01:20:26 +0000
@@ -0,0 +1,16 @@
+import os
+
+# Get compiler from pkg-config
+compiler = commands.getoutput('pkg-config --variable=compiler dolfin')
+
+# Create a SCons Environment based on the main os environment
+env = Environment(ENV=os.environ, CXX=compiler)
+
+# Get compiler flags from pkg-config
+env.ParseConfig('pkg-config --cflags --libs dolfin')
+
+env.Program(['primitive_intersection_test.cpp'])
+env.Program(['tree_intersection_test.cpp'])
+env.Program(['tree_intersection_fail_example.cpp'])
+env.Program(['failed_tetrahedron_triangle_intersection.cpp'])
+env.Program(['failed_tetrahedron_triangle_intersection_2.cpp'])
=== added file 'sandbox/cgal/primitive_intersection/failed_tetrahedron_triangle_intersection.cpp'
--- sandbox/cgal/primitive_intersection/failed_tetrahedron_triangle_intersection.cpp 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/primitive_intersection/failed_tetrahedron_triangle_intersection.cpp 2009-11-30 01:20:26 +0000
@@ -0,0 +1,44 @@
+#include <CGAL/Bbox_3.h>
+//#include <CGAL/Simple_cartesian.h>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <iostream>
+
+//typedef CGAL::Simple_cartesian<double> K;
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+typedef K::Point_3 Point_3;
+typedef K::Segment_3 Segment_3;
+typedef K::Triangle_3 Triangle_3;
+typedef K::Tetrahedron_3 Tetrahedron_3;
+typedef K::Iso_cuboid_3 Iso_cuboid_3;
+typedef K::Ray_3 Ray_3;
+typedef K::Line_3 Line_3;
+
+typedef CGAL::Bbox_3 Bbox_3;
+
+using std::cout;
+using std::endl;
+
+int main()
+{
+ Point_3 A(0.9037749551350623, 0.9037749551350623, 0.9037749551350623);
+ Point_3 B(1.096225044864938, 0.9037749551350623, 0.9037749551350623);
+ Point_3 C(1.096225044864938, 0.9037749551350623,1.096225044864938);
+ Point_3 D(1.096225044864938, 1.096225044864938, 1.096225044864938);
+
+ Point_3 a(0.6666666666666666,0.6666666666666666,0.5);
+ Point_3 b(1,0.6666666666666666,0.5);
+ Point_3 c(1.0,1.0,1.0);
+
+ Triangle_3 tri_1(a,b,c);
+ Tetrahedron_3 tet_2(A,B,C,D);
+
+ if (CGAL::do_intersect(tet_2,tri_1))
+ cout <<"Intersection of Tetraeder with face triangle" << endl;
+ else
+ cout <<"NO Intersection." << endl;
+
+ return 0;
+}
=== added file 'sandbox/cgal/primitive_intersection/failed_tetrahedron_triangle_intersection_2.cpp'
--- sandbox/cgal/primitive_intersection/failed_tetrahedron_triangle_intersection_2.cpp 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/primitive_intersection/failed_tetrahedron_triangle_intersection_2.cpp 2009-11-30 01:20:26 +0000
@@ -0,0 +1,44 @@
+#include <CGAL/Bbox_3.h>
+//#include <CGAL/Simple_cartesian.h>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <iostream>
+
+typedef CGAL::Simple_cartesian<double> K;
+//typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+typedef K::Point_3 Point_3;
+typedef K::Segment_3 Segment_3;
+typedef K::Triangle_3 Triangle_3;
+typedef K::Tetrahedron_3 Tetrahedron_3;
+typedef K::Iso_cuboid_3 Iso_cuboid_3;
+typedef K::Ray_3 Ray_3;
+typedef K::Line_3 Line_3;
+
+typedef CGAL::Bbox_3 Bbox_3;
+
+using std::cout;
+using std::endl;
+
+int main()
+{
+ Point_3 A(0.3,0.0,0.0);
+ Point_3 B(0.3,0.0,0.1);
+ Point_3 C(0.3,0.1,0.1);
+ Point_3 D(0.4,0.1,0.1);
+
+ Point_3 a(0.2,0.0,0.0);
+ Point_3 b(0.2,0.1,0.0);
+ Point_3 c(0.3,0.1,0.1);
+
+ Triangle_3 tri_1(a,b,c);
+ Tetrahedron_3 tet_2(A,B,C,D);
+
+ if (CGAL::do_intersect(tet_2,tri_1))
+ cout <<"Intersection of Tetraeder with face triangle" << endl;
+ else
+ cout <<"NO Intersection." << endl;
+
+ return 0;
+}
=== added file 'sandbox/cgal/primitive_intersection/point_3_primitive_3_intersection_test.cpp'
--- sandbox/cgal/primitive_intersection/point_3_primitive_3_intersection_test.cpp 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/primitive_intersection/point_3_primitive_3_intersection_test.cpp 2009-11-30 01:20:26 +0000
@@ -0,0 +1,19 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/intersections.h>
+#include <CGAL/point_generators_3.h>
+#include <CGAL/Bbox_3.h>
+#include <CGAL/box_intersection_d.h>
+#include <CGAL/function_objects.h>
+#include <CGAL/Join_input_iterator.h>
+#include <CGAL/algorithm.h>
+#include <vector>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef Kernel::Point_3 Point_3;
+typedef Kernel::Triangle_3 Triangle_3;
+typedef std::vector<Triangle_3> Triangles;
+typedef std::vector<Point_3> Points;
+typedef Triangles::iterator Iterator;
+typedef CGAL::Box_intersection_d::Box_with_handle_d<double,3,Iterator> Box;
+
+
=== added file 'sandbox/cgal/primitive_intersection/primitive_intersection_test.cpp'
--- sandbox/cgal/primitive_intersection/primitive_intersection_test.cpp 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/primitive_intersection/primitive_intersection_test.cpp 2009-11-30 01:20:26 +0000
@@ -0,0 +1,185 @@
+#include <dolfin/mesh/CGAL_includes.h>
+
+#include <CGAL/AABB_triangle_primitive.h>
+
+#include <vector>
+#include <iostream>
+#include <string>
+
+//using dolfin::Point;
+
+//typedef CGAL::Simple_cartesian<double> K;
+
+///Definition of the used geometric primitives types.
+//typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+//typedef CGAL::Simple_cartesian<double> Kernel;
+//typedef K::Point_3 Point_3;
+//typedef K::Segment_3 Segment_3;
+//typedef K::Triangle_3 Triangle_3;
+//typedef K::Tetrahedron_3 Tetrahedron_3;
+//typedef K::Iso_cuboid_3 Iso_cuboid_3;
+//typedef K::Ray_3 Ray_3;
+//typedef K::Line_3 Line_3;
+
+//typedef CGAL::Bbox_3 Bbox_3;
+
+
+typedef std::vector<Point_3> Point_3s;
+//typedef std::vector<Point> Points;
+typedef std::vector<Segment_3> Segments;
+typedef std::vector<Triangle_3> Triangles;
+typedef std::vector<Tetrahedron_3> Tetrahedrons;
+
+typedef Point_3s::iterator Poi3_Iterator;
+//typedef Points::iterator Poi_Iterator;
+typedef Segments::iterator Seg_Iterator;
+typedef Triangles::iterator Tri_Iterator;
+typedef Tetrahedrons::iterator Tet_Iterator;
+
+//using namespace dolfin;
+//typedef dolfin::MeshPrimitive<TriangleCellPrimitive> CellPrimitive;
+typedef CGAL::AABB_triangle_primitive<K,Tri_Iterator> CellPrimitive;
+typedef CGAL::AABB_traits<K, CellPrimitive> AABB_triangle_traits;
+typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+
+
+using std::cout;
+using std::endl;
+using std::string;
+
+void print_message(const string & s1, const string & s2, bool result)
+{
+ if (result)
+ cout << s1 << " and " << s2 << " intersect" << endl;
+ else
+ cout << s1 << " and " << s2 << " DO NOT intersect" << endl;
+}
+
+int main()
+{
+ //Building points
+ //Unit points
+ Point_3 e0(0.0,0.0,0.0);
+ Point_3 e1(1.0,0.0,0.0);
+ Point_3 e2(0.0,1.0,0.0);
+ Point_3 e3(0.0,0.0,1.0);
+ Point_3 e4(1.0,1.0,1.0);
+
+
+ double frac = 1.0/3.0;
+ Point_3 a(frac, frac, 0.0);
+ Point_3 b(0.0, frac, 0.0);
+ Point_3 c(frac, 2*frac, 0.0);
+ Point_3 d(0.0, 2*frac, 0.0);
+
+ Point_3 f(0.0, frac, 0.0);
+ Point_3 g(-1.0,0.0,0.0);
+ Point_3 h(-1.0,1.0,0.0);
+
+ Point_3 x(-1.0, 0.0, 0.0);
+ Point_3 y(-1.0, 1.0, 0.0);
+
+ //points shifted by (2,0,0)
+ Point_3 A(3.0, 0.0, 0.0);
+ Point_3 B(2.0, 1.0, 0.0);
+ Point_3 C(2.0, 0.0, 1.0);
+ Point_3 D(2.0, 0.0, 0.0);
+
+// Point de0(0.0,0.0,0.0);
+// Point de1(1.0,0.0,0.0);
+// Point de2(0.0,1.0,0.0);
+// Point de3(0.0,0.0,1.0);
+// Point de4(1.0,1.0,1.0);
+
+ //vector of dolfin points to test intersection incoporation the conversion
+ //operator.
+
+ Segment_3 seg_1(e0,e1);
+ Triangle_3 tri_1(e0,e1,e2);
+ Tetrahedron_3 tet_1(e0,e1,e2,e3);
+ Tetrahedron_3 tet_2(A,B,C,D);
+ Iso_cuboid_3 cub_1(e0,e4);
+ Bbox_3 box_1(cub_1.bbox());
+
+ string poistring("Point");
+ string segstring("Segment");
+ string tristring("Triangle");
+ string tetstring("Tetrahedron");
+ string cubstring("Cuboid");
+ string boxstring("Bbox");
+
+ std::cout << "Principial intersection capabilities test: " << std::endl;
+ //point primitive intersection test
+ print_message(poistring,poistring,CGAL::do_intersect(e0,e0)); //FIXED
+ print_message(poistring,segstring,CGAL::do_intersect(e0,seg_1)); //FIXED
+ print_message(poistring,tristring,CGAL::do_intersect(e0,tri_1));
+ print_message(poistring,tetstring,CGAL::do_intersect(e0,tet_1)); //FIXED
+ print_message(poistring,cubstring,CGAL::do_intersect(e0,cub_1)); //FIXED
+ print_message(poistring,boxstring,CGAL::do_intersect(e0,box_1)); //FIXED
+
+ //segment primitive intersection test
+// print_message(segstring,segstring,CGAL::do_intersect(seg_1,seg_1)); //FAILED
+ print_message(segstring,tristring,CGAL::do_intersect(seg_1,tri_1));
+// print_message(segstring,tetstring,CGAL::do_intersect(seg_1,tet_1)); //FAILED
+ print_message(segstring,boxstring,CGAL::do_intersect(seg_1,box_1)); //FIXED
+
+ //triangles primitive intersection test
+ print_message(tristring,tristring,CGAL::do_intersect(tri_1,tri_1));
+ print_message(tristring,tetstring,CGAL::do_intersect(tri_1,tet_1));
+ print_message(tristring,boxstring,CGAL::do_intersect(tri_1,box_1)); //FAILED
+
+ //tetrahedron primitive intersection test
+ print_message(tetstring,tetstring,CGAL::do_intersect(tet_1,tet_1)); //FIXED
+ print_message(tetstring,boxstring,CGAL::do_intersect(box_1,tet_1)); //FIXED
+ print_message(tetstring,boxstring,CGAL::do_intersect(box_1,tet_1)); //FIXED
+ print_message(tetstring,boxstring,CGAL::do_intersect(box_1,tet_2)); //FIXED
+
+ std::cout << "Additional intersection test: " << std::endl;
+ Triangles triangles;
+ triangles.push_back(Triangle_3(e0,a,b));
+ triangles.push_back(Triangle_3(b,a,c));
+ triangles.push_back(Triangle_3(b,c,d));
+
+ Triangle_3 test_tri(b,x,y);
+
+ int trij = 0;
+ for (Tri_Iterator i = triangles.begin(); i != triangles.end(); ++i)
+ {
+ ++trij;
+ std::cout <<"Point " <<std::endl;
+ if (CGAL::do_intersect(b,*i))
+ std::cout <<"intersects with triangle " << trij <<std::endl;
+ else
+ std::cout <<"DOES NOT intersect with triangle " << trij <<std::endl;
+ std::cout <<"Triangle " <<std::endl;
+ if (CGAL::do_intersect(test_tri,*i))
+ std::cout <<"intersects with triangle " << trij <<std::endl;
+ else
+ std::cout <<"DOES NOT intersect with triangle " << trij <<std::endl;
+ }
+
+ //Now intersection using CGAL AABB_tree
+ Tree tree(triangles.begin(), triangles.end());
+ std::cout <<"Intersection detection using AABB tre" << std::endl;
+ if(tree.do_intersect(test_tri))
+ std::cout << std::endl <<"intersection(s) with triangle test_tri";
+ std::cout << std::endl << tree.number_of_intersected_primitives(test_tri)
+ << " intersection(s) with triangle test_tri " << trij << std::endl;
+
+
+// triangles.push_back(Triangle_3(e0,a,b)); //
+// triangles.push_back(Triangle_3(e0,e1,b)); //
+// triangles.push_back(Triangle_3(e0,e1,c)); //
+// triangles.push_back(Triangle_3(x,y,z)); //
+// triangles.push_back(Triangle_3(e0,x,y)); //
+// triangles.push_back(Triangle_3(d,x,y)); //
+// triangles.push_back(Triangle_3(e,f,g)); //
+// triangles.push_back(Triangle_3(A,B,C)); //
+// triangles.push_back(Triangle_3(x,y,d)); //
+// triangles.push_back(Triangle_3(x,y,D)); //
+
+// std::list<Tetrahedron> tetrahedrons;
+// tetrahedrons.push_back(Tetrahedron(a,b,c,d));
+// tetrahedrons.push_back(Tetrahedron(A,B,C,D));
+ return 0;
+}
=== added file 'sandbox/cgal/primitive_intersection/tetrahedron_triangle_intersection_test.cpp'
--- sandbox/cgal/primitive_intersection/tetrahedron_triangle_intersection_test.cpp 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/primitive_intersection/tetrahedron_triangle_intersection_test.cpp 2009-11-30 01:20:26 +0000
@@ -0,0 +1,128 @@
+#include <iostream>
+#include <list>
+
+#include <CGAL/AABB_tree.h> // must be inserted before kernel
+#include <CGAL/AABB_traits.h>
+#include "AABB_tetrahedron_primitive.h"
+#include <CGAL/AABB_triangle_primitive.h>
+
+#include <CGAL/Simple_cartesian.h>
+//#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/intersections.h>
+#include <CGAL/Bbox_3.h>
+#include <CGAL/box_intersection_d.h>
+#include <CGAL/function_objects.h>
+#include <CGAL/Join_input_iterator.h>
+#include <CGAL/algorithm.h>
+#include <vector>
+
+
+typedef CGAL::Simple_cartesian<double> K;
+//typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+typedef K::FT FT;
+typedef K::Ray_3 Ray;
+typedef K::Line_3 Line;
+typedef K::Point_3 Point;
+typedef K::Triangle_3 Triangle;
+typedef K::Tetrahedron_3 Tetrahedron;
+
+//types
+typedef std::list<Tetrahedron>::iterator Tet_iterator;
+typedef std::list<Triangle>::iterator Tri_iterator;
+typedef std::list<Point>::iterator P_iterator;
+
+//typedef CGAL::AABB_tetrahedron_primitive<K,Tet_iterator> Primitive;
+typedef CGAL::AABB_triangle_primitive<K,Tri_iterator> Primitive;
+typedef CGAL::AABB_traits<K, Primitive> AABB_traits;
+typedef CGAL::AABB_tree<AABB_traits> Tree;
+
+struct My_Point {
+ double x,y,z;
+ operator Point() const { return Point(x,y,z); }
+};
+
+Point shift_point(const My_Point & p)
+{
+ return p;
+}
+
+
+int main()
+{
+ Point a(1.0, 0.0, 0.0);
+ Point b(0.0, 1.0, 0.0);
+ Point c(0.0, 0.0, 1.0);
+ Point d(0.0, 0.0, 0.0);
+
+ //points shifted by (2,0,0)
+ Point A(3.0, 0.0, 0.0);
+ Point B(2.0, 1.0, 0.0);
+ Point C(2.0, 0.0, 1.0);
+ Point D(2.0, 0.0, 0.0);
+
+ //points shifted by (-1,0,0)
+ Point x(-1.0, 0.0, 0.0);
+ Point y(-1.0, 1.0, 0.0);
+ Point z(-1.0, 0.0, 1.0);
+
+ std::list<Tetrahedron> tetrahedrons;
+ std::list<Triangle> triangles;
+
+ tetrahedrons.push_back(Tetrahedron(a,b,c,d));
+ tetrahedrons.push_back(Tetrahedron(A,B,C,D));
+
+ triangles.push_back(Triangle(a,b,c)); //intersects 1st
+ triangles.push_back(Triangle(A,B,C)); //intersects 2nd
+ triangles.push_back(Triangle(x,y,z)); //intersects none
+ triangles.push_back(Triangle(x,y,d)); //intersects 1
+ triangles.push_back(Triangle(x,y,D)); //intersects 1st and 2nd
+
+ std::list<Point> points;
+
+ points.push_back(a);
+ points.push_back(b);
+ points.push_back(c);
+ points.push_back(d);
+ points.push_back(A);
+ points.push_back(B);
+ points.push_back(C);
+ points.push_back(D);
+
+ //construct AABB tree
+// Tree tree(tetrahedrons.begin(),tetrahedrons.end());
+ Tree tree(triangles.begin(),triangles.end());
+
+ //counts intersections for each triangle.
+// for (P_iterator i = points.begin(); i != points.end(); ++i)
+// {
+// ++j;
+// if(tree.do_intersect(*i))
+// std::cout << "intersection(s) with point " << j << std::endl;
+// else
+// std::cout << "no intersection with point" << j << std::endl;
+
+// std::cout << tree.number_of_intersected_primitives(*i)
+// << " intersection(s) with point number " << j << std::endl;
+// }
+
+ //counts intersections for each triangle.
+ int j = 0;
+ for (Tri_iterator i = triangles.begin(); i != triangles.end(); ++i)
+ {
+ ++j;
+ if(tree.do_intersect(*i))
+ std::cout << "intersection(s) with triangle " << j << std::endl;
+ else
+ std::cout << "no intersection with triangle" << j << std::endl;
+
+// //computes intersections with segment query
+ std::cout << tree.number_of_intersected_primitives(*i)
+ << " intersection(s) with triangle number " << j << std::endl;
+ }
+ std::cout << tree.number_of_intersected_primitives(triangles.front())
+ << " intersections(s) with triangle query" << std::endl;
+// std::cout << A;
+ return 0;
+}
=== added file 'sandbox/cgal/primitive_intersection/tetrahedron_triangle_intersection_test_2.cpp'
--- sandbox/cgal/primitive_intersection/tetrahedron_triangle_intersection_test_2.cpp 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/primitive_intersection/tetrahedron_triangle_intersection_test_2.cpp 2009-11-30 01:20:26 +0000
@@ -0,0 +1,128 @@
+#include <iostream>
+#include <list>
+
+#include <CGAL/AABB_tree.h> // must be inserted before kernel
+#include <CGAL/AABB_traits.h>
+#include "AABB_tetrahedron_primitive.h"
+#include <CGAL/AABB_triangle_primitive.h>
+
+#include <CGAL/Simple_cartesian.h>
+//#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/intersections.h>
+#include <CGAL/Bbox_3.h>
+#include <CGAL/box_intersection_d.h>
+#include <CGAL/function_objects.h>
+#include <CGAL/Join_input_iterator.h>
+#include <CGAL/algorithm.h>
+#include <vector>
+
+
+typedef CGAL::Simple_cartesian<double> K;
+//typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+typedef K::FT FT;
+typedef K::Ray_3 Ray;
+typedef K::Line_3 Line;
+typedef K::Point_3 Point;
+typedef K::Triangle_3 Triangle;
+typedef K::Tetrahedron_3 Tetrahedron;
+
+//types
+typedef std::list<Tetrahedron>::iterator Tet_iterator;
+typedef std::list<Triangle>::iterator Tri_iterator;
+typedef std::list<Point>::iterator P_iterator;
+
+//typedef CGAL::AABB_tetrahedron_primitive<K,Tet_iterator> Primitive;
+typedef CGAL::AABB_triangle_primitive<K,Tri_iterator> Primitive;
+typedef CGAL::AABB_traits<K, Primitive> AABB_traits;
+typedef CGAL::AABB_tree<AABB_traits> Tree;
+
+struct My_Point {
+ double x,y,z;
+ operator Point() const { return Point(x,y,z); }
+};
+
+Point shift_point(const My_Point & p)
+{
+ return p;
+}
+
+
+int main()
+{
+ Point a(1.0, 0.0, 0.0);
+ Point b(0.0, 1.0, 0.0);
+ Point c(0.0, 0.0, 1.0);
+ Point d(0.0, 0.0, 0.0);
+
+ //points shifted by (2,0,0)
+ Point A(3.0, 0.0, 0.0);
+ Point B(2.0, 1.0, 0.0);
+ Point C(2.0, 0.0, 1.0);
+ Point D(2.0, 0.0, 0.0);
+
+ //points shifted by (-1,0,0)
+ Point x(-1.0, 0.0, 0.0);
+ Point y(-1.0, 1.0, 0.0);
+ Point z(-1.0, 0.0, 1.0);
+
+ std::list<Tetrahedron> tetrahedrons;
+ std::list<Triangle> triangles;
+
+ tetrahedrons.push_back(Tetrahedron(a,b,c,d));
+ tetrahedrons.push_back(Tetrahedron(A,B,C,D));
+
+ triangles.push_back(Triangle(a,b,c)); //intersects 1st
+ triangles.push_back(Triangle(A,B,C)); //intersects 2nd
+ triangles.push_back(Triangle(x,y,z)); //intersects none
+ triangles.push_back(Triangle(x,y,d)); //intersects 1
+ triangles.push_back(Triangle(x,y,D)); //intersects 1st and 2nd
+
+ std::list<Point> points;
+
+ points.push_back(a);
+ points.push_back(b);
+ points.push_back(c);
+ points.push_back(d);
+ points.push_back(A);
+ points.push_back(B);
+ points.push_back(C);
+ points.push_back(D);
+
+ //construct AABB tree
+// Tree tree(tetrahedrons.begin(),tetrahedrons.end());
+ Tree tree(triangles.begin(),triangles.end());
+
+ //counts intersections for each triangle.
+// for (P_iterator i = points.begin(); i != points.end(); ++i)
+// {
+// ++j;
+// if(tree.do_intersect(*i))
+// std::cout << "intersection(s) with point " << j << std::endl;
+// else
+// std::cout << "no intersection with point" << j << std::endl;
+
+// std::cout << tree.number_of_intersected_primitives(*i)
+// << " intersection(s) with point number " << j << std::endl;
+// }
+
+ //counts intersections for each triangle.
+ int j = 0;
+ for (Tri_iterator i = triangles.begin(); i != triangles.end(); ++i)
+ {
+ ++j;
+ if(tree.do_intersect(*i))
+ std::cout << "intersection(s) with triangle " << j << std::endl;
+ else
+ std::cout << "no intersection with triangle" << j << std::endl;
+
+// //computes intersections with segment query
+ std::cout << tree.number_of_intersected_primitives(*i)
+ << " intersection(s) with triangle number " << j << std::endl;
+ }
+ std::cout << tree.number_of_intersected_primitives(triangles.front())
+ << " intersections(s) with triangle query" << std::endl;
+// std::cout << A;
+ return 0;
+}
=== added file 'sandbox/cgal/primitive_intersection/tree_intersection_fail_example.cpp'
--- sandbox/cgal/primitive_intersection/tree_intersection_fail_example.cpp 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/primitive_intersection/tree_intersection_fail_example.cpp 2009-11-30 01:20:26 +0000
@@ -0,0 +1,80 @@
+#include <CGAL/AABB_triangle_primitive.h>
+#include <CGAL/AABB_tree.h> // *Must* be inserted before kernel!
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+
+#include <CGAL/Bbox_3.h>
+
+#include <CGAL/Simple_cartesian.h> //used kernel
+typedef CGAL::Simple_cartesian<double> K;
+
+///Definition of the used geometric primitives types.
+typedef K::Point_3 Point_3;
+typedef K::Segment_3 Segment_3;
+typedef K::Triangle_3 Triangle_3;
+typedef K::Tetrahedron_3 Tetrahedron_3;
+typedef K::Iso_cuboid_3 Iso_cuboid_3;
+typedef K::Ray_3 Ray_3;
+typedef K::Line_3 Line_3;
+
+typedef CGAL::Bbox_3 Bbox_3;
+
+typedef std::vector<Triangle_3> Triangles;
+typedef Triangles::iterator Tri_Iterator;
+
+typedef CGAL::AABB_triangle_primitive<K,Tri_Iterator> CellPrimitive;
+typedef CGAL::AABB_traits<K, CellPrimitive> AABB_triangle_traits;
+typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+typedef Tree::Primitive_id Primitive_id;
+
+int main()
+{
+
+ double frac = 1.0/3.0;
+ Point_3 e0(0.0,0.0,0.0);
+ Point_3 a(frac, frac, 0.0);
+ Point_3 b(0.0, frac, 0.0);
+ Point_3 c(frac, 2*frac, 0.0);
+ Point_3 d(0.0, 2*frac, 0.0);
+ Point_3 x(-1.0, 0.0, 0.0);
+ Point_3 y(-1.0, 1.0, 0.0);
+
+ Triangles triangles;
+ triangles.push_back(Triangle_3(e0,a,b));
+ triangles.push_back(Triangle_3(b,a,c));
+ triangles.push_back(Triangle_3(b,c,d));
+
+ Triangle_3 test_tri(b,x,y);
+
+ std::cout <<"Intersection detection via CGAL::do_intersect:" << std::endl;
+ int trij = 0;
+ int counter = 0;
+ for (Tri_Iterator i = triangles.begin(); i != triangles.end(); ++i)
+ {
+ ++trij;
+ if (CGAL::do_intersect(test_tri,*i))
+ {
+ std::cout <<"Test triangle intersects with triangle " << trij <<std::endl;
+ ++counter;
+ }
+ else
+ std::cout <<"DOES NOT intersect with triangle " << trij <<std::endl;
+ }
+ std::cout <<"------------------------------" <<std::endl;
+ std::cout <<counter << " intersection(s) found." <<std::endl;
+
+
+ //Now intersection using CGAL AABB_tree
+ Tree tree(triangles.begin(), triangles.end());
+ std::cout << std::endl
+ <<"Intersection detection using AABB tree" << std::endl;
+ std::list<Primitive_id> primitives;
+ tree.all_intersected_primitives(test_tri, std::back_inserter(primitives));
+ for (std::list<Primitive_id>::const_iterator i = primitives.begin(); i != primitives.end(); ++i)
+ std::cout <<"Test triangle intersects with triangle " << (*i - triangles.begin()) << std::endl;
+ std::cout <<"------------------------------" <<std::endl;
+ std::cout << tree.number_of_intersected_primitives(test_tri)
+ << " intersection(s) found." << std::endl;
+
+ return 0;
+}
=== added file 'sandbox/cgal/primitive_intersection/tree_intersection_test.cpp'
--- sandbox/cgal/primitive_intersection/tree_intersection_test.cpp 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/primitive_intersection/tree_intersection_test.cpp 2009-11-30 01:20:26 +0000
@@ -0,0 +1,97 @@
+#include <CGAL/AABB_tree.h> // must be inserted before kernel
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+
+#include <CGAL/Simple_cartesian.h>
+
+//#include <dolfin/mesh/MeshPrimitives.h>
+
+//using namespace dolfin;
+//typedef Triangle_3 Triangle_3;
+
+//typedef Primitive<2> Triangle_Primitive;
+//typedef CGAL::AABB_traits<K,Primitive<2> > AABB_PrimitiveTraits;
+//typedef CGAL::AABB_tree<AABB_PrimitiveTraits> Tree;
+//typedef Point_3 Point;
+
+typedef CGAL::Simple_cartesian<double> K;
+typedef K::Tetrahedron_3 Tetrahedron;
+typedef K::Triangle_3 Triangle_3;
+typedef K::Point_3 Point_3;
+
+typedef std::list<Triangle_3>::iterator Tri_iterator;
+typedef CGAL::AABB_triangle_primitive<K,Tri_iterator> Primitive;
+typedef CGAL::AABB_traits<K, Primitive> AABB_triangle_traits;
+typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+
+typedef std::list<Point_3>::iterator P_iterator;
+
+int main()
+{
+ //Building points
+ Point_3 a(1.0, 0.0, 0.0);
+ Point_3 b(0.0, 1.0, 0.0);
+ Point_3 c(0.0, 0.0, 1.0);
+ Point_3 d(0.0, 0.0, 0.0);
+
+ //points shifted by (2,0,0)
+ Point_3 A(3.0, 0.0, 0.0);
+ Point_3 B(2.0, 1.0, 0.0);
+ Point_3 C(2.0, 0.0, 1.0);
+ Point_3 D(2.0, 0.0, 0.0);
+
+ //points shifted by (-1,0,0)
+ Point_3 x(-1.0, 0.0, 0.0);
+ Point_3 y(-1.0, 1.0, 0.0);
+ Point_3 z(-1.0, 0.0, 1.0);
+
+ std::list<Point_3> points;
+ points.push_back(a);
+ points.push_back(b);
+ points.push_back(c);
+ points.push_back(d);
+ points.push_back(A);
+ points.push_back(B);
+ points.push_back(C);
+ points.push_back(D);
+
+ std::list<Triangle_3> triangles;
+ triangles.push_back(Triangle_3(a,b,c)); //intersects 1st
+ triangles.push_back(Triangle_3(A,B,C)); //intersects 2nd
+ triangles.push_back(Triangle_3(x,y,z)); //intersects none
+ triangles.push_back(Triangle_3(x,y,d)); //intersects 1
+ triangles.push_back(Triangle_3(x,y,D)); //intersects 1st and 2nd
+
+// std::list<Tetrahedron> tetrahedrons;
+// tetrahedrons.push_back(Tetrahedron(a,b,c,d));
+// tetrahedrons.push_back(Tetrahedron(A,B,C,D));
+
+ Tree tree(triangles.begin(),triangles.end());
+
+// for (P_iterator i = points.begin(); i != points.end(); ++i)
+// {
+// ++j;
+// if(tree.do_intersect(*i))
+// std::cout << "intersection(s) with point " << j << std::endl;
+// else
+// std::cout << "no intersection with point" << j << std::endl;
+
+// std::cout << tree.number_of_intersected_primitives(*i)
+// << " intersection(s) with triangle number " << j << std::endl;
+// }
+
+// counts intersections for each triangle.
+ int j = 0;
+ for (Tri_iterator i = triangles.begin(); i != triangles.end(); ++i)
+ {
+ ++j;
+ if(tree.do_intersect(*i))
+ std::cout << "intersection(s) with triangle " << j << std::endl;
+ else
+ std::cout << "no intersection with triangle" << j << std::endl;
+
+ std::cout << tree.number_of_intersected_primitives(*i)
+ << " intersection(s) with triangle number " << j << std::endl;
+ }
+ return 0;
+}
=== added directory 'sandbox/cgal/test'
=== added file 'sandbox/cgal/test/P1.h'
--- sandbox/cgal/test/P1.h 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/test/P1.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,817 @@
+// This code conforms with the UFC specification version 1.0
+// and was automatically generated by FFC version 0.7.0.
+//
+// Warning: This code was generated with the option '-l dolfin'
+// and contains DOLFIN-specific wrappers that depend on DOLFIN.
+
+#ifndef __P1_H
+#define __P1_H
+
+#include <cmath>
+#include <stdexcept>
+#include <fstream>
+#include <ufc.h>
+
+/// This class defines the interface for a finite element.
+
+class p1_0_finite_element_0: public ufc::finite_element
+{
+public:
+
+ /// Constructor
+ p1_0_finite_element_0() : ufc::finite_element()
+ {
+ // Do nothing
+ }
+
+ /// Destructor
+ virtual ~p1_0_finite_element_0()
+ {
+ // Do nothing
+ }
+
+ /// Return a string identifying the finite element
+ virtual const char* signature() const
+ {
+ return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
+ }
+
+ /// Return the cell shape
+ virtual ufc::shape cell_shape() const
+ {
+ return ufc::triangle;
+ }
+
+ /// Return the dimension of the finite element function space
+ virtual unsigned int space_dimension() const
+ {
+ return 3;
+ }
+
+ /// Return the rank of the value space
+ virtual unsigned int value_rank() const
+ {
+ return 0;
+ }
+
+ /// Return the dimension of the value space for axis i
+ virtual unsigned int value_dimension(unsigned int i) const
+ {
+ return 1;
+ }
+
+ /// Evaluate basis function i at given point in cell
+ virtual void evaluate_basis(unsigned int i,
+ double* values,
+ const double* coordinates,
+ const ufc::cell& c) const
+ {
+ // Extract vertex coordinates
+ const double * const * element_coordinates = c.coordinates;
+
+ // Compute Jacobian of affine map from reference cell
+ const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
+ const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
+ const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
+ const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
+
+ // Compute determinant of Jacobian
+ const double detJ = J_00*J_11 - J_01*J_10;
+
+ // Compute inverse of Jacobian
+
+ // Get coordinates and map to the reference (UFC) element
+ double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
+ element_coordinates[0][0]*element_coordinates[2][1] +\
+ J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
+ double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
+ element_coordinates[1][0]*element_coordinates[0][1] -\
+ J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
+
+ // Map coordinates to the reference square
+ if (std::abs(y - 1.0) < 1e-14)
+ x = -1.0;
+ else
+ x = 2.0 *x/(1.0 - y) - 1.0;
+ y = 2.0*y - 1.0;
+
+ // Reset values
+ *values = 0;
+
+ // Map degree of freedom to element degree of freedom
+ const unsigned int dof = i;
+
+ // Generate scalings
+ const double scalings_y_0 = 1;
+ const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
+
+ // Compute psitilde_a
+ const double psitilde_a_0 = 1;
+ const double psitilde_a_1 = x;
+
+ // Compute psitilde_bs
+ const double psitilde_bs_0_0 = 1;
+ const double psitilde_bs_0_1 = 1.5*y + 0.5;
+ const double psitilde_bs_1_0 = 1;
+
+ // Compute basisvalues
+ const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
+ const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
+ const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
+
+ // Table(s) of coefficients
+ static const double coefficients0[3][3] = \
+ {{0.471404520791032, -0.288675134594813, -0.166666666666667},
+ {0.471404520791032, 0.288675134594813, -0.166666666666667},
+ {0.471404520791032, 0, 0.333333333333333}};
+
+ // Extract relevant coefficients
+ const double coeff0_0 = coefficients0[dof][0];
+ const double coeff0_1 = coefficients0[dof][1];
+ const double coeff0_2 = coefficients0[dof][2];
+
+ // Compute value(s)
+ *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;
+ }
+
+ /// Evaluate all basis functions at given point in cell
+ virtual void evaluate_basis_all(double* values,
+ const double* coordinates,
+ const ufc::cell& c) const
+ {
+ throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
+ }
+
+ /// Evaluate order n derivatives of basis function i at given point in cell
+ virtual void evaluate_basis_derivatives(unsigned int i,
+ unsigned int n,
+ double* values,
+ const double* coordinates,
+ const ufc::cell& c) const
+ {
+ // Extract vertex coordinates
+ const double * const * element_coordinates = c.coordinates;
+
+ // Compute Jacobian of affine map from reference cell
+ const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
+ const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
+ const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
+ const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
+
+ // Compute determinant of Jacobian
+ const double detJ = J_00*J_11 - J_01*J_10;
+
+ // Compute inverse of Jacobian
+
+ // Get coordinates and map to the reference (UFC) element
+ double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
+ element_coordinates[0][0]*element_coordinates[2][1] +\
+ J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
+ double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
+ element_coordinates[1][0]*element_coordinates[0][1] -\
+ J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
+
+ // Map coordinates to the reference square
+ if (std::abs(y - 1.0) < 1e-14)
+ x = -1.0;
+ else
+ x = 2.0 *x/(1.0 - y) - 1.0;
+ y = 2.0*y - 1.0;
+
+ // Compute number of derivatives
+ unsigned int num_derivatives = 1;
+
+ for (unsigned int j = 0; j < n; j++)
+ num_derivatives *= 2;
+
+
+ // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
+ unsigned int **combinations = new unsigned int *[num_derivatives];
+
+ for (unsigned int j = 0; j < num_derivatives; j++)
+ {
+ combinations[j] = new unsigned int [n];
+ for (unsigned int k = 0; k < n; k++)
+ combinations[j][k] = 0;
+ }
+
+ // Generate combinations of derivatives
+ for (unsigned int row = 1; row < num_derivatives; row++)
+ {
+ for (unsigned int num = 0; num < row; num++)
+ {
+ for (unsigned int col = n-1; col+1 > 0; col--)
+ {
+ if (combinations[row][col] + 1 > 1)
+ combinations[row][col] = 0;
+ else
+ {
+ combinations[row][col] += 1;
+ break;
+ }
+ }
+ }
+ }
+
+ // Compute inverse of Jacobian
+ const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
+
+ // Declare transformation matrix
+ // Declare pointer to two dimensional array and initialise
+ double **transform = new double *[num_derivatives];
+
+ for (unsigned int j = 0; j < num_derivatives; j++)
+ {
+ transform[j] = new double [num_derivatives];
+ for (unsigned int k = 0; k < num_derivatives; k++)
+ transform[j][k] = 1;
+ }
+
+ // Construct transformation matrix
+ for (unsigned int row = 0; row < num_derivatives; row++)
+ {
+ for (unsigned int col = 0; col < num_derivatives; col++)
+ {
+ for (unsigned int k = 0; k < n; k++)
+ transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
+ }
+ }
+
+ // Reset values
+ for (unsigned int j = 0; j < 1*num_derivatives; j++)
+ values[j] = 0;
+
+ // Map degree of freedom to element degree of freedom
+ const unsigned int dof = i;
+
+ // Generate scalings
+ const double scalings_y_0 = 1;
+ const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
+
+ // Compute psitilde_a
+ const double psitilde_a_0 = 1;
+ const double psitilde_a_1 = x;
+
+ // Compute psitilde_bs
+ const double psitilde_bs_0_0 = 1;
+ const double psitilde_bs_0_1 = 1.5*y + 0.5;
+ const double psitilde_bs_1_0 = 1;
+
+ // Compute basisvalues
+ const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
+ const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
+ const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
+
+ // Table(s) of coefficients
+ static const double coefficients0[3][3] = \
+ {{0.471404520791032, -0.288675134594813, -0.166666666666667},
+ {0.471404520791032, 0.288675134594813, -0.166666666666667},
+ {0.471404520791032, 0, 0.333333333333333}};
+
+ // Interesting (new) part
+ // Tables of derivatives of the polynomial base (transpose)
+ static const double dmats0[3][3] = \
+ {{0, 0, 0},
+ {4.89897948556636, 0, 0},
+ {0, 0, 0}};
+
+ static const double dmats1[3][3] = \
+ {{0, 0, 0},
+ {2.44948974278318, 0, 0},
+ {4.24264068711928, 0, 0}};
+
+ // Compute reference derivatives
+ // Declare pointer to array of derivatives on FIAT element
+ double *derivatives = new double [num_derivatives];
+
+ // Declare coefficients
+ double coeff0_0 = 0;
+ double coeff0_1 = 0;
+ double coeff0_2 = 0;
+
+ // Declare new coefficients
+ double new_coeff0_0 = 0;
+ double new_coeff0_1 = 0;
+ double new_coeff0_2 = 0;
+
+ // Loop possible derivatives
+ for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
+ {
+ // Get values from coefficients array
+ new_coeff0_0 = coefficients0[dof][0];
+ new_coeff0_1 = coefficients0[dof][1];
+ new_coeff0_2 = coefficients0[dof][2];
+
+ // Loop derivative order
+ for (unsigned int j = 0; j < n; j++)
+ {
+ // Update old coefficients
+ coeff0_0 = new_coeff0_0;
+ coeff0_1 = new_coeff0_1;
+ coeff0_2 = new_coeff0_2;
+
+ if(combinations[deriv_num][j] == 0)
+ {
+ new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];
+ new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];
+ new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];
+ }
+ if(combinations[deriv_num][j] == 1)
+ {
+ new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];
+ new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];
+ new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];
+ }
+
+ }
+ // Compute derivatives on reference element as dot product of coefficients and basisvalues
+ derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;
+ }
+
+ // Transform derivatives back to physical element
+ for (unsigned int row = 0; row < num_derivatives; row++)
+ {
+ for (unsigned int col = 0; col < num_derivatives; col++)
+ {
+ values[row] += transform[row][col]*derivatives[col];
+ }
+ }
+ // Delete pointer to array of derivatives on FIAT element
+ delete [] derivatives;
+
+ // Delete pointer to array of combinations of derivatives and transform
+ for (unsigned int row = 0; row < num_derivatives; row++)
+ {
+ delete [] combinations[row];
+ delete [] transform[row];
+ }
+
+ delete [] combinations;
+ delete [] transform;
+ }
+
+ /// Evaluate order n derivatives of all basis functions at given point in cell
+ virtual void evaluate_basis_derivatives_all(unsigned int n,
+ double* values,
+ const double* coordinates,
+ const ufc::cell& c) const
+ {
+ throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
+ }
+
+ /// Evaluate linear functional for dof i on the function f
+ virtual double evaluate_dof(unsigned int i,
+ const ufc::function& f,
+ const ufc::cell& c) const
+ {
+ // The reference points, direction and weights:
+ static const double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};
+ static const double W[3][1] = {{1}, {1}, {1}};
+ static const double D[3][1][1] = {{{1}}, {{1}}, {{1}}};
+
+ const double * const * x = c.coordinates;
+ double result = 0.0;
+ // Iterate over the points:
+ // Evaluate basis functions for affine mapping
+ const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
+ const double w1 = X[i][0][0];
+ const double w2 = X[i][0][1];
+
+ // Compute affine mapping y = F(X)
+ double y[2];
+ y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
+ y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
+
+ // Evaluate function at physical points
+ double values[1];
+ f.evaluate(values, y, c);
+
+ // Map function values using appropriate mapping
+ // Affine map: Do nothing
+
+ // Note that we do not map the weights (yet).
+
+ // Take directional components
+ for(int k = 0; k < 1; k++)
+ result += values[k]*D[i][0][k];
+ // Multiply by weights
+ result *= W[i][0];
+
+ return result;
+ }
+
+ /// Evaluate linear functionals for all dofs on the function f
+ virtual void evaluate_dofs(double* values,
+ const ufc::function& f,
+ const ufc::cell& c) const
+ {
+ throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
+ }
+
+ /// Interpolate vertex values from dof values
+ virtual void interpolate_vertex_values(double* vertex_values,
+ const double* dof_values,
+ const ufc::cell& c) const
+ {
+ // Evaluate at vertices and use affine mapping
+ vertex_values[0] = dof_values[0];
+ vertex_values[1] = dof_values[1];
+ vertex_values[2] = dof_values[2];
+ }
+
+ /// Return the number of sub elements (for a mixed element)
+ virtual unsigned int num_sub_elements() const
+ {
+ return 1;
+ }
+
+ /// Create a new finite element for sub element i (for a mixed element)
+ virtual ufc::finite_element* create_sub_element(unsigned int i) const
+ {
+ return new p1_0_finite_element_0();
+ }
+
+};
+
+/// This class defines the interface for a local-to-global mapping of
+/// degrees of freedom (dofs).
+
+class p1_0_dof_map_0: public ufc::dof_map
+{
+private:
+
+ unsigned int __global_dimension;
+
+public:
+
+ /// Constructor
+ p1_0_dof_map_0() : ufc::dof_map()
+ {
+ __global_dimension = 0;
+ }
+
+ /// Destructor
+ virtual ~p1_0_dof_map_0()
+ {
+ // Do nothing
+ }
+
+ /// Return a string identifying the dof map
+ virtual const char* signature() const
+ {
+ return "FFC dof map for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
+ }
+
+ /// Return true iff mesh entities of topological dimension d are needed
+ virtual bool needs_mesh_entities(unsigned int d) const
+ {
+ switch ( d )
+ {
+ case 0:
+ return true;
+ break;
+ case 1:
+ return false;
+ break;
+ case 2:
+ return false;
+ break;
+ }
+ return false;
+ }
+
+ /// Initialize dof map for mesh (return true iff init_cell() is needed)
+ virtual bool init_mesh(const ufc::mesh& m)
+ {
+ __global_dimension = m.num_entities[0];
+ return false;
+ }
+
+ /// Initialize dof map for given cell
+ virtual void init_cell(const ufc::mesh& m,
+ const ufc::cell& c)
+ {
+ // Do nothing
+ }
+
+ /// Finish initialization of dof map for cells
+ virtual void init_cell_finalize()
+ {
+ // Do nothing
+ }
+
+ /// Return the dimension of the global finite element function space
+ virtual unsigned int global_dimension() const
+ {
+ return __global_dimension;
+ }
+
+ /// Return the dimension of the local finite element function space for a cell
+ virtual unsigned int local_dimension(const ufc::cell& c) const
+ {
+ return 3;
+ }
+
+ /// Return the maximum dimension of the local finite element function space
+ virtual unsigned int max_local_dimension() const
+ {
+ return 3;
+ }
+
+ // Return the geometric dimension of the coordinates this dof map provides
+ virtual unsigned int geometric_dimension() const
+ {
+ return 2;
+ }
+
+ /// Return the number of dofs on each cell facet
+ virtual unsigned int num_facet_dofs() const
+ {
+ return 2;
+ }
+
+ /// Return the number of dofs associated with each cell entity of dimension d
+ virtual unsigned int num_entity_dofs(unsigned int d) const
+ {
+ throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
+ }
+
+ /// Tabulate the local-to-global mapping of dofs on a cell
+ virtual void tabulate_dofs(unsigned int* dofs,
+ const ufc::mesh& m,
+ const ufc::cell& c) const
+ {
+ dofs[0] = c.entity_indices[0][0];
+ dofs[1] = c.entity_indices[0][1];
+ dofs[2] = c.entity_indices[0][2];
+ }
+
+ /// Tabulate the local-to-local mapping from facet dofs to cell dofs
+ virtual void tabulate_facet_dofs(unsigned int* dofs,
+ unsigned int facet) const
+ {
+ switch ( facet )
+ {
+ case 0:
+ dofs[0] = 1;
+ dofs[1] = 2;
+ break;
+ case 1:
+ dofs[0] = 0;
+ dofs[1] = 2;
+ break;
+ case 2:
+ dofs[0] = 0;
+ dofs[1] = 1;
+ break;
+ }
+ }
+
+ /// Tabulate the local-to-local mapping of dofs on entity (d, i)
+ virtual void tabulate_entity_dofs(unsigned int* dofs,
+ unsigned int d, unsigned int i) const
+ {
+ throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
+ }
+
+ /// Tabulate the coordinates of all dofs on a cell
+ virtual void tabulate_coordinates(double** coordinates,
+ const ufc::cell& c) const
+ {
+ const double * const * x = c.coordinates;
+ coordinates[0][0] = x[0][0];
+ coordinates[0][1] = x[0][1];
+ coordinates[1][0] = x[1][0];
+ coordinates[1][1] = x[1][1];
+ coordinates[2][0] = x[2][0];
+ coordinates[2][1] = x[2][1];
+ }
+
+ /// Return the number of sub dof maps (for a mixed element)
+ virtual unsigned int num_sub_dof_maps() const
+ {
+ return 1;
+ }
+
+ /// Create a new dof_map for sub dof map i (for a mixed element)
+ virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
+ {
+ return new p1_0_dof_map_0();
+ }
+
+};
+
+/// This class defines the interface for the assembly of the global
+/// tensor corresponding to a form with r + n arguments, that is, a
+/// mapping
+///
+/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
+///
+/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
+/// global tensor A is defined by
+///
+/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
+///
+/// where each argument Vj represents the application to the
+/// sequence of basis functions of Vj and w1, w2, ..., wn are given
+/// fixed functions (coefficients).
+
+class p1_form_0: public ufc::form
+{
+public:
+
+ /// Constructor
+ p1_form_0() : ufc::form()
+ {
+ // Do nothing
+ }
+
+ /// Destructor
+ virtual ~p1_form_0()
+ {
+ // Do nothing
+ }
+
+ /// Return a string identifying the form
+ virtual const char* signature() const
+ {
+ return "None";
+ }
+
+ /// Return the rank of the global tensor (r)
+ virtual unsigned int rank() const
+ {
+ return -1;
+ }
+
+ /// Return the number of coefficients (n)
+ virtual unsigned int num_coefficients() const
+ {
+ return 0;
+ }
+
+ /// Return the number of cell integrals
+ virtual unsigned int num_cell_integrals() const
+ {
+ return 0;
+ }
+
+ /// Return the number of exterior facet integrals
+ virtual unsigned int num_exterior_facet_integrals() const
+ {
+ return 0;
+ }
+
+ /// Return the number of interior facet integrals
+ virtual unsigned int num_interior_facet_integrals() const
+ {
+ return 0;
+ }
+
+ /// Create a new finite element for argument function i
+ virtual ufc::finite_element* create_finite_element(unsigned int i) const
+ {
+ return 0;
+ }
+
+ /// Create a new dof map for argument function i
+ virtual ufc::dof_map* create_dof_map(unsigned int i) const
+ {
+ return 0;
+ }
+
+ /// Create a new cell integral on sub domain i
+ virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
+ {
+ return 0;
+ }
+
+ /// Create a new exterior facet integral on sub domain i
+ virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
+ {
+ return 0;
+ }
+
+ /// Create a new interior facet integral on sub domain i
+ virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
+ {
+ return 0;
+ }
+
+};
+
+// DOLFIN wrappers
+
+// Standard library includes
+#include <string>
+
+// DOLFIN includes
+#include <dolfin/common/NoDeleter.h>
+#include <dolfin/fem/FiniteElement.h>
+#include <dolfin/fem/DofMap.h>
+#include <dolfin/fem/Form.h>
+#include <dolfin/function/FunctionSpace.h>
+#include <dolfin/function/GenericFunction.h>
+#include <dolfin/function/CoefficientAssigner.h>
+
+namespace P1
+{
+
+class Form_0_FunctionSpace_0: public dolfin::FunctionSpace
+{
+public:
+
+ Form_0_FunctionSpace_0(const dolfin::Mesh& mesh):
+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new p1_0_finite_element_0()))),
+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new p1_0_dof_map_0()), mesh)))
+ {
+ // Do nothing
+ }
+
+ Form_0_FunctionSpace_0(dolfin::Mesh& mesh):
+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new p1_0_finite_element_0()))),
+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new p1_0_dof_map_0()), mesh)))
+ {
+ // Do nothing
+ }
+
+ Form_0_FunctionSpace_0(boost::shared_ptr<dolfin::Mesh> mesh):
+ dolfin::FunctionSpace(mesh,
+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new p1_0_finite_element_0()))),
+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new p1_0_dof_map_0()), *mesh)))
+ {
+ // Do nothing
+ }
+
+ Form_0_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh):
+ dolfin::FunctionSpace(mesh,
+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new p1_0_finite_element_0()))),
+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new p1_0_dof_map_0()), *mesh)))
+ {
+ // Do nothing
+ }
+
+
+ ~Form_0_FunctionSpace_0()
+ {
+ }
+
+};
+
+class Form_0: public dolfin::Form
+{
+public:
+
+ // Constructor
+ Form_0(const dolfin::FunctionSpace& V0):
+ dolfin::Form(1, 0)
+ {
+ _function_spaces[0] = reference_to_no_delete_pointer(V0);
+
+ _ufc_form = boost::shared_ptr<const ufc::form>(new p1_form_0());
+ }
+
+ // Constructor
+ Form_0(boost::shared_ptr<const dolfin::FunctionSpace> V0):
+ dolfin::Form(1, 0)
+ {
+ _function_spaces[0] = V0;
+
+ _ufc_form = boost::shared_ptr<const ufc::form>(new p1_form_0());
+ }
+
+ // Destructor
+ ~Form_0()
+ {}
+
+ /// Return the number of the coefficient with this name
+ virtual dolfin::uint coefficient_number(const std::string& name) const
+ {
+
+ dolfin::error("No coefficients.");
+ return 0;
+ }
+
+ /// Return the name of the coefficient with this number
+ virtual std::string coefficient_name(dolfin::uint i) const
+ {
+
+ dolfin::error("No coefficients.");
+ return "unnamed";
+ }
+
+ // Typedefs
+ typedef Form_0_FunctionSpace_0 TestSpace;
+
+ // Coefficients
+};
+
+// Class typedefs
+typedef Form_0 LinearForm;
+typedef Form_0::TestSpace FunctionSpace;
+
+} // namespace P1
+
+#endif
=== added file 'sandbox/cgal/test/P1.ufl'
--- sandbox/cgal/test/P1.ufl 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/test/P1.ufl 2009-11-30 01:20:26 +0000
@@ -0,0 +1,13 @@
+# Copyright (C) 2009 Garth N. Wells.
+# Licensed under the GNU LGPL Version 2.1.
+#
+# First added: 2009-06-18
+# Last changed:
+#
+# The bilinear form a(v, u) and linear form L(v) for
+# projection onto piecewise quadratics.
+#
+# Compile this form with FFC: ffc -l dolfin P1.ufl
+
+element = FiniteElement("Lagrange", "triangle", 1)
+
=== added file 'sandbox/cgal/test/P3.h'
--- sandbox/cgal/test/P3.h 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/test/P3.h 2009-11-30 01:20:26 +0000
@@ -0,0 +1,959 @@
+// This code conforms with the UFC specification version 1.0
+// and was automatically generated by FFC version 0.7.0.
+//
+// Warning: This code was generated with the option '-l dolfin'
+// and contains DOLFIN-specific wrappers that depend on DOLFIN.
+
+#ifndef __P3_H
+#define __P3_H
+
+#include <cmath>
+#include <stdexcept>
+#include <fstream>
+#include <ufc.h>
+
+/// This class defines the interface for a finite element.
+
+class p3_0_finite_element_0: public ufc::finite_element
+{
+public:
+
+ /// Constructor
+ p3_0_finite_element_0() : ufc::finite_element()
+ {
+ // Do nothing
+ }
+
+ /// Destructor
+ virtual ~p3_0_finite_element_0()
+ {
+ // Do nothing
+ }
+
+ /// Return a string identifying the finite element
+ virtual const char* signature() const
+ {
+ return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 3)";
+ }
+
+ /// Return the cell shape
+ virtual ufc::shape cell_shape() const
+ {
+ return ufc::triangle;
+ }
+
+ /// Return the dimension of the finite element function space
+ virtual unsigned int space_dimension() const
+ {
+ return 10;
+ }
+
+ /// Return the rank of the value space
+ virtual unsigned int value_rank() const
+ {
+ return 0;
+ }
+
+ /// Return the dimension of the value space for axis i
+ virtual unsigned int value_dimension(unsigned int i) const
+ {
+ return 1;
+ }
+
+ /// Evaluate basis function i at given point in cell
+ virtual void evaluate_basis(unsigned int i,
+ double* values,
+ const double* coordinates,
+ const ufc::cell& c) const
+ {
+ // Extract vertex coordinates
+ const double * const * element_coordinates = c.coordinates;
+
+ // Compute Jacobian of affine map from reference cell
+ const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
+ const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
+ const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
+ const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
+
+ // Compute determinant of Jacobian
+ const double detJ = J_00*J_11 - J_01*J_10;
+
+ // Compute inverse of Jacobian
+
+ // Get coordinates and map to the reference (UFC) element
+ double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
+ element_coordinates[0][0]*element_coordinates[2][1] +\
+ J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
+ double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
+ element_coordinates[1][0]*element_coordinates[0][1] -\
+ J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
+
+ // Map coordinates to the reference square
+ if (std::abs(y - 1.0) < 1e-14)
+ x = -1.0;
+ else
+ x = 2.0 *x/(1.0 - y) - 1.0;
+ y = 2.0*y - 1.0;
+
+ // Reset values
+ *values = 0;
+
+ // Map degree of freedom to element degree of freedom
+ const unsigned int dof = i;
+
+ // Generate scalings
+ const double scalings_y_0 = 1;
+ const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
+ const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
+ const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
+
+ // Compute psitilde_a
+ const double psitilde_a_0 = 1;
+ const double psitilde_a_1 = x;
+ const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
+ const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
+
+ // Compute psitilde_bs
+ const double psitilde_bs_0_0 = 1;
+ const double psitilde_bs_0_1 = 1.5*y + 0.5;
+ const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;
+ const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
+ const double psitilde_bs_1_0 = 1;
+ const double psitilde_bs_1_1 = 2.5*y + 1.5;
+ const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
+ const double psitilde_bs_2_0 = 1;
+ const double psitilde_bs_2_1 = 3.5*y + 2.5;
+ const double psitilde_bs_3_0 = 1;
+
+ // Compute basisvalues
+ const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
+ const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
+ const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
+ const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
+ const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
+ const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
+ const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
+ const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
+ const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
+ const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
+
+ // Table(s) of coefficients
+ static const double coefficients0[10][10] = \
+ {{0.047140452079103, -0.0288675134594813, -0.0166666666666666, 0.0782460796435951, 0.0606091526731326, 0.0349927106111883, -0.0601337794302955, -0.0508223195384204, -0.0393667994375868, -0.0227284322524248},
+ {0.0471404520791031, 0.0288675134594813, -0.0166666666666666, 0.0782460796435952, -0.0606091526731326, 0.0349927106111883, 0.0601337794302955, -0.0508223195384204, 0.0393667994375868, -0.0227284322524248},
+ {0.0471404520791032, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},
+ {0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},
+ {0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},
+ {0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},
+ {0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},
+ {0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524247},
+ {0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384204, 0.0131222664791956, -0.0227284322524247},
+ {0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153682, 0, 0.090913729009699}};
+
+ // Extract relevant coefficients
+ const double coeff0_0 = coefficients0[dof][0];
+ const double coeff0_1 = coefficients0[dof][1];
+ const double coeff0_2 = coefficients0[dof][2];
+ const double coeff0_3 = coefficients0[dof][3];
+ const double coeff0_4 = coefficients0[dof][4];
+ const double coeff0_5 = coefficients0[dof][5];
+ const double coeff0_6 = coefficients0[dof][6];
+ const double coeff0_7 = coefficients0[dof][7];
+ const double coeff0_8 = coefficients0[dof][8];
+ const double coeff0_9 = coefficients0[dof][9];
+
+ // Compute value(s)
+ *values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5 + coeff0_6*basisvalue6 + coeff0_7*basisvalue7 + coeff0_8*basisvalue8 + coeff0_9*basisvalue9;
+ }
+
+ /// Evaluate all basis functions at given point in cell
+ virtual void evaluate_basis_all(double* values,
+ const double* coordinates,
+ const ufc::cell& c) const
+ {
+ throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
+ }
+
+ /// Evaluate order n derivatives of basis function i at given point in cell
+ virtual void evaluate_basis_derivatives(unsigned int i,
+ unsigned int n,
+ double* values,
+ const double* coordinates,
+ const ufc::cell& c) const
+ {
+ // Extract vertex coordinates
+ const double * const * element_coordinates = c.coordinates;
+
+ // Compute Jacobian of affine map from reference cell
+ const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
+ const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
+ const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
+ const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
+
+ // Compute determinant of Jacobian
+ const double detJ = J_00*J_11 - J_01*J_10;
+
+ // Compute inverse of Jacobian
+
+ // Get coordinates and map to the reference (UFC) element
+ double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
+ element_coordinates[0][0]*element_coordinates[2][1] +\
+ J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
+ double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
+ element_coordinates[1][0]*element_coordinates[0][1] -\
+ J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
+
+ // Map coordinates to the reference square
+ if (std::abs(y - 1.0) < 1e-14)
+ x = -1.0;
+ else
+ x = 2.0 *x/(1.0 - y) - 1.0;
+ y = 2.0*y - 1.0;
+
+ // Compute number of derivatives
+ unsigned int num_derivatives = 1;
+
+ for (unsigned int j = 0; j < n; j++)
+ num_derivatives *= 2;
+
+
+ // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
+ unsigned int **combinations = new unsigned int *[num_derivatives];
+
+ for (unsigned int j = 0; j < num_derivatives; j++)
+ {
+ combinations[j] = new unsigned int [n];
+ for (unsigned int k = 0; k < n; k++)
+ combinations[j][k] = 0;
+ }
+
+ // Generate combinations of derivatives
+ for (unsigned int row = 1; row < num_derivatives; row++)
+ {
+ for (unsigned int num = 0; num < row; num++)
+ {
+ for (unsigned int col = n-1; col+1 > 0; col--)
+ {
+ if (combinations[row][col] + 1 > 1)
+ combinations[row][col] = 0;
+ else
+ {
+ combinations[row][col] += 1;
+ break;
+ }
+ }
+ }
+ }
+
+ // Compute inverse of Jacobian
+ const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
+
+ // Declare transformation matrix
+ // Declare pointer to two dimensional array and initialise
+ double **transform = new double *[num_derivatives];
+
+ for (unsigned int j = 0; j < num_derivatives; j++)
+ {
+ transform[j] = new double [num_derivatives];
+ for (unsigned int k = 0; k < num_derivatives; k++)
+ transform[j][k] = 1;
+ }
+
+ // Construct transformation matrix
+ for (unsigned int row = 0; row < num_derivatives; row++)
+ {
+ for (unsigned int col = 0; col < num_derivatives; col++)
+ {
+ for (unsigned int k = 0; k < n; k++)
+ transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
+ }
+ }
+
+ // Reset values
+ for (unsigned int j = 0; j < 1*num_derivatives; j++)
+ values[j] = 0;
+
+ // Map degree of freedom to element degree of freedom
+ const unsigned int dof = i;
+
+ // Generate scalings
+ const double scalings_y_0 = 1;
+ const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);
+ const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);
+ const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);
+
+ // Compute psitilde_a
+ const double psitilde_a_0 = 1;
+ const double psitilde_a_1 = x;
+ const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;
+ const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;
+
+ // Compute psitilde_bs
+ const double psitilde_bs_0_0 = 1;
+ const double psitilde_bs_0_1 = 1.5*y + 0.5;
+ const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;
+ const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;
+ const double psitilde_bs_1_0 = 1;
+ const double psitilde_bs_1_1 = 2.5*y + 1.5;
+ const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;
+ const double psitilde_bs_2_0 = 1;
+ const double psitilde_bs_2_1 = 3.5*y + 2.5;
+ const double psitilde_bs_3_0 = 1;
+
+ // Compute basisvalues
+ const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
+ const double basisvalue1 = 1.73205080756888*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;
+ const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;
+ const double basisvalue3 = 2.73861278752583*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;
+ const double basisvalue4 = 2.12132034355964*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;
+ const double basisvalue5 = 1.22474487139159*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;
+ const double basisvalue6 = 3.74165738677394*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;
+ const double basisvalue7 = 3.16227766016838*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;
+ const double basisvalue8 = 2.44948974278318*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;
+ const double basisvalue9 = 1.4142135623731*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;
+
+ // Table(s) of coefficients
+ static const double coefficients0[10][10] = \
+ {{0.047140452079103, -0.0288675134594813, -0.0166666666666666, 0.0782460796435951, 0.0606091526731326, 0.0349927106111883, -0.0601337794302955, -0.0508223195384204, -0.0393667994375868, -0.0227284322524248},
+ {0.0471404520791031, 0.0288675134594813, -0.0166666666666666, 0.0782460796435952, -0.0606091526731326, 0.0349927106111883, 0.0601337794302955, -0.0508223195384204, 0.0393667994375868, -0.0227284322524248},
+ {0.0471404520791032, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},
+ {0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},
+ {0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},
+ {0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},
+ {0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},
+ {0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524247},
+ {0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384204, 0.0131222664791956, -0.0227284322524247},
+ {0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153682, 0, 0.090913729009699}};
+
+ // Interesting (new) part
+ // Tables of derivatives of the polynomial base (transpose)
+ static const double dmats0[10][10] = \
+ {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {4.89897948556636, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {0, 9.48683298050514, 0, 0, 0, 0, 0, 0, 0, 0},
+ {4, 0, 7.07106781186547, 0, 0, 0, 0, 0, 0, 0},
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {5.29150262212918, 0, -2.99332590941915, 13.6626010212795, 0, 0.61101009266078, 0, 0, 0, 0},
+ {0, 4.38178046004133, 0, 0, 12.5219806739988, 0, 0, 0, 0, 0},
+ {3.46410161513775, 0, 7.83836717690617, 0, 0, 8.4, 0, 0, 0, 0},
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
+
+ static const double dmats1[10][10] = \
+ {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {2.44948974278318, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {4.24264068711929, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {2.58198889747161, 4.74341649025257, -0.912870929175274, 0, 0, 0, 0, 0, 0, 0},
+ {2, 6.12372435695795, 3.53553390593274, 0, 0, 0, 0, 0, 0, 0},
+ {-2.30940107675851, 0, 8.16496580927726, 0, 0, 0, 0, 0, 0, 0},
+ {2.64575131106459, 5.18459255872629, -1.49666295470957, 6.83130051063973, -1.05830052442584, 0.305505046330391, 0, 0, 0, 0},
+ {2.23606797749979, 2.19089023002067, 2.5298221281347, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0, 0, 0, 0},
+ {1.73205080756887, -5.09116882454314, 3.91918358845308, 0, 9.69948452238571, 4.2, 0, 0, 0, 0},
+ {5.00000000000001, 0, -2.82842712474619, 0, 0, 12.1243556529821, 0, 0, 0, 0}};
+
+ // Compute reference derivatives
+ // Declare pointer to array of derivatives on FIAT element
+ double *derivatives = new double [num_derivatives];
+
+ // Declare coefficients
+ double coeff0_0 = 0;
+ double coeff0_1 = 0;
+ double coeff0_2 = 0;
+ double coeff0_3 = 0;
+ double coeff0_4 = 0;
+ double coeff0_5 = 0;
+ double coeff0_6 = 0;
+ double coeff0_7 = 0;
+ double coeff0_8 = 0;
+ double coeff0_9 = 0;
+
+ // Declare new coefficients
+ double new_coeff0_0 = 0;
+ double new_coeff0_1 = 0;
+ double new_coeff0_2 = 0;
+ double new_coeff0_3 = 0;
+ double new_coeff0_4 = 0;
+ double new_coeff0_5 = 0;
+ double new_coeff0_6 = 0;
+ double new_coeff0_7 = 0;
+ double new_coeff0_8 = 0;
+ double new_coeff0_9 = 0;
+
+ // Loop possible derivatives
+ for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
+ {
+ // Get values from coefficients array
+ new_coeff0_0 = coefficients0[dof][0];
+ new_coeff0_1 = coefficients0[dof][1];
+ new_coeff0_2 = coefficients0[dof][2];
+ new_coeff0_3 = coefficients0[dof][3];
+ new_coeff0_4 = coefficients0[dof][4];
+ new_coeff0_5 = coefficients0[dof][5];
+ new_coeff0_6 = coefficients0[dof][6];
+ new_coeff0_7 = coefficients0[dof][7];
+ new_coeff0_8 = coefficients0[dof][8];
+ new_coeff0_9 = coefficients0[dof][9];
+
+ // Loop derivative order
+ for (unsigned int j = 0; j < n; j++)
+ {
+ // Update old coefficients
+ coeff0_0 = new_coeff0_0;
+ coeff0_1 = new_coeff0_1;
+ coeff0_2 = new_coeff0_2;
+ coeff0_3 = new_coeff0_3;
+ coeff0_4 = new_coeff0_4;
+ coeff0_5 = new_coeff0_5;
+ coeff0_6 = new_coeff0_6;
+ coeff0_7 = new_coeff0_7;
+ coeff0_8 = new_coeff0_8;
+ coeff0_9 = new_coeff0_9;
+
+ if(combinations[deriv_num][j] == 0)
+ {
+ new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0] + coeff0_6*dmats0[6][0] + coeff0_7*dmats0[7][0] + coeff0_8*dmats0[8][0] + coeff0_9*dmats0[9][0];
+ new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1] + coeff0_6*dmats0[6][1] + coeff0_7*dmats0[7][1] + coeff0_8*dmats0[8][1] + coeff0_9*dmats0[9][1];
+ new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2] + coeff0_6*dmats0[6][2] + coeff0_7*dmats0[7][2] + coeff0_8*dmats0[8][2] + coeff0_9*dmats0[9][2];
+ new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3] + coeff0_6*dmats0[6][3] + coeff0_7*dmats0[7][3] + coeff0_8*dmats0[8][3] + coeff0_9*dmats0[9][3];
+ new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4] + coeff0_6*dmats0[6][4] + coeff0_7*dmats0[7][4] + coeff0_8*dmats0[8][4] + coeff0_9*dmats0[9][4];
+ new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5] + coeff0_6*dmats0[6][5] + coeff0_7*dmats0[7][5] + coeff0_8*dmats0[8][5] + coeff0_9*dmats0[9][5];
+ new_coeff0_6 = coeff0_0*dmats0[0][6] + coeff0_1*dmats0[1][6] + coeff0_2*dmats0[2][6] + coeff0_3*dmats0[3][6] + coeff0_4*dmats0[4][6] + coeff0_5*dmats0[5][6] + coeff0_6*dmats0[6][6] + coeff0_7*dmats0[7][6] + coeff0_8*dmats0[8][6] + coeff0_9*dmats0[9][6];
+ new_coeff0_7 = coeff0_0*dmats0[0][7] + coeff0_1*dmats0[1][7] + coeff0_2*dmats0[2][7] + coeff0_3*dmats0[3][7] + coeff0_4*dmats0[4][7] + coeff0_5*dmats0[5][7] + coeff0_6*dmats0[6][7] + coeff0_7*dmats0[7][7] + coeff0_8*dmats0[8][7] + coeff0_9*dmats0[9][7];
+ new_coeff0_8 = coeff0_0*dmats0[0][8] + coeff0_1*dmats0[1][8] + coeff0_2*dmats0[2][8] + coeff0_3*dmats0[3][8] + coeff0_4*dmats0[4][8] + coeff0_5*dmats0[5][8] + coeff0_6*dmats0[6][8] + coeff0_7*dmats0[7][8] + coeff0_8*dmats0[8][8] + coeff0_9*dmats0[9][8];
+ new_coeff0_9 = coeff0_0*dmats0[0][9] + coeff0_1*dmats0[1][9] + coeff0_2*dmats0[2][9] + coeff0_3*dmats0[3][9] + coeff0_4*dmats0[4][9] + coeff0_5*dmats0[5][9] + coeff0_6*dmats0[6][9] + coeff0_7*dmats0[7][9] + coeff0_8*dmats0[8][9] + coeff0_9*dmats0[9][9];
+ }
+ if(combinations[deriv_num][j] == 1)
+ {
+ new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0] + coeff0_6*dmats1[6][0] + coeff0_7*dmats1[7][0] + coeff0_8*dmats1[8][0] + coeff0_9*dmats1[9][0];
+ new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1] + coeff0_6*dmats1[6][1] + coeff0_7*dmats1[7][1] + coeff0_8*dmats1[8][1] + coeff0_9*dmats1[9][1];
+ new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2] + coeff0_6*dmats1[6][2] + coeff0_7*dmats1[7][2] + coeff0_8*dmats1[8][2] + coeff0_9*dmats1[9][2];
+ new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3] + coeff0_6*dmats1[6][3] + coeff0_7*dmats1[7][3] + coeff0_8*dmats1[8][3] + coeff0_9*dmats1[9][3];
+ new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4] + coeff0_6*dmats1[6][4] + coeff0_7*dmats1[7][4] + coeff0_8*dmats1[8][4] + coeff0_9*dmats1[9][4];
+ new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5] + coeff0_6*dmats1[6][5] + coeff0_7*dmats1[7][5] + coeff0_8*dmats1[8][5] + coeff0_9*dmats1[9][5];
+ new_coeff0_6 = coeff0_0*dmats1[0][6] + coeff0_1*dmats1[1][6] + coeff0_2*dmats1[2][6] + coeff0_3*dmats1[3][6] + coeff0_4*dmats1[4][6] + coeff0_5*dmats1[5][6] + coeff0_6*dmats1[6][6] + coeff0_7*dmats1[7][6] + coeff0_8*dmats1[8][6] + coeff0_9*dmats1[9][6];
+ new_coeff0_7 = coeff0_0*dmats1[0][7] + coeff0_1*dmats1[1][7] + coeff0_2*dmats1[2][7] + coeff0_3*dmats1[3][7] + coeff0_4*dmats1[4][7] + coeff0_5*dmats1[5][7] + coeff0_6*dmats1[6][7] + coeff0_7*dmats1[7][7] + coeff0_8*dmats1[8][7] + coeff0_9*dmats1[9][7];
+ new_coeff0_8 = coeff0_0*dmats1[0][8] + coeff0_1*dmats1[1][8] + coeff0_2*dmats1[2][8] + coeff0_3*dmats1[3][8] + coeff0_4*dmats1[4][8] + coeff0_5*dmats1[5][8] + coeff0_6*dmats1[6][8] + coeff0_7*dmats1[7][8] + coeff0_8*dmats1[8][8] + coeff0_9*dmats1[9][8];
+ new_coeff0_9 = coeff0_0*dmats1[0][9] + coeff0_1*dmats1[1][9] + coeff0_2*dmats1[2][9] + coeff0_3*dmats1[3][9] + coeff0_4*dmats1[4][9] + coeff0_5*dmats1[5][9] + coeff0_6*dmats1[6][9] + coeff0_7*dmats1[7][9] + coeff0_8*dmats1[8][9] + coeff0_9*dmats1[9][9];
+ }
+
+ }
+ // Compute derivatives on reference element as dot product of coefficients and basisvalues
+ derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5 + new_coeff0_6*basisvalue6 + new_coeff0_7*basisvalue7 + new_coeff0_8*basisvalue8 + new_coeff0_9*basisvalue9;
+ }
+
+ // Transform derivatives back to physical element
+ for (unsigned int row = 0; row < num_derivatives; row++)
+ {
+ for (unsigned int col = 0; col < num_derivatives; col++)
+ {
+ values[row] += transform[row][col]*derivatives[col];
+ }
+ }
+ // Delete pointer to array of derivatives on FIAT element
+ delete [] derivatives;
+
+ // Delete pointer to array of combinations of derivatives and transform
+ for (unsigned int row = 0; row < num_derivatives; row++)
+ {
+ delete [] combinations[row];
+ delete [] transform[row];
+ }
+
+ delete [] combinations;
+ delete [] transform;
+ }
+
+ /// Evaluate order n derivatives of all basis functions at given point in cell
+ virtual void evaluate_basis_derivatives_all(unsigned int n,
+ double* values,
+ const double* coordinates,
+ const ufc::cell& c) const
+ {
+ throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
+ }
+
+ /// Evaluate linear functional for dof i on the function f
+ virtual double evaluate_dof(unsigned int i,
+ const ufc::function& f,
+ const ufc::cell& c) const
+ {
+ // The reference points, direction and weights:
+ static const double X[10][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.666666666666667, 0.333333333333333}}, {{0.333333333333333, 0.666666666666667}}, {{0, 0.333333333333333}}, {{0, 0.666666666666667}}, {{0.333333333333333, 0}}, {{0.666666666666667, 0}}, {{0.333333333333333, 0.333333333333333}}};
+ static const double W[10][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};
+ static const double D[10][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};
+
+ const double * const * x = c.coordinates;
+ double result = 0.0;
+ // Iterate over the points:
+ // Evaluate basis functions for affine mapping
+ const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
+ const double w1 = X[i][0][0];
+ const double w2 = X[i][0][1];
+
+ // Compute affine mapping y = F(X)
+ double y[2];
+ y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
+ y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
+
+ // Evaluate function at physical points
+ double values[1];
+ f.evaluate(values, y, c);
+
+ // Map function values using appropriate mapping
+ // Affine map: Do nothing
+
+ // Note that we do not map the weights (yet).
+
+ // Take directional components
+ for(int k = 0; k < 1; k++)
+ result += values[k]*D[i][0][k];
+ // Multiply by weights
+ result *= W[i][0];
+
+ return result;
+ }
+
+ /// Evaluate linear functionals for all dofs on the function f
+ virtual void evaluate_dofs(double* values,
+ const ufc::function& f,
+ const ufc::cell& c) const
+ {
+ throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
+ }
+
+ /// Interpolate vertex values from dof values
+ virtual void interpolate_vertex_values(double* vertex_values,
+ const double* dof_values,
+ const ufc::cell& c) const
+ {
+ // Evaluate at vertices and use affine mapping
+ vertex_values[0] = dof_values[0];
+ vertex_values[1] = dof_values[1];
+ vertex_values[2] = dof_values[2];
+ }
+
+ /// Return the number of sub elements (for a mixed element)
+ virtual unsigned int num_sub_elements() const
+ {
+ return 1;
+ }
+
+ /// Create a new finite element for sub element i (for a mixed element)
+ virtual ufc::finite_element* create_sub_element(unsigned int i) const
+ {
+ return new p3_0_finite_element_0();
+ }
+
+};
+
+/// This class defines the interface for a local-to-global mapping of
+/// degrees of freedom (dofs).
+
+class p3_0_dof_map_0: public ufc::dof_map
+{
+private:
+
+ unsigned int __global_dimension;
+
+public:
+
+ /// Constructor
+ p3_0_dof_map_0() : ufc::dof_map()
+ {
+ __global_dimension = 0;
+ }
+
+ /// Destructor
+ virtual ~p3_0_dof_map_0()
+ {
+ // Do nothing
+ }
+
+ /// Return a string identifying the dof map
+ virtual const char* signature() const
+ {
+ return "FFC dof map for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 3)";
+ }
+
+ /// Return true iff mesh entities of topological dimension d are needed
+ virtual bool needs_mesh_entities(unsigned int d) const
+ {
+ switch ( d )
+ {
+ case 0:
+ return true;
+ break;
+ case 1:
+ return true;
+ break;
+ case 2:
+ return true;
+ break;
+ }
+ return false;
+ }
+
+ /// Initialize dof map for mesh (return true iff init_cell() is needed)
+ virtual bool init_mesh(const ufc::mesh& m)
+ {
+ __global_dimension = m.num_entities[0] + 2*m.num_entities[1] + m.num_entities[2];
+ return false;
+ }
+
+ /// Initialize dof map for given cell
+ virtual void init_cell(const ufc::mesh& m,
+ const ufc::cell& c)
+ {
+ // Do nothing
+ }
+
+ /// Finish initialization of dof map for cells
+ virtual void init_cell_finalize()
+ {
+ // Do nothing
+ }
+
+ /// Return the dimension of the global finite element function space
+ virtual unsigned int global_dimension() const
+ {
+ return __global_dimension;
+ }
+
+ /// Return the dimension of the local finite element function space for a cell
+ virtual unsigned int local_dimension(const ufc::cell& c) const
+ {
+ return 10;
+ }
+
+ /// Return the maximum dimension of the local finite element function space
+ virtual unsigned int max_local_dimension() const
+ {
+ return 10;
+ }
+
+ // Return the geometric dimension of the coordinates this dof map provides
+ virtual unsigned int geometric_dimension() const
+ {
+ return 2;
+ }
+
+ /// Return the number of dofs on each cell facet
+ virtual unsigned int num_facet_dofs() const
+ {
+ return 4;
+ }
+
+ /// Return the number of dofs associated with each cell entity of dimension d
+ virtual unsigned int num_entity_dofs(unsigned int d) const
+ {
+ throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
+ }
+
+ /// Tabulate the local-to-global mapping of dofs on a cell
+ virtual void tabulate_dofs(unsigned int* dofs,
+ const ufc::mesh& m,
+ const ufc::cell& c) const
+ {
+ dofs[0] = c.entity_indices[0][0];
+ dofs[1] = c.entity_indices[0][1];
+ dofs[2] = c.entity_indices[0][2];
+ unsigned int offset = m.num_entities[0];
+ dofs[3] = offset + 2*c.entity_indices[1][0];
+ dofs[4] = offset + 2*c.entity_indices[1][0] + 1;
+ dofs[5] = offset + 2*c.entity_indices[1][1];
+ dofs[6] = offset + 2*c.entity_indices[1][1] + 1;
+ dofs[7] = offset + 2*c.entity_indices[1][2];
+ dofs[8] = offset + 2*c.entity_indices[1][2] + 1;
+ offset = offset + 2*m.num_entities[1];
+ dofs[9] = offset + c.entity_indices[2][0];
+ }
+
+ /// Tabulate the local-to-local mapping from facet dofs to cell dofs
+ virtual void tabulate_facet_dofs(unsigned int* dofs,
+ unsigned int facet) const
+ {
+ switch ( facet )
+ {
+ case 0:
+ dofs[0] = 1;
+ dofs[1] = 2;
+ dofs[2] = 3;
+ dofs[3] = 4;
+ break;
+ case 1:
+ dofs[0] = 0;
+ dofs[1] = 2;
+ dofs[2] = 5;
+ dofs[3] = 6;
+ break;
+ case 2:
+ dofs[0] = 0;
+ dofs[1] = 1;
+ dofs[2] = 7;
+ dofs[3] = 8;
+ break;
+ }
+ }
+
+ /// Tabulate the local-to-local mapping of dofs on entity (d, i)
+ virtual void tabulate_entity_dofs(unsigned int* dofs,
+ unsigned int d, unsigned int i) const
+ {
+ throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
+ }
+
+ /// Tabulate the coordinates of all dofs on a cell
+ virtual void tabulate_coordinates(double** coordinates,
+ const ufc::cell& c) const
+ {
+ const double * const * x = c.coordinates;
+ coordinates[0][0] = x[0][0];
+ coordinates[0][1] = x[0][1];
+ coordinates[1][0] = x[1][0];
+ coordinates[1][1] = x[1][1];
+ coordinates[2][0] = x[2][0];
+ coordinates[2][1] = x[2][1];
+ coordinates[3][0] = 0.666666666666667*x[1][0] + 0.333333333333333*x[2][0];
+ coordinates[3][1] = 0.666666666666667*x[1][1] + 0.333333333333333*x[2][1];
+ coordinates[4][0] = 0.333333333333333*x[1][0] + 0.666666666666667*x[2][0];
+ coordinates[4][1] = 0.333333333333333*x[1][1] + 0.666666666666667*x[2][1];
+ coordinates[5][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[2][0];
+ coordinates[5][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[2][1];
+ coordinates[6][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[2][0];
+ coordinates[6][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[2][1];
+ coordinates[7][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[1][0];
+ coordinates[7][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[1][1];
+ coordinates[8][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[1][0];
+ coordinates[8][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[1][1];
+ coordinates[9][0] = 0.333333333333333*x[0][0] + 0.333333333333333*x[1][0] + 0.333333333333333*x[2][0];
+ coordinates[9][1] = 0.333333333333333*x[0][1] + 0.333333333333333*x[1][1] + 0.333333333333333*x[2][1];
+ }
+
+ /// Return the number of sub dof maps (for a mixed element)
+ virtual unsigned int num_sub_dof_maps() const
+ {
+ return 1;
+ }
+
+ /// Create a new dof_map for sub dof map i (for a mixed element)
+ virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
+ {
+ return new p3_0_dof_map_0();
+ }
+
+};
+
+/// This class defines the interface for the assembly of the global
+/// tensor corresponding to a form with r + n arguments, that is, a
+/// mapping
+///
+/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
+///
+/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
+/// global tensor A is defined by
+///
+/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
+///
+/// where each argument Vj represents the application to the
+/// sequence of basis functions of Vj and w1, w2, ..., wn are given
+/// fixed functions (coefficients).
+
+class p3_form_0: public ufc::form
+{
+public:
+
+ /// Constructor
+ p3_form_0() : ufc::form()
+ {
+ // Do nothing
+ }
+
+ /// Destructor
+ virtual ~p3_form_0()
+ {
+ // Do nothing
+ }
+
+ /// Return a string identifying the form
+ virtual const char* signature() const
+ {
+ return "None";
+ }
+
+ /// Return the rank of the global tensor (r)
+ virtual unsigned int rank() const
+ {
+ return -1;
+ }
+
+ /// Return the number of coefficients (n)
+ virtual unsigned int num_coefficients() const
+ {
+ return 0;
+ }
+
+ /// Return the number of cell integrals
+ virtual unsigned int num_cell_integrals() const
+ {
+ return 0;
+ }
+
+ /// Return the number of exterior facet integrals
+ virtual unsigned int num_exterior_facet_integrals() const
+ {
+ return 0;
+ }
+
+ /// Return the number of interior facet integrals
+ virtual unsigned int num_interior_facet_integrals() const
+ {
+ return 0;
+ }
+
+ /// Create a new finite element for argument function i
+ virtual ufc::finite_element* create_finite_element(unsigned int i) const
+ {
+ return 0;
+ }
+
+ /// Create a new dof map for argument function i
+ virtual ufc::dof_map* create_dof_map(unsigned int i) const
+ {
+ return 0;
+ }
+
+ /// Create a new cell integral on sub domain i
+ virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
+ {
+ return 0;
+ }
+
+ /// Create a new exterior facet integral on sub domain i
+ virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
+ {
+ return 0;
+ }
+
+ /// Create a new interior facet integral on sub domain i
+ virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
+ {
+ return 0;
+ }
+
+};
+
+// DOLFIN wrappers
+
+// Standard library includes
+#include <string>
+
+// DOLFIN includes
+#include <dolfin/common/NoDeleter.h>
+#include <dolfin/fem/FiniteElement.h>
+#include <dolfin/fem/DofMap.h>
+#include <dolfin/fem/Form.h>
+#include <dolfin/function/FunctionSpace.h>
+#include <dolfin/function/GenericFunction.h>
+#include <dolfin/function/CoefficientAssigner.h>
+
+namespace P3
+{
+
+class Form_0_FunctionSpace_0: public dolfin::FunctionSpace
+{
+public:
+
+ Form_0_FunctionSpace_0(const dolfin::Mesh& mesh):
+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new p3_0_finite_element_0()))),
+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new p3_0_dof_map_0()), mesh)))
+ {
+ // Do nothing
+ }
+
+ Form_0_FunctionSpace_0(dolfin::Mesh& mesh):
+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new p3_0_finite_element_0()))),
+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new p3_0_dof_map_0()), mesh)))
+ {
+ // Do nothing
+ }
+
+ Form_0_FunctionSpace_0(boost::shared_ptr<dolfin::Mesh> mesh):
+ dolfin::FunctionSpace(mesh,
+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new p3_0_finite_element_0()))),
+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new p3_0_dof_map_0()), *mesh)))
+ {
+ // Do nothing
+ }
+
+ Form_0_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh):
+ dolfin::FunctionSpace(mesh,
+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new p3_0_finite_element_0()))),
+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new p3_0_dof_map_0()), *mesh)))
+ {
+ // Do nothing
+ }
+
+
+ ~Form_0_FunctionSpace_0()
+ {
+ }
+
+};
+
+class Form_0: public dolfin::Form
+{
+public:
+
+ // Constructor
+ Form_0(const dolfin::FunctionSpace& V0):
+ dolfin::Form(1, 0)
+ {
+ _function_spaces[0] = reference_to_no_delete_pointer(V0);
+
+ _ufc_form = boost::shared_ptr<const ufc::form>(new p3_form_0());
+ }
+
+ // Constructor
+ Form_0(boost::shared_ptr<const dolfin::FunctionSpace> V0):
+ dolfin::Form(1, 0)
+ {
+ _function_spaces[0] = V0;
+
+ _ufc_form = boost::shared_ptr<const ufc::form>(new p3_form_0());
+ }
+
+ // Destructor
+ ~Form_0()
+ {}
+
+ /// Return the number of the coefficient with this name
+ virtual dolfin::uint coefficient_number(const std::string& name) const
+ {
+
+ dolfin::error("No coefficients.");
+ return 0;
+ }
+
+ /// Return the name of the coefficient with this number
+ virtual std::string coefficient_name(dolfin::uint i) const
+ {
+
+ dolfin::error("No coefficients.");
+ return "unnamed";
+ }
+
+ // Typedefs
+ typedef Form_0_FunctionSpace_0 TestSpace;
+
+ // Coefficients
+};
+
+// Class typedefs
+typedef Form_0 LinearForm;
+typedef Form_0::TestSpace FunctionSpace;
+
+} // namespace P3
+
+#endif
=== added file 'sandbox/cgal/test/P3.ufl'
--- sandbox/cgal/test/P3.ufl 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/test/P3.ufl 2009-11-30 01:20:26 +0000
@@ -0,0 +1,13 @@
+# Copyright (C) 2009 Garth N. Wells.
+# Licensed under the GNU LGPL Version 2.1.
+#
+# First added: 2009-06-18
+# Last changed:
+#
+# The bilinear form a(v, u) and linear form L(v) for
+# projection onto piecewise quadratics.
+#
+# Compile this form with FFC: ffc -l dolfin P3.ufl
+
+element = FiniteElement("Lagrange", "triangle", 3)
+
=== added file 'sandbox/cgal/test/test_function_eval.cpp'
--- sandbox/cgal/test/test_function_eval.cpp 1970-01-01 00:00:00 +0000
+++ sandbox/cgal/test/test_function_eval.cpp 2009-11-30 01:20:26 +0000
@@ -0,0 +1,79 @@
+// Copyright (C) 2009 Andre Massing
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added: 2009-11-07
+// Last changed: 2009-11-10
+//
+// Added to generate wiht GTS a list of reference points and function values
+// to check function evaluation with CGAL against it.
+
+#include <dolfin.h>
+#include "P1.h"
+#include "P3.h"
+
+using namespace dolfin;
+using dolfin::uint;
+
+#include "P3.h"
+
+#ifdef HAS_GTS
+
+class F : public Expression
+{
+public:
+
+ F() : Expression(3) {}
+
+ void eval(double* values, const double* x) const
+ {
+ values[0] = sin(3.0*x[0])*sin(3.0*x[1])*sin(3.0*x[2]);
+ }
+
+};
+
+#elif HAS_GCAL
+
+class F : public Expression
+{
+public:
+
+ F() : Expression(3) {}
+
+ void eval(double* values, const double* x) const
+ {
+ values[0] = sin(3.0*x[0])*sin(3.0*x[1])*sin(3.0*x[2]);
+ }
+
+};
+
+#endif
+
+int main()
+{
+ not_working_in_parallel("This demo");
+ UnitCube mesh(20, 20, 20);
+
+ // Create function spaces
+ P3::FunctionSpace V0(mesh);
+
+ // Create functions
+ Function f0(V0);
+
+ F f;
+ f0.interpolate(f);
+
+ double value = 0.0;
+ double X[3] = {0.0, 0.0, 0.0};
+ srand(1);
+ uint num_points = 100000;
+ for (uint i = 1; i <= num_points; ++i)
+ {
+ X[0] = std::rand()/static_cast<double>(RAND_MAX);
+ X[1] = std::rand()/static_cast<double>(RAND_MAX);
+ X[2] = std::rand()/static_cast<double>(RAND_MAX);
+ f.eval(&value, X);
+ info("x = %.12e\ty = %.12e\tz = %.12e\tf(x) = %.12e",X[0],X[1],X[2],value);
+ }
+
+return 0;
+}
=== modified file 'site-packages/dolfin/__init__.py'
--- site-packages/dolfin/__init__.py 2009-11-11 19:59:27 +0000
+++ site-packages/dolfin/__init__.py 2009-11-30 01:20:26 +0000
@@ -66,7 +66,7 @@
from cpp import SubDomain, DomainBoundary
from cpp import UnitInterval, Interval, UnitSquare, Rectangle, UnitCircle, UnitCube, Box, UnitSphere
#from cpp import DofMap, DofMapSet
-from cpp import IntersectionDetector
+from cpp import IntersectionOperator
#from cpp import BoundaryCondition
#from cpp import DirichletBC, PeriodicBC
#from cpp import Form, Assembler
@@ -92,7 +92,7 @@
from cpp import lagrange, hermite, harmonic
# Linear algebra backend related imports
-from cpp import (has_la_backend, has_mpi, has_slepc, has_scotch, has_gts,
+from cpp import (has_la_backend, has_mpi, has_slepc, has_scotch, has_gts, has_cgal,
has_umfpack, has_cholmod, has_parmetis, has_gmp, has_zlib)
from cpp import DefaultFactory
=== modified file 'test/unit/mesh/python/test.py'
--- test/unit/mesh/python/test.py 2009-10-08 19:11:44 +0000
+++ test/unit/mesh/python/test.py 2009-11-30 01:20:26 +0000
@@ -200,5 +200,19 @@
mesh = UnitSquare(5, 5)
self.assertEqual(len(mesh.cells()), 50)
+class IntersectionOperator(unittest.TestCase):
+ def testIntersectPoints(self):
+ pass
+
+ def testIntersectPoints(self):
+ pass
+
+ def testIntersectMesh2D(self):
+ pass
+
+ def testIntersectMesh3D(self):
+ pass
+
+
if __name__ == "__main__":
unittest.main()
