Commit 9f44e925 authored by Praetorius, Simon's avatar Praetorius, Simon

ConstantExpr with reference_wrapper argument, LocalAssembler accepts operator by reference_wrapper

parent 7e44b707
Pipeline #940 failed with stage
in 1 minute and 31 seconds
......@@ -8,7 +8,7 @@
#include <dune/amdis/DirichletBC.hpp>
#include <dune/amdis/LinearAlgebra.hpp>
#include <dune/amdis/Operator.hpp>
#include <dune/amdis/LocalAssemblerBase.hpp>
#include <dune/amdis/common/Mpl.hpp>
#include <dune/amdis/common/TypeDefs.hpp>
......
......@@ -42,7 +42,7 @@ void Assembler<Traits>::assemble(
{
auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
auto rowLocalView = rowBasis.localView();
rowLocalView.bind(element); // NOTE: Is this necessary
rowLocalView.bind(element); // NOTE: Is this necessary?
auto& rhsOp = rhsOperators_[rowNode];
if (rhsOp.assemble(asmVector) && !rhsOp.empty())
......@@ -54,7 +54,7 @@ void Assembler<Traits>::assemble(
if (matOp.assemble(asmMatrix) && !matOp.empty()) {
auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
auto colLocalView = colBasis.localView();
colLocalView.bind(element); // NOTE: Is this necessary
colLocalView.bind(element); // NOTE: Is this necessary?
assembleElementOperators(elementMatrix, matrix, matOp, geometry, rowLocalView, colLocalView);
}
......@@ -159,9 +159,7 @@ void Assembler<Traits>::initMatrixVector(
{
auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
for (auto bc : constraints_[rowNode][colNode].scalar)
bc->init(matrix, solution, rhs, rowBasis, colBasis);
for (auto bc : constraints_[rowNode][colNode].vector)
for (auto bc : constraints_[rowNode][colNode])
bc->init(matrix, solution, rhs, rowBasis, colBasis);
});
});
......@@ -195,9 +193,7 @@ std::size_t Assembler<Traits>::finishMatrixVector(
matOp.assembled = true;
// finish boundary condition
for (auto bc : constraints_[rowNode][colNode].scalar)
bc->finish(matrix, solution, rhs, rowBasis, colBasis);
for (auto bc : constraints_[rowNode][colNode].vector)
for (auto bc : constraints_[rowNode][colNode])
bc->finish(matrix, solution, rhs, rowBasis, colBasis);
});
});
......
......@@ -8,7 +8,7 @@ dune_add_library("duneamdis" NO_EXPORT
AMDiS.cpp
Initfile.cpp
ProblemInstatBase.cpp
ProblemInstat.cpp
# ProblemInstat.cpp
ProblemStat.cpp
StandardProblemIteration.cpp
#linear_algebra/istl/SystemMatrix.cpp
......
......@@ -7,6 +7,7 @@
#include <dune/amdis/Output.hpp>
#include <dune/amdis/common/Concepts.hpp>
#include <dune/amdis/utility/RangeType.hpp>
#include <dune/amdis/utility/TreeData.hpp>
namespace AMDiS
......@@ -34,8 +35,8 @@ namespace AMDiS
REQUIRES(Concepts::Functor<Predicate, bool(Domain)> &&
Concepts::Functor<Values, Range(Domain)>) >
DirichletBC(Predicate&& predicate, Values&& values)
: predicate(std::forward<Predicate>(predicate))
, values(std::forward<Values>(values))
: predicate_(std::forward<Predicate>(predicate))
, values_(std::forward<Values>(values))
{}
......@@ -78,11 +79,11 @@ namespace AMDiS
void initImpl(RB const&, RowNodeTag) {}
private:
std::function<bool(Domain)> predicate;
std::function<Range(Domain)> values;
std::function<bool(Domain)> predicate_;
std::function<Range(Domain)> values_;
bool initialized = false;
std::vector<char> dirichletNodes;
bool initialized_ = false;
std::vector<char> dirichletNodes_;
};
......@@ -92,14 +93,14 @@ namespace AMDiS
using WorldVector = typename GlobalBasis::GridView::template Codim<0>::Geometry::GlobalCoordinate;
template <class RowNode, class ColNode>
struct type
{
std::list<std::shared_ptr<DirichletBC<WorldVector, double>>> scalar;
std::list<std::shared_ptr<DirichletBC<WorldVector, WorldVector>>> vector;
void push_back(std::shared_ptr<DirichletBC<WorldVector, double>> const& bc) { scalar.push_back(bc); }
void push_back(std::shared_ptr<DirichletBC<WorldVector, WorldVector>> const& bc) { vector.push_back(bc); }
};
using type = std::list< std::shared_ptr<DirichletBC<WorldVector, RangeType<RowNode>>> >;
// {
// using Range = RangeType<RowNode>;
// std::list<std::shared_ptr<DirichletBC<WorldVector, Range>>> bc;
// void push_back(std::shared_ptr<DirichletBC<WorldVector, double>> const& bc) { scalar.push_back(bc); }
// void push_back(std::shared_ptr<DirichletBC<WorldVector, WorldVector>> const& bc) { vector.push_back(bc); }
// };
};
template <class GlobalBasis>
......
......@@ -14,9 +14,9 @@ namespace AMDiS
{
using Dune::Functions::interpolate;
if (!initialized) {
interpolate(rowBasis, hierarchicVectorWrapper(dirichletNodes), predicate);
initialized = true;
if (!initialized_) {
interpolate(rowBasis, hierarchicVectorWrapper(dirichletNodes_), predicate_);
initialized_ = true;
}
}
......@@ -27,13 +27,13 @@ namespace AMDiS
{
using Dune::Functions::interpolate;
if (!initialized) {
if (!initialized_) {
auto tp = rowBasis.prefixPath();
auto const& basis = rowBasis.rootBasis();
for (std::size_t i = 0; i < degree(rowBasis.localView().tree()); ++i)
interpolate(Dune::Functions::subspaceBasis(basis, push_back(tp,i)),
hierarchicVectorWrapper(dirichletNodes), predicate);
initialized = true;
hierarchicVectorWrapper(dirichletNodes_), predicate_);
initialized_ = true;
}
}
......@@ -47,11 +47,13 @@ namespace AMDiS
{
using Dune::Functions::interpolate;
test_exit_dbg(initialized, "Boundary condition not initialized!");
auto columns = matrix.applyDirichletBC(dirichletNodes);
test_exit_dbg(initialized_, "Boundary condition not initialized!");
auto columns = matrix.applyDirichletBC(dirichletNodes_);
interpolate(rowBasis, hierarchicVectorWrapper(rhs.getVector()), values, hierarchicVectorWrapper(dirichletNodes));
interpolate(colBasis, hierarchicVectorWrapper(solution.getVector()), values, hierarchicVectorWrapper(dirichletNodes));
interpolate(rowBasis, hierarchicVectorWrapper(rhs.getVector()), values_,
hierarchicVectorWrapper(dirichletNodes_));
interpolate(colBasis, hierarchicVectorWrapper(solution.getVector()), values_,
hierarchicVectorWrapper(dirichletNodes_));
// subtract columns of dirichlet nodes from rhs
for (auto const& triplet : columns)
......
......@@ -34,8 +34,8 @@ namespace AMDiS
/// Constructor. Stores a copy of operator `op`.
explicit LocalAssembler(Operator const& op)
: opStorage_(std::make_unique<Operator>(op))
, op_(*opStorage_)
: storage_(std::make_unique<Operator>(op))
, op_(*storage_)
{}
/// Constructor. Stores the reference to the operator.
......@@ -189,7 +189,7 @@ namespace AMDiS
private:
std::unique_ptr<Operator> opStorage_; //< the stored operator, implementing \ref ExpressionOperatorBase
std::unique_ptr<Operator> storage_; //< the stored operator, implementing \ref ExpressionOperatorBase
Operator& op_;
Element const* element_ = nullptr;
......
......@@ -193,18 +193,18 @@ namespace AMDiS
/// Generate an \ref ExpressionOperator with a given tag `t` and pre-expression `epxr`.
template <class Tag, class Expr_>
auto makeOperator(Tag t, Expr_ const& expr)
template <class Tag, class PreExpr>
auto makeOperator(Tag t, PreExpr const& expr)
{
using Expr = ToTerm_t<Expr_>;
using Expr = ToTerm_t<PreExpr>;
return ExpressionOperator<Tag, Expr>{t, toTerm(expr)};
}
/// Generate a shared_ptr to \ref ExpressionOperator with a given tag `t` and pre-expression `epxr`.
template <class Tag, class Expr_>
auto makeOperatorPtr(Tag t, Expr_ const& expr)
template <class Tag, class PreExpr>
auto makeOperatorPtr(Tag t, PreExpr const& expr)
{
using Expr = ToTerm_t<Expr_>;
using Expr = ToTerm_t<PreExpr>;
return std::make_shared<ExpressionOperator<Tag, Expr>>(t, toTerm(expr));
}
......
......@@ -6,5 +6,5 @@
namespace AMDiS
{
// explicit template instatiation
template class ProblemInstat<YaspGridBasis<2,1>>;
// template class ProblemInstat<YaspGridBasis<2,1>>;
} // end namespace AMDiS
......@@ -79,11 +79,6 @@ namespace AMDiS
std::unique_ptr<SystemVector> oldSolution;
};
#ifndef AMDIS_NO_EXTERN_PROBLEMINSTAT
extern template class ProblemInstat<YaspGridBasis<2,1>>;
#endif
} // end namespace AMDiS
#include "ProblemInstat.inc.hpp"
......@@ -21,7 +21,6 @@
#include <dune/amdis/Initfile.hpp>
#include <dune/amdis/LinearAlgebra.hpp>
#include <dune/amdis/Mesh.hpp>
#include <dune/amdis/Operator.hpp>
#include <dune/amdis/ProblemStatBase.hpp>
#include <dune/amdis/ProblemStatTraits.hpp>
#include <dune/amdis/StandardProblemIteration.hpp>
......@@ -143,14 +142,22 @@ namespace AMDiS
RowTreePath row, ColTreePath col,
Values const& values);
template <class Predicate, class RowTreePath, class ColTreePath>
void addDirichletBC(Predicate const& predicate,
RowTreePath row, ColTreePath col,
double constant)
{
addDirichletBC(predicate, row, col, [constant](auto const&) { return constant; });
}
private: // implementation detail
template <class Predicate, class RowNode, class ColNode>
void addDirichletBCImpl(Predicate const& predicate,
RowNode const& rowNode, ColNode const& colNode,
RangeType<RowNode> const& value);
template <class Predicate, class RowNode, class ColNode, class Values>
std::enable_if_t<not std::is_convertible<Values, RangeType<RowNode>>::value>
addDirichletBCImpl(Predicate const& predicate,
RowNode const& rowNode, ColNode const& colNode,
Values const& values);
public:
/// Implementation of \ref ProblemStatBase::solve
virtual void solve(AdaptInfo& adaptInfo,
......
......@@ -9,6 +9,7 @@
#include <dune/amdis/AdaptInfo.hpp>
#include <dune/amdis/Assembler.hpp>
#include <dune/amdis/FileWriter.hpp>
#include <dune/amdis/LocalAssembler.hpp>
#include <dune/amdis/common/Loops.hpp>
#include <dune/amdis/common/Timer.hpp>
......@@ -205,19 +206,38 @@ void ProblemStat<Traits>::
addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Values const& values)
{
static_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
"Function passed to addDirichletBC for predicate does not model the Functor<bool(WorldVector)> concept");
static_assert( Concepts::Callable<Values, WorldVector>,
"Function passed to addDirichletBC for values does not model the Callable<WorldVector> concept");
"Function passed to addDirichletBC for `predicate` does not model the Functor<bool(WorldVector)> concept");
auto i = child(globalBasis->localView().tree(), makeTreePath(row));
auto j = child(globalBasis->localView().tree(), makeTreePath(col));
using Range = decltype(values(std::declval<WorldVector>()));
addDirichletBCImpl(predicate, i, j, values);
}
template <class Traits>
template <class Predicate, class RowNode, class ColNode>
void ProblemStat<Traits>::
addDirichletBCImpl(Predicate const& predicate, RowNode const& row, ColNode const& col, RangeType<RowNode> const& value)
{
using Range = RangeType<RowNode>;
using BcType = DirichletBC<WorldVector,Range>;
auto bc = std::make_shared<BcType>(predicate, [value](auto const&) { return value; });
constraints[row][col].push_back(bc);
}
template <class Traits>
template <class Predicate, class RowNode, class ColNode, class Values>
std::enable_if_t<not std::is_convertible<Values, RangeType<RowNode>>::value>
ProblemStat<Traits>::
addDirichletBCImpl(Predicate const& predicate, RowNode const& row, ColNode const& col, Values const& values)
{
using Range = RangeType<RowNode>;
static_assert( Concepts::Functor<Values, Range(WorldVector)>,
"Function passed to addDirichletBC for `values` does not model the Functor<Range(WorldVector)> concept");
using BcType = DirichletBC<WorldVector,Range>;
auto bc = std::make_shared<BcType>(predicate, values);
constraints[i][i].push_back(bc);
constraints[row][col].push_back(bc);
}
......@@ -233,10 +253,6 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
solution->compress();
msg("matrix: ",num_rows(systemMatrix->getMatrix()), " x ", num_cols(systemMatrix->getMatrix()));
msg("x: ",num_rows(solution->getVector()));
msg("rhs: ",num_rows(rhs->getVector()));
linearSolver->solve(systemMatrix->getMatrix(),
solution->getVector(), rhs->getVector(),
solverInfo);
......
#pragma once
#include <functional>
#include <type_traits>
#include <dune/common/fmatrix.hh>
......@@ -34,6 +35,12 @@ namespace AMDiS
using type = tag::scalar;
};
template <class T>
struct ValueCategory<std::reference_wrapper<T>, std::enable_if_t< std::is_arithmetic<T>::value >>
{
using type = tag::scalar;
};
template <class K, int SIZE>
struct ValueCategory< Dune::FieldVector<K, SIZE> >
{
......
......@@ -34,8 +34,9 @@ namespace AMDiS
std::integral_constant<bool, sameFE> flagSameFE,
std::integral_constant<bool, sameNode> flagSameNode)
{
auto elementMatrixTransposed = trans(elementMatrix);
Transposed::calculateElementMatrix(
context, quad, elementMatrix, colNode, rowNode, flagSameFE, flagSameNode);
context, quad, elementMatrixTransposed, colNode, rowNode, flagSameFE, flagSameNode);
}
};
......
......@@ -34,8 +34,9 @@ namespace AMDiS
std::integral_constant<bool, sameFE> flagSameFE,
std::integral_constant<bool, sameNode> flagSameNode)
{
auto elementMatrixTransposed = trans(elementMatrix);
Transposed::calculateElementMatrix(
context, quad, elementMatrix, colNode, rowNode, flagSameFE, flagSameNode);
context, quad, elementMatrixTransposed, colNode, rowNode, flagSameFE, flagSameNode);
}
};
......
......@@ -34,8 +34,9 @@ namespace AMDiS
std::integral_constant<bool, sameFE> flagSameFE,
std::integral_constant<bool, sameNode> flagSameNode)
{
auto elementMatrixTransposed = trans(elementMatrix);
Transposed::calculateElementMatrix(
context, quad, elementMatrix, colNode, rowNode, flagSameFE, flagSameNode);
context, quad, elementMatrixTransposed, colNode, rowNode, flagSameFE, flagSameNode);
}
};
......
......@@ -37,8 +37,9 @@ namespace AMDiS
std::integral_constant<bool, sameFE> flagSameFE,
std::integral_constant<bool, sameNode> flagSameNode)
{
auto elementMatrixTransposed = trans(elementMatrix);
Transposed::calculateElementMatrix(
context, quad, elementMatrix, colNode, rowNode, flagSameFE, flagSameNode);
context, quad, elementMatrixTransposed, colNode, rowNode, flagSameFE, flagSameNode);
}
};
......
......@@ -34,8 +34,9 @@ namespace AMDiS
std::integral_constant<bool, sameFE> flagSameFE,
std::integral_constant<bool, sameNode> flagSameNode)
{
auto elementMatrixTransposed = trans(elementMatrix);
Transposed::calculateElementMatrix(
context, quad, elementMatrix, colNode, rowNode, flagSameFE, flagSameNode);
context, quad, elementMatrixTransposed, colNode, rowNode, flagSameFE, flagSameNode);
}
};
......
......@@ -2,7 +2,6 @@
#include <type_traits>
#include <dune/amdis/OperatorTermBase.hpp>
#include <dune/amdis/common/ScalarTypes.hpp>
#include <dune/amdis/terms/ExpressionBase.hpp>
......@@ -14,15 +13,15 @@ namespace AMDiS
**/
/// An expression representing a constant (arithmetic) value, \see constant
template <class ValueType>
template <class T>
class ConstantTerm
: public ExpressionBase
{
public:
// TODO: possibly convert to plain type, in case of expression templates.
using value_type = ValueType;
using value_type = T;
explicit ConstantTerm(value_type value)
explicit ConstantTerm(value_type const& value)
: ExpressionBase(0)
, value_(value)
{}
......@@ -37,9 +36,33 @@ namespace AMDiS
};
/// An expression representing a reference to a constant (arithmetic) value, \see constant
template <class T>
class ConstantTerm<std::reference_wrapper<T>>
: public ExpressionBase
{
public:
// TODO: possibly convert to plain type, in case of expression templates.
using value_type = T;
explicit ConstantTerm(std::reference_wrapper<value_type> const& ref)
: ExpressionBase(0)
, ref_(ref)
{}
value_type const& operator[](std::size_t iq) const
{
return ref_;
}
private:
std::reference_wrapper<value_type> ref_;
};
/// generator function for \ref ConstantTerm expressions
template <class T>
auto constant(T value) { return ConstantTerm<T>{value}; }
auto constant(T const& value) { return ConstantTerm<T>{value}; }
/** @} **/
......
......@@ -5,7 +5,6 @@
#include <dune/istl/bvector.hh>
#include <dune/amdis/OperatorTermBase.hpp>
#include <dune/amdis/common/Concepts.hpp>
#include <dune/amdis/terms/ExpressionBase.hpp>
......
......@@ -7,7 +7,6 @@
#include <dune/istl/bvector.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/amdis/OperatorTermBase.hpp>
#include <dune/amdis/terms/ExpressionBase.hpp>
#include <dune/amdis/utility/GetDegree.hpp>
......
......@@ -6,10 +6,10 @@
#include <dune/istl/bvector.hh>
#include <dune/amdis/OperatorTermBase.hpp>
#include <dune/amdis/common/IndexSeq.hpp>
#include <dune/amdis/common/Concepts.hpp>
#include <dune/amdis/terms/ExpressionBase.hpp>
#include <dune/amdis/terms/TermGenerator.hpp>
namespace AMDiS
{
......@@ -28,11 +28,10 @@ namespace AMDiS
public:
using value_type = typename std::result_of<Functor(Value_t<Terms>...)>::type;
template <class F, class... Ts,
REQUIRES(Concepts::Compatible<Functor, F>) >
explicit FunctorTerm(F&& f, Ts&&... ts)
template <class... Ts>
explicit FunctorTerm(Functor const& fct, Ts&&... ts)
: ExpressionBase(2) // TODO: provide functor degree
, fct_(std::forward<F>(f))
, fct_(fct)
, terms_(std::forward<Ts>(ts)...)
{}
......@@ -80,10 +79,10 @@ namespace AMDiS
/// Generator function for \ref FunctorTerm expressions
template <class F, class... Ts>
auto func(F&& f, Ts&&... terms)
template <class Functor, class... Ts>
auto func(Functor const& f, Ts&&... terms)
{
return FunctorTerm<std::decay_t<F>, std::decay_t<Ts>...>{std::forward<F>(f), std::forward<Ts>(terms)...};
return FunctorTerm<Functor, ToTerm_t<Ts>...>{f, toTerm(std::forward<Ts>(terms))...};
}
/** @} **/
......
/** \file TermGenerator.hpp */
#pragma once
// std c++ includes
......@@ -30,17 +28,16 @@ namespace AMDiS
!std::is_same<ValueCategory_t<T>, tag::unknown>::value;
// TODO: generalize to functors with arbitrary return values
#ifndef AMDIS_DOW
template <class T>
constexpr bool FunctorCategory =
Concepts::Functor<T, double(Dune::FieldVector<double, 1>)> ||
Concepts::Functor<T, double(Dune::FieldVector<double, 2>)> ||
Concepts::Functor<T, double(Dune::FieldVector<double, 3>)>;
Concepts::Callable<T, Dune::FieldVector<double, 1>> ||
Concepts::Callable<T, Dune::FieldVector<double, 2>> ||
Concepts::Callable<T, Dune::FieldVector<double, 3>>;
#else
template <class T>
constexpr bool FunctorCategory =
Concepts::Functor<T, double(Dune::FieldVector<double, AMDIS_DOW>)>;
Concepts::Callable<T, Dune::FieldVector<double, AMDIS_DOW>>;
#endif
......
......@@ -14,7 +14,9 @@ namespace AMDiS
struct RangeType<Node, true, false, false>
{
using LocalBasis = typename Node::FiniteElement::Traits::LocalBasisType;
using type = typename LocalBasis::Traits::RangeType;
using type_ = typename LocalBasis::Traits::RangeType;
using type = std::conditional_t<std::is_same<type_, Dune::FieldVector<double,1>>::value,
double, type_>;
};
template <class Node>
......
......@@ -4,9 +4,8 @@
#include <iostream>
#include <dune/amdis/AMDiS.hpp>
#include <dune/amdis/AdaptInstationary.hpp>
#include <dune/amdis/ProblemInstat.hpp>
#include <dune/amdis/ProblemStat.hpp>
#include <dune/amdis/Operators.hpp>
#include <dune/amdis/Terms.hpp>
#include <dune/amdis/common/Literals.hpp>
......@@ -27,7 +26,7 @@ int main(int argc, char** argv)
AdaptInfo adaptInfo("adapt");
auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
prob.addMatrixOperator(std::ref(opL), _0, _0);
prob.addMatrixOperator(opL, _0, _0);
auto opForce = makeOperator(tag::test{}, [](auto const& x) { return -1.0; });
prob.addVectorOperator(opForce, _0);
......
......@@ -32,7 +32,7 @@ int main(int argc, char** argv)
double* invTau = probInstat.getInvTau();
auto opTimeLhs = makeOperator(tag::test_trial{}, [invTau](auto const&) { return *invTau; });
auto opTimeLhs = makeOperator(tag::test_trial{}, std::ref(*invTau));
prob.addMatrixOperator(opTimeLhs, 0, 0);
auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
......
......@@ -4,6 +4,7 @@
#include <dune/amdis/AMDiS.hpp>
#include <dune/amdis/ProblemStat.hpp>
#include <dune/amdis/Operators.hpp>
#include <dune/amdis/Terms.hpp>
#ifdef DOW
......@@ -33,9 +34,8 @@ int main(int argc, char** argv)
auto _v = 0_c;
auto _p = 1_c;
// define the differential operators
// <viscosity*grad(u_i), grad(v_i)>
for (std::size_t i = 0; i < DOW; ++i) {
// <viscosity*grad(u_i), grad(v_i)>
auto opL = makeOperator(tag::gradtest_gradtrial{}, viscosity);
prob.addMatrixOperator(opL, treepath(_v,i), treepath(_v,i));
}
......@@ -71,12 +71,23 @@ int main(int argc, char** argv)
prob.addDirichletBC(left, _v, _v, parabolic_y);
prob.addDirichletBC(not_left, _v, _v, zero);
prob.addDirichletBC([](auto const& x) { return x[0] < 1.e-8 && x[1] < 1.e-8; }, _p, _p, zero);
prob.addDirichletBC([](auto const& x) { return x[0] < 1.e-8 && x[1] < 1.e-8; }, _p, _p, 0.0);