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

merged remote and cleanup

parents 8394ba2d cd2a5c1d
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_CRVSRF_BACKUPRESTORE_HH
#define DUNE_CRVSRF_BACKUPRESTORE_HH
#ifndef DUNE_CURVED_SURFACE_GRID_BACKUPRESTORE_HH
#define DUNE_CURVED_SURFACE_GRID_BACKUPRESTORE_HH
#include <dune/grid/common/backuprestore.hh>
......@@ -9,43 +9,42 @@
namespace Dune
{
namespace crvsrf
namespace CGeo
{
// BackupRestoreFacilities
// -----------------------
template< class Grid, bool hasBackupRestoreFacilities = Capabilities::hasBackupRestoreFacilities< Grid > ::v >
template< class Grid, bool hasBackupRestoreFacilities = Capabilities::hasBackupRestoreFacilities<Grid> ::v >
class BackupRestoreFacilities
{};
template< class Grid >
class BackupRestoreFacilities< Grid, true >
class BackupRestoreFacilities<Grid, true>
{
typedef BackupRestoreFacilities< Grid, true > This;
using Self = BackupRestoreFacilities<Grid, true>;
protected:
BackupRestoreFacilities ()
{}
private:
BackupRestoreFacilities ( const This & );
This &operator= ( const This & );
BackupRestoreFacilities (const Self&);
Self& operator= (const Self&);
protected:
const Grid &asImp () const
const Grid& asImp () const
{
return static_cast< const Grid & >( *this );
return static_cast<const Grid&>(*this);
}
Grid &asImp ()
Grid& asImp ()
{
return static_cast< Grid & >( *this );
return static_cast<Grid&>(*this);
}
};
} // namespace crvsrf
} // namespace CGeo
......@@ -53,40 +52,39 @@ namespace Dune
// -------------------------------------------
template< class HostGrid, class CoordFunction, int order, class Allocator >
struct BackupRestoreFacility< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > >
struct BackupRestoreFacility<CurvedSurfaceGrid<HostGrid, CoordFunction, order, Allocator>>
{
typedef CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > Grid;
typedef BackupRestoreFacility< HostGrid > HostBackupRestoreFacility;
using Grid = CurvedSurfaceGrid<HostGrid, CoordFunction, order, Allocator>;
using HostBackupRestoreFacility = BackupRestoreFacility<HostGrid>;
static void backup ( const Grid &grid, const std::string &filename )
static void backup (const Grid& grid, const std::string& filename)
{
// notice: We should also backup the coordinate function
HostBackupRestoreFacility::backup( grid.hostGrid(), filename );
}
static void backup ( const Grid &grid, const std::ostream &stream )
static void backup (const Grid& grid, const std::ostream& stream)
{
// notice: We should also backup the coordinate function
HostBackupRestoreFacility::backup( grid.hostGrid(), stream );
HostBackupRestoreFacility::backup(grid.hostGrid(), stream);
}
static Grid *restore ( const std::string &filename )
static Grid* restore (const std::string& filename)
{
// notice: We should also restore the coordinate function
HostGrid *hostGrid = HostBackupRestoreFacility::restore( filename );
CoordFunction *coordFunction = new CoordFunction();
return new Grid( hostGrid, coordFunction );
HostGrid* hostGrid = HostBackupRestoreFacility::restore(filename);
CoordFunction* coordFunction = new CoordFunction();
return new Grid(hostGrid, coordFunction);
}
static Grid *restore ( const std::istream &stream )
static Grid* restore (const std::istream& stream)
{
// notice: We should also restore the coordinate function
HostGrid *hostGrid = HostBackupRestoreFacility::restore( stream );
CoordFunction *coordFunction = new CoordFunction();
return new Grid( hostGrid, coordFunction );
HostGrid* hostGrid = HostBackupRestoreFacility::restore(stream);
CoordFunction* coordFunction = new CoordFunction();
return new Grid(hostGrid, coordFunction);
}
};
} // namespace Dune
#endif // #ifndef DUNE_CRVSRF_BACKUPRESTORE_HH
#endif // DUNE_CURVED_SURFACE_GRID_BACKUPRESTORE_HH
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_CRVSRF_CAPABILITIES_HH
#define DUNE_CRVSRF_CAPABILITIES_HH
#ifndef DUNE_CURVED_SURFACE_GRID_CAPABILITIES_HH
#define DUNE_CURVED_SURFACE_GRID_CAPABILITIES_HH
#include <cassert>
......@@ -24,61 +24,61 @@ namespace Dune
// Capabilities from dune-grid
// ---------------------------
template< class HostGrid, class CoordFunction, int order, class Allocator >
struct hasSingleGeometryType< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > >
template< class HostGrid, class CoordFunction, int order, bool geoCaching, class Allocator >
struct hasSingleGeometryType< CurvedSurfaceGrid< HostGrid, CoordFunction, order, geoCaching, Allocator > >
{
static const bool v = hasSingleGeometryType< HostGrid > :: v;
static const unsigned int topologyId = hasSingleGeometryType< HostGrid > :: topologyId;
};
template< class HostGrid, class CoordFunction, int order, class Allocator, int codim >
struct hasEntity< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator >, codim >
template< class HostGrid, class CoordFunction, int order, bool geoCaching, class Allocator, int codim >
struct hasEntity< CurvedSurfaceGrid< HostGrid, CoordFunction, order, geoCaching, Allocator >, codim >
{
static const bool v = true;
};
template< class HostGrid, class CoordFunction, int order, class Allocator, int codim >
struct hasEntityIterator< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator >, codim >
template< class HostGrid, class CoordFunction, int order, bool geoCaching, class Allocator, int codim >
struct hasEntityIterator< CurvedSurfaceGrid< HostGrid, CoordFunction, order, geoCaching, Allocator >, codim >
{
static const bool v = true;
};
template< class HostGrid, class CoordFunction, int order, class Allocator, int codim >
struct canCommunicate< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator >, codim >
template< class HostGrid, class CoordFunction, int order, bool geoCaching, class Allocator, int codim >
struct canCommunicate< CurvedSurfaceGrid< HostGrid, CoordFunction, order, geoCaching, Allocator >, codim >
{
static const bool v = canCommunicate< HostGrid, codim >::v && hasEntity< HostGrid, codim >::v;
};
template< class HostGrid, class CoordFunction, int order, class Allocator >
struct hasBackupRestoreFacilities< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > >
template< class HostGrid, class CoordFunction, int order, bool geoCaching, class Allocator >
struct hasBackupRestoreFacilities< CurvedSurfaceGrid< HostGrid, CoordFunction, order, geoCaching, Allocator > >
{
static const bool v = hasBackupRestoreFacilities< HostGrid >::v;
};
template< class HostGrid, class CoordFunction, int order, class Allocator >
struct isLevelwiseConforming< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > >
template< class HostGrid, class CoordFunction, int order, bool geoCaching, class Allocator >
struct isLevelwiseConforming< CurvedSurfaceGrid< HostGrid, CoordFunction, order, geoCaching, Allocator > >
{
static const bool v = isLevelwiseConforming< HostGrid >::v;
};
template< class HostGrid, class CoordFunction, int order, class Allocator >
struct isLeafwiseConforming< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > >
template< class HostGrid, class CoordFunction, int order, bool geoCaching, class Allocator >
struct isLeafwiseConforming< CurvedSurfaceGrid< HostGrid, CoordFunction, order, geoCaching, Allocator > >
{
static const bool v = isLeafwiseConforming< HostGrid >::v;
};
template< class HostGrid, class CoordFunction, int order, class Allocator >
struct threadSafe< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > >
template< class HostGrid, class CoordFunction, int order, bool geoCaching, class Allocator >
struct threadSafe< CurvedSurfaceGrid< HostGrid, CoordFunction, order, geoCaching, Allocator > >
{
static const bool v = false;
};
template< class HostGrid, class CoordFunction, int order, class Allocator >
struct viewThreadSafe< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > >
template< class HostGrid, class CoordFunction, int order, bool geoCaching, class Allocator >
struct viewThreadSafe< CurvedSurfaceGrid< HostGrid, CoordFunction, order, geoCaching, Allocator > >
{
static const bool v = false;
};
......@@ -88,8 +88,8 @@ namespace Dune
// hasHostEntity
// -------------
template< class HostGrid, class CoordFunction, int order, class Allocator, int codim >
struct hasHostEntity< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator >, codim >
template< class HostGrid, class CoordFunction, int order, bool geoCaching, class Allocator, int codim >
struct hasHostEntity< CurvedSurfaceGrid< HostGrid, CoordFunction, order, geoCaching, Allocator >, codim >
{
static const bool v = hasEntity< HostGrid, codim >::v;
};
......@@ -98,4 +98,4 @@ namespace Dune
} // namespace Dune
#endif // #ifndef DUNE_CRVSRF_CAPABILITIES_HH
#endif // DUNE_CURVED_SURFACE_GRID_CAPABILITIES_HH
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_CRVSRF_COORDFUNCTIONIMPLEMENTATIONS_HH
#define DUNE_CRVSRF_COORDFUNCTIONIMPLEMENTATIONS_HH
#ifndef DUNE_CURVED_SURFACE_GRID_COORDFUNCTIONIMPLEMENTATIONS_HH
#define DUNE_CURVED_SURFACE_GRID_COORDFUNCTIONIMPLEMENTATIONS_HH
#include <dune/curvedsurfacegrid/surfacedistance/surfacedistance.hh>
#include <dune/curvedsurfacegrid/surfacedistance/vertexmap.hh>
#include <dune/curvedsurfacegrid/surfacedistance/vtureader.hh>
#include <dune/grid/geometrygrid/coordfunction.hh>
#include <cfloat>
#define COLLECTOR_BENCHMARK 0
......@@ -27,6 +28,48 @@ double dbl_time(){
namespace Dune
{
//convertDuneGridToCrvsrfMesh
template<class GridPtr>
CGeo::Mesh *convertDuneGridToCrvsrfMesh(GridPtr &grid){
auto gridView = grid->leafGridView();
//determine coordinate-extrema
CGeo::Vertex minc(DBL_MAX, DBL_MAX, DBL_MAX), maxc(-DBL_MAX, -DBL_MAX, -DBL_MAX);
for(const auto &vertexEntity : vertices(gridView)){
const auto &vertex = vertexEntity.geometry().corner(0);
for(int i = 0; i < grid->dimensionworld; ++i){
if(vertex[i] < minc[i]) minc[i] = vertex[i];
if(vertex[i] > maxc[i]) maxc[i] = vertex[i];
}
}
CGeo::Mesh *surface_mesh = new CGeo::Mesh(grid->dimension, grid->dimensionworld);
CGeo::VertexMap<int> vmap(minc, maxc, grid->dimensionworld);
//build mesh
int next_index = 0;
for(const auto &element : elements(gridView)){
if(element.geometry().corners() != 3){
CGeo::error("Element is not a triangle in convertDuneGridToCrvsrfMesh.");
}
int32_t el[grid->dimension + 1];
for(int i = 0; i <= grid->dimension; ++i){
auto coords = element.geometry().corner(i);
CGeo::Vertex v(coords[0], coords[1], coords[2]);
if(!vmap.find_3d(v, el[i])){
surface_mesh->add_vertex(v);
el[i] = next_index;
vmap.insert(v, next_index);
++next_index;
}
}
surface_mesh->add_element(el);
}
return surface_mesh;
}
//CenterSphereCoordFunction
struct CenterSphereCoordFunction
: public AnalyticalCoordFunction<double, 3, 3, CenterSphereCoordFunction>
......@@ -44,13 +87,13 @@ namespace Dune
}
#if COLLECTOR_BENCHMARK != 0
void collector_benchmark(bool fresh_cache){
void collectorBenchmark(bool fresh_cache){
double start_time = dbl_time();
for(auto &input : inputs){
input *= radius_ / input.two_norm();
}
double run_time = dbl_time() - start_time;
std::cout << "csphere time: " << crvsrf::ft(run_time, 11, 9)
std::cout << "csphere time: " << CGeo::ft(run_time, 11, 9)
<< ", count: " << inputs.size() << std::endl;
}
#endif
......@@ -84,7 +127,7 @@ namespace Dune
}
#if COLLECTOR_BENCHMARK != 0
void collector_benchmark(bool fresh_cache){
void collectorBenchmark(bool fresh_cache){
double start_time = dbl_time();
for(auto &input : inputs){
FieldVector<double, 3> direction = input - center_;
......@@ -92,7 +135,7 @@ namespace Dune
for(int i = 0; i < 3; ++i) input[i] = center_[i] + direction[i] * factor;
}
double run_time = dbl_time() - start_time;
std::cout << "sphere time: " << crvsrf::ft(run_time, 11, 9)
std::cout << "sphere time: " << CGeo::ft(run_time, 11, 9)
<< ", count: " << inputs.size() << std::endl;
}
#endif
......@@ -108,96 +151,162 @@ namespace Dune
//SurfaceDistanceCoordFunction
template<bool cached=true>
struct SurfaceDistanceCoordFunction
: public AnalyticalCoordFunction<double, 3, 3, SurfaceDistanceCoordFunction>
: public AnalyticalCoordFunction<double, 3, 3, SurfaceDistanceCoordFunction<cached> >
{
const bool coordCaching = cached;
//constructor using a pre-loaded Dune-surface-grid
template<class GridPtr>
SurfaceDistanceCoordFunction(const GridPtr &gridPtr,
double expansion_width=-1.0, bool no_scaling=false) : cache(nullptr){
if(gridPtr->dimension != 2 || gridPtr->dimensionworld != 3){
CGeo::error("SurfaceDistanceCoordFunction needs a triangle-surface-mesh!");
}
surface_mesh = convertDuneGridToCrvsrfMesh(gridPtr);
backgrid = new CGeo::BackGrid(*surface_mesh, expansion_width, no_scaling);
if(cached){
CGeo::Vertex c_begin, c_end;
surface_mesh->determine_coordinate_extrema(c_begin, c_end);
CGeo::Vertex c_expansion = (c_end - c_begin) * 0.1;
cache = new CGeo::VertexMap<CGeo::Vertex>(c_begin - c_expansion, c_end + c_expansion, 3);
}
}
//constructor using vtu-reader to read "file"
SurfaceDistanceCoordFunction(const std::string &file, double expansion_width=-1.0, bool no_scaling=false){
surface_mesh = new crvsrf::Mesh(2, 3);
SurfaceDistanceCoordFunction(const char *file,
double expansion_width=-1.0, bool no_scaling=false) : cache(nullptr){
surface_mesh = new CGeo::Mesh(2, 3);
read_vtu_file(file, *surface_mesh);
backgrid = new crvsrf::BackGrid(*surface_mesh, expansion_width, no_scaling);
crvsrf::Vertex c_begin, c_end;
surface_mesh->determine_coordinate_extrema(c_begin, c_end);
crvsrf::Vertex c_expansion = (c_end - c_begin) * 0.1;
cache = new crvsrf::VertexMap<crvsrf::Vertex>(c_begin - c_expansion, c_end + c_expansion, 3);
backgrid = new CGeo::BackGrid(*surface_mesh, expansion_width, no_scaling);
if(cached){
CGeo::Vertex c_begin, c_end;
surface_mesh->determine_coordinate_extrema(c_begin, c_end);
CGeo::Vertex c_expansion = (c_end - c_begin) * 0.1;
cache = new CGeo::VertexMap<CGeo::Vertex>(c_begin - c_expansion, c_end + c_expansion, 3);
}
}
//constructor using a pre-loaded Dune-surface-mesh
/*!SurfaceDistanceCoordFunction(Grid &mesh, double expansion_width=-1.0, bool no_scaling=false){
surface_mesh = convert_dune_mesh_to_crvsrf::mesh...
backgrid = new BackGrid(*surface_mesh, expansion_width, no_scaling);
//! implement cache-init
}*/
//constructor using vtu-reader to read "file"
SurfaceDistanceCoordFunction(const std::string &file,
double expansion_width=-1.0, bool no_scaling=false)
: SurfaceDistanceCoordFunction(file.c_str(), expansion_width, no_scaling)
{}
//destructor
~SurfaceDistanceCoordFunction(){
delete cache;
if(cached) delete cache;
delete backgrid;
delete surface_mesh;
}
//operator()
FieldVector<double, 3> operator()(const FieldVector<double, 3> &x) const{
#if COLLECTOR_BENCHMARK != 0
inputs.push_back(x);
#endif
crvsrf::Vertex hit_vertex;
crvsrf::Vertex v(x[0], x[1], x[2]);
//uncached
//backgrid->distance_to_surface_3d(v, hit_vertex);
//cached
if(!cache->find_3d(v, hit_vertex)){
backgrid->distance_to_surface_3d(v, hit_vertex);
cache->insert(v, hit_vertex);
//resetCache
void resetCache(CGeo::Vertex begin = CGeo::Vertex(),
CGeo::Vertex end = CGeo::Vertex(),
double expansion_factor = 0.1){
if(cached){
delete cache;
if(begin.is_origin() && end.is_origin()){
surface_mesh->determine_coordinate_extrema(begin, end);
}
CGeo::Vertex c_expansion = (end - begin) * expansion_factor;
cache = new CGeo::VertexMap<CGeo::Vertex>(begin - c_expansion, end + c_expansion, 3);
}else{
CGeo::error("Cannot reset cache in uncached SurfaceDistanceCoordFunction!");
}
return {hit_vertex[0], hit_vertex[1], hit_vertex[2]};
}
//operator()
FieldVector<double, 3> operator()(const FieldVector<double, 3> &x) const;
//output_backgrid_dimensions
void output_backgrid_dimensions(){
backgrid->debug_output_dimensions();
}
#if COLLECTOR_BENCHMARK != 0
void collector_benchmark(bool fresh_cache){
if(fresh_cache) cache->clear();
int uncached_cnt = 0;
double start_time = dbl_time();
for(auto &input : inputs){
crvsrf::Vertex hit_vertex;
crvsrf::Vertex v(input[0], input[1], input[2]);
//uncached
/*backgrid->distance_to_surface_3d(v, hit_vertex);
++uncached_cnt;*/
//cached
if(!cache->find_3d(v, hit_vertex)){
void collectorBenchmark(bool fresh_cache){
if(cached){
if(fresh_cache) cache->clear();
int uncached_cnt = 0;
double start_time = dbl_time();
for(auto &input : inputs){
CGeo::Vertex hit_vertex;
CGeo::Vertex v(input[0], input[1], input[2]);
if(!cache->find_3d(v, hit_vertex)){
backgrid->distance_to_surface_3d(v, hit_vertex);
cache->insert(v, hit_vertex);
++uncached_cnt;
}
}
double run_time = dbl_time() - start_time;
std::cout << "surfdist time: " << CGeo::ft(run_time, 11, 9) << ", count: "<< inputs.size()
<< ", uncached: " << uncached_cnt << std::endl;
}else{ //uncached
double start_time = dbl_time();
for(auto &input : inputs){
CGeo::Vertex hit_vertex;
CGeo::Vertex v(input[0], input[1], input[2]);
backgrid->distance_to_surface_3d(v, hit_vertex);
cache->insert(v, hit_vertex);
++uncached_cnt;
}
double run_time = dbl_time() - start_time;
std::cout << "surfdist time: " << CGeo::ft(run_time, 11, 9)
<< ", count: "<< inputs.size() << std::endl;
}
}
#endif
input = {hit_vertex[0], hit_vertex[1], hit_vertex[2]};
#if VERTEXMAP_DEBUG_MODE != 0
void output_vertexmap_statistics(){
if(cached){
cache->debug_output_statistics();
}else{
CGeo::error("Cannot output vertexmap-statistics in "
"uncached SurfaceDistanceCoordFunction!");
}
double run_time = dbl_time() - start_time;
std::cout << "surfdist time: " << crvsrf::ft(run_time, 11, 9) << ", count: "<< inputs.size()
<< ", uncached: " << uncached_cnt << std::endl;
}
#endif
private:
crvsrf::Mesh *surface_mesh;
crvsrf::BackGrid *backgrid;
crvsrf::VertexMap<crvsrf::Vertex> *cache;
CGeo::Mesh *surface_mesh;
CGeo::BackGrid *backgrid;
CGeo::VertexMap<CGeo::Vertex> *cache;
#if COLLECTOR_BENCHMARK != 0
mutable std::vector<FieldVector<double, 3> > inputs;
#endif
};
//operator() (cached)
template<>
inline FieldVector<double, 3> SurfaceDistanceCoordFunction<true>::operator()(
const FieldVector<double, 3> &x) const{
#if COLLECTOR_BENCHMARK != 0
inputs.push_back(x);
#endif
CGeo::Vertex hit_vertex;
CGeo::Vertex v(x[0], x[1], x[2]);
if(!cache->find_3d(v, hit_vertex)){
backgrid->distance_to_surface_3d(v, hit_vertex);
cache->insert(v, hit_vertex);
}
return {hit_vertex[0], hit_vertex[1], hit_vertex[2]};
}
//operator() (uncached)
template<>
inline FieldVector<double, 3> SurfaceDistanceCoordFunction<false>::operator()(
const FieldVector<double, 3> &x) const{
#if COLLECTOR_BENCHMARK != 0
inputs.push_back(x);
#endif
CGeo::Vertex hit_vertex;
CGeo::Vertex v(x[0], x[1], x[2]);
backgrid->distance_to_surface_3d(v, hit_vertex);
return {hit_vertex[0], hit_vertex[1], hit_vertex[2]};
}
} // namespace Dune
#endif // #ifndef DUNE_CRVSRF_COORDFUNCTIONIMPLEMENTATIONS_HH
#endif // DUNE_CURVED_SURFACE_GRID_COORDFUNCTIONIMPLEMENTATIONS_HH
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_CRVSRF_COORDPROVIDER_HH
#define DUNE_CRVSRF_COORDPROVIDER_HH
#include <array>
#include <dune/localfunctions/lagrange/equidistantpoints.hh>
#include <dune/localfunctions/lagrange/lagrangelfecache.hh>
namespace Dune
{
namespace crvsrf
{
template< class Coordinate>
Coordinate edgeCenter(Coordinate gc1, const Coordinate &gc2)
{
for(typename Coordinate::size_type i = 0, end = gc1.size(); i < end; ++i)
{
gc1[i] = (gc1[i] + gc2[i]) / 2.0;
}
return gc1;
}
template< class Coordinate>
Coordinate edgeVertex(Coordinate gc1, const Coordinate &gc2, double factor)
{
for(typename Coordinate::size_type i = 0, end = gc1.size(); i < end; ++i)
{
gc1[i] = gc1[i] + (gc2[i] - gc1[i]) * factor;
}
return gc1;
}
template< class