Skip to content
Snippets Groups Projects
Commit ef59d5d4 authored by Dr. Felix Tobias Schindler's avatar Dr. Felix Tobias Schindler
Browse files

[spaces.cg] implement new interface

parent e4c75df0
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/gdt/spaces/cg/default.hh
+ 55
− 92
View file @ ef59d5d4
......@@ -15,14 +15,9 @@
#include <dune/common/typetraits.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>
#include <dune/grid/common/capabilities.hh>
#include <dune/grid/common/rangegenerators.hh>
#include <dune/xt/common/exceptions.hh>
#include <dune/xt/common/numeric_cast.hh>
#include <dune/gdt/local/finite-elements/lagrange.hh>
#include <dune/gdt/spaces/basis/default.hh>
......@@ -33,40 +28,6 @@ namespace Dune {
namespace GDT {
template <class GL, int p, class R = double>
class ContinuousLagrangeSpace;
namespace internal {
template <class GL, int p, class R>
class ContinuousLagrangeSpaceTraits
{
static_assert(XT::Grid::is_layer<GL>::value, "");
static_assert(1 <= p && p <= 2, "Not implemented yet!");
using G = XT::Grid::extract_grid_t<GL>;
public:
using derived_type = ContinuousLagrangeSpace<GL, p, R>;
static const constexpr int polOrder = p;
static const constexpr size_t dimDomain = GL::dimension;
static const constexpr size_t dimRange = 1;
static const constexpr size_t dimRangeCols = 1;
static const constexpr bool continuous = false;
using GridLayerType = GL;
using RangeFieldType = R;
using BackendType = double;
static const constexpr XT::Grid::Backends layer_backend = XT::Grid::Backends::view;
static const constexpr Backends backend_type{Backends::gdt};
typedef DofCommunicationChooser<GridLayerType> DofCommunicationChooserType;
typedef typename DofCommunicationChooserType::Type DofCommunicatorType;
}; // class ContinuousLagrangeSpaceTraits
} // namespace internal
/**
* The following dimensions/orders/elements are tested to work:
*
......@@ -81,44 +42,33 @@ public:
*
* \sa make_lagrange_local_finite_element
*/
template <class GL, int p, class R>
class ContinuousLagrangeSpace
: public SpaceInterface<internal::ContinuousLagrangeSpaceTraits<GL, p, R>, GL::dimension, 1>
template <class GV, int p, class R = double>
class ContinuousLagrangeSpace : public SpaceInterface<GV, 1, 1, R>
{
public:
using Traits = internal::ContinuousLagrangeSpaceTraits<GL, p, R>;
private:
using BaseType = SpaceInterface<Traits, GL::dimension, 1>;
using ThisType = ContinuousLagrangeSpace<GL, p, R>;
using D = typename GL::ctype;
static const constexpr size_t d = BaseType::dimDomain;
using DofCommunicationChooserType = typename Traits::DofCommunicationChooserType;
using ThisType = ContinuousLagrangeSpace<GV, p, R>;
using BaseType = SpaceInterface<GV, 1, 1, R>;
public:
using typename BaseType::GridLayerType;
using typename BaseType::EntityType;
using typename BaseType::MapperType;
using typename BaseType::D;
using BaseType::d;
using typename BaseType::GridViewType;
using typename BaseType::GlobalBasisType;
using typename BaseType::MapperType;
using typename BaseType::FiniteElementType;
using DomainType = FieldVector<D, d>;
using DofCommunicatorType = typename Traits::DofCommunicatorType;
private:
using MapperImplementation = ContinuousMapper<GridLayerType, FiniteElementType>;
using GlobalBasisImplementation = DefaultGlobalBasis<GridLayerType, 1, 1, R>;
using MapperImplementation = ContinuousMapper<GridViewType, FiniteElementType>;
using GlobalBasisImplementation = DefaultGlobalBasis<GridViewType, 1, 1, R>;
public:
ContinuousLagrangeSpace(GridLayerType grd_lr)
: grid_layer_(grd_lr)
, communicator_(DofCommunicationChooserType::create(grid_layer_))
, backend_(0)
ContinuousLagrangeSpace(GridViewType grd_vw)
: grid_view_(grd_vw)
, finite_elements_(new std::map<GeometryType, std::shared_ptr<FiniteElementType>>())
, mapper_(nullptr)
, basis_(nullptr)
{
// create finite elements
for (auto&& geometry_type : grid_layer_.indexSet().types(0))
for (auto&& geometry_type : grid_view_.indexSet().types(0))
finite_elements_->insert(
std::make_pair(geometry_type, make_lagrange_local_finite_element<D, d, R>(geometry_type, p)));
// check
......@@ -126,9 +76,10 @@ public:
DUNE_THROW(Exceptions::space_error,
"when creating a ContinuousLagrangeSpace: non-conforming intersections are not (yet) "
"supported, and more than one element type in 3d leads to non-conforming intersections!");
// create mapper and basis
mapper_ = std::make_shared<MapperImplementation>(grid_layer_, finite_elements_);
basis_ = std::make_shared<GlobalBasisImplementation>(grid_layer_, finite_elements_);
// create mapper, basis and communicator
mapper_ = std::make_shared<MapperImplementation>(grid_view_, finite_elements_);
basis_ = std::make_shared<GlobalBasisImplementation>(grid_view_, finite_elements_);
this->create_communicator();
} // ContinuousLagrangeSpace(...)
ContinuousLagrangeSpace(const ThisType&) = default;
......@@ -137,39 +88,22 @@ public:
ThisType& operator=(const ThisType&) = delete;
ThisType& operator=(ThisType&&) = delete;
const GridLayerType& grid_layer() const
{
return grid_layer_;
}
const double& backend() const
const GridViewType& grid_view() const override final
{
return backend_;
return grid_view_;
}
const MapperType& mapper() const
const MapperType& mapper() const override final
{
return *mapper_;
}
const GlobalBasisType& basis() const
const GlobalBasisType& basis() const override final
{
return *basis_;
}
DofCommunicatorType& dof_communicator() const
{
DofCommunicationChooserType::prepare(*this, *communicator_);
return *communicator_;
}
const std::vector<DomainType>& lagrange_points(const EntityType& entity) const
{
return get_finite_element(entity.geometry().type()).lagrange_points();
}
private:
const FiniteElementType& get_finite_element(const GeometryType& geometry_type) const
const FiniteElementType& finite_element(const GeometryType& geometry_type) const override final
{
const auto finite_element_search_result = finite_elements_->find(geometry_type);
if (finite_element_search_result == finite_elements_->end())
......@@ -178,11 +112,40 @@ private:
<< "\n entity.geometry().type() = "
<< geometry_type);
return *finite_element_search_result->second;
} // ... get_finite_element(...)
}
SpaceType type() const override final
{
return SpaceType::continuous_lagrange;
}
int min_polorder() const override final
{
return p;
}
const GridLayerType grid_layer_;
mutable std::shared_ptr<DofCommunicatorType> communicator_;
const double backend_;
int max_polorder() const override final
{
return p;
}
bool continuous(const int diff_order) const override final
{
return diff_order == 0;
}
bool continuous_normal_components() const override final
{
return true;
}
bool is_lagrangian() const override final
{
return true;
}
private:
const GridViewType grid_view_;
std::shared_ptr<std::map<GeometryType, std::shared_ptr<FiniteElementType>>> finite_elements_;
std::shared_ptr<MapperImplementation> mapper_;
std::shared_ptr<GlobalBasisImplementation> basis_;
......
dune/gdt/test/spaces/spaces_cg.cc
+ 78
− 31
View file @ ef59d5d4
......@@ -23,17 +23,17 @@
#include <dune/xt/common/numeric_cast.hh>
#include <dune/xt/grid/gridprovider/cube.hh>
#include <dune/gdt/spaces/cg.hh>
#include <dune/gdt/spaces/cg/default.hh>
template <class GridLayerType, int p>
template <class GridViewType, int p>
struct ContinuousLagrangeSpace : public ::testing::Test
{
using SpaceType = Dune::GDT::ContinuousLagrangeSpace<GridLayerType, p>;
using D = typename SpaceType::DomainFieldType;
static const constexpr size_t d = SpaceType::dimDomain;
using SpaceType = Dune::GDT::ContinuousLagrangeSpace<GridViewType, p>;
using D = typename SpaceType::D;
static const constexpr size_t d = SpaceType::d;
virtual std::shared_ptr<GridLayerType> grid_layer() = 0;
virtual std::shared_ptr<GridViewType> grid_view() = 0;
std::shared_ptr<SpaceType> space;
......@@ -41,8 +41,8 @@ struct ContinuousLagrangeSpace : public ::testing::Test
void SetUp() override final
{
ASSERT_NE(grid_layer(), nullptr);
space = std::shared_ptr<SpaceType>(new SpaceType(*grid_layer()));
ASSERT_NE(grid_view(), nullptr);
space = std::shared_ptr<SpaceType>(new SpaceType(*grid_view()));
}
void TearDown() override final
......@@ -50,52 +50,66 @@ struct ContinuousLagrangeSpace : public ::testing::Test
space.reset();
}
void gives_correct_identification()
{
ASSERT_NE(grid_view(), nullptr);
ASSERT_NE(space, nullptr);
EXPECT_EQ(Dune::GDT::SpaceType::continuous_lagrange, space->type());
EXPECT_EQ(p, space->min_polorder());
EXPECT_EQ(p, space->max_polorder());
EXPECT_TRUE(space->continuous(0));
for (int diff_order : {1, 2, 3, 4, 5})
EXPECT_FALSE(space->continuous(diff_order));
EXPECT_TRUE(space->continuous_normal_components());
EXPECT_TRUE(space->is_lagrangian());
}
void basis_exists_on_each_element_with_correct_size()
{
ASSERT_NE(grid_layer(), nullptr);
ASSERT_NE(grid_view(), nullptr);
ASSERT_NE(space, nullptr);
for (auto&& element : elements(*grid_layer()))
for (auto&& element : elements(*grid_view()))
EXPECT_EQ(Dune::numLagrangePoints(element.geometry().type().id(), d, p),
space->basis().localize(element)->size());
}
void basis_exists_on_each_element_with_correct_order()
{
ASSERT_NE(grid_layer(), nullptr);
ASSERT_NE(grid_view(), nullptr);
ASSERT_NE(space, nullptr);
for (auto&& element : elements(*grid_layer()))
for (auto&& element : elements(*grid_view()))
EXPECT_EQ(p, space->basis().localize(element)->order());
}
void mapper_reports_correct_num_DoFs_on_each_element()
{
ASSERT_NE(grid_layer(), nullptr);
ASSERT_NE(grid_view(), nullptr);
ASSERT_NE(space, nullptr);
for (auto&& element : elements(*grid_layer()))
for (auto&& element : elements(*grid_view()))
EXPECT_EQ(Dune::numLagrangePoints(element.geometry().type().id(), d, p), space->mapper().local_size(element));
}
void mapper_reports_correct_max_num_DoFs()
{
ASSERT_NE(grid_layer(), nullptr);
ASSERT_NE(grid_view(), nullptr);
ASSERT_NE(space, nullptr);
size_t max_num_dofs = 0;
for (auto&& element : elements(*grid_layer()))
for (auto&& element : elements(*grid_view()))
max_num_dofs = std::max(max_num_dofs, space->mapper().local_size(element));
EXPECT_LE(max_num_dofs, space->mapper().max_local_size());
}
void mapper_maps_correctly()
{
ASSERT_NE(grid_layer(), nullptr);
ASSERT_NE(grid_view(), nullptr);
ASSERT_NE(space, nullptr);
// collect all global ids that are associated with a global lagrange point
std::map<Dune::FieldVector<D, d>, std::set<size_t>, Dune::XT::Common::FieldVectorLess>
global_lagrange_point_to_global_indices_map;
for (auto&& element : elements(*grid_layer())) {
for (auto&& element : elements(*grid_view())) {
const auto global_indices = space->mapper().global_indices(element);
EXPECT_LE(space->mapper().local_size(element), global_indices.size());
const auto lagrange_points = space->lagrange_points(element);
const auto lagrange_points = space->finite_element(element.geometry().type()).lagrange_points();
EXPECT_EQ(lagrange_points.size(), space->mapper().local_size(element));
for (size_t ii = 0; ii < lagrange_points.size(); ++ii) {
const auto global_lagrange_point = element.geometry().global(lagrange_points[ii]);
......@@ -123,19 +137,20 @@ struct ContinuousLagrangeSpace : public ::testing::Test
void lagrange_points_exist_on_each_element_with_correct_size()
{
ASSERT_NE(grid_layer(), nullptr);
ASSERT_NE(grid_view(), nullptr);
ASSERT_NE(space, nullptr);
for (auto&& element : elements(*grid_layer()))
EXPECT_EQ(Dune::numLagrangePoints(element.geometry().type().id(), d, p), space->lagrange_points(element).size());
for (auto&& element : elements(*grid_view()))
EXPECT_EQ(Dune::numLagrangePoints(element.geometry().type().id(), d, p),
space->finite_element(element.geometry().type()).lagrange_points().size());
}
void basis_is_lagrange_basis()
{
ASSERT_NE(grid_layer(), nullptr);
ASSERT_NE(grid_view(), nullptr);
ASSERT_NE(space, nullptr);
for (auto&& element : elements(*grid_layer())) {
for (auto&& element : elements(*grid_view())) {
const auto basis = space->basis().localize(element);
const auto lagrange_points = space->lagrange_points(element);
const auto lagrange_points = space->finite_element(element.geometry().type()).lagrange_points();
EXPECT_EQ(lagrange_points.size(), basis->size());
for (size_t ii = 0; ii < lagrange_points.size(); ++ii) {
const auto values = basis->evaluate(lagrange_points[ii]);
......@@ -150,9 +165,9 @@ struct ContinuousLagrangeSpace : public ::testing::Test
void basis_jacobians_seem_to_be_correct()
{
ASSERT_NE(grid_layer(), nullptr);
ASSERT_NE(grid_view(), nullptr);
ASSERT_NE(space, nullptr);
for (auto&& element : elements(*grid_layer())) {
for (auto&& element : elements(*grid_view())) {
const auto& reference_element = Dune::ReferenceElements<D, d>::general(element.geometry().type());
const auto basis = space->basis().localize(element);
const double h = 1e-6;
......@@ -222,7 +237,7 @@ struct ContinuousLagrangeSpaceOnSimplicialLeafView
~ContinuousLagrangeSpaceOnSimplicialLeafView() = default;
std::shared_ptr<LeafGridViewType> grid_layer() override final
std::shared_ptr<LeafGridViewType> grid_view() override final
{
return leaf_view;
}
......@@ -255,6 +270,10 @@ using SimplicialGrids = ::testing::Types<ONED_1D,
template <class G>
using Order1SimplicialContinuousLagrangeSpace = ContinuousLagrangeSpaceOnSimplicialLeafView<G, 1>;
TYPED_TEST_CASE(Order1SimplicialContinuousLagrangeSpace, SimplicialGrids);
TYPED_TEST(Order1SimplicialContinuousLagrangeSpace, gives_correct_identification)
{
this->gives_correct_identification();
}
TYPED_TEST(Order1SimplicialContinuousLagrangeSpace, basis_exists_on_each_element_with_correct_size)
{
this->basis_exists_on_each_element_with_correct_size();
......@@ -292,6 +311,10 @@ TYPED_TEST(Order1SimplicialContinuousLagrangeSpace, basis_jacobians_seem_to_be_c
template <class G>
using Order2SimplicialContinuousLagrangeSpace = ContinuousLagrangeSpaceOnSimplicialLeafView<G, 2>;
TYPED_TEST_CASE(Order2SimplicialContinuousLagrangeSpace, SimplicialGrids);
TYPED_TEST(Order2SimplicialContinuousLagrangeSpace, gives_correct_identification)
{
this->gives_correct_identification();
}
TYPED_TEST(Order2SimplicialContinuousLagrangeSpace, basis_exists_on_each_element_with_correct_size)
{
this->basis_exists_on_each_element_with_correct_size();
......@@ -351,7 +374,7 @@ struct ContinuousLagrangeSpaceOnCubicLeafView
~ContinuousLagrangeSpaceOnCubicLeafView() = default;
std::shared_ptr<LeafGridViewType> grid_layer() override final
std::shared_ptr<LeafGridViewType> grid_view() override final
{
return leaf_view;
}
......@@ -382,6 +405,10 @@ using CubicGrids = ::testing::Types<YASP_2D_EQUIDISTANT_OFFSET
template <class G>
using Order1CubicContinuousLagrangeSpace = ContinuousLagrangeSpaceOnCubicLeafView<G, 1>;
TYPED_TEST_CASE(Order1CubicContinuousLagrangeSpace, CubicGrids);
TYPED_TEST(Order1CubicContinuousLagrangeSpace, gives_correct_identification)
{
this->gives_correct_identification();
}
TYPED_TEST(Order1CubicContinuousLagrangeSpace, basis_exists_on_each_element_with_correct_size)
{
this->basis_exists_on_each_element_with_correct_size();
......@@ -419,6 +446,10 @@ TYPED_TEST(Order1CubicContinuousLagrangeSpace, basis_jacobians_seem_to_be_correc
template <class G>
using Order2CubicContinuousLagrangeSpace = ContinuousLagrangeSpaceOnCubicLeafView<G, 2>;
TYPED_TEST_CASE(Order2CubicContinuousLagrangeSpace, CubicGrids);
TYPED_TEST(Order2CubicContinuousLagrangeSpace, gives_correct_identification)
{
this->gives_correct_identification();
}
TYPED_TEST(Order2CubicContinuousLagrangeSpace, basis_exists_on_each_element_with_correct_size)
{
this->basis_exists_on_each_element_with_correct_size();
......@@ -492,7 +523,7 @@ struct ContinuousLagrangeSpaceOnPrismLeafView
~ContinuousLagrangeSpaceOnPrismLeafView() = default;
std::shared_ptr<LeafGridViewType> grid_layer() override final
std::shared_ptr<LeafGridViewType> grid_view() override final
{
return leaf_view;
}
......@@ -509,6 +540,10 @@ using PrismGrids = ::testing::Types<
template <class G>
using Order1PrismContinuousLagrangeSpace = ContinuousLagrangeSpaceOnPrismLeafView<G, 1>;
TYPED_TEST_CASE(Order1PrismContinuousLagrangeSpace, PrismGrids);
TYPED_TEST(Order1PrismContinuousLagrangeSpace, gives_correct_identification)
{
this->gives_correct_identification();
}
TYPED_TEST(Order1PrismContinuousLagrangeSpace, basis_exists_on_each_element_with_correct_size)
{
this->basis_exists_on_each_element_with_correct_size();
......@@ -546,6 +581,10 @@ TYPED_TEST(Order1PrismContinuousLagrangeSpace, basis_jacobians_seem_to_be_correc
template <class G>
using Order2PrismContinuousLagrangeSpace = ContinuousLagrangeSpaceOnPrismLeafView<G, 2>;
TYPED_TEST_CASE(Order2PrismContinuousLagrangeSpace, PrismGrids);
TYPED_TEST(Order2PrismContinuousLagrangeSpace, gives_correct_identification)
{
this->gives_correct_identification();
}
TYPED_TEST(Order2PrismContinuousLagrangeSpace, basis_exists_on_each_element_with_correct_size)
{
this->basis_exists_on_each_element_with_correct_size();
......@@ -661,7 +700,7 @@ struct ContinuousLagrangeSpaceOnMixedLeafView
~ContinuousLagrangeSpaceOnMixedLeafView() = default;
std::shared_ptr<LeafGridViewType> grid_layer() override final
std::shared_ptr<LeafGridViewType> grid_view() override final
{
return leaf_view;
}
......@@ -679,6 +718,10 @@ using MixedGrids = ::testing::Types<
template <class G>
using Order1MixedContinuousLagrangeSpace = ContinuousLagrangeSpaceOnMixedLeafView<G, 1>;
TYPED_TEST_CASE(Order1MixedContinuousLagrangeSpace, MixedGrids);
TYPED_TEST(Order1MixedContinuousLagrangeSpace, gives_correct_identification)
{
this->gives_correct_identification();
}
TYPED_TEST(Order1MixedContinuousLagrangeSpace, basis_exists_on_each_element_with_correct_size)
{
this->basis_exists_on_each_element_with_correct_size();
......@@ -716,6 +759,10 @@ TYPED_TEST(Order1MixedContinuousLagrangeSpace, basis_jacobians_seem_to_be_correc
template <class G>
using Order2MixedContinuousLagrangeSpace = ContinuousLagrangeSpaceOnMixedLeafView<G, 2>;
TYPED_TEST_CASE(Order2MixedContinuousLagrangeSpace, MixedGrids);
TYPED_TEST(Order2MixedContinuousLagrangeSpace, gives_correct_identification)
{
this->gives_correct_identification();
}
TYPED_TEST(Order2MixedContinuousLagrangeSpace, basis_exists_on_each_element_with_correct_size)
{
this->basis_exists_on_each_element_with_correct_size();
......
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