Commit 10aab891 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Merge branch 'issue/global_id_set_cleanup' into 'master'

added twist utilities to flip the order of edge DOFs

See merge request amdis/amdis-core!85
parents d7705dd7 449811fe
...@@ -11,11 +11,12 @@ ...@@ -11,11 +11,12 @@
#include <dune/grid/common/gridenums.hh> #include <dune/grid/common/gridenums.hh>
#include <dune/localfunctions/common/localkey.hh> #include <dune/localfunctions/common/localkey.hh>
#include <amdis/Output.hpp>
#include <amdis/common/Apply.hpp> #include <amdis/common/Apply.hpp>
#include <amdis/common/ForEach.hpp> #include <amdis/common/ForEach.hpp>
#include <amdis/common/TupleUtility.hpp> #include <amdis/common/TupleUtility.hpp>
#include <amdis/functions/Nodes.hpp> #include <amdis/functions/Nodes.hpp>
#include <amdis/Output.hpp> #include <amdis/utility/Twist.hpp>
namespace Dune namespace Dune
{ {
...@@ -88,9 +89,20 @@ namespace AMDiS ...@@ -88,9 +89,20 @@ namespace AMDiS
{ {
using Super = std::pair<EntityIdType, std::size_t>; using Super = std::pair<EntityIdType, std::size_t>;
IdType(int i = 0) : Super() {}; IdType(std::size_t i = 0) : Super(EntityIdType{}, i) {};
using Super::Super; using Super::Super;
IdType& operator++()
{
++this->second;
return *this;
}
IdType operator+(std::size_t shift) const
{
return IdType{this->first, this->second + shift};
}
friend std::ostream& operator<<(std::ostream& os, IdType const& id) friend std::ostream& operator<<(std::ostream& os, IdType const& id)
{ {
os << "(" << id.first << "," << id.second << ")"; os << "(" << id.first << "," << id.second << ")";
...@@ -101,11 +113,13 @@ namespace AMDiS ...@@ -101,11 +113,13 @@ namespace AMDiS
using PreBasis = typename GlobalBasis::PreBasis; using PreBasis = typename GlobalBasis::PreBasis;
using TreePath = typename GlobalBasis::PrefixPath; using TreePath = typename GlobalBasis::PrefixPath;
using NodeIdSet = AMDiS::NodeIdSet<PreBasis, TreePath>; using NodeIdSet = AMDiS::NodeIdSet<PreBasis, TreePath>;
using Twist = AMDiS::Twist<typename Grid::GlobalIdSet, GridView::dimension>;
public: public:
GlobalBasisIdSet(GlobalBasis const& globalBasis) GlobalBasisIdSet(GlobalBasis const& globalBasis)
: tree_(makeNode(globalBasis.preBasis(), TreePath{})) : tree_(makeNode(globalBasis.preBasis(), TreePath{}))
, nodeIdSet_(globalBasis.gridView()) , nodeIdSet_(globalBasis.gridView())
, twist_(globalBasis.gridView().grid().globalIdSet())
{ {
Dune::Functions::initializeTree(tree_); Dune::Functions::initializeTree(tree_);
} }
...@@ -120,8 +134,9 @@ namespace AMDiS ...@@ -120,8 +134,9 @@ namespace AMDiS
{ {
Dune::Functions::bindTree(tree_, element); Dune::Functions::bindTree(tree_, element);
nodeIdSet_.bind(tree_); nodeIdSet_.bind(tree_);
twist_.bind(element);
data_.resize(size()); data_.resize(size());
nodeIdSet_.fillIn(data_.begin()); nodeIdSet_.fillIn(twist_, data_.begin());
} }
/// \brief unbind from the element /// \brief unbind from the element
...@@ -162,6 +177,7 @@ namespace AMDiS ...@@ -162,6 +177,7 @@ namespace AMDiS
protected: protected:
Tree tree_; Tree tree_;
NodeIdSet nodeIdSet_; NodeIdSet nodeIdSet_;
Twist twist_;
using Data = std::pair<IdType, PartitionType>; using Data = std::pair<IdType, PartitionType>;
std::vector<Data> data_; std::vector<Data> data_;
}; };
...@@ -195,6 +211,7 @@ namespace AMDiS ...@@ -195,6 +211,7 @@ namespace AMDiS
public: public:
NodeIdSet(GridView const& gridView) NodeIdSet(GridView const& gridView)
: gridView_(gridView) : gridView_(gridView)
, sizes_{}
{} {}
/// \brief Bind the idset to a tree node /// \brief Bind the idset to a tree node
...@@ -202,6 +219,15 @@ namespace AMDiS ...@@ -202,6 +219,15 @@ namespace AMDiS
{ {
node_ = &node; node_ = &node;
size_ = node_->finiteElement().size(); size_ = node_->finiteElement().size();
std::fill(std::begin(sizes_), std::end(sizes_), 0u);
for (size_type i = 0; i < size_ ; ++i) {
Dune::LocalKey const& localKey = node_->finiteElement().localCoefficients().localKey(i);
sizes_[localKey.codim()]++;
}
auto refElem = Dune::referenceElement<double,GridView::dimension>(node_->element().type());
for (size_type c = 0; c <= GridView::dimension ; ++c)
sizes_[c] /= refElem.size(c);
} }
/// \brief Unbind the idset /// \brief Unbind the idset
...@@ -218,8 +244,8 @@ namespace AMDiS ...@@ -218,8 +244,8 @@ namespace AMDiS
/// \brief Maps from subtree index set [0..size-1] to a globally unique id in global basis /// \brief Maps from subtree index set [0..size-1] to a globally unique id in global basis
// [[expects: node_ != nullptr]] // [[expects: node_ != nullptr]]
template <class It> template <class Twist, class It>
It fillIn(It it, size_type shift = 0) const It fillIn(Twist const& twist, It it, size_type shift = 0) const
{ {
const auto& gridIdSet = gridView_.grid().globalIdSet(); const auto& gridIdSet = gridView_.grid().globalIdSet();
...@@ -227,11 +253,7 @@ namespace AMDiS ...@@ -227,11 +253,7 @@ namespace AMDiS
Dune::LocalKey localKey = node_->finiteElement().localCoefficients().localKey(i); Dune::LocalKey localKey = node_->finiteElement().localCoefficients().localKey(i);
unsigned int s = localKey.subEntity(); unsigned int s = localKey.subEntity();
unsigned int c = localKey.codim(); unsigned int c = localKey.codim();
unsigned int idx = localKey.index(); it->first = {gridIdSet.subId(node_->element(), s, c), shift + twist.get(localKey,sizes_[c])};
if (!(c == GridView::dimension || c == 0 || idx == 0))
DUNE_THROW(Dune::NotImplemented, "Bases with more then one DoF per subentity are not supported.");
it->first = {gridIdSet.subId(node_->element(), s, c), shift + idx};
it->second = Dune::Hybrid::switchCases(std::make_index_sequence<dim+1>{}, c, it->second = Dune::Hybrid::switchCases(std::make_index_sequence<dim+1>{}, c,
[&](auto codim) { return node_->element().template subEntity<codim>(s).partitionType(); }, [&](auto codim) { return node_->element().template subEntity<codim>(s).partitionType(); },
...@@ -248,6 +270,7 @@ namespace AMDiS ...@@ -248,6 +270,7 @@ namespace AMDiS
GridView gridView_; GridView gridView_;
const Node* node_ = nullptr; const Node* node_ = nullptr;
size_type size_ = 0; size_type size_ = 0;
std::array<unsigned int,GridView::dimension+1> sizes_;
}; };
...@@ -293,13 +316,13 @@ namespace AMDiS ...@@ -293,13 +316,13 @@ namespace AMDiS
/// \brief Maps from subtree index set [0..size-1] to a globally unique id in global basis /// \brief Maps from subtree index set [0..size-1] to a globally unique id in global basis
// [[expects: node_ != nullptr]] // [[expects: node_ != nullptr]]
template <class It> template <class Twist, class It>
It fillIn(It it, size_type shift = 0) const It fillIn(Twist const& twist, It it, size_type shift = 0) const
{ {
for (std::size_t child = 0; child < children; ++child) for (std::size_t child = 0; child < children; ++child)
{ {
size_type subTreeSize = subIds_.size(); size_type subTreeSize = subIds_.size();
it = subIds_.fillIn(it, shift); it = subIds_.fillIn(twist, it, shift);
shift += subTreeSize; shift += subTreeSize;
} }
return it; return it;
...@@ -370,13 +393,13 @@ namespace AMDiS ...@@ -370,13 +393,13 @@ namespace AMDiS
/// \brief Maps from subtree index set [0..size-1] to a globally unique id in global basis /// \brief Maps from subtree index set [0..size-1] to a globally unique id in global basis
// [[expects: node_ != nullptr]] // [[expects: node_ != nullptr]]
template <class It> template <class Twist, class It>
It fillIn(It it, size_type shift = 0) const It fillIn(Twist const& twist, It it, size_type shift = 0) const
{ {
Tools::for_each(idsTuple_, [&](auto const& ids) Tools::for_each(idsTuple_, [&](auto const& ids)
{ {
size_type subTreeSize = ids.size(); size_type subTreeSize = ids.size();
it = ids.fillIn(it, shift); it = ids.fillIn(twist, it, shift);
shift += subTreeSize; shift += subTreeSize;
}); });
return it; return it;
...@@ -443,15 +466,15 @@ namespace AMDiS ...@@ -443,15 +466,15 @@ namespace AMDiS
/// \brief Maps from subtree index set [0..size-1] to a globally unique id in global basis /// \brief Maps from subtree index set [0..size-1] to a globally unique id in global basis
// [[expects: node_ != nullptr]] // [[expects: node_ != nullptr]]
template <class It> template <class Twist, class It>
It fillIn(It it, size_type shift = 0) const It fillIn(Twist const& twist, It it, size_type shift = 0) const
{ {
for (int child = 0; child < dow; ++child) { for (int child = 0; child < dow; ++child) {
size_type subTreeSize = pq2NodeIdSet_.size(); size_type subTreeSize = pq2NodeIdSet_.size();
it = pq2NodeIdSet_.fillIn(it, shift); it = pq2NodeIdSet_.fillIn(twist, it, shift);
shift += subTreeSize; shift += subTreeSize;
} }
it = pq1NodeIdSet_.fillIn(it, shift); it = pq1NodeIdSet_.fillIn(twist, it, shift);
return it; return it;
} }
...@@ -496,8 +519,8 @@ namespace AMDiS ...@@ -496,8 +519,8 @@ namespace AMDiS
/// \brief Maps from subtree index set [0..size-1] to a globally unique id in global basis /// \brief Maps from subtree index set [0..size-1] to a globally unique id in global basis
// [[expects: node_ != nullptr]] // [[expects: node_ != nullptr]]
template <class It> template <class Twist, class It>
It fillIn(It it, size_type shift = 0) const It fillIn(Twist const& /*twist*/, It it, size_type shift = 0) const
{ {
const auto& gridIdSet = gridView_.grid().globalIdSet(); const auto& gridIdSet = gridView_.grid().globalIdSet();
auto elementId = gridIdSet.id(node_->element()); auto elementId = gridIdSet.id(node_->element());
......
...@@ -79,7 +79,6 @@ namespace AMDiS ...@@ -79,7 +79,6 @@ namespace AMDiS
/// Dummy implementation for ISTL-specific communication when no MPI is found /// Dummy implementation for ISTL-specific communication when no MPI is found
template <class Basis> template <class Basis>
class ISTLCommunication class ISTLCommunication
: public DefaultCommunication<Basis>
{ {
using Self = ISTLCommunication; using Self = ISTLCommunication;
......
...@@ -456,7 +456,7 @@ namespace AMDiS ...@@ -456,7 +456,7 @@ namespace AMDiS
// default iterative solver // default iterative solver
Map::addCreator("default", gmres); Map::addCreator("default", gmres);
init_direct(std::is_same<typename Dune::FieldTraits<T>::real_type, double>{}); init_direct(std::is_same<typename Dune::FieldTraits<typename Matrix::value_type>::real_type, double>{});
} }
static void init_direct(std::false_type) {} static void init_direct(std::false_type) {}
......
...@@ -48,7 +48,6 @@ namespace AMDiS ...@@ -48,7 +48,6 @@ namespace AMDiS
/// Implementation of \ref RunnerInterface::init() /// Implementation of \ref RunnerInterface::init()
void init(Matrix const& matrix, Comm&) override void init(Matrix const& matrix, Comm&) override
{ {
DUNE_UNUSED_PARAMETER(comm);
try { try {
solver_.reset(new SolverType(matrix, symmetricStrategy_, allocInit_)); solver_.reset(new SolverType(matrix, symmetricStrategy_, allocInit_));
} catch (mtl::mat::umfpack::error const& e) { } catch (mtl::mat::umfpack::error const& e) {
......
...@@ -4,5 +4,6 @@ install(FILES ...@@ -4,5 +4,6 @@ install(FILES
LocalToGlobalAdapter.hpp LocalToGlobalAdapter.hpp
MacroGridFactory.hpp MacroGridFactory.hpp
QuadratureFactory.hpp QuadratureFactory.hpp
Twist.hpp
UniqueBorderPartition.hpp UniqueBorderPartition.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/utility) DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/utility)
#pragma once
#include <vector>
#include <dune/geometry/referenceelements.hh>
#include <dune/localfunctions/common/localkey.hh>
#include <amdis/Output.hpp>
namespace AMDiS
{
/// \brief Permutate the dof indices w.r.t. a global entity orientation
template <class IdSet, int dim>
class Twist
{
using IdType = typename IdSet::IdType;
using RefElements = typename Dune::ReferenceElements<double,dim>;
using RefElement = typename RefElements::ReferenceElement;
public:
/// Constructor. Stores a reference to the passed idSet
Twist(IdSet const& idSet)
: idSet_(idSet)
{}
/// Bind the twist to a codim-0 element. Calculates the global ids of all the element vertices.
template <class Element>
void bind(Element const& element)
{
static_assert(dim == Element::mydimension, "");
refElem_ = &RefElements::general(element.type());
ids_.resize(refElem_->size(dim));
for (int i = 0; i < refElem_->size(dim); ++i)
ids_[i] = idSet_.subId(element, i, dim);
}
/// Get the permutated local dof index, living on a subEntity of the bound element.
/**
* If the global entity orientation differs from the local orientation, perform a
* permutation of the local dof indices on that entity. This is currently implemented
* for edge entities only.
*
* \param localKey The Dune::LocalKey of the local dof, identifying the entity and local index
* \param size The number of local dofs on that entity.
*/
unsigned int get(Dune::LocalKey const& localKey, unsigned int size) const
{
int subDim = dim - localKey.codim();
if (subDim == 1 && 1 < dim) // facet == edge
return id(localKey,0) < id(localKey,1) ? localKey.index() : size - localKey.index() - 1;
else if (subDim == 2 && 2 < dim) { // facet == face
test_exit(size == 1, "Bases with more then one DoF per subentity are not yet supported.");
return localKey.index();
}
return localKey.index();
}
private:
IdType const& id(Dune::LocalKey const& localKey, int ii) const
{
return ids_[refElem_->subEntity(localKey.subEntity(), localKey.codim(), ii, dim)];
}
private:
IdSet const& idSet_;
RefElement const* refElem_ = nullptr;
std::vector<IdType> ids_;
};
} // end namespace AMDiS
...@@ -40,12 +40,15 @@ bool test(Basis& basis, std::string const& prefix) ...@@ -40,12 +40,15 @@ bool test(Basis& basis, std::string const& prefix)
// Make communicator and exchange dofs // Make communicator and exchange dofs
auto comm = Comm::create(basis, prefix); auto comm = Comm::create(basis, prefix);
#if HAVE_MPI
comm->copyOwnerToAll(dofs, dofs); comm->copyOwnerToAll(dofs, dofs);
// Compare post-exchange data with reference // Compare post-exchange data with reference
for (std::size_t i = 0; i < dofs.size(); ++i) for (std::size_t i = 0; i < dofs.size(); ++i)
if (std::abs(dofs[i]-ref[i]) > AMDIS_TEST_TOL) if (std::abs(dofs[i]-ref[i]) > AMDIS_TEST_TOL)
passed = false; passed = false;
#endif
return passed; return passed;
} }
......
Supports Markdown
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