Commit 686fbf8d authored by Oliver Sander's avatar Oliver Sander

Merge branch 'improve-raviart-thomas-elements' into 'master'

Improve Raviart-Thomas-Basis

See merge request staging/dune-functions!270
parents 75fd4b6d f5ad40e2
Pipeline #30184 failed with stage
in 27 minutes and 34 seconds
......@@ -15,6 +15,8 @@ corresponding version of the Dune core modules.
transformation. Domain-space transformations still have to be done
by the calling code. The `GlobalValuedLocalFiniteElement` still
implements the `LocalFiniteElement` interface of `dune-localfunctions`.
- The `RaviartThomasBasis` class now supports tetrahedral grids for `order=0`,
quadrilateral grids for `order=2`, and hexahedral grids for `order=1`.
## Release 2.7
......
......@@ -15,6 +15,7 @@
#include <dune/localfunctions/raviartthomas/raviartthomas0cube2d.hh>
#include <dune/localfunctions/raviartthomas/raviartthomas0cube3d.hh>
#include <dune/localfunctions/raviartthomas/raviartthomas02d.hh>
#include <dune/localfunctions/raviartthomas/raviartthomas03d.hh>
#include <dune/localfunctions/raviartthomas/raviartthomas1cube2d.hh>
#include <dune/localfunctions/raviartthomas/raviartthomas1cube3d.hh>
#include <dune/localfunctions/raviartthomas/raviartthomas12d.hh>
......@@ -33,62 +34,63 @@ namespace Impl {
template<int dim, typename D, typename R, std::size_t k>
struct RaviartThomasSimplexLocalInfo
{
static_assert((AlwaysFalse<D>::value),"The requested type of Raviart-Thomas element is not implemented, sorry!");
// Dummy type, must be something that we can have a std::unique_ptr to
using FiniteElement = void*;
};
template<typename D, typename R>
struct RaviartThomasSimplexLocalInfo<2,D,R,0>
{
using FiniteElement = RT02DLocalFiniteElement<D,R>;
static const std::size_t Variants = 8;
};
template<typename D, typename R>
struct RaviartThomasSimplexLocalInfo<2,D,R,1>
{
using FiniteElement = RT12DLocalFiniteElement<D,R>;
static const std::size_t Variants = 8;
};
template<typename D, typename R>
struct RaviartThomasSimplexLocalInfo<3,D,R,0>
{
using FiniteElement = RT03DLocalFiniteElement<D,R>;
};
template<int dim, typename D, typename R, std::size_t k>
struct RaviartThomasCubeLocalInfo
{
static_assert((AlwaysFalse<D>::value),"The requested type of Raviart-Thomas element is not implemented, sorry!");
// Dummy type, must be something that we can have a std::unique_ptr to
using FiniteElement = void*;
};
template<typename D, typename R>
struct RaviartThomasCubeLocalInfo<2,D,R,0>
{
using FiniteElement = RT0Cube2DLocalFiniteElement<D,R>;
static const std::size_t Variants = 16;
};
template<typename D, typename R>
struct RaviartThomasCubeLocalInfo<2,D,R,1>
{
using FiniteElement = RT1Cube2DLocalFiniteElement<D,R>;
static const std::size_t Variants = 16;
};
template<typename D, typename R>
struct RaviartThomasCubeLocalInfo<2,D,R,2>
{
using FiniteElement = RT2Cube2DLocalFiniteElement<D,R>;
static const std::size_t Variants = 16;
};
template<typename D, typename R>
struct RaviartThomasCubeLocalInfo<3,D,R,0>
{
using FiniteElement = RT0Cube3DLocalFiniteElement<D,R>;
static const std::size_t Variants = 64;
};
template<typename D, typename R>
struct RaviartThomasCubeLocalInfo<3,D,R,1>
{
using FiniteElement = RT1Cube3DLocalFiniteElement<D,R>;
static const std::size_t Variants = 64;
};
template<typename GV, int dim, typename R, std::size_t k>
......@@ -117,18 +119,34 @@ namespace Impl {
typename std::conditional<type.isCube(),CubeFiniteElement,SimplexFiniteElement>::type,
LocalFiniteElementVirtualInterface<T> >::type;
static_assert(!std::is_same_v<FiniteElement,void*>,"The requested type of Raviart-Thomas element is not implemented, sorry!");
// Each element facet can have its orientation reversed, hence there are
// 2^#facets different variants.
static std::size_t numVariants(GeometryType type)
{
auto numFacets = referenceElement<D,dim>(type).size(1);
return power(2,numFacets);
}
RaviartThomasLocalFiniteElementMap(const GV& gv)
: is_(&(gv.indexSet())), orient_(gv.size(0))
{
cubeVariant_.resize(RaviartThomasCubeLocalInfo<dim, D, R, k>::Variants);
simplexVariant_.resize(RaviartThomasSimplexLocalInfo<dim, D, R, k>::Variants);
cubeVariant_.resize(numVariants(GeometryTypes::cube(dim)));
simplexVariant_.resize(numVariants(GeometryTypes::simplex(dim)));
// create all variants
for (size_t i = 0; i < cubeVariant_.size(); i++)
cubeVariant_[i] = std::make_unique<CubeFiniteElementImp>(CubeFiniteElement(i));
// create all variants, if they exist
if constexpr (!std::is_same_v<CubeFiniteElement,void*>)
{
for (size_t i = 0; i < cubeVariant_.size(); i++)
cubeVariant_[i] = std::make_unique<CubeFiniteElementImp>(CubeFiniteElement(i));
}
if constexpr (!std::is_same_v<SimplexFiniteElement,void*>)
{
for (size_t i = 0; i < simplexVariant_.size(); i++)
simplexVariant_[i] = std::make_unique<SimplexFiniteElementImp>(SimplexFiniteElement(i));
}
// compute orientation for all elements
// loop once over the grid
......@@ -248,8 +266,8 @@ public:
DUNE_THROW(Dune::NotImplemented, "Raviart-Thomas basis is only implemented for grids with a single element type");
GeometryType type = gv.template begin<0>()->type();
const static int dofsPerElement = (dim == 2) ? (type.isCube() ? k*(k+1)*dim : k*dim) : k*(k+1)*(k+1)*dim;
constexpr int dofsPerFace = (dim == 2) ? k+1 : 3*k+1;
const static int dofsPerElement = type.isCube() ? ((dim == 2) ? k*(k+1)*dim : k*(k+1)*(k+1)*dim) : k*dim;
const static int dofsPerFace = type.isCube() ? (dim-2)*2*k+k+1 : (dim-1)*k+1 ;
dofsPerCodim_ = {{dofsPerElement, dofsPerFace}};
}
......@@ -316,9 +334,9 @@ public:
for (auto&& type : gridView_.indexSet().types(0))
{
size_t numFaces = ReferenceElements<double,dim>::general(type).size(1);
const static int dofsPerCodim0 = (dim == 2) ? (type.isCube() ? k*(k+1)*dim : k*dim) : k*(k+1)*(k+1)*dim;
constexpr int dofsPerCodim1 = (dim == 2) ? k+1 : 3*k+1;
result = std::max(result, dofsPerCodim0 + dofsPerCodim1 * numFaces);
const static int dofsPerElement = type.isCube() ? ((dim == 2) ? k*(k+1)*dim : k*(k+1)*(k+1)*dim) : k*dim;
const static int dofsPerFace = type.isCube() ? (dim-2)*2*k+k+1 : (dim-1)*k+1 ;
result = std::max(result, dofsPerElement + dofsPerFace * numFaces);
}
return result;
......
......@@ -7,6 +7,7 @@
#include <dune/grid/uggrid.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/functions/functionspacebases/raviartthomasbasis.hh>
......@@ -14,48 +15,71 @@
using namespace Dune;
int main (int argc, char* argv[])
template<int k, class GridView>
void testRaviartThomasBasis(TestSuite& test, const GridView& gridView)
{
MPIHelper::instance(argc, argv);
TestSuite test;
std::cout<<" Testing order: "<< k <<std::endl;
// Generate grid for testing
const int dim = 2;
typedef YaspGrid<dim> GridType;
FieldVector<double,dim> l(1);
std::array<int,dim> elements = {{10, 10}};
GridType grid(l,elements);
// check RaviartThomasBasis created 'manually'
// Check RaviartThomasBasis created 'manually'
{
typedef GridType::LeafGridView GridView;
const GridView& gridView = grid.leafGridView();
Functions::RaviartThomasBasis<GridView,0> basis(gridView);
Functions::RaviartThomasBasis<GridView,k> basis(gridView);
test.subTest(checkBasis(basis, EnableNormalContinuityCheck()));
}
// check RaviartThomasBasis created using basis builder mechanism
// Check RaviartThomasBasis created using basis builder mechanism
{
using namespace Functions::BasisFactory;
auto basis = makeBasis(grid.leafGridView(), raviartThomas<0>());
auto basis = makeBasis(gridView, raviartThomas<k>());
test.subTest(checkBasis(basis, EnableNormalContinuityCheck()));
}
}
// check RaviartThomasBasis on a grid without a compile-time-fixed element type
{
using Grid = UGGrid<dim>;
std::shared_ptr<Grid> grid = StructuredGridFactory<Grid>::createCubeGrid({0.0,0.0}, l, {{10,10}});
Functions::RaviartThomasBasis<Grid::LeafGridView,0> basis(grid->leafGridView());
test.subTest(checkBasis(basis, EnableNormalContinuityCheck()));
}
int main (int argc, char* argv[])
{
MPIHelper::instance(argc, argv);
TestSuite test;
// Test with grid that only supports cube elements
// (Grids with only a single element type receive special treatment by RaviartThomasBasis)
std::cout<<"Testing RaviartThomasBasis in 2D with cube grids\n";
YaspGrid<2> quadGrid({1.0, 1.0}, {5,5});
auto quadGridView = quadGrid.leafGridView();
testRaviartThomasBasis<0>(test, quadGridView);
testRaviartThomasBasis<1>(test, quadGridView);
testRaviartThomasBasis<2>(test, quadGridView);
std::cout<<"Testing RaviartThomasBasis in 3D with cube grids\n";
YaspGrid<3> hexaGrid({1.0, 1.0, 1.0}, {4,4,4});
auto hexaGridView = hexaGrid.leafGridView();
testRaviartThomasBasis<0>(test, hexaGridView);
testRaviartThomasBasis<1>(test, hexaGridView);
// Test with pure simplex grid
// (Unfortunately there is no grid implementation available that only supports simplices.)
std::cout<<"Testing RaviartThomasBasis in 2D with simplex grid\n";
using Mixed2dGrid = UGGrid<2>;
const std::string path = std::string(DUNE_GRID_EXAMPLE_GRIDS_PATH) + "gmsh/";
auto triangleGrid = GmshReader<Mixed2dGrid>::read(path + "curved2d.msh");
auto triangleGridView = triangleGrid->leafGridView();
testRaviartThomasBasis<0>(test, triangleGridView);
testRaviartThomasBasis<1>(test, triangleGridView);
std::cout<<"Testing RaviartThomasBasis in 3D with simplex grid\n";
using Mixed3dGrid = UGGrid<3>;
auto tetraGrid = GmshReader<Mixed3dGrid>::read(path + "telescope.msh");
auto tetraGridView = tetraGrid->leafGridView();
testRaviartThomasBasis<0>(test, tetraGridView);
// Test with mixed-element 2d grid
// The mixed-element case needs additional fixes in RaviartThomasBasis.
// std::cout<<"Testing RaviartThomasBasis in 2D with mixed-element grid\n";
// auto mixed2dGrid = GmshReader<Mixed2dGrid>::read(path + "hybrid-testgrid-2d.msh");
// auto mixed2dGridView = mixed2dGrid->leafGridView();
// testRaviartThomasBasis<0>(test, mixed2dGridView);
// testRaviartThomasBasis<1>(test, mixed2dGridView);
return test.exit();
}
Markdown is supported
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