Skip to content
GitLab
Commit 9c87c77c authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

implementation of curved geometries based on grid-function evaluation

parent 6c3f2da0
Hide whitespace changes
Inline Side-by-side
dune/curvedsurfacegrid/backuprestore.hh
View file @ 9c87c77c
......@@ -51,10 +51,11 @@ namespace Dune
// BackupRestoreFacility for CurvedSurfaceGrid
// -------------------------------------------
template< class HostGrid, class CoordFunction, int order, class Allocator >
struct BackupRestoreFacility<CurvedSurfaceGrid<HostGrid, CoordFunction, order, Allocator>>
template< class GridFunction, int order >
struct BackupRestoreFacility<CurvedSurfaceGrid<GridFunction, order>>
{
using Grid = CurvedSurfaceGrid<HostGrid, CoordFunction, order, Allocator>;
using Grid = CurvedSurfaceGrid<GridFunction, order>;
using HostGrid = typename Grid::HostGrid;
using HostBackupRestoreFacility = BackupRestoreFacility<HostGrid>;
static void backup (const Grid& grid, const std::string& filename)
......@@ -71,15 +72,14 @@ namespace Dune
static Grid* restore (const std::string& filename)
{
// notice: We should also restore the coordinate function
HostGrid* hostGrid = HostBackupRestoreFacility::restore(filename);
return new Grid(hostGrid, CoordFunction{});
assert(false && "Restore not yet implemented for CurvedSurfaceGrid");
return nullptr;
}
static Grid* restore (const std::istream& stream)
{
HostGrid* hostGrid = HostBackupRestoreFacility::restore(stream);
return new Grid(hostGrid, CoordFunction{});
assert(false && "Restore not yet implemented for CurvedSurfaceGrid");
return nullptr;
}
};
......
dune/curvedsurfacegrid/capabilities.hh
View file @ 9c87c77c
......@@ -24,61 +24,66 @@ namespace Dune
// Capabilities from dune-grid
// ---------------------------
template< class HostGrid, class CoordFunction, int order, class Allocator >
struct hasSingleGeometryType< CurvedSurfaceGrid<HostGrid, CoordFunction, order, Allocator> >
template< class GridFunction, int order >
struct hasSingleGeometryType< CurvedSurfaceGrid<GridFunction,order> >
{
using HostGrid = typename GridFunction::GridView::Grid;
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 GridFunction, int order, int codim >
struct hasEntity< CurvedSurfaceGrid<GridFunction,order>, 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 GridFunction, int order, int codim >
struct hasEntityIterator< CurvedSurfaceGrid<GridFunction,order>, 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 GridFunction, int order, int codim >
struct canCommunicate< CurvedSurfaceGrid<GridFunction,order>, codim >
{
using HostGrid = typename GridFunction::GridView::Grid;
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 GridFunction, int order >
struct hasBackupRestoreFacilities< CurvedSurfaceGrid<GridFunction,order> >
{
using HostGrid = typename GridFunction::GridView::Grid;
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 GridFunction, int order >
struct isLevelwiseConforming< CurvedSurfaceGrid<GridFunction,order> >
{
using HostGrid = typename GridFunction::GridView::Grid;
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 GridFunction, int order >
struct isLeafwiseConforming< CurvedSurfaceGrid<GridFunction,order> >
{
using HostGrid = typename GridFunction::GridView::Grid;
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 GridFunction, int order >
struct threadSafe< CurvedSurfaceGrid<GridFunction,order> >
{
static const bool v = false;
};
template< class HostGrid, class CoordFunction, int order, class Allocator >
struct viewThreadSafe< CurvedSurfaceGrid<HostGrid, CoordFunction, order, Allocator> >
template< class GridFunction, int order >
struct viewThreadSafe< CurvedSurfaceGrid<GridFunction,order> >
{
static const bool v = false;
};
......@@ -88,9 +93,10 @@ namespace Dune
// hasHostEntity
// -------------
template< class HostGrid, class CoordFunction, int order, class Allocator, int codim >
struct hasHostEntity< CurvedSurfaceGrid<HostGrid, CoordFunction, order, Allocator>, codim >
template< class GridFunction, int order, int codim >
struct hasHostEntity< CurvedSurfaceGrid<GridFunction,order>, codim >
{
using HostGrid = typename GridFunction::GridView::Grid;
static const bool v = hasEntity< HostGrid, codim >::v;
};
......
dune/curvedsurfacegrid/datahandle.hh
View file @ 9c87c77c
......@@ -47,7 +47,7 @@ namespace Dune
{
using Entity = typename Grid::Traits::template Codim<HostEntity::codimension>::Entity;
using EntityImpl = typename Grid::Traits::template Codim<HostEntity::codimension>::EntityImpl;
Entity entity(EntityImpl(grid_, hostEntity));
Entity entity(EntityImpl(grid_.gridFunction(), hostEntity));
return wrappedHandle_.size(entity);
}
......@@ -56,7 +56,7 @@ namespace Dune
{
using Entity = typename Grid::Traits::template Codim<HostEntity::codimension>::Entity;
using EntityImpl = typename Grid::Traits::template Codim<HostEntity::codimension>::EntityImpl;
Entity entity(EntityImpl(grid_, hostEntity));
Entity entity(EntityImpl(grid_.gridFunction(), hostEntity));
wrappedHandle_.gather(buffer, entity);
}
......@@ -65,7 +65,7 @@ namespace Dune
{
using Entity = typename Grid::Traits::template Codim< HostEntity::codimension >::Entity;
using EntityImpl = typename Grid::Traits::template Codim< HostEntity::codimension >::EntityImpl;
Entity entity(EntityImpl(grid_, hostEntity));
Entity entity(EntityImpl(grid_.gridFunction(), hostEntity));
wrappedHandle_.scatter(buffer, entity, size);
}
......
dune/curvedsurfacegrid/declaration.hh
View file @ 9c87c77c
......@@ -6,7 +6,7 @@
namespace Dune
{
template< class HostGrid, class CoordFunction, int lagrangeOrder, class Allocator >
template< class GridFunction, int order >
class CurvedSurfaceGrid;
} // namespace Dune
......
dune/curvedsurfacegrid/entity.hh
View file @ 9c87c77c
......@@ -7,6 +7,7 @@
#include <dune/geometry/referenceelements.hh>
#include <dune/grid/common/grid.hh>
#include <dune/curvedsurfacegrid/geometry.hh>
namespace Dune
{
......@@ -53,21 +54,10 @@ namespace Dune
class IntersectionIterator;
// EntityBase (real)
// -----------------
/** \copydoc EntityBase
*
* This specialization implements the case, where the host grid provides
* the entity for this codimension, i.e., \em fake = \b false.
*
* \nosubgrouping
*/
template< int codim, class Grid >
template< int codim, class G >
class EntityBase
{
using Traits = typename std::remove_const_t<Grid>::Traits;
using Traits = typename std::remove_const_t<G>::Traits;
public:
/** \name Attributes
......@@ -90,65 +80,57 @@ namespace Dune
//! coordinate type of the grid
using ctype = typename Traits::ctype;
//! type of corresponding geometry
using Geometry = typename Traits::template Codim<codimension>::Geometry;
/** \} */
private:
using HostGrid = typename Traits::HostGrid;
public:
/** \name Host Types
* \{ */
//! type of corresponding host entity
using HostEntity = typename HostGrid::template Codim<codimension>::Entity;
//! type of corresponding entity seed
using EntitySeed = typename Traits::template Codim<codimension>::EntitySeed;
//! type of host elements, i.e., of host entities of codimension 0
using HostElement = typename HostGrid::template Codim<0>::Entity;
/** \} */
private:
using GeometryImpl = typename Traits::template Codim<codimension>::GeometryImpl;
using HostGeometry = typename HostGrid::template Codim<codimension>::Geometry;
public:
/** \name Construction, Initialization and Destruction
* \{ */
using Grid = typename Traits::Grid;
using GridFunction = typename Traits::GridFunction;
//! type of the host grid
using HostGrid = typename Traits::HostGrid;
//! type of corresponding host entity
using HostEntity = typename HostGrid::template Codim<codim>::Entity;
//! type of host elements, i.e., of host entities of codimension 0
using HostElement = typename HostGrid::template Codim<0>::Entity;
public:
EntityBase () = default;
// Construct the entity from an entity seed
EntityBase (const Grid& grid, const EntitySeed& seed)
: hostEntity_(grid.hostGrid().entity(seed.impl().hostEntitySeed()))
, grid_(&grid)
, gridFunction_(&grid.gridFunction())
{}
// construct the entity from a subentity of a host-entity
EntityBase (const Grid& grid, const HostElement& hostElement, int i)
EntityBase (const GridFunction& gridFunction, const HostElement& hostElement, int i)
: hostEntity_(hostElement.template subEntity<codim>(i))
, grid_(&grid)
, gridFunction_(&gridFunction)
{}
// construct the entity from a host-entity
EntityBase (const Grid& grid, const HostEntity& hostEntity)
EntityBase (const GridFunction& gridFunction, const HostEntity& hostEntity)
: hostEntity_(hostEntity)
, grid_(&grid)
, gridFunction_(&gridFunction)
{}
// construct the entity from a host-entity
EntityBase (const Grid& grid, HostEntity&& hostEntity)
EntityBase (const GridFunction& gridFunction, HostEntity&& hostEntity)
: hostEntity_(std::move(hostEntity))
, grid_(&grid)
, gridFunction_(&gridFunction)
{}
/** \} */
//! compare two entities
bool equals (const EntityBase& other) const
{
......@@ -180,42 +162,6 @@ namespace Dune
return hostEntity().partitionType();
}
//! obtain the geometry of this entity
/**
* Each DUNE entity encapsulates a geometry object, representing the map
* from the reference element to world coordinates. Wrapping the geometry
* is the main objective of the CurvedSurfaceGrid.
*
* The CurvedSurfaceGrid provides geometries by parametrization with local basis
* functions, using the CurvedGeometry.
*
* \returns a new curvedgeometry object
*/
Geometry geometry () const
{
if (!geo_) {
// mapping from local to curved global coordinates
auto ff = [f=grid().coordFunction(),geo=hostEntity().geometry()](const auto& local) {
return f(geo.global(local));
};
if (grid_->useGeometryCaching()) {
auto const& idSet = grid_->hostGrid().localIdSet();
auto id = idSet.id(hostEntity());
auto& cache = std::get<codim>(grid_->geometryCache_);
auto [it,inserted] = cache.try_emplace(id, type(), ff);
geo_ = it->second;
}
else {
geo_ = GeometryImpl(type(), ff);
}
}
return Geometry(*geo_);
}
//! obtain number of sub-entities of the current entity
unsigned int subEntities (unsigned int cc) const
{
......@@ -228,16 +174,13 @@ namespace Dune
return typename EntitySeed::Implementation(hostEntity().seed());
}
/** \} */
/** \name Methods Supporting the Grid Implementation
* \{ */
const Grid& grid () const
const GridFunction& gridFunction () const
{
assert(grid_);
return *grid_;
return *gridFunction_;
}
//! return the wrapped host-entity
......@@ -249,20 +192,99 @@ namespace Dune
/** \} */
private:
HostEntity hostEntity_ = {};
const Grid* grid_ = nullptr;
mutable Std::optional<GeometryImpl> geo_;
HostEntity hostEntity_;
const GridFunction* gridFunction_ = nullptr;
};
template< int codim, class G >
class EntityBaseGeometry
: public EntityBase<codim, G>
{
using Super = EntityBase<codim, G>;
using Traits = typename std::remove_const_t<G>::Traits;
public:
//! type of corresponding geometry
using Geometry = typename Traits::template Codim<codim>::Geometry;
public:
using Super::Super;
//! obtain the geometry of this entity
Geometry geometry () const
{
if (!geo_) {
auto localFct = localFunction(Super::gridFunction());
// localFct.bind(Super::hostEntity());
// construct a fake local geometry for transformations.
// auto refElem = referenceElement<typename Geometry::ctype, Super::dimension>(GeometryTypes::simplex(Super::dimension));
geo_ = std::make_shared<GeometryImpl>(Super::type(), localFct);
}
return Geometry(*geo_);
//return Super::hostEntity().geometry();
}
private:
using GeometryImpl = typename Traits::template Codim<codim>::GeometryImpl;
mutable std::shared_ptr<GeometryImpl> geo_;
};
template< class G >
class EntityBaseGeometry<0, G>
: public EntityBase<0, G>
{
using Super = EntityBase<0, G>;
using Traits = typename std::remove_const_t<G>::Traits;
public:
//! type of corresponding geometry
using Geometry = typename Traits::template Codim<0>::Geometry;
public:
using Super::Super;
//! obtain the geometry of this entity
/**
* Each DUNE entity encapsulates a geometry object, representing the map
* from the reference element to world coordinates. Wrapping the geometry
* is the main objective of the CurvedSurfaceGrid.
*
* The CurvedSurfaceGrid provides geometries by parametrization with local basis
* functions, using the CurvedGeometry.
*
* \returns a new curvedgeometry object
*/
Geometry geometry () const
{
if (!geo_) {
auto localFct = localFunction(Super::gridFunction());
localFct.bind(Super::hostEntity());
geo_ = std::make_shared<GeometryImpl>(Super::type(), localFct, Dune::CGeo::Impl::DefaultLocalGeometry{});
}
return Geometry(*geo_);
}
private:
using GeometryImpl = typename Traits::template Codim<0>::GeometryImpl;
mutable std::shared_ptr<GeometryImpl> geo_;
};
// Entity
// ------
template< int codim, int dim, class Grid >
template< int codim, int dim, class G >
class Entity
: public EntityBase<codim, Grid>
: public EntityBaseGeometry<codim, G>
{
using Super = EntityBase<codim, Grid>;
using Super = EntityBaseGeometry<codim, G>;
public:
// import constructors from base class
......@@ -273,13 +295,15 @@ namespace Dune
// Entity for codimension 0
// ------------------------
template< int dim, class Grid >
class Entity<0, dim, Grid>
: public EntityBase<0, Grid>
template< int dim, class G >
class Entity<0, dim, G>
: public EntityBaseGeometry<0, G>
{
using Super = EntityBase<0, Grid>;
using Super = EntityBaseGeometry<0, G>;
using Traits = typename std::remove_const_t<Grid>::Traits;
using Traits = typename std::remove_const_t<G>::Traits;
using Grid = typename Traits::Grid;
using GridFunction = typename Grid::GridFunction;
using HostGrid = typename Traits::HostGrid;
public:
......@@ -291,7 +315,7 @@ namespace Dune
using LocalGeometry = typename Traits::template Codim<0>::LocalGeometry;
//! facade type for entities
using EntityFacade = Dune::Entity<0, dim, Grid, Dune::CGeo::Entity>;
using EntityFacade = Dune::Entity<0, dim, G, Dune::CGeo::Entity>;
//! type of hierarchic iterator
using HierarchicIterator = typename Traits::HierarchicIterator;
......@@ -311,7 +335,7 @@ namespace Dune
typename Grid::template Codim<codim>::Entity subEntity (int i) const
{
using EntityImpl = typename Traits::template Codim<codim>::EntityImpl;
return EntityImpl(Super::grid(), Super::hostEntity(), i);
return EntityImpl(Super::gridFunction(), Super::hostEntity(), i);
}
LevelIntersectionIterator ilevelbegin () const
......@@ -350,7 +374,7 @@ namespace Dune
EntityFacade father () const
{
return Entity(Super::grid(), Super::hostEntity().father());
return Entity(Super::gridFunction(), Super::hostEntity().father());
}
bool hasFather () const
......@@ -365,14 +389,14 @@ namespace Dune
HierarchicIterator hbegin (int maxLevel) const
{
using IteratorImpl = CGeo::HierarchicIterator<Grid>;
return IteratorImpl(Super::grid(), Super::hostEntity().hbegin(maxLevel));
using IteratorImpl = CGeo::HierarchicIterator<G>;
return IteratorImpl(Super::gridFunction(), Super::hostEntity().hbegin(maxLevel));
}
HierarchicIterator hend (int maxLevel) const
{
using IteratorImpl = CGeo::HierarchicIterator<Grid>;
return IteratorImpl(Super::grid(), Super::hostEntity().hend(maxLevel));
using IteratorImpl = CGeo::HierarchicIterator<G>;
return IteratorImpl(Super::gridFunction(), Super::hostEntity().hend(maxLevel));
}
bool isRegular () const
......
dune/curvedsurfacegrid/entityseed.hh
View file @ 9c87c77c
......@@ -16,10 +16,10 @@ namespace Dune
// EntitySeed
// ----------
template< int codim, class Grd >
template< int codim, class G >
class EntitySeed
{
using Traits = typename std::remove_const_t<Grd>::Traits;
using Traits = typename std::remove_const_t<G>::Traits;
using HostGrid = typename Traits::HostGrid;
using HostEntitySeed = typename HostGrid::template Codim<codim>::EntitySeed;
......
dune/curvedsurfacegrid/geometry.hh 0 → 100644
View file @ 9c87c77c
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_CURVED_SURFACE_GRID_GEOMETRY_HH
#define DUNE_CURVED_SURFACE_GRID_GEOMETRY_HH
#include <cassert>
#include <functional>
#include <iterator>
#include <limits>
#include <vector>
#include <dune/common/diagonalmatrix.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/std/optional.hh>
#include <dune/common/std/type_traits.hh>
#include <dune/geometry/affinegeometry.hh>
#include <dune/geometry/multilineargeometry.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>
namespace Dune
{
namespace CGeo
{
namespace Impl
{
struct IdentityMatrix {};
template <class K, int n, int m, int l>
FieldMatrix<K,n,m> ABt(const FieldMatrix<K,n,l>& A, const FieldMatrix<K,m,l>& B)
{
FieldMatrix<K,n,m> ABt;
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j) {
ABt[i][j] = 0;
for (int k = 0; k < l; ++k)
ABt[i][j] += A[i][k] * B[j][k];
}
return ABt;
}
template <class K, int n, int m>
FieldMatrix<K,n,m> ABt(const DiagonalMatrix<K,n>& A, const FieldMatrix<K,m,n>& B)
{
FieldMatrix<K,n,m> ABt;
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j) {
ABt[i][j] = A[i][i] * B[j][i];
}
return ABt;
}
template <class K, int n, int m, int l>
FieldMatrix<K,n,m> AB(const FieldMatrix<K,n,l>& A, const FieldMatrix<K,l,m>& B)
{