Skip to content
GitLab
Commit 8912c2e0 authored by Stenger, Florian's avatar Stenger, Florian
Browse files

surfdist-caching, CoordFunction is now applied to intersections

parent fc7e16ac
Expand all Hide whitespace changes
Inline Side-by-side
dune/curvedsurfacegrid/backuprestore.hh
View file @ 8912c2e0
......@@ -52,10 +52,10 @@ namespace Dune
// BackupRestoreFacility for CurvedSurfaceGrid
// -------------------------------------------
template< class HostGrid, class CoordFunction, int order, class Allocator > //! $flo
struct BackupRestoreFacility< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > > //! $flo
template< class HostGrid, class CoordFunction, int order, class Allocator >
struct BackupRestoreFacility< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > >
{
typedef CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > Grid; //! $flo
typedef CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > Grid;
typedef BackupRestoreFacility< HostGrid > HostBackupRestoreFacility;
static void backup ( const Grid &grid, const std::string &filename )
......
dune/curvedsurfacegrid/capabilities.hh
View file @ 8912c2e0
......@@ -24,7 +24,7 @@ namespace Dune
// Capabilities from dune-grid
// ---------------------------
template< class HostGrid, class CoordFunction, int order, class Allocator > //! $flo
template< class HostGrid, class CoordFunction, int order, class Allocator >
struct hasSingleGeometryType< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > >
{
static const bool v = hasSingleGeometryType< HostGrid > :: v;
......@@ -32,52 +32,52 @@ namespace Dune
};
template< class HostGrid, class CoordFunction, int order, class Allocator, int codim > //! $flo
template< class HostGrid, class CoordFunction, int order, class Allocator, int codim >
struct hasEntity< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator >, codim >
{
static const bool v = true;
};
template< class HostGrid, class CoordFunction, int order, class Allocator, int codim > //! $flo
template< class HostGrid, class CoordFunction, int order, class Allocator, int codim >
struct hasEntityIterator< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator >, codim >
{
static const bool v = true;
};
template< class HostGrid, class CoordFunction, int order, class Allocator, int codim > //! $flo
template< class HostGrid, class CoordFunction, int order, class Allocator, int codim >
struct canCommunicate< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator >, codim >
{
static const bool v = canCommunicate< HostGrid, codim >::v && hasEntity< HostGrid, codim >::v;
};
template< class HostGrid, class CoordFunction, int order, class Allocator > //! $flo
template< class HostGrid, class CoordFunction, int order, class Allocator >
struct hasBackupRestoreFacilities< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > >
{
static const bool v = hasBackupRestoreFacilities< HostGrid >::v;
};
template< class HostGrid, class CoordFunction, int order, class Allocator > //! $flo
template< class HostGrid, class CoordFunction, int order, class Allocator >
struct isLevelwiseConforming< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > >
{
static const bool v = isLevelwiseConforming< HostGrid >::v;
};
template< class HostGrid, class CoordFunction, int order, class Allocator > //! $flo
template< class HostGrid, class CoordFunction, int order, class Allocator >
struct isLeafwiseConforming< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > >
{
static const bool v = isLeafwiseConforming< HostGrid >::v;
};
template< class HostGrid, class CoordFunction, int order, class Allocator > //! $flo
template< class HostGrid, class CoordFunction, int order, class Allocator >
struct threadSafe< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > >
{
static const bool v = false;
};
template< class HostGrid, class CoordFunction, int order, class Allocator > //! $flo
template< class HostGrid, class CoordFunction, int order, class Allocator >
struct viewThreadSafe< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator > >
{
static const bool v = false;
......@@ -85,59 +85,15 @@ namespace Dune
// hasHostEntity
// -------------
//! obsolete (already defined in geometrygrid/capabilities.hh)
// template< class Grid, int codim >
// struct hasHostEntity;
//
// template< class Grid, int codim >
// struct hasHostEntity< const Grid, codim >
// {
// static const bool v = hasHostEntity< Grid, codim >::v;
// };
template< class HostGrid, class CoordFunction, int order, class Allocator, int codim > //! $flo
template< class HostGrid, class CoordFunction, int order, class Allocator, int codim >
struct hasHostEntity< CurvedSurfaceGrid< HostGrid, CoordFunction, order, Allocator >, codim >
{
static const bool v = hasEntity< HostGrid, codim >::v;
};
// CodimCache
// ----------
//! obsolete (already defined in geometrygrid/capabilities.hh)
// template< class Grid >
// class CodimCache
// {
// static const int dimension = Grid::dimension;
//
// bool hasHostEntity_[ Grid::dimension + 1 ];
//
// CodimCache ()
// {
// Hybrid::forEach( Std::make_index_sequence< dimension+1 >{},
// [ & ]( auto i ){ hasHostEntity_[ i ] = Capabilities::hasHostEntity< Grid, i >::v; } );
// }
//
// static CodimCache &instance ()
// {
// static CodimCache singleton;
// return singleton;
// }
//
// public:
// static bool hasHostEntity ( int codim )
// {
// assert( (codim >= 0) && (codim <= dimension) );
// return instance().hasHostEntity_[ codim ];
// }
// };
} // namespace Capabilities
} // namespace Dune
......
dune/curvedsurfacegrid/coordfunctionimplementations.hh
View file @ 8912c2e0
......@@ -4,12 +4,66 @@
#define DUNE_CRVSRF_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>
#define COLLECTOR_BENCHMARK 0
#if COLLECTOR_BENCHMARK != 0
#include <ctime>
//wallclock-time in seconds
double dbl_time(){
#ifdef NO_CLOCK_GETTIME
return (double)time(NULL);
#else
struct timespec now;
clock_gettime(CLOCK_MONOTONIC, &now);
return (double)now.tv_nsec / 1000000000.0 + now.tv_sec;
#endif
}
#endif
namespace Dune
{
//CenterSphereCoordFunction
struct CenterSphereCoordFunction
: public AnalyticalCoordFunction<double, 3, 3, CenterSphereCoordFunction>
{
//constrcutor
CenterSphereCoordFunction(double radius) : radius_(radius){}
//operator()
FieldVector<double, 3> operator()(FieldVector<double, 3> x) const{
#if COLLECTOR_BENCHMARK != 0
inputs.push_back(x);
#endif
x *= radius_ / x.two_norm();
return x;
}
#if COLLECTOR_BENCHMARK != 0
void collector_benchmark(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)
<< ", count: " << inputs.size() << std::endl;
}
#endif
private:
double radius_;
#if COLLECTOR_BENCHMARK != 0
mutable std::vector<FieldVector<double, 3> > inputs;
#endif
};
//SphereCoordFunction
struct SphereCoordFunction
: public AnalyticalCoordFunction<double, 3, 3, SphereCoordFunction>
......@@ -20,18 +74,39 @@ namespace Dune
//operator()
FieldVector<double, 3> operator()(FieldVector<double, 3> x) const{
#if COLLECTOR_BENCHMARK != 0
inputs.push_back(x);
#endif
FieldVector<double, 3> direction = x - center_;
double len = direction.two_norm();
for(int i = 0; i < 3; ++i) x[i] = center_[i] + direction[i] * radius_ / len;
double factor = radius_ / direction.two_norm();
for(int i = 0; i < 3; ++i) x[i] = center_[i] + direction[i] * factor;
return x;
}
#if COLLECTOR_BENCHMARK != 0
void collector_benchmark(bool fresh_cache){
double start_time = dbl_time();
for(auto &input : inputs){
FieldVector<double, 3> direction = input - center_;
double factor = radius_ / direction.two_norm();
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)
<< ", count: " << inputs.size() << std::endl;
}
#endif
private:
FieldVector<double, 3> center_;
double radius_;
#if COLLECTOR_BENCHMARK != 0
mutable std::vector<FieldVector<double, 3> > inputs;
#endif
};
//SurfaceDistanceCoordFunction
struct SurfaceDistanceCoordFunction
: public AnalyticalCoordFunction<double, 3, 3, SurfaceDistanceCoordFunction>
......@@ -41,25 +116,43 @@ namespace Dune
surface_mesh = new crvsrf::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);
}
//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
}*/
//destructor
~SurfaceDistanceCoordFunction(){
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]);
backgrid->distance_to_surface_3d(v, hit_vertex);
//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);
}
return {hit_vertex[0], hit_vertex[1], hit_vertex[2]};
}
......@@ -68,9 +161,41 @@ namespace Dune
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)){
backgrid->distance_to_surface_3d(v, hit_vertex);
cache->insert(v, hit_vertex);
++uncached_cnt;
}
input = {hit_vertex[0], hit_vertex[1], hit_vertex[2]};
}
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;
#if COLLECTOR_BENCHMARK != 0
mutable std::vector<FieldVector<double, 3> > inputs;
#endif
};
} // namespace Dune
......
dune/curvedsurfacegrid/coordprovider.hh
View file @ 8912c2e0
......@@ -3,7 +3,6 @@
#ifndef DUNE_CRVSRF_COORDPROVIDER_HH
#define DUNE_CRVSRF_COORDPROVIDER_HH
//!#include <dune/grid/geometrygrid/coordfunctioncaller.hh>
#include <array>
namespace Dune
......@@ -12,7 +11,6 @@ namespace Dune
namespace crvsrf
{
//! $flo
template< class Coordinate>
Coordinate edgeCenter(Coordinate gc1, const Coordinate &gc2)
{
......@@ -44,7 +42,8 @@ namespace Dune
return gc1;
}
//! $flo
// InterpolatoryVerticesGenerator
// ------------------------------
......@@ -192,41 +191,7 @@ namespace Dune
}
};
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//! $flo: $continue
/*template< class CoordFunction >
struct CoordFunctionCaller;
template<> //! more arguments! see CoordFunctionCaller from GeometryGrid!
struct CoordFunctionCaller<AnalyticalCoordFunctionInterface< ct, dimD, dimR, Impl > >
{
typedef AnalyticalCoordFunctionInterface< ct, dimD, dimR, Impl > CoordFunctionInterface;
CoordFunctionCaller(Entity &entity, vector< FieldVector< ctype, dim > > &vertices,
CoordFunctionInterface coordFunction)
: entity_(entity), vertices_(vertices), coordFunction_(coordFunction){}
void evaluate(int index, FieldVector< ctype, dim > &result)
{
coordFunction.evaluate(vertices_[index], result);
}
private:
Entity &entity_
vector< FieldVector< ctype, dim > > &vertices_;
CoordFunctionInterface &coordFunction_;
};
template<>
struct CoordFunctionCaller<DiscreteCoordFunction>
{
void evaluate(const vector<FieldVector< ctype, dim > > &vertices,
int index, FieldVector< ctype, dim > &outp)
{
//!coordFunction_.evaluate(new_local_entity?, index, outp); //! how is this going to work?
}
};*/
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
// CoordProvider (adapted CoordVector from GeometryGrid)
// -----------------------------------------------------
......@@ -262,14 +227,14 @@ namespace Dune
hostEntity_(hostEntity)
{}
void calculate ( std::vector< Coordinate > &vertices ) const //! $flo
void calculate ( std::vector< Coordinate > &vertices ) const
{
const std::size_t numCorners = hostEntity_.geometry().corners();
std::array<Coordinate, (1 << mydim)> corners; //! $flo
for( std::size_t i = 0; i < numCorners; ++i ) //! $flo
std::array<Coordinate, (1 << mydim)> corners;
for( std::size_t i = 0; i < numCorners; ++i )
corners[i] = hostEntity_.geometry().corner(i);
InterpolatoryVerticesGenerator<Coordinate, mydim, Grid::order>::generate(corners, vertices); //! $flo
for(auto &v : vertices) coordFunction_.evaluate(v, v); //! $flo
InterpolatoryVerticesGenerator<Coordinate, mydim, Grid::order>::generate(corners, vertices);
for(auto &v : vertices) coordFunction_.evaluate(v, v);
}
private:
......@@ -302,28 +267,28 @@ namespace Dune
const unsigned int subEntity,
const CoordFunction &coordFunction )
: coordFunction_(coordFunction),
hostElement_(hostElement), //! $flo
hostElement_(hostElement),
subEntity_( subEntity )
{}
void calculate ( std::vector< Coordinate > &vertices ) const //! $flo
void calculate ( std::vector< Coordinate > &vertices ) const
{
const GeometryType type = hostElement_.geometry().type(); //! $flo
const GeometryType type = hostElement_.geometry().type();
auto refElement = referenceElement< ctype, dimension >( type );
const std::size_t numCorners = refElement.size( subEntity_, codimension, dimension );
std::array<Coordinate, (1 << mydim)> corners; //! $flo
std::array<Coordinate, (1 << mydim)> corners;
for( std::size_t i = 0; i < numCorners; ++i )
{
const std::size_t j = refElement.subEntity( subEntity_, codimension, i, dimension );
corners[i] = hostElement_.geometry().corner(j);
}
InterpolatoryVerticesGenerator<Coordinate, mydim, Grid::order>::generate(corners, vertices); //! $flo
for(auto &v : vertices) coordFunction_.evaluate(v, v); //! $flo
InterpolatoryVerticesGenerator<Coordinate, mydim, Grid::order>::generate(corners, vertices);
for(auto &v : vertices) coordFunction_.evaluate(v, v);
}
private:
const CoordFunction &coordFunction_;
const HostElement &hostElement_; //! $flo
const HostElement &hostElement_;
const unsigned int subEntity_;
};
......@@ -355,20 +320,21 @@ namespace Dune
public:
IntersectionCoordProvider ( const ElementGeometryImpl &elementGeometry,
const HostLocalGeometry &hostLocalGeometry,
const CoordFunction & /*coordFunction*/ )
: //coordFunction_(coordFunction), //! $flo (obsolete)
const CoordFunction &coordFunction )
: coordFunction_(coordFunction),
elementGeometry_( elementGeometry ),
hostLocalGeometry_( hostLocalGeometry )
{}
void calculate ( std::vector< Coordinate > &vertices ) const //! $flo
void calculate ( std::vector< Coordinate > &vertices ) const
{
const std::size_t numCorners = hostLocalGeometry_.corners();
std::array<Coordinate, (1 << mydimension)> corners; //! $flo
//! $flo: this does not apply coordFunction on the corners!!!
std::array<Coordinate, (1 << mydimension)> corners;
// $flo: coordFunction is applied here which is in contrast to GeomentryGrid's behaviour!
for( std::size_t i = 0; i < numCorners; ++i )
corners[ i ] = elementGeometry_.global( hostLocalGeometry_.corner( i ) );
InterpolatoryVerticesGenerator<Coordinate, mydimension, Grid::order>::generate(corners, vertices); //! $flo
InterpolatoryVerticesGenerator<Coordinate, mydimension, Grid::order>::generate(corners, vertices);
for(auto &v : vertices) coordFunction_.evaluate(v, v);
}
template< unsigned int numCorners >
......@@ -378,53 +344,11 @@ namespace Dune
}
private:
//!const CoordFunction &coordFunction_; //! $flo (obsolete)
const CoordFunction &coordFunction_;
const ElementGeometryImpl &elementGeometry_;
HostLocalGeometry hostLocalGeometry_;
};
//! $flo (obsolete)
//! CornerStorage
//! -------------
/*!template< int mydim, int cdim, class Grid >
class CornerStorage
{
typedef typename std::remove_const< Grid >::type::Traits Traits;
typedef typename Traits::ctype ctype;
typedef FieldVector< ctype, cdim > Coordinate;
//! $flo: need more entries in array to accommodate for all dofs!
typedef std::array< Coordinate, (1 << mydim) > Coords;
public:
typedef typename Coords::const_iterator const_iterator;
template< bool fake >
explicit CornerStorage ( const CoordVector< mydim, Grid, fake > &coords )
{
coords.calculate( coords_ );
}
explicit CornerStorage ( const IntersectionCoordVector< Grid > &coords )
{
coords.calculate( coords_ );
}
const Coordinate &operator[] ( unsigned int i ) const
{
return coords_[ i ];
}
const_iterator begin () const { return coords_.begin(); }
const_iterator end () const { return coords_.end(); }
private:
Coords coords_;
};*/
} // namespace crvsrf
} // namespace Dune
......
dune/curvedsurfacegrid/entity.hh
View file @ 8912c2e0
......@@ -126,7 +126,7 @@ namespace Dune
private:
typedef typename HostGrid::template Codim< codimension >::Geometry HostGeometry;
typedef crvsrf::CoordProvider< mydimension, Grid, fake > CoordProvider; //! $flo
typedef crvsrf::CoordProvider< mydimension, Grid, fake > CoordProvider;
public:
/** \name Construction, Initialization and Destruction
......@@ -252,11 +252,9 @@ namespace Dune
{
if( !geo_ )
{
//!std::cout << "create geometry " << codimension << std::endl; //! $flo (debug)
CoordProvider coords( hostEntity(), grid().coordFunction() ); //! $flo
CoordProvider coords( hostEntity(), grid().coordFunction() );
geo_ = GeometryImpl( grid(), type(), coords );
}
//!else std::cout << "load cached geometry " << codimension << std::endl; //! $flo (debug)
return Geometry( geo_ );
}
......@@ -417,7 +415,7 @@ namespace Dune
private:
typedef typename HostGrid::template Codim< 0 >::Geometry HostGeometry;
typedef crvsrf::CoordProvider< mydimension, Grid, fake > CoordProvider; //! $flo
typedef crvsrf::CoordProvider< mydimension, Grid, fake > CoordProvider;
public:
/** \name Construction, Initialization and Destruction
......@@ -571,7 +569,7 @@ namespace Dune
{
if( !geo_ )
{
CoordProvider coords( hostElement(), subEntity_, grid().coordFunction() ); //! $flo
CoordProvider coords( hostElement(), subEntity_, grid().coordFunction() );
geo_ = GeometryImpl( grid(), type(), coords );
}
return Geometry( geo_ );
......
dune/curvedsurfacegrid/geometry.hh
View file @ 8912c2e0
......@@ -7,7 +7,7 @@
#include <dune/common/typetraits.hh>
#include <dune/curvilineargeometry/curvilineargeometry.hh> //! $flo
#include <dune/curvilineargeometry/curvilineargeometry.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>
......@@ -63,12 +63,6 @@ namespace Dune
static ctype tolerance () { return 16 * std::numeric_limits< ctype >::epsilon(); }
//!template< int mydim, int cdim > //! $flo: (obsolete)
//!struct CornerStorage