Commit be3f0e91 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

remove required dependency on dune-functions and dune-localfunctions

parent ca0597ef
...@@ -4,8 +4,7 @@ ...@@ -4,8 +4,7 @@
#Name of the module #Name of the module
Module: dune-curvedsurfacegrid Module: dune-curvedsurfacegrid
Version: 0.1.0 Version: 2.8
Maintainer: florian.stenger@tu-dresden.de Maintainer: florian.stenger@tu-dresden.de
#depending on Depends: dune-common dune-geometry dune-grid dune-curvedgeometry
Depends: dune-common dune-geometry dune-grid dune-localfunctions dune-curvedgeometry dune-functions Suggests: dune-localfunctions dune-functions dune-foamgrid
Suggests: dune-alugrid dune-foamgrid
...@@ -280,6 +280,12 @@ public: ...@@ -280,6 +280,12 @@ public:
* Internally, the local function is used to construct this parametrization * Internally, the local function is used to construct this parametrization
* either implicitly, by first interpolating the local function into a local basis * either implicitly, by first interpolating the local function into a local basis
* or explicitly as reinterpreting the local function directly as geometry. * or explicitly as reinterpreting the local function directly as geometry.
*
* NOTE: The geometry expects a LocalGeometry parametrization for the entity in the
* grid element. Since this specialization is for grid elements, this parametrization
* mapping is just the identity. This is implemented by using the DefaultLocalGeometry
* that simply forwards the input arguments of all the functions or returns an identity
* object.
*/ */
Geometry geometry () const Geometry geometry () const
{ {
......
...@@ -4,17 +4,18 @@ ...@@ -4,17 +4,18 @@
#include <cmath> #include <cmath>
#include <dune/curvedsurfacegrid/gridfunctions/analyticgridfunction.hh> #include <dune/curvedsurfacegrid/gridfunctions/analyticgridfunction.hh>
#include <dune/curvedsurfacegrid/utility/kdtree.hh>
#include <dune/curvedsurfacegrid/utility/vertexcache.hh>
#include "kdtree.hh" #ifndef COLLECTOR_BENCHMARK
#include "vertexcache.hh" #define COLLECTOR_BENCHMARK 1 // $benchmark
#endif
#define COLLECTOR_BENCHMARK 1 /// \brief $benchmark
#if COLLECTOR_BENCHMARK != 0 #if COLLECTOR_BENCHMARK != 0
#include <dune/common/timer.hh> #include <dune/common/timer.hh>
#endif #endif
namespace Dune namespace Dune {
{
/// \brief Closest-point projection to a surface given by a dune-grid /// \brief Closest-point projection to a surface given by a dune-grid
/** /**
...@@ -29,14 +30,20 @@ class ExplicitSurfaceProjection ...@@ -29,14 +30,20 @@ class ExplicitSurfaceProjection
using Self = ExplicitSurfaceProjection<T>; using Self = ExplicitSurfaceProjection<T>;
using Domain = FieldVector<T, 3>; using Domain = FieldVector<T, 3>;
struct SurfaceElement{ struct SurfaceElement
{
unsigned v0, v1, v2; unsigned v0, v1, v2;
SurfaceElement(unsigned v0, unsigned v1, unsigned v2) : v0(v0), v1(v1), v2(v2){} SurfaceElement (unsigned v0, unsigned v1, unsigned v2)
: v0(v0), v1(v1), v2(v2)
{}
}; };
//calculates the distance between a given point and the line segment between 'lin0' and 'lin1'. //calculates the distance between a given point and the line segment between 'lin0' and 'lin1'.
static T distance_point_line_3d(const Domain &point, const Domain &lin0, const Domain &lin1, static T distance_point_line_3d (const Domain& point, const Domain& lin0, const Domain& lin1,
Domain &hit_vertex){ Domain &hit_vertex)
{
using std::sqrt;
//the (complete) line through 'lin0' and 'lin1' is given by 'l = lin0 + t*g' with g as follows: //the (complete) line through 'lin0' and 'lin1' is given by 'l = lin0 + t*g' with g as follows:
T gx = lin1[0] - lin0[0]; T gx = lin1[0] - lin0[0];
T gy = lin1[1] - lin0[1]; T gy = lin1[1] - lin0[1];
...@@ -49,17 +56,17 @@ class ExplicitSurfaceProjection ...@@ -49,17 +56,17 @@ class ExplicitSurfaceProjection
//if the orthogonal projection of 'point' onto l is not on the given line segment then one of the //if the orthogonal projection of 'point' onto l is not on the given line segment then one of the
//vertices 'lin0' or 'lin1' is the nearest point to 'point' //vertices 'lin0' or 'lin1' is the nearest point to 'point'
if(t < 0){ if (t < 0) {
//'lin0' is nearest //'lin0' is nearest
hit_vertex[0] = lin0[0]; hit_vertex[0] = lin0[0];
hit_vertex[1] = lin0[1]; hit_vertex[1] = lin0[1];
hit_vertex[2] = lin0[2]; hit_vertex[2] = lin0[2];
}else if(t > 1){ } else if(t > 1) {
//'lin1' is nearest //'lin1' is nearest
hit_vertex[0] = lin1[0]; hit_vertex[0] = lin1[0];
hit_vertex[1] = lin1[1]; hit_vertex[1] = lin1[1];
hit_vertex[2] = lin1[2]; hit_vertex[2] = lin1[2];
}else{ } else {
//orthogonal projection is nearest //orthogonal projection is nearest
hit_vertex[0] = lin0[0] + t * gx; hit_vertex[0] = lin0[0] + t * gx;
hit_vertex[1] = lin0[1] + t * gy; hit_vertex[1] = lin0[1] + t * gy;
...@@ -75,9 +82,12 @@ class ExplicitSurfaceProjection ...@@ -75,9 +82,12 @@ class ExplicitSurfaceProjection
} }
//calculates the distance between a given 'point' and the triangle given by 'tri0', 'tri1', 'tri2'. //calculates the distance between a given 'point' and the triangle given by 'tri0', 'tri1', 'tri2'.
static T distance_point_triangle_3d(const Domain &point, static T distance_point_triangle_3d (const Domain& point,
const Domain &tri0, const Domain &tri1, const Domain &tri2, Domain &hit_vertex){ const Domain& tri0, const Domain& tri1, const Domain& tri2, Domain& hit_vertex)
const T epsilon = std::sqrt(std::numeric_limits<T>::epsilon()); {
using std::abs;
using std::sqrt;
const T epsilon = sqrt(std::numeric_limits<T>::epsilon());
//triangle lies in plane 'p = tri0 + u*g + v*h', where g and h are defined as follows: //triangle lies in plane 'p = tri0 + u*g + v*h', where g and h are defined as follows:
T gx = tri1[0] - tri0[0]; T gx = tri1[0] - tri0[0];
...@@ -101,17 +111,17 @@ class ExplicitSurfaceProjection ...@@ -101,17 +111,17 @@ class ExplicitSurfaceProjection
//solve the system: //solve the system:
T u, v; T u, v;
if(fabs(b1) < epsilon){ if (abs(b1) < epsilon) {
v = c2 / b2; v = c2 / b2;
u = c1 / a1; u = c1 / a1;
}else{ } else {
v = (c1 * b1 - c2 * a1) / (b1 * b1 - b2 * a1); v = (c1 * b1 - c2 * a1) / (b1 * b1 - b2 * a1);
u = (c1 - v * b1) / a1; u = (c1 - v * b1) / a1;
} }
//now test whether the orthogonal projection of 'point' onto p lies inside the triangle //now test whether the orthogonal projection of 'point' onto p lies inside the triangle
T distance; T distance;
if(u >= 0 && v >= 0 && u + v <= 1){ if (u >= 0 && v >= 0 && u + v <= 1) {
//yes, inside. The length of the normal-vector is the desired distance //yes, inside. The length of the normal-vector is the desired distance
T fx = hit_vertex[0] = tri0[0] + u * gx + v * hx; T fx = hit_vertex[0] = tri0[0] + u * gx + v * hx;
T fy = hit_vertex[1] = tri0[1] + u * gy + v * hy; T fy = hit_vertex[1] = tri0[1] + u * gy + v * hy;
...@@ -120,17 +130,17 @@ class ExplicitSurfaceProjection ...@@ -120,17 +130,17 @@ class ExplicitSurfaceProjection
fy -= point[1]; fy -= point[1];
fz -= point[2]; fz -= point[2];
distance = sqrt(fx * fx + fy * fy + fz * fz); distance = sqrt(fx * fx + fy * fy + fz * fz);
}else{ } else {
//no, not inside. The desired distance is the distance to one of the triangle-edges //no, not inside. The desired distance is the distance to one of the triangle-edges
distance = distance_point_line_3d(point, tri0, tri1, hit_vertex); distance = distance_point_line_3d(point, tri0, tri1, hit_vertex);
Domain current_hit_vertex; Domain current_hit_vertex;
T current_distance = distance_point_line_3d(point, tri1, tri2, current_hit_vertex); T current_distance = distance_point_line_3d(point, tri1, tri2, current_hit_vertex);
if(current_distance < distance){ if (current_distance < distance) {
distance = current_distance; distance = current_distance;
hit_vertex = current_hit_vertex; hit_vertex = current_hit_vertex;
} }
current_distance = distance_point_line_3d(point, tri2, tri0, current_hit_vertex); current_distance = distance_point_line_3d(point, tri2, tri0, current_hit_vertex);
if(current_distance < distance){ if (current_distance < distance) {
distance = current_distance; distance = current_distance;
hit_vertex = current_hit_vertex; hit_vertex = current_hit_vertex;
} }
...@@ -141,7 +151,7 @@ class ExplicitSurfaceProjection ...@@ -141,7 +151,7 @@ class ExplicitSurfaceProjection
public: public:
//constructor using a pre-loaded Dune-surface-grid //constructor using a pre-loaded Dune-surface-grid
template <class Grid> template <class Grid>
ExplicitSurfaceProjection(const Grid& grid, bool cached=true) ExplicitSurfaceProjection (const Grid& grid, bool cached = true)
: vertices_(grid.size(Grid::dimension)) : vertices_(grid.size(Grid::dimension))
, adjacent_elements_(grid.size(Grid::dimension)) , adjacent_elements_(grid.size(Grid::dimension))
, minCorner_(std::numeric_limits<T>::max()) , minCorner_(std::numeric_limits<T>::max())
...@@ -180,13 +190,13 @@ public: ...@@ -180,13 +190,13 @@ public:
kdtree_.emplace(vertices_, minCorner_, maxCorner_); kdtree_.emplace(vertices_, minCorner_, maxCorner_);
kdtree_->update(); kdtree_->update();
if(cached_){ if (cached_) {
Domain c_expansion = (maxCorner_ - minCorner_) * 0.1; Domain c_expansion = (maxCorner_ - minCorner_) * 0.1;
cache_.emplace(minCorner_ - c_expansion, maxCorner_ + c_expansion); cache_.emplace(minCorner_ - c_expansion, maxCorner_ + c_expansion);
} }
} }
ExplicitSurfaceProjection(Self const& other) ExplicitSurfaceProjection (const Self& other)
: vertices_(other.vertices_) : vertices_(other.vertices_)
, adjacent_elements_(other.adjacent_elements_) , adjacent_elements_(other.adjacent_elements_)
, elements_(other.elements_) , elements_(other.elements_)
...@@ -203,7 +213,7 @@ public: ...@@ -203,7 +213,7 @@ public:
kdtree_->update(); kdtree_->update();
} }
ExplicitSurfaceProjection(Self&& other) ExplicitSurfaceProjection (Self&& other)
: vertices_(std::move(other.vertices_)) : vertices_(std::move(other.vertices_))
, adjacent_elements_(std::move(other.adjacent_elements_)) , adjacent_elements_(std::move(other.adjacent_elements_))
, elements_(std::move(other.elements_)) , elements_(std::move(other.elements_))
...@@ -223,15 +233,15 @@ public: ...@@ -223,15 +233,15 @@ public:
kdtree_->update(); /// \brief $crash: slow but works kdtree_->update(); /// \brief $crash: slow but works
} }
Self& operator=(Self const&) = delete; Self& operator= (const Self&) = delete;
Self& operator=(Self&&) = delete; Self& operator= (Self&&) = delete;
void resetCache() void resetCache ()
{ {
if(cached_) cache_->clear(); if(cached_) cache_->clear();
} }
Domain closestPoint(const Domain& x, std::size_t closest_vertex) const Domain closestPoint (const Domain& x, std::size_t closest_vertex) const
{ {
T min_dist = std::numeric_limits<T>::max(); T min_dist = std::numeric_limits<T>::max();
Domain hit_vertex, best_hit_vertex; Domain hit_vertex, best_hit_vertex;
...@@ -248,14 +258,14 @@ public: ...@@ -248,14 +258,14 @@ public:
} }
//evaluation of the closest-point projection //evaluation of the closest-point projection
Domain operator()(const Domain& x) const Domain operator() (const Domain& x) const
{ {
#if COLLECTOR_BENCHMARK != 0 #if COLLECTOR_BENCHMARK != 0
inputs_.push_back(x); inputs_.push_back(x);
#endif #endif
Domain result; Domain result;
if(!cached_ || !cache_->find(x, result)){ if (!cached_ || !cache_->find(x, result)){
//find closest surface-vertex //find closest surface-vertex
std::vector<std::size_t> outIndices(1); std::vector<std::size_t> outIndices(1);
std::vector<T> outDistSq(1); std::vector<T> outDistSq(1);
...@@ -265,20 +275,21 @@ public: ...@@ -265,20 +275,21 @@ public:
result = closestPoint(x, outIndices[0]); result = closestPoint(x, outIndices[0]);
//cache result //cache result
if(cached_) cache_->insert(x, result); if (cached_) cache_->insert(x, result);
} }
return result; return result;
} }
#if COLLECTOR_BENCHMARK != 0 #if COLLECTOR_BENCHMARK != 0
void collectorBenchmark(bool fresh_cache=true) const{ void collectorBenchmark (bool fresh_cache=true) const
if(cached_ && fresh_cache) cache_->clear(); {
if (cached_ && fresh_cache) cache_->clear();
int uncached_cnt = 0; int uncached_cnt = 0;
Timer t; Timer t;
for(auto &input : inputs_){ for (auto &input : inputs_){
Domain result; Domain result;
if(!cached_ || !cache_->find(input, result)){ if (!cached_ || !cache_->find(input, result)) {
//find closest surface-vertex //find closest surface-vertex
std::vector<std::size_t> outIndices(1); std::vector<std::size_t> outIndices(1);
std::vector<T> outDistSq(1); std::vector<T> outDistSq(1);
...@@ -288,7 +299,7 @@ public: ...@@ -288,7 +299,7 @@ public:
result = closestPoint(input, outIndices[0]); result = closestPoint(input, outIndices[0]);
//cache result //cache result
if(cached_) cache_->insert(input, result); if (cached_) cache_->insert(input, result);
++uncached_cnt; ++uncached_cnt;
} }
input = result; input = result;
...@@ -296,7 +307,7 @@ public: ...@@ -296,7 +307,7 @@ public:
std::cout << "surfGF time: " << t.elapsed() << ", count: "<< inputs_.size() std::cout << "surfGF time: " << t.elapsed() << ", count: "<< inputs_.size()
<< ", uncached: " << uncached_cnt << std::endl; << ", uncached: " << uncached_cnt << std::endl;
} }
#endif #endif
private: private:
std::vector<Domain> vertices_; std::vector<Domain> vertices_;
...@@ -304,17 +315,17 @@ private: ...@@ -304,17 +315,17 @@ private:
std::vector<SurfaceElement> elements_; std::vector<SurfaceElement> elements_;
Domain minCorner_; Domain minCorner_;
Domain maxCorner_; Domain maxCorner_;
Std::optional<KDTree<T>> kdtree_; std::optional<KDTree<T>> kdtree_;
bool cached_; bool cached_;
mutable Std::optional<VertexCache<Domain, T>> cache_; mutable std::optional<VertexCache<Domain, T>> cache_;
#if COLLECTOR_BENCHMARK != 0 #if COLLECTOR_BENCHMARK != 0
mutable std::vector<Domain> inputs_; mutable std::vector<Domain> inputs_;
#endif #endif
}; };
/// \brief construct a grid function representing an explicit surface parametrization /// \brief construct a grid function representing an explicit surface parametrization
template <class Grid, class T > template <class Grid, class T>
auto explicitSurfaceGridFunction (Grid &grid, bool cached=true) auto explicitSurfaceGridFunction (Grid& grid, bool cached = true)
{ {
return analyticGridFunction<Grid>(ExplicitSurfaceProjection<T>{grid, cached}); return analyticGridFunction<Grid>(ExplicitSurfaceProjection<T>{grid, cached});
} }
......
...@@ -5,6 +5,10 @@ ...@@ -5,6 +5,10 @@
#include <type_traits> #include <type_traits>
#include <utility> #include <utility>
#if !HAVE_DUNE_LOCALFUNCTIONS
#error "Need dune-localfunctions for the definition of the AnalyticDiscreteFunction"
#endif
#include <dune/common/typeutilities.hh> #include <dune/common/typeutilities.hh>
#include <dune/curvedsurfacegrid/gridfunctions/discretegridviewfunction.hh> #include <dune/curvedsurfacegrid/gridfunctions/discretegridviewfunction.hh>
#include <dune/curvedsurfacegrid/gridfunctions/gridentityset.hh> #include <dune/curvedsurfacegrid/gridfunctions/gridentityset.hh>
......
if (dune-foamgrid_FOUND) foreach (geometry RANGE 1 3)
foreach (geometry RANGE 1 3) dune_add_test(NAME convergence-${geometry}
dune_add_test(NAME convergence-${geometry} SOURCES convergence.cc
SOURCES convergence.cc
COMPILE_DEFINITIONS
GEOMETRY_TYPE=${geometry}
DUNE_GRID_PATH=\"${PROJECT_SOURCE_DIR}/doc/grids/\"
CMAKE_GUARD
dune-foamgrid_FOUND)
endforeach (geometry)
dune_add_test(SOURCES discretegridviewfunction.cc
COMPILE_DEFINITIONS COMPILE_DEFINITIONS
GEOMETRY_TYPE=${geometry}
DUNE_GRID_PATH=\"${PROJECT_SOURCE_DIR}/doc/grids/\" DUNE_GRID_PATH=\"${PROJECT_SOURCE_DIR}/doc/grids/\"
USE_FLOAT128=0
CMAKE_GUARD CMAKE_GUARD
dune-functions_FOUND
dune-foamgrid_FOUND) dune-foamgrid_FOUND)
endif () endforeach (geometry)
dune_add_test(SOURCES discretegridviewfunction.cc
COMPILE_DEFINITIONS
DUNE_GRID_PATH=\"${PROJECT_SOURCE_DIR}/doc/grids/\"
USE_FLOAT128=0
CMAKE_GUARD
dune-functions_FOUND
dune-foamgrid_FOUND)
...@@ -10,7 +10,6 @@ ...@@ -10,7 +10,6 @@
#include <dune/curvedsurfacegrid/geometries/ellipsoid.hh> #include <dune/curvedsurfacegrid/geometries/ellipsoid.hh>
#include <dune/curvedsurfacegrid/geometries/sphere.hh> #include <dune/curvedsurfacegrid/geometries/sphere.hh>
#include <dune/curvedsurfacegrid/geometries/torus.hh> #include <dune/curvedsurfacegrid/geometries/torus.hh>
#include <dune/curvedsurfacegrid/gridfunctions/discretegridviewfunction.hh>
#include <dune/foamgrid/foamgrid.hh> #include <dune/foamgrid/foamgrid.hh>
#include <dune/geometry/quadraturerules.hh> #include <dune/geometry/quadraturerules.hh>
#include <dune/grid/io/file/gmshreader.hh> #include <dune/grid/io/file/gmshreader.hh>
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment