Commit 2b0d9042 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Merge branch 'issue/bugfixes' into 'develop'

added virtual destructors, wrong extension string in mesh reader, forward…

See merge request !8
parents d355c059 98c0ebbe
Pipeline #1031 failed with stage
in 1 minute and 24 seconds
......@@ -51,6 +51,8 @@ namespace AMDiS
public:
virtual ~LocalAssemblerBase() {}
virtual void bind(Element const& element, Geometry const& geometry) = 0;
virtual void unbind() = 0;
......@@ -97,6 +99,21 @@ namespace AMDiS
{
return flag && (!assembled || changing);
}
template <class Geo>
void bind(Element const& elem, Geo const& geo)
{
for (auto& scaled : element) scaled.op->bind(elem,geo);
for (auto& scaled : boundary) scaled.op->bind(elem,geo);
for (auto& scaled : intersection) scaled.op->bind(elem,geo);
}
void unbind()
{
for (auto& scaled : element) scaled.op->unbind();
for (auto& scaled : boundary) scaled.op->unbind();
for (auto& scaled : intersection) scaled.op->unbind();
}
};
public:
......
......@@ -113,14 +113,14 @@ namespace AMDiS
auto ext = fn.extension();
#if HAVE_ALBERTA
if (ext == "1d" || ext == "2d" || ext == "3d") {
if (ext == ".1d" || ext == ".2d" || ext == ".3d") {
Dune::GridFactory<Grid> factory;
Dune::AlbertaReader<Grid> reader;
reader.readGrid(filename, factory);
return std::unique_ptr<Grid>{factory.createGrid()};
}
#endif
if (ext == "msh") {
if (ext == ".msh") {
Dune::GmshReader<Grid> reader;
return std::unique_ptr<Grid>{reader.read(filename)};
}
......
......@@ -14,5 +14,6 @@
#include <amdis/operations/Basic.hpp>
#include <amdis/operations/CMath.hpp>
#include <amdis/operations/Composer.hpp>
#include <amdis/operations/Composer.impl.hpp>
#include <amdis/operations/FieldMatVec.hpp>
#include <amdis/operations/MaxMin.hpp>
......@@ -79,6 +79,23 @@ namespace AMDiS
std::unique_ptr<SystemVector> oldSolution;
};
#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
// Deduction rule
template <class Traits>
ProblemInstat(std::string name, ProblemStat<Traits>& prob)
-> ProblemInstat<Traits>;
#endif
// Generator for ProblemInstat with given ProblemStat
template <class Traits>
ProblemInstat<Traits>
makeProblemInstat(std::string name, ProblemStat<Traits>& prob)
{
return {std::move(name), prob};
}
} // end namespace AMDiS
#include "ProblemInstat.inc.hpp"
......@@ -65,11 +65,11 @@ namespace AMDiS
const auto exprValue = Super::coefficient(local);
// The gradients of the shape functions on the reference element
std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
thread_local std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
rowLocalFE.localBasis().evaluateJacobian(local, referenceGradients);
// Compute the shape function gradients on the real element
std::vector<Dune::FieldVector<double,Context::dow> > gradients(referenceGradients.size());
thread_local std::vector<Dune::FieldVector<double,Context::dow> > gradients(referenceGradients.size());
for (std::size_t i = 0; i < gradients.size(); ++i)
jacobian.mv(referenceGradients[i][0], gradients[i]);
......@@ -113,7 +113,7 @@ namespace AMDiS
tag::scalar)
{
auto const& localFE = node.finiteElement();
std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
thread_local std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
......@@ -129,7 +129,7 @@ namespace AMDiS
localFE.localBasis().evaluateJacobian(local, referenceGradients);
// Compute the shape function gradients on the real element
std::vector<Dune::FieldVector<double,Context::dow> > gradients(referenceGradients.size());
thread_local std::vector<Dune::FieldVector<double,Context::dow> > gradients(referenceGradients.size());
for (std::size_t i = 0; i < gradients.size(); ++i)
jacobian.mv(referenceGradients[i][0], gradients[i]);
......@@ -160,7 +160,7 @@ namespace AMDiS
tag::matrix)
{
auto const& localFE = node.finiteElement();
std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
thread_local std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
......@@ -177,7 +177,7 @@ namespace AMDiS
localFE.localBasis().evaluateJacobian(local, referenceGradients);
// Compute the shape function gradients on the real element
std::vector<Dune::FieldVector<double,Context::dow> > gradients(referenceGradients.size());
thread_local std::vector<Dune::FieldVector<double,Context::dow> > gradients(referenceGradients.size());
for (std::size_t i = 0; i < gradients.size(); ++i)
jacobian.mv(referenceGradients[i][0], gradients[i]);
......
......@@ -139,7 +139,7 @@ namespace AMDiS
&& all_of_v<Concepts::HasOrder<LFs>...>)>
int order(FunctorLocalFunction<Sig,F,LFs...> const& lf)
{
return Dune::Std::apply([&lf](auto const&... gfs) { return order(lf.fct(), order(gfs)...); },
return Dune::Std::apply([&lf](auto const&... lgfs) { return order(lf.fct(), order(lgfs)...); },
lf.localFcts());
}
......
......@@ -24,6 +24,8 @@ namespace AMDiS
Parameters::get(base + "->name", name_);
}
virtual ~FileWriterInterface() {}
// pure virtual method to be implemented by derived classes
virtual void writeFiles(AdaptInfo& adaptInfo, bool force) = 0;
......
......@@ -67,34 +67,10 @@ namespace AMDiS
return Dune::Std::apply([&](auto const&... gs) { return order(c.f_, deg(gs)...); }, c.gs_);
}
// forward declaration
struct Multiplies;
struct Plus;
/// Partial derivative of composed function:
/// Implements: // sum_i [ d_i(f)[g...] * d_j(g_i) ]
template <int J, class F, class... Gs>
auto partial(Composer<F,Gs...> const& c, index_t<J> _j)
{
auto index_seq = MakeSeq_t<sizeof...(Gs)>{};
// d_i(f)[g...] * d_j(g_i)
auto term_i = [&](auto const _i)
{
auto di_f = Dune::Std::apply([&](auto const&... gs) {
return compose(partial(c.f_, _i), gs...);
}, c.gs_);
auto const& g_i = std::get<_i>(c.gs_);
return compose(Multiplies{}, di_f, partial(g_i, _j));
};
// sum_i [ d_i(f)[g...] * d_j(g_i) ]
return Dune::Std::apply([&](auto const... _i)
{
return compose(Plus{}, term_i(_i)...);
}, index_seq);
}
auto partial(Composer<F,Gs...> const& c, index_t<J> _j);
#ifndef DOXYGEN
......
......@@ -4,219 +4,114 @@
#include <dune/typetree/treepath.hh>
#include <dune/typetree/visitor.hh>
namespace AMDiS {
#ifndef DOXYGEN // these are all internals and not public API. Only access is using applyToTree().
namespace AMDiS
{
// forward declaration of main engine struct
template <typename tag = Dune::TypeTree::StartTag, bool doApply = true>
template <typename NodeTag, bool visit = true>
struct TraverseTree;
// This struct is the core of the algorithm. While this specialization simply serves as the starting point
// of the traversal and takes care of some setup work, the struct has to be specialized for each TreeType node type it
// should support.
// The first parameter is the tag of the node type
// and the second one must always be specialized as true, as a value of false means the node should in fact not be visited.
// That case is already handled by a specialization of the struct.
template <bool doApply>
struct TraverseTree<Dune::TypeTree::StartTag, doApply>
{
template <typename Node, typename Visitor>
static void apply(Node&& node, Visitor&& visitor)
{
TraverseTree<Dune::TypeTree::NodeTag<Node>>::apply(std::forward<Node>(node),
std::forward<Visitor>(visitor),
Dune::TypeTree::hybridTreePath());
}
};
// Do not visit nodes the visitor is not interested in
template <typename NodeTag>
struct TraverseTree<NodeTag, false>
{
// we won't do anything with the objects, so having them all const
// works fine.
template <typename Node, typename Visitor, typename TreePath>
static void apply(const Node& node, const Visitor& visitor, TreePath treePath)
static void apply(const Node& node, const Visitor& visitor, TreePath const& tp)
{}
};
#ifndef DOXYGEN
// some implementation details
// ********************************************************************************
// LeafNode
// ********************************************************************************
template <class Node, class Index>
struct HybridChildType
: HybridChildType<std::remove_const_t<Node>, std::remove_const_t<Index>> {};
// LeafNode - just call the leaf() callback
template <>
struct TraverseTree<Dune::TypeTree::LeafNodeTag, true>
template <class Node>
struct HybridChildType<Node, std::size_t>
{
template <typename N, typename V, typename TreePath>
static void apply(N&& n, V&& v, TreePath tp)
{
v.leaf(std::forward<N>(n),tp);
}
using type = typename Node::template Child<0>::Type;
};
template <class Node, std::size_t K>
struct HybridChildType<Node, Dune::index_constant<K>>
{
using type = typename Node::template Child<K>::Type;
};
// ********************************************************************************
// PwerNode
// ********************************************************************************
template <class NodeTag, class Node>
constexpr std::size_t hybridDegree(NodeTag, Node const& node)
{
return Dune::TypeTree::degree(node);
}
template <>
struct TraverseTree<Dune::TypeTree::PowerNodeTag, true>
template <class Node>
constexpr auto hybridDegree(Dune::TypeTree::CompositeNodeTag, Node const& node)
{
template <typename N, typename V, typename TreePath>
static void apply(N&& n, V&& v, TreePath tp)
{
v.pre(std::forward<N>(n),tp);
return Dune::index_constant<Node::CHILDREN>{};
}
typedef typename std::remove_reference<N>::type Node;
typedef typename std::remove_reference<V>::type Visitor;
// get child type
typedef typename Node::template Child<0>::Type C;
template <std::size_t k, std::size_t n>
constexpr bool notLastChild(Dune::index_constant<k> const&, Dune::index_constant<n> const&)
{
return k < n-1;
}
constexpr bool notLastChild(std::size_t k, std::size_t n)
{
return k < n-1;
}
// Do we have to visit the children? As the TreePath is dynamic, it does not
// contain any information that could be evaluated at compile time, so we only
// have to query the visitor once.
const bool visit = Visitor::template VisitChild<Node,C,TreePath>::value;
#endif
auto child_tp = push_back(tp, std::size_t(0));
static const std::size_t last = treePathSize(tp);
// iterate over children
for (std::size_t k = 0; k < degree(n); ++k)
template <class NodeTag>
struct TraverseTree<NodeTag, true>
{
template <typename N, typename V, typename TreePath>
static void apply(N&& n, V&& v, TreePath const& tp)
{
using Node = std::remove_reference_t<N>;
using Visitor = std::remove_reference_t<V>;
v.pre(std::forward<N>(n),tp);
auto const deg = hybridDegree(NodeTag{}, n);
forEach(Dune::range(deg), [&](auto const _k)
{
// always call beforeChild(), regardless of the value of visit
v.beforeChild(std::forward<N>(n),n.child(k),tp,k);
// update TreePath
std::get<last>(child_tp._data) = k;
v.beforeChild(std::forward<N>(n),n.child(_k),tp,_k);
// descend to child
TraverseTree<Dune::TypeTree::NodeTag<C>,visit>::apply(n.child(k),std::forward<V>(v),child_tp);
using C = typename HybridChildType<Node, decltype(_k)>::type;
const bool visit = Visitor::template VisitChild<Node,C,TreePath>::value;
TraverseTree<Dune::TypeTree::NodeTag<C>,visit>::apply(n.child(_k),std::forward<V>(v),push_back(tp, _k));
// always call afterChild(), regardless of the value of visit
v.afterChild(std::forward<N>(n),n.child(k),tp,k);
v.afterChild(std::forward<N>(n),n.child(_k),tp,_k);
// if this is not the last child, call infix callback
if (k < degree(n) - 1)
if (notLastChild(_k, deg))
v.in(std::forward<N>(n),tp);
}
});
v.post(std::forward<N>(n),tp);
}
};
namespace Impl
{
// ********************************************************************************
// Static Version
// ********************************************************************************
template <std::size_t inverse_k, std::size_t count>
struct TraverseTreeStatic
{
template <typename N, typename V, typename TreePath>
static void apply(N&& n, V&& v, TreePath tp)
{
// make sure we do not try to work with references to the actual types
typedef typename std::remove_reference<N>::type Node;
typedef typename std::remove_reference<V>::type Visitor;
// get child type
typedef typename Node::template Child<count-inverse_k>::Type C;
// extend TreePath by child index
auto k = std::integral_constant<std::size_t, count-inverse_k>{};
auto child_tp = push_back(tp, k);
// is the visitor interested in this child?
const bool visit = Visitor::template VisitChild<Node,C,decltype(child_tp)>::value;
// beforeChild() gets called regardless of the value of visit
v.beforeChild(std::forward<N>(n),n.template child<count-inverse_k>(),tp,k);
// traverse to child
TraverseTree<Dune::TypeTree::NodeTag<C>,visit>::apply(n.template child<count-inverse_k>(),
std::forward<V>(v),
child_tp);
// afterChild() gets called regardless of the value of visit
v.afterChild(std::forward<N>(n),n.template child<count-inverse_k>(),tp,k);
// we are not at the last child (that is specialized), so call infix visitor callback
v.in(std::forward<N>(n),tp);
// continue with next child
TraverseTreeStatic<inverse_k-1,count>::apply(std::forward<N>(n), std::forward<V>(v), tp);
}
};
// Specialization for last child. This specialization stops the recursion and
// does not call the infix visitor on the CompositeNode.
template <std::size_t count>
struct TraverseTreeStatic<1,count>
{
template<typename N, typename V, typename TreePath>
static void apply(N&& n, V&& v, TreePath tp)
{
typedef typename std::remove_reference<N>::type Node;
typedef typename std::remove_reference<V>::type Visitor;
typedef typename Node::template Child<count-1>::Type C;
auto k = std::integral_constant<std::size_t, count-1>{};
auto child_tp = push_back(tp, k);
const bool visit = Visitor::template VisitChild<Node,C,decltype(child_tp)>::value;
v.beforeChild(std::forward<N>(n),n.template child<count-1>(),tp,k);
TraverseTree<Dune::TypeTree::NodeTag<C>,visit>::apply(n.template child<count-1>(),
std::forward<V>(v),
child_tp);
v.afterChild(std::forward<N>(n),n.template child<count-1>(),tp,k);
}
};
// Specialization for CompositeNode without any children.
template <>
struct TraverseTreeStatic<0,0>
{
template <typename N, typename V, typename TreePath>
static void apply(N&& n, V&& v, TreePath tp) {}
};
} // end namespace Impl
// ********************************************************************************
// CompositeNode
// ********************************************************************************
// Traverse CompositeNode - just forward to the generic algorithm
// LeafNode - just call the leaf() callback
template <>
struct TraverseTree<Dune::TypeTree::CompositeNodeTag,true>
struct TraverseTree<Dune::TypeTree::LeafNodeTag, true>
{
template <typename N, typename V, typename TreePath>
static void apply(N&& n, V&& v, TreePath tp)
static void apply(N&& n, V&& v, TreePath const& tp)
{
v.pre(std::forward<N>(n),tp);
typedef typename std::remove_reference<N>::type Node;
static const std::size_t C = Dune::TypeTree::StaticDegree<Node>::value;
Impl::TraverseTreeStatic<C,C>::apply(std::forward<N>(n), std::forward<V>(v), tp);
v.post(std::forward<N>(n),tp);
v.leaf(std::forward<N>(n),tp);
}
};
#endif // DOXYGEN
// ********************************************************************************
// Public Interface
// ********************************************************************************
//! Apply visitor to TypeTree.
/**
* This function applies the given visitor to the given tree. Both visitor and tree may be const
......@@ -231,7 +126,10 @@ namespace AMDiS {
template <typename Tree, typename Visitor>
void traverseTree(Tree&& tree, Visitor&& visitor)
{
TraverseTree<>::apply(std::forward<Tree>(tree), std::forward<Visitor>(visitor));
using Node = std::remove_reference_t<Tree>;
using NodeTag = Dune::TypeTree::NodeTag<Node>;
using TreePath = Dune::TypeTree::HybridTreePath<>;
TraverseTree<NodeTag>::apply(std::forward<Tree>(tree), std::forward<Visitor>(visitor), TreePath{});
}
} // end namespace AMDiS
......@@ -145,6 +145,26 @@ namespace AMDiS
return "";
}
namespace Impl
{
template <class TreePath, std::size_t... I>
std::array<std::size_t, sizeof...(I)> toArrayImpl(TreePath const& tp, std::index_sequence<I...>)
{
return {{std::size_t(Dune::TypeTree::treePathEntry<I>(tp))...}};
}
}
template <class T0, class... T>
auto to_array(Dune::TypeTree::HybridTreePath<T0,T...> const& tp)
{
return Impl::toArrayImpl(tp, std::make_index_sequence<1+sizeof...(T)>{});
}
std::array<std::size_t,1> to_array(Dune::TypeTree::HybridTreePath<> const& tp)
{
return {{0u}};
}
/// \brief Generate a TreePath from a sequence of integers and integral-constants
/**
* Converts a sequence of arguments to a \ref Dune::TypeTree::HybridTreePath.
......@@ -158,4 +178,33 @@ namespace AMDiS
return Dune::TypeTree::HybridTreePath<T...>(t...);
}
namespace Impl
{
template <class TreePath, std::size_t... I>
auto popFrontImpl(TreePath const& tp, std::index_sequence<I...>)
{
return Dune::TypeTree::hybridTreePath(Dune::TypeTree::treePathEntry<I+1>(tp)...);
}
template <class TreePath, std::size_t... I>
auto popBackImpl(TreePath const& tp, std::index_sequence<I...>)
{
return Dune::TypeTree::hybridTreePath(Dune::TypeTree::treePathEntry<I>(tp)...);
}
}
template <class T0, class... T>
Dune::TypeTree::HybridTreePath<T...> pop_front(Dune::TypeTree::HybridTreePath<T0,T...> const& tp)
{
return Impl::popFrontImpl(tp, std::make_index_sequence<sizeof...(T)>{});
}
template <class... T, class TN>
Dune::TypeTree::HybridTreePath<T...> pop_front(Dune::TypeTree::HybridTreePath<T...,TN> const& tp)
{
return Impl::popBackImpl(tp, std::make_index_sequence<sizeof...(T)>{});
}
} // end namespace AMDiS
......@@ -9,22 +9,9 @@ namespace AMDiS
// from dune-typetree merge-request !2
namespace Impl
{
struct CallInnerNode
{
template <class Node, class /*TreePath*/>
using Visit = bool_t<Node::isPower && Node::template Child<0>::type::isLeaf>;
};
struct CallAnyNode
{
template <class /*Node*/, class /*TreePath*/>
using Visit = std::true_type;
};
template <class PreFunc, class LeafFunc, class PostFunc>
class CallbackVisitor
: public Dune::TypeTree::TreeVisitor
, public CallAnyNode
{
public:
CallbackVisitor(PreFunc& preFunc, LeafFunc& leafFunc, PostFunc& postFunc)
......@@ -36,9 +23,7 @@ namespace AMDiS
template <typename Node, typename TreePath>
void pre(Node&& node, TreePath treePath)
{
static_assert(std::decay_t<Node>::isPower || std::decay_t<Node>::isComposite, "");
Dune::Hybrid::ifElse(Visit<std::decay_t<Node>, TreePath>{},
[&](auto id) { id(preFunc_)(node, treePath); });
preFunc_(node, treePath);
}
template <typename Node, typename TreePath>
......@@ -50,9 +35,7 @@ namespace AMDiS
template <typename Node, typename TreePath>
void post(Node&& node, TreePath treePath)
{
static_assert(std::decay_t<Node>::isPower || std::decay_t<Node>::isComposite, "");
Dune::Hybrid::ifElse(Visit<std::decay_t<Node>, TreePath>{},
[&](auto id) { id(postFunc_)(node, treePath); });
postFunc_(node, treePath);
}
private:
......