Skip to content

GitLab

Commit 210d43c1 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

added poisson example from wiki, made vtk-write compatible grid

parent 412e05e3
Hide whitespace changes
Inline Side-by-side
dune/dec/FlatSharp.hpp
View file @ 210d43c1
......@@ -71,7 +71,7 @@ namespace Dec
static auto impl(GridView const& gv, Vector const& field, int_t<k>)
{
static_assert( dim == GridView::dimension, "" );
return impl(gv, field, int_<k>, valid_);
return Flat::impl(gv, field, int_<k>, valid_);
}
};
......@@ -80,7 +80,7 @@ namespace Dec
struct Flat<1, dim>
{
// define the transformation for vector fields living on the vertices
template <class GridView, class WorldVector = typename GridView::GlobalCoords>
template <class GridView, class WorldVector = typename GridView::GlobalCoordinate>
static DOFVector<float_type> impl(GridView const& gv, DOFVector<WorldVector> const& field, int_t<0>)
{
static_assert( dim == GridView::dimension, "" );
......@@ -146,7 +146,7 @@ namespace Dec
}
template <class GridView, class Vector, int k>
static auto impl(GridView const& gv, Vector const& form, int_t<k>) { return impl(gv, form, int_<k>, valid_); }
static auto impl(GridView const& gv, Vector const& form, int_t<k>) { return Sharp::impl(gv, form, int_<k>, valid_); }
};
......@@ -155,7 +155,7 @@ namespace Dec
{
// sharp a 2-form
template <class GridView, class Vector>
static Vector impl(GridView const& gv, Vector const& form, int_t<dim>, valid)
static Vector impl(GridView const& gv, Vector const& form, int_t<dim>, bool_t<false>)
{
static_assert( dim == GridView::dimension, "" );
......@@ -171,8 +171,8 @@ namespace Dec
// sharp a 1-form to a vector-form by averaging over all edges of a cell
// => cell-vector
template <class GridView, class V, class WorldVector = typename GridView::GlobalCoords>
static DOFVector<WorldVector> impl(GridView const& gv, DOFVector<V> const& form, int_t<1>, requires_t<(dim>1)>)
template <class GridView, class V, class WorldVector = typename GridView::GlobalCoordinate>
static DOFVector<WorldVector> impl(GridView const& gv, DOFVector<V> const& form, int_t<1>, bool_t<true>)
{
constexpr int nVertices = dim+1;
constexpr int nEdges = kombinations(nVertices, 2);
......@@ -210,7 +210,7 @@ namespace Dec
}
template <class GridView, class Vector, int k>
static auto impl(GridView const& gv, Vector const& form, int_t<k>) { return impl(gv, form, int_<k>, valid_); }
static auto impl(GridView const& gv, Vector const& form, int_t<k>) { return Sharp::impl(gv, form, int_<k>, bool_<(k==1 && dim > 1)>); }
};
......
dune/dec/halfedgegrid/Entity.hpp
View file @ 210d43c1
......@@ -29,7 +29,7 @@ namespace Dec
using Geometry = CircumcenterGeometry<mydimension, dimensionworld, Grid>;
using RefElements = Dune::ReferenceElements<Dec::float_type, mydimension>;
using GlobalCoords = typename Geometry::GlobalCoords;
using GlobalCoordinate = typename Geometry::GlobalCoordinate;
template <class> friend class HalfEdgeGrid;
template <int> friend class HalfEdgeIndexSet;
......@@ -39,7 +39,7 @@ namespace Dec
HalfEdgeEntity() = default;
/// Constructor, stores a pointer to `indexSet` and `coordinates`.
HalfEdgeEntity(IndexSet const& indexSet, std::vector<GlobalCoords> const& coordinates, IndexType he)
HalfEdgeEntity(IndexSet const& indexSet, std::vector<GlobalCoordinate> const& coordinates, IndexType he)
: indexSet_(&indexSet)
, coordinates_(&coordinates)
, he_(he)
......@@ -54,6 +54,9 @@ namespace Dec
public:
/// Partition type of this entity
Dune::PartitionType partitionType() const { return Dune::InteriorEntity; }
/// Returns the index of this entity
IndexType index() const
{
......@@ -206,7 +209,7 @@ namespace Dec
IndexSet const* indexSet_ = nullptr;
/// a pointer to the vector of all grid coordinates
std::vector<GlobalCoords> const* coordinates_ = nullptr;
std::vector<GlobalCoordinate> const* coordinates_ = nullptr;
/// an half-edge index representing this entity
IndexType he_ = HalfEdge::invalid;
......
dune/dec/halfedgegrid/Geometry.hpp
View file @ 210d43c1
......@@ -32,10 +32,10 @@ namespace Dec
using ctype = float_type;
/// Type for local coordinate vector
using LocalCoords = Dune::FieldVector<float_type, mydimension>;
using LocalCoordinate = Dune::FieldVector<float_type, mydimension>;
/// Type for coordinate vector in world space
using GlobalCoords = Dune::FieldVector<float_type, dimensionworld>;
using GlobalCoordinate = Dune::FieldVector<float_type, dimensionworld>;
using DuneGeometry = Dune::AffineGeometry<float_type, mydimension, dimensionworld>;
using JacobianTransposed = typename DuneGeometry::JacobianTransposed;
......@@ -54,7 +54,7 @@ namespace Dec
/// Constructor, stores reference to `entity`, `indexSet` and `coordinates`.
CircumcenterGeometry(Entity const& entity,
IndexSet const& indexSet,
std::vector<GlobalCoords> const& coordinates)
std::vector<GlobalCoordinate> const& coordinates)
: entity_(entity)
, indexSet_(indexSet)
, coordinates_(coordinates)
......@@ -62,14 +62,14 @@ namespace Dec
/// Return the center of this entity
GlobalCoords center() const
GlobalCoordinate center() const
{
return center(entity_, Dim_<mydimension>);
}
/// Return the barycenter of this entity
/// NOTE: Redirects to a Dune::AffineGeometry implementation
GlobalCoords barycenter() const
GlobalCoordinate barycenter() const
{
return duneGeometry().center();
}
......@@ -97,7 +97,7 @@ namespace Dec
/// \brief Obtain coordinates of the i-th corner.
/// NOTE: Not efficient, use ranges instead, \see coordinates.
GlobalCoords const& corner(int i) const
GlobalCoordinate const& corner(int i) const
{
auto vertices = indexSet_.vertexIndices(entity_);
auto it = vertices.begin();
......@@ -113,28 +113,28 @@ namespace Dec
/// \brief Return the factor appearing in the integral transformation formula
/// NOTE: Redirects to a Dune::AffineGeometry implementation
float_type integrationElement(LocalCoords const& local) const
float_type integrationElement(LocalCoordinate const& local) const
{
return duneGeometry().integrationElement(local);
}
/// \brief Obtain the transposed of the Jacobian.
/// NOTE: Redirects to a Dune::AffineGeometry implementation
JacobianTransposed const& jacobianTransposed(LocalCoords const& local) const
JacobianTransposed const& jacobianTransposed(LocalCoordinate const& local) const
{
return duneGeometry().jacobianTransposed(local);
}
/// \brief Obtain the transposed of the Jacobian's inverse.
/// NOTE: Redirects to a Dune::AffineGeometry implementation
JacobianInverseTransposed const& jacobianInverseTransposed(LocalCoords const& local) const
JacobianInverseTransposed const& jacobianInverseTransposed(LocalCoordinate const& local) const
{
return duneGeometry().jacobianInverseTransposed(local);
}
/// \brief Return a vector pointing in the direction of the oriented edge.
/// NOTE: Available for edge-entities only.
GlobalCoords vector() const
GlobalCoordinate vector() const
{
static_assert( mydimension == 1, "An edge-vector is defined for edges only!" );
......@@ -148,7 +148,7 @@ namespace Dec
/// \brief Return a vector pointing in the direction of the oriented dual_edge, i.e.
/// pointing from center to center of the ancident cells.
/// NOTE: Available for edge-entities only.
GlobalCoords dual_vector() const
GlobalCoordinate dual_vector() const
{
static_assert( mydimension == 1, "An dual-edge-vector is defined for edges only!" );
......@@ -171,7 +171,7 @@ namespace Dec
// vertex coordinate
template <class EntityType>
GlobalCoords center(EntityType const& entity, Dim_t<0>) const
GlobalCoordinate center(EntityType const& entity, Dim_t<0>) const
{
assert( EntityType::mydimension == 0 );
assert( entity.index() < coordinates_.size() );
......@@ -180,7 +180,7 @@ namespace Dec
// Center of an edge
template <class EntityType>
GlobalCoords center(EntityType const& entity, Dim_t<1>) const
GlobalCoordinate center(EntityType const& entity, Dim_t<1>) const
{
assert( EntityType::mydimension == 1 );
auto vertices = indexSet_.vertexIndices(entity);
......@@ -191,7 +191,7 @@ namespace Dec
// Center of a triangle, TODO: use more efficient formula based on angles
template <class EntityType>
GlobalCoords center(EntityType const& entity, Dim_t<2>) const
GlobalCoordinate center(EntityType const& entity, Dim_t<2>) const
{
assert( EntityType::mydimension == 2 );
auto vertices = indexSet_.vertexIndices(entity);
......@@ -199,9 +199,9 @@ namespace Dec
auto v1 = std::next(v0);
auto v2 = std::next(v1);
GlobalCoords e01 = coordinates_[*v0] - coordinates_[*v1];
GlobalCoords e02 = coordinates_[*v0] - coordinates_[*v2];
GlobalCoords e21 = coordinates_[*v2] - coordinates_[*v1];
GlobalCoordinate e01 = coordinates_[*v0] - coordinates_[*v1];
GlobalCoordinate e02 = coordinates_[*v0] - coordinates_[*v2];
GlobalCoordinate e21 = coordinates_[*v2] - coordinates_[*v1];
auto sqr_e01 = unary_dot(e01);
auto sqr_e02 = unary_dot(e02);
......@@ -270,7 +270,7 @@ namespace Dec
template <class EntityType, int cd>
float_type dual_volume(EntityType const& entity, Codim_t<cd>) const
{
std::vector<GlobalCoords> corners(cd+1);
std::vector<GlobalCoordinate> corners(cd+1);
return fill_dual_volume(entity, corners, 0, Codim_<cd>, Codim_<cd>);
}
......@@ -329,7 +329,7 @@ namespace Dec
Entity const& entity_; //< a reference to the entity, the geometry is related to.
IndexSet const& indexSet_; //< a reference to the global indexSet
std::vector<GlobalCoords> const& coordinates_; //< a reference to the vector of all grid coordinates
std::vector<GlobalCoordinate> const& coordinates_; //< a reference to the vector of all grid coordinates
ReferenceElement const& refElem_ = ReferenceElements::simplex(); //< the corresponding reference element
float_type const tol = std::numeric_limits<float_type>::epsilon();
......
dune/dec/halfedgegrid/GeometryCache.hpp
View file @ 210d43c1
......@@ -26,7 +26,7 @@ namespace Dec
using IndexSet = HalfEdgeIndexSet<dim>;
using IndexType = typename IndexSet::IndexType;
using GlobalCoords = Dune::FieldVector<float_type, Grid::dimensionworld>;
using GlobalCoordinate = Dune::FieldVector<float_type, Grid::dimensionworld>;
template <int cd>
struct Codim
......@@ -61,16 +61,16 @@ namespace Dec
/// Returns a vector of grid vertex coordinates
std::vector<GlobalCoords> const& coordinates() const { return std::get<0>(centers_); }
std::vector<GlobalCoordinate> const& coordinates() const { return std::get<0>(centers_); }
/// Returns the coordinate of the i'th grid vertices
GlobalCoords const& coordinate(std::size_t i) const { return coordinates()[i]; }
GlobalCoordinate const& coordinate(std::size_t i) const { return coordinates()[i]; }
public: // Methods for calculating the center and volume:
template <int codim>
GlobalCoords const& center(std::size_t i) const
GlobalCoordinate const& center(std::size_t i) const
{
return center_impl(i, Dim_<dim-codim>);
}
......@@ -92,10 +92,10 @@ namespace Dec
private: // Implementation details:
GlobalCoords const& center_impl(std::size_t i, Dim_t<0>) const { return coordinate(i); }
GlobalCoordinate const& center_impl(std::size_t i, Dim_t<0>) const { return coordinate(i); }
template <int d>
GlobalCoords const& center_impl(std::size_t i, Dim_t<d>) const
GlobalCoordinate const& center_impl(std::size_t i, Dim_t<d>) const
{
assert( i < std::get<d>(centers_).size() );
return std::get<d>(centers_)[i];
......@@ -171,7 +171,7 @@ namespace Dec
private:
std::array<std::vector<GlobalCoords>, dim+1> centers_; // center coordinates of dim d
std::array<std::vector<GlobalCoordinate>, dim+1> centers_; // center coordinates of dim d
std::array<std::vector<float_type>, dim+1> volumes_; // element volumes of dim d
std::array<std::vector<float_type>, dim+1> dual_volumes_; // element dual-volumes of codim cd
};
......
dune/dec/halfedgegrid/Grid.hpp
View file @ 210d43c1
......@@ -38,7 +38,7 @@ namespace Dec
using IndexType = typename HalfEdge::IndexType;
using Geometry = GeometryCache<HalfEdgeGridFamily<GridBase>>;
using GlobalCoords = typename Geometry::GlobalCoords;
using GlobalCoordinate = typename Geometry::GlobalCoordinate;
using GridView = HalfEdgeGridView<Self>;
using IndexSet = HalfEdgeIndexSet<dimension>;
using LeafGridView = GridView;
......@@ -162,6 +162,9 @@ namespace Dec
return {leafIndexSet(), coordinates(), seed.he_};
}
/// return const reference to a collective communication object.
auto const& comm() const { return grid_.comm(); }
public: // access to geometry information
......@@ -193,7 +196,7 @@ namespace Dec
/// Returns the (circum-)center of the `entity`.
template <class Entity>
GlobalCoords const& center(Entity const& entity) const
GlobalCoordinate const& center(Entity const& entity) const
{
return leafGeometry().template center<Entity::codimension>(entity.index());
}
......
dune/dec/halfedgegrid/GridView.hpp
View file @ 210d43c1
......@@ -8,10 +8,12 @@
namespace Dec
{
template <class Grid>
template <class GridType>
class HalfEdgeGridView
{
public:
using Grid = GridType;
static constexpr int dimension = Grid::dimension;
static constexpr int dimensionworld = Grid::dimensionworld;
......@@ -19,7 +21,7 @@ namespace Dec
using IndexSet = typename Grid::IndexSet;
using Geometry = typename Grid::Geometry;
using GlobalCoords = typename Geometry::GlobalCoords;
using GlobalCoordinate = typename Geometry::GlobalCoordinate;
using ctype = float_type;
......@@ -33,6 +35,7 @@ namespace Dec
{
using value_type = Entity;
using difference_type = std::ptrdiff_t;
using Reference = Entity;
Iterator(HalfEdgeGridView const& gv, IndexType index)
: gv_(gv)
......@@ -62,6 +65,10 @@ namespace Dec
bool operator>=(Iterator const& other) const { return index_ >= other.index_; }
bool operator<=(Iterator const& other) const { return index_ <= other.index_; }
std::shared_ptr<Entity> operator->() const
{
return std::make_shared<Entity>(gv_.indexSet(), gv_.coordinates(), halfEdge());
}
Entity operator*() const { return {gv_.indexSet(), gv_.coordinates(), halfEdge()}; }
Entity operator[](difference_type idx)
{
......@@ -76,6 +83,12 @@ namespace Dec
HalfEdgeGridView gv_;
IndexType index_;
};
template <Dune::PartitionIteratorType /*pit*/>
struct Partition
{
using Iterator = Codim::Iterator;
};
};
......@@ -115,14 +128,14 @@ namespace Dec
}
/// Obtain begin iterator for this view.
template <int cd>
template <int cd, Dune::PartitionIteratorType pitype = Dune::All_Partition>
typename Codim<cd>::Iterator begin() const
{
return {*this, 0u};
}
/// Obtain end iterator for this view.
template <int cd>
template <int cd, Dune::PartitionIteratorType pitype = Dune::All_Partition>
typename Codim<cd>::Iterator end() const
{
return {*this, size(cd)};
......@@ -130,6 +143,8 @@ namespace Dec
int level() const { return level_; }
auto const& comm() const { return grid_.comm(); }
public: // access to geometry information
/// Return the Geomery-Cache for the current level
......@@ -143,7 +158,7 @@ namespace Dec
/// Returns the (circum-)center of the `entity`.
template <class Entity>
GlobalCoords const& center(Entity const& entity) const
GlobalCoordinate const& center(Entity const& entity) const
{
return geometry().template center<Entity::codimension>(entity.index());
}
......
dune/dec/linear_algebra/block/BlockMatrix.hpp
View file @ 210d43c1
......@@ -13,12 +13,12 @@ namespace Dec
{
template <class T, std::size_t N, std::size_t M>
class BlockMatrix
: public FixMat<DOFMatrix<T>, N, M> // std::array<std::array<DOFMatrix<T>, M>, N>
: public FixMat<DOFMatrix<T>, N, M>
{
using Super = FixMat<DOFMatrix<T>, N, M>; //std::array<std::array<DOFMatrix<T>, M>, N>;
using Super = FixMat<DOFMatrix<T>, N, M>;
public:
using value_type = T;
using size_type = std::size_t;
using Matrix = DOFMatrix<T>;
......@@ -29,10 +29,6 @@ namespace Dec
Expects( n == N && m == M );
}
/// Access the (i,j)'th block
// DOFMatrix<T>& operator()(size_type i, size_type j) { return (*this)[i][j]; }
// DOFMatrix<T> const& operator()(size_type i, size_type j) const { return (*this)[i][j]; }
/// Return number of non-zeros
std::size_t nonZeros() const
{
......
dune/dec/linear_algebra/block/BlockVector.hpp
View file @ 210d43c1
......@@ -39,9 +39,9 @@ namespace Dec
template <class T, std::size_t N>
class BlockVector
: public FixVec<DOFVector<T>, N> //std::array<DOFVector<T>, N>
: public FixVec<DOFVector<T>, N>
{
using Super = FixVec<DOFVector<T>, N>; //std::array<DOFVector<T>, N>;
using Super = FixVec<DOFVector<T>, N>;
public:
using value_type = T;
......@@ -51,7 +51,6 @@ namespace Dec
BlockVector(std::vector<std::size_t> sizes)
: mapper_(sizes)
{
msg("constructor(sizes)");
Expects( N == sizes.size() );
for (std::size_t i = 0; i < N; ++i)
(*this)[i].resize(sizes[i]);
......@@ -61,7 +60,6 @@ namespace Dec
REQUIRES( concepts::GridView<GV> )>
BlockVector(GV const& gv, std::vector<int> Ks)
{
msg("constructor(gv,sizes)");
Expects( N == Ks.size() );
std::vector<std::size_t> sizes(N);
for (std::size_t i = 0; i < N; ++i) {
......@@ -77,7 +75,6 @@ namespace Dec
BlockVector(Splitter<Vector_> const& splitter)
: mapper_(splitter.mapper())
{
msg("constructor(splitter)");
auto const& input = splitter.vector();
for (std::size_t i = 0; i < N; ++i) {
(*this)[i].resize(mapper_.size(i));
......@@ -90,7 +87,6 @@ namespace Dec
: mapper_(that.mapper_)
, vector_(mapper_.size())
{
msg("constructor(that,ressource)");
for (std::size_t i = 0; i < N; ++i)
(*this)[i].resize(mapper_.size(i));
}
......@@ -98,7 +94,6 @@ namespace Dec
BlockVector& operator=(DOFVector<T> const& that)
{
msg("operator=(DOFVector)");
Expects( that.size() == size_type(mapper_.size()) );
for (std::size_t i = 0; i < N; ++i)
(*this)[i] = that.segment(mapper_.shift(i),mapper_.size(i));
......@@ -109,7 +104,6 @@ namespace Dec
template <class Vector_>
BlockVector& operator=(Splitter<Vector_> const& splitter)
{
msg("operator=(splitter)");
mapper_ = splitter.mapper();
auto const& input = splitter.vector();
for (std::size_t i = 0; i < N; ++i) {
......@@ -122,48 +116,6 @@ namespace Dec
using Super::operator=;
// BlockVector& operator+=(BlockVector const& that)
// {
// for (std::size_t i = 0; i < N; ++i)
// (*this)[i] += that[i];
// return *this;
// }
//
// BlockVector& operator-=(BlockVector const& that)
// {
// for (std::size_t i = 0; i < N; ++i)
// (*this)[i] -= that[i];
// return *this;
// }
//
// BlockVector& operator*=(float_type factor)
// {
// for (std::size_t i = 0; i < N; ++i)
// (*this)[i] *= factor;
// return *this;
// }
//
// template <class M>
// BlockVector& operator=(LinearAlgebraExpression<M> const& expr)
// {
// expr.assign(*this);
// return *this;
// }
//
// template <class M>
// BlockVector& operator+=(LinearAlgebraExpression<M> const& expr)
// {
// expr.add_assign(*this);
// return *this;
// }
//
// template <class M>
// BlockVector& operator-=(LinearAlgebraExpression<M> const& expr)
// {
// expr.add_assign(*this, -1.0);
// return *this;
// }
/// Sets all vectors to zero
void setZero()
{
......
dune/dec/linear_algebra/block/Operations.hpp
View file @ 210d43c1
......@@ -18,32 +18,6 @@ namespace Dec
return vec;
}
// template <class T, std::size_t N>
// T unary_dot(BlockVector<T,N> const& a)
// {
// T result{};
// for (std::size_t i = 0; i < N; ++i)
// result += unary_dot(a[i]);
// return result;
// }
//
// template <class T, std::size_t N>
// inline T two_norm(BlockVector<T,N> const& vec)
// {
// using std::sqrt;
// return sqrt(unary_dot(vec));
// }
// template <class T, std::size_t N>
// T dot(BlockVector<T,N> const& a, BlockVector<T,N> const& b)
// {
// T result{};
// for (std::size_t i = 0; i < N; ++i)
// result += dot(a[i], b[i]);
// return result;
// }
template <class T, std::size_t N, std::size_t M>
BlockVector<T, N> operator*(BlockMatrix<T, N, M> const& mat, BlockVector<T, M> const& vec)
{
......@@ -69,91 +43,4 @@ namespace Dec
return two_norm(r);
}
// //< lhs + rhs
// template <class T, std::size_t N>
// inline VecAddExpression<T,N> operator+(BlockVector<T,N> const& lhs, BlockVector<T,N> const& rhs)
// {
// return {lhs, rhs.self()};
// }
//
// //< lhs - rhs
// template <class T, std::size_t N>
// inline VecAddExpression<T,N> operator-(BlockVector<T,N> lhs, BlockVector<T,N> const& rhs)
// {
// return {lhs, rhs.self(), -1};
// }
//
//
// //< lhs + rhs
// template <class T, std::size_t N, class M>
// inline BlockVector<T,N> operator+(BlockVector<T,N> lhs, LinearAlgebraExpression<M> const& rhs)
// {
// return lhs += rhs;
// }
//
// //< lhs - rhs
// template <class T, std::size_t N, class M>
// inline BlockVector<T,N> operator-(BlockVector<T,N> lhs, LinearAlgebraExpression<M> const& rhs)
// {
// return lhs -= rhs;
// }
//
//
// //< lhs + rhs
// template <class T, std::size_t N, class M>
// inline BlockVector<T,N> operator+(VecScaleExpression<T,N> const& lhs, LinearAlgebraExpression<M> const& rhs)
// {
// BlockVector<T,N> result{lhs.vec_, tag::ressource{}}; result = lhs;
// return result += rhs;
// }
//