Commit 58acd3f3 authored by Martin Nolte's avatar Martin Nolte

Merge branch 'cherry-pick-058d0fc0' into 'releases/2.6'

Merge branch 'feature/yaspgrid-constexpr-tables' into 'master'

See merge request !245
parents a483815f 8c6daef7
Pipeline #5887 passed with stage
in 76 minutes and 49 seconds
......@@ -17,6 +17,8 @@ dune_add_test(NAME test-yaspgrid-backuprestore-tensor
TIMEOUT 666
)
dune_add_test(SOURCES test-yaspgrid-entityshifttable.cc)
dune_add_test(NAME test-yaspgrid-tensorgridfactory
SOURCES test-yaspgrid-tensorgridfactory.cc
MPI_RANKS 1 2
......
#include "config.h"
#include <iostream>
#include <dune/common/test/testsuite.hh>
#include <dune/grid/yaspgrid.hh>
using Dune::TestSuite;
template<int n>
TestSuite testBinomialTable(std::vector<int> const* reference = nullptr, bool dump = false)
{
TestSuite t;
using Table = Dune::Yasp::BinomialTable<n>;
if (!reference)
std::cerr << "W: No reference given!\n";
if (dump)
std::cout << " std::vector<int> reference{";
int i = 0;
for (int d = 0; d <= n; ++d) {
for (int c = 0; c <= d; ++c, ++i) {
const auto value = Table::evaluate(d, c);
t.check(value == Table::binomial(d, c));
if (reference)
t.check(Table::evaluate(d, c) == reference->at(i));
if (dump)
std::cout << value << ", ";
}
}
if (dump)
std::cout << "};\n";
return t;
}
template<class F, int dim>
TestSuite testEntityShiftTable(std::vector<unsigned long long> const* reference = nullptr, bool dump = false)
{
TestSuite t;
using Table = Dune::Yasp::EntityShiftTable<F, dim>;
if (!reference)
std::cerr << "W: No reference given!\n";
if (dump)
std::cout << " std::vector<unsigned long long> reference{";
int i = 0;
for (int codim = 0; codim <= dim; ++codim) {
for (int j = 0; j < Dune::Yasp::subEnt<dim>(dim, codim); ++j, ++i) {
const auto value = Table::evaluate(j, codim);
t.check(value == F::evaluate(j, codim));
if (reference)
t.check(value.to_ullong() == reference->at(i));
if (dump)
std::cout << "0b" << value.to_string() << "ull, ";
}
}
if (dump)
std::cout << "};\n";
return t;
}
int main(int argc, char** argv)
{
bool dump = false;
if (argc > 1)
dump = true;
TestSuite t;
{
std::vector<int> reference{1, 1, 1};
t.subTest(testBinomialTable<1>(&reference, dump));
}
{
std::vector<int> reference{1, 1, 1, 1, 2, 1};
t.subTest(testBinomialTable<2>(&reference, dump));
}
{
std::vector<int> reference{1, 1, 1, 1, 2, 1, 1, 3, 3, 1};
t.subTest(testBinomialTable<3>(&reference, dump));
}
{
std::vector<int> reference{1, 1, 1, 1, 2, 1, 1, 3, 3, 1, 1, 4, 6, 4, 1};
t.subTest(testBinomialTable<4>(&reference, dump));
}
{
std::vector<unsigned long long> reference{0b1ull, 0b0ull, 0b0ull};
t.subTest(testEntityShiftTable<Dune::Yasp::calculate_entity_shift<1>, 1>(&reference, dump));
}
{
std::vector<unsigned long long> reference{
0b11ull, 0b10ull, 0b10ull, 0b01ull, 0b01ull,
0b00ull, 0b00ull, 0b00ull, 0b00ull
};
t.subTest(testEntityShiftTable<Dune::Yasp::calculate_entity_shift<2>, 2>(&reference, dump));
}
{
std::vector<unsigned long long> reference{
0b111ull, 0b110ull, 0b110ull, 0b101ull, 0b101ull, 0b011ull, 0b011ull,
0b100ull, 0b100ull, 0b100ull, 0b100ull, 0b010ull, 0b010ull, 0b001ull,
0b001ull, 0b010ull, 0b010ull, 0b001ull, 0b001ull, 0b000ull, 0b000ull,
0b000ull, 0b000ull, 0b000ull, 0b000ull, 0b000ull, 0b000ull
};
t.subTest(testEntityShiftTable<Dune::Yasp::calculate_entity_shift<3>, 3>(&reference, dump));
}
{
std::vector<unsigned long long> reference{
0b1111ull, 0b1110ull, 0b1110ull, 0b1101ull, 0b1101ull, 0b1011ull,
0b1011ull, 0b0111ull, 0b0111ull, 0b1100ull, 0b1100ull, 0b1100ull,
0b1100ull, 0b1010ull, 0b1010ull, 0b1001ull, 0b1001ull, 0b1010ull,
0b1010ull, 0b1001ull, 0b1001ull, 0b0110ull, 0b0110ull, 0b0101ull,
0b0101ull, 0b0011ull, 0b0011ull, 0b0110ull, 0b0110ull, 0b0101ull,
0b0101ull, 0b0011ull, 0b0011ull, 0b1000ull, 0b1000ull, 0b1000ull,
0b1000ull, 0b1000ull, 0b1000ull, 0b1000ull, 0b1000ull, 0b0100ull,
0b0100ull, 0b0100ull, 0b0100ull, 0b0010ull, 0b0010ull, 0b0001ull,
0b0001ull, 0b0010ull, 0b0010ull, 0b0001ull, 0b0001ull, 0b0100ull,
0b0100ull, 0b0100ull, 0b0100ull, 0b0010ull, 0b0010ull, 0b0001ull,
0b0001ull, 0b0010ull, 0b0010ull, 0b0001ull, 0b0001ull, 0b0000ull,
0b0000ull, 0b0000ull, 0b0000ull, 0b0000ull, 0b0000ull, 0b0000ull,
0b0000ull, 0b0000ull, 0b0000ull, 0b0000ull, 0b0000ull, 0b0000ull,
0b0000ull, 0b0000ull, 0b0000ull
};
t.subTest(testEntityShiftTable<Dune::Yasp::calculate_entity_shift<4>, 4>(&reference, dump));
}
{
std::vector<unsigned long long> reference{0b0ull, 0b0ull, 0b1ull};
t.subTest(testEntityShiftTable<Dune::Yasp::calculate_entity_move<1>, 1>(&reference, dump));
}
{
std::vector<unsigned long long> reference{
0b00ull, 0b00ull, 0b01ull, 0b00ull, 0b10ull,
0b00ull, 0b01ull, 0b10ull, 0b11ull
};
t.subTest(testEntityShiftTable<Dune::Yasp::calculate_entity_move<2>, 2>(&reference, dump));
}
{
std::vector<unsigned long long> reference{
0b000ull, 0b000ull, 0b001ull, 0b000ull, 0b010ull, 0b000ull, 0b100ull,
0b000ull, 0b001ull, 0b010ull, 0b011ull, 0b000ull, 0b001ull, 0b000ull,
0b010ull, 0b100ull, 0b101ull, 0b100ull, 0b110ull, 0b000ull, 0b001ull,
0b010ull, 0b011ull, 0b100ull, 0b101ull, 0b110ull, 0b111ull
};
t.subTest(testEntityShiftTable<Dune::Yasp::calculate_entity_move<3>, 3>(&reference, dump));
}
{
std::vector<unsigned long long> reference{
0b0000ull, 0b0000ull, 0b0001ull, 0b0000ull, 0b0010ull, 0b0000ull,
0b0100ull, 0b0000ull, 0b1000ull, 0b0000ull, 0b0001ull, 0b0010ull,
0b0011ull, 0b0000ull, 0b0001ull, 0b0000ull, 0b0010ull, 0b0100ull,
0b0101ull, 0b0100ull, 0b0110ull, 0b0000ull, 0b0001ull, 0b0000ull,
0b0010ull, 0b0000ull, 0b0100ull, 0b1000ull, 0b1001ull, 0b1000ull,
0b1010ull, 0b1000ull, 0b1100ull, 0b0000ull, 0b0001ull, 0b0010ull,
0b0011ull, 0b0100ull, 0b0101ull, 0b0110ull, 0b0111ull, 0b0000ull,
0b0001ull, 0b0010ull, 0b0011ull, 0b0000ull, 0b0001ull, 0b0000ull,
0b0010ull, 0b0100ull, 0b0101ull, 0b0100ull, 0b0110ull, 0b1000ull,
0b1001ull, 0b1010ull, 0b1011ull, 0b1000ull, 0b1001ull, 0b1000ull,
0b1010ull, 0b1100ull, 0b1101ull, 0b1100ull, 0b1110ull, 0b0000ull,
0b0001ull, 0b0010ull, 0b0011ull, 0b0100ull, 0b0101ull, 0b0110ull,
0b0111ull, 0b1000ull, 0b1001ull, 0b1010ull, 0b1011ull, 0b1100ull,
0b1101ull, 0b1110ull, 0b1111ull
};
t.subTest(testEntityShiftTable<Dune::Yasp::calculate_entity_move<4>, 4>(&reference, dump));
}
return t.exit();
}
......@@ -678,11 +678,6 @@ namespace Dune {
void init()
{
Hybrid::forEach( std::make_integer_sequence<int,dim+1>(), [] ( auto d ) {
Yasp::BinomialTable<d>::init();
Yasp::EntityShiftTable<Yasp::calculate_entity_shift<d>,d>::init();
Yasp::EntityShiftTable<Yasp::calculate_entity_move<d>,d>::init();
} );
indexsets.push_back( std::make_shared< YaspIndexSet<const YaspGrid<dim, Coordinates>, false > >(*this,0) );
boundarysegmentssize();
}
......
......@@ -31,39 +31,15 @@ namespace Dune {
template<int n>
struct BinomialTable
{
static void init()
{
if (_initialized)
return;
int offset = 0;
for (int d = 0; d <= n; ++d)
{
_offsets[d] = offset;
for (int c = 0; c <= d; ++c, ++offset)
_values[offset] = binomial(d,c);
}
_initialized = true;
}
// evaluation - note that in general d!=n, n is only the
// maximum value of d (in our case dimworld)
static int evaluate(int d, int c)
static constexpr int evaluate(int d, int c)
{
return _values[_offsets[d] + c];
}
private:
// prevent construction
BinomialTable();
static bool _initialized;
static std::array<int,(n+1)*(n+2)/2> _values;
static std::array<int,n+1> _offsets;
public:
// the actual implementation
static int binomial(int d, int c)
static constexpr int binomial(int d, int c)
{
long binomial=1;
for (int i=d-c+1; i<=d; i++)
......@@ -72,14 +48,46 @@ namespace Dune {
binomial /= i;
return binomial;
}
private:
// prevent construction
BinomialTable() = delete;
// compute binomial(r, c) and advance row `r` and column `c`
static constexpr int nextValue(int& r, int& c)
{
const auto result = binomial(r, c);
c += 1;
if (c > r) {
r += 1;
c = 0;
}
return result;
}
template<std::size_t... I>
static constexpr std::array<int, sizeof...(I)> computeValues(std::index_sequence<I...>)
{
int r = 0, c = 0;
return { ((void)I, nextValue(r, c))... };
}
template<std::size_t... I>
static constexpr std::array<int, sizeof...(I)> computeOffsets(std::index_sequence<I...>)
{ return { (I*(I+1)/2)... }; }
static constexpr std::array<int,(n+1)*(n+2)/2> _values = computeValues(std::make_index_sequence<(n+1)*(n+2)/2>{});
static constexpr std::array<int,n+1> _offsets = computeOffsets(std::make_index_sequence<n+1>{});
};
#if __cplusplus < 201703L
template<int n>
bool BinomialTable<n>::_initialized = false;
template<int n>
std::array<int,(n+1)*(n+2)/2> BinomialTable<n>::_values;
constexpr std::array<int,(n+1)*(n+2)/2> BinomialTable<n>::_values;
template<int n>
std::array<int,n+1> BinomialTable<n>::_offsets;
constexpr std::array<int,n+1> BinomialTable<n>::_offsets;
#endif
/** \returns number of subentities of given codim in a cube of dimension dim
* \tparam dimworld the maximum dimension the table holds entries for
......@@ -88,7 +96,7 @@ namespace Dune {
* That number is d choose c times 2^c.
*/
template<int dimworld>
int subEnt(int d, int c)
constexpr int subEnt(int d, int c)
{
return (d < c ? 0 : BinomialTable<dimworld>::evaluate(d,c) << c);
}
......@@ -98,24 +106,8 @@ namespace Dune {
template<typename F, int dim>
struct EntityShiftTable
{
typedef std::bitset<dim> value_type;
static void init()
{
if (_initialized)
return;
F f;
int offset = 0;
for (int codim = 0; codim <= dim; ++codim)
{
_offsets[codim] = offset;
for (int i = 0; i < subEnt<dim>(dim,codim); ++i, ++offset)
_values[offset] = static_cast<unsigned char>(f(i,codim).to_ulong());
}
_initialized = true;
}
static value_type evaluate(int i, int codim)
{
return {_values[_offsets[codim] + i]};
......@@ -124,39 +116,73 @@ namespace Dune {
private:
// prevent construction
EntityShiftTable();
EntityShiftTable() = delete;
// compute offset of codimension `codim` entities and advance `offset`
static constexpr int nextOffset(int& offset, int codim)
{
if (codim == 0) {
offset = 0;
return 0;
}
static bool _initialized;
static std::array<int,dim+1> _offsets;
static std::array<unsigned char,StaticPower<3,dim>::power> _values;
offset += subEnt<dim>(dim, codim-1);
return offset;
}
template<std::size_t... I>
static constexpr std::array<int, sizeof...(I)> computeOffsets(std::index_sequence<I...>)
{
int offset = 0;
return { (nextOffset(offset, I))... };
}
// compute shift table entry for (`codim`, `i`) and advance `codim`, `i`
static constexpr unsigned char nextValue(int& codim, int& i)
{
const auto result = F::evaluate(i, codim);
i += 1;
if (i >= subEnt<dim>(dim, codim)) {
codim += 1;
i = 0;
}
return result;
}
template<std::size_t... I>
static constexpr std::array<unsigned char, sizeof...(I)> computeValues(std::index_sequence<I...>)
{
int codim = 0, i = 0;
return { ((void)I, nextValue(codim, i))... };
}
static constexpr std::array<int,dim+1> _offsets = computeOffsets(std::make_index_sequence<dim+1>{});
static constexpr std::array<unsigned char,StaticPower<3,dim>::power> _values = computeValues(std::make_index_sequence<StaticPower<3,dim>::power>{});
};
#if __cplusplus < 201703L
template<typename F, int dim>
bool EntityShiftTable<F,dim>::_initialized = false;
template<typename F, int dim>
std::array<int,dim+1> EntityShiftTable<F,dim>::_offsets;
constexpr std::array<int,dim+1> EntityShiftTable<F, dim>::_offsets;
template<typename F, int dim>
std::array<unsigned char,StaticPower<3,dim>::power> EntityShiftTable<F,dim>::_values;
constexpr std::array<unsigned char,StaticPower<3,dim>::power> EntityShiftTable<F, dim>::_values;
#endif
// functor for doing the actual entity shift calculation
template<int dim>
struct calculate_entity_shift
{
calculate_entity_shift()
{
BinomialTable<dim>::init();
}
std::bitset<dim> operator()(int index, int cc) const
static constexpr unsigned long long evaluate(int index, int cc)
{
std::bitset<dim> result(0ull);
auto result = 0ull;
for (int d = dim; d>0; d--)
{
if (cc == d)
return result;
if (index < subEnt<dim>(d-1,cc))
result[d-1]=true;
result |= 1ull << (d-1);
else
{
index = (index - subEnt<dim>(d-1, cc)) % subEnt<dim>(d-1,cc-1);
......@@ -185,27 +211,24 @@ namespace Dune {
template<int dim>
struct calculate_entity_move
{
calculate_entity_move()
static constexpr unsigned long long evaluate(int index, int cc)
{
BinomialTable<dim>::init();
}
std::bitset<dim> operator()(int index, int cc) const
{
std::bitset<dim> result(0ull);
auto result = 0ull;
for (int d = dim; d>0; d--)
{
if (d == cc)
{
result[d-1] = index & (1<<(d-1));
// result[d-1] = index & (1<<(d-1));
result &= ~(1ull << (d-1));
result |= index & (1ull << (d-1));
index &= ~(1<<(d-1));
}
if (index >= subEnt<dim>(d-1,cc))
{
if ((index - subEnt<dim>(d-1,cc)) / subEnt<dim>(d-1,cc-1) == 1)
{
result[d-1] = true;
result |= 1ull << (d-1);
}
index = (index - subEnt<dim>(d-1, cc)) % subEnt<dim>(d-1,cc-1);
cc--;
......
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