Commit 289b2b9a authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

working version of dune-functions assembling

parent 7668c6b7
......@@ -6,7 +6,9 @@
#ifdef HAVE_UMFPACK
#include <algorithm>
#include <iostream>
#include <string>
#include <boost/numeric/mtl/operation/two_norm.hpp>
#include <boost/numeric/mtl/interface/umfpack_solve.hpp>
......@@ -31,7 +33,8 @@ namespace AMDiS
class UmfpackRunner;
template <class Matrix, class Vector>
class UmfpackRunnerBase : public RunnerInterface<Matrix, Vector, Vector>
class UmfpackRunnerBase
: public RunnerInterface<Matrix, Vector, Vector>
{
protected:
using Super = RunnerInterface<Matrix, Vector>;
......@@ -43,10 +46,6 @@ namespace AMDiS
public:
/// Constructor. Reads UMFPACK parameters from initfile
UmfpackRunnerBase(std::string prefix)
: solver(NULL)
, store_symbolic(0)
, symmetric_strategy(0)
, alloc_init(0.7)
{
Parameters::get(prefix + "->store symbolic", store_symbolic); // ?
Parameters::get(prefix + "->symmetric strategy", symmetric_strategy);
......@@ -86,25 +85,21 @@ namespace AMDiS
protected:
std::shared_ptr<SolverType> solver;
int store_symbolic;
int symmetric_strategy;
double alloc_init;
int store_symbolic = 0;
int symmetric_strategy = 0;
double alloc_init = 0.7;
};
// specialization for block-matrices
template <class SubMatrix, std::size_t _N, std::size_t _M, class Vector>
class UmfpackRunner<BlockMTLMatrix<SubMatrix, _N, _M>, Vector>
: public UmfpackRunnerBase<BlockMTLMatrix<SubMatrix, _N, _M>, Vector>
: public UmfpackRunnerBase<BlockMTLMatrix<SubMatrix, _N, _M>, Vector>
{
using Matrix = BlockMTLMatrix<SubMatrix, _N, _M>;
using Super = UmfpackRunnerBase<Matrix, Vector>;
using SolverType = typename Super::SolverType;
using Super::symmetric_strategy;
using Super::alloc_init;
using Super::solver;
public:
UmfpackRunner(std::string prefix)
: Super(prefix)
......@@ -116,8 +111,8 @@ namespace AMDiS
copy(A, fullMatrix);
try {
Super::solver.reset(new SolverType(fullMatrix, symmetric_strategy, alloc_init));
} catch (mtl::mat::umfpack::error& e) {
Super::solver.reset(new SolverType(fullMatrix, Super::symmetric_strategy, Super::alloc_init));
} catch (mtl::mat::umfpack::error const& e) {
error_exit("UMFPACK_ERROR(factorize, ", e.code, ") = ", e.what());
}
}
......@@ -129,17 +124,13 @@ namespace AMDiS
// specialization for full-matrices
template <class T, class Param, class Vector>
class UmfpackRunner<mtl::compressed2D<T, Param>, Vector>
: public UmfpackRunnerBase<mtl::compressed2D<T, Param>, Vector>
: public UmfpackRunnerBase<mtl::compressed2D<T, Param>, Vector>
{
using FullMatrix = mtl::compressed2D<T, Param>;
using Super = UmfpackRunnerBase<FullMatrix, Vector>;
using SolverType = typename Super::SolverType;
using Super::symmetric_strategy;
using Super::alloc_init;
using Super::solver;
public:
UmfpackRunner(std::string prefix)
: Super(prefix)
......@@ -149,8 +140,8 @@ namespace AMDiS
virtual void init(FullMatrix const& fullMatrix) override
{
try {
solver.reset(new SolverType(fullMatrix, symmetric_strategy, alloc_init));
} catch (mtl::mat::umfpack::error& e) {
Super::solver.reset(new SolverType(fullMatrix, Super::symmetric_strategy, Super::alloc_init));
} catch (mtl::mat::umfpack::error const& e) {
error_exit("UMFPACK_ERROR(factorize, ", e.code, ") = ", e.what());
}
}
......
#pragma once
#include <vector>
#include <string>
#include <memory>
#include <tuple>
#include <dune/amdis/Output.hpp>
#include <dune/amdis/common/Loops.hpp>
#include <dune/amdis/common/ScalarTypes.hpp>
#include <dune/amdis/common/TupleUtility.hpp>
#include <dune/amdis/common/Utility.hpp>
#include <dune/amdis/linear_algebra/mtl/BlockMTLVector.hpp>
#include <dune/amdis/linear_algebra/mtl/DOFVector.hpp>
// for explicit instantiation
#include <dune/amdis/ProblemStatTraits.hpp>
namespace AMDiS
{
namespace Impl
{
// DOFVectors = std::tuple< DOFVector<FeSpace, ValueType>... >
template <class DOFVectors>
struct BuildDOFVectors
{
template <std::size_t I>
using DOFVectorIdx = std::tuple_element_t<I, DOFVectors>;
template <std::size_t... I, class FeSpaces, class MultiVector>
static DOFVectors make(Indices<I...>,
FeSpaces&& feSpaces,
std::vector<std::string> const& names,
MultiVector&& multiVector)
{
return DOFVectors{DOFVectorIdx<I>(
std::get<I>(std::forward<FeSpaces>(feSpaces)),
names[I],
multiVector[I]
)...};
}
};
// construction method to construct a tuple of DOFVectors
template <class DOFVectors, class FeSpaces, class MultiVector>
DOFVectors buildDOFVectors(FeSpaces&& feSpaces,
std::vector<std::string> const& names,
MultiVector&& multiVector)
{
return BuildDOFVectors<DOFVectors>::make(
MakeSeq_t<std::tuple_size<std::decay_t<FeSpaces>>::value>(),
std::forward<FeSpaces>(feSpaces),
names,
std::forward<MultiVector>(multiVector));
}
} // end namespace Impl
/// \brief Container that repesents multiple data-Vectors of different value types.
/**
* Represents a std::array of \ref mtl::dense_vector + a tuple of corresponding
* feSpaces. This I'th matrix combined with the I'th FeSpace
* builds a \ref DOFVector and can be return by the \ref operator().
**/
template <class GlobalBasis, class ValueType = double>
class SystemVector
{
using Self = SystemVector;
using Vector = mtl::dense_vector<ValueType>;
// NOTE: Vector must have a hierarchic structure, e.g. MultiTypeBlockVector< VelocityVectorType, PressureVectorType >
using Wrapper = Dune::HierarchicVectorWrapper<Vector>;
using DiscreteFunction = Dune::DiscreteGlobalBasisFunction<
GlobalBasis, Dune::TypeTree::HybridTreePath<>, Wrapper>;
template <std::size_t I>
using DOFVectorType = DOFVector<
GlobalBasis, Dune::TypeTree::HybridTreePath<index_t<I>>, ValueType>;
public:
/// Constructor that constructs new base-vectors
SystemVector(GlobalBasis const& globalBasis, std::vector<std::string> const& names)
: vector()
, wrapper(vector)
, globalFunction(Dune::makeDiscreteGlobalBasisFunction(globalBasis, wrapper))
, names(names)
{}
/// Return the name of the i'th DOFVector
std::string getName(std::size_t i) const
{
return names[i];
}
template <class MultiIndex>
auto const& operator[](MultiIndex const& index) const
{
return wrapper[index];
}
template <class MultiIndex>
auto& operator[](MultiIndex const& index)
{
return wrapper[index];
}
/// Return the I'th DOFVector
template <class... Indices>
DOFVectorType<I> getDOFVector(Indices const... indices)
{
auto treePath = Dune::TypeTree::HybridTreePath(indices...);
return {globalFunction.basis(), treePath, wrapper, getName(I)};
}
/// Return the I'th DOFVector
template <class... Indices>
DOFVectorType<I> const getDOFVector(Indices const... indices) const
{
auto treePath = Dune::TypeTree::HybridTreePath(indices...);
return {globalFunction.basis(), treePath, wrapper, getName(I)};
}
private:
Vector vector;
Wrapper wrapper;
DiscreteFunction globalFunction;
/// The names of the DOFVectors
std::vector<std::string> names;
};
#ifndef AMDIS_NO_EXTERN_SYSTEMVECTOR
// explicit instantiation in SystemVector.cpp
extern template class SystemVector<typename TestTraits<2,1,2>::FeSpaces>;
#endif
} // end namespace AMDiS
......@@ -20,10 +20,13 @@ namespace AMDiS
// TODO: possibly convert to plain type, in case of expression templates.
using value_type = ValueType;
ConstantTerm(value_type value)
explicit ConstantTerm(value_type value)
: value(value)
{}
template <class Element>
void bind(Element const&) {}
template <class Element, class PointList>
void init(Element const& element,
PointList const& points) { /* do nothing */ }
......@@ -33,7 +36,7 @@ namespace AMDiS
return value;
}
int getDegree(Dune::GeometryType const&) const
int getDegree() const
{
return 0;
}
......@@ -45,7 +48,7 @@ namespace AMDiS
/// generator function for \ref ConstantTerm expressions
template <class T>
ConstantTerm<T> constant(T value) { return {value}; }
auto constant(T value) { return ConstantTerm<T>{value}; }
/** @} **/
......
#pragma once
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <vector>
#include <type_traits>
......@@ -55,11 +51,14 @@ namespace AMDiS
template <class F,
class = std::enable_if_t<Concepts::Compatible<Functor, F>> >
CoordsFunctorTerm(F&& f, int degree = 1)
: fct(std::forward<F>(f))
, degree(degree)
explicit CoordsFunctorTerm(F&& f, int degree = 1)
: fct_(std::forward<F>(f))
, degree_(degree)
{}
template <class Element>
void bind(Element const&) {}
/// \brief Evaluates the functor \ref fct at all quadrature points and
/// stores the result in a vector \ref values.
template <class Element, class PointList>
......@@ -69,7 +68,7 @@ namespace AMDiS
auto geo = element.geometry();
values.resize(points.size());
for (std::size_t iq = 0; iq < points.size(); ++iq)
values[iq] = fct(geo.global(points[iq].position()));
values[iq] = fct_(geo.global(points[iq].position()));
}
value_type const& operator[](std::size_t iq) const
......@@ -77,14 +76,14 @@ namespace AMDiS
return values[iq];
}
int getDegree(Dune::GeometryType const&) const
int getDegree() const
{
return degree;
return degree_;
}
private:
Functor fct;
int degree;
Functor fct_;
int degree_;
std::vector<value_type> values;
};
......@@ -92,7 +91,7 @@ namespace AMDiS
/// Generator function for \ref CoordsFunctorTerm expressions
template <class F>
CoordsFunctorTerm< std::decay_t<F> > eval(F&& f) { return {std::forward<F>(f)}; }
auto eval(F&& f) { return CoordsFunctorTerm<std::decay_t<F>>{std::forward<F>(f)}; }
......@@ -107,6 +106,9 @@ namespace AMDiS
public:
using value_type = Dune::FieldVector<double, dow>;
template <class Element>
void bind(Element const&) {}
/// \brief Evaluates the functor \ref fct at all quadrature points and
/// stores the result in a vector \ref values.
template <class Entity, class PointList>
......@@ -124,7 +126,7 @@ namespace AMDiS
return values[iq];
}
int getDegree(Dune::GeometryType const&) const
int getDegree() const
{
return 1;
}
......
#pragma once
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <algorithm>
#include <vector>
#include <type_traits>
......@@ -21,77 +17,6 @@ namespace AMDiS
* @{
**/
#if 0
/// \brief An expression that evalues a DOFVector at a given element, either on the
/// dofs or on the quadrature points. \see valueOf.
template <class DOFVectorType>
class DOFVectorTerm
{
public:
using value_type = Value_t<DOFVectorType>;
DOFVectorTerm(DOFVectorType const& dofvector, double factor = 1.0)
: vector(dofvector.getVector())
, factor(factor)
, localView(dofvector.getFeSpace().localView())
, localIndexSet(dofvector.getFeSpace().localIndexSet())
{}
template <class Element, class PointList>
void init(Element const& element,
PointList const& points)
{
localView.bind(element);
localIndexSet.bind(localView);
const auto& localFiniteElem = localView.tree().finiteElement();
const std::size_t nBasisFct = localFiniteElem.size();
// update degree of basis functions
degree = localFiniteElem.localBasis().order();
std::vector<Dune::FieldVector<double,1> > shapeValues(nBasisFct);
std::vector<value_type> localVec(nBasisFct);
for (std::size_t j = 0; j < nBasisFct; ++j) {
const auto global_idx = localIndexSet.index(j);
localVec[j] = factor * vector[global_idx];
}
values.resize(points.size());
for (std::size_t iq = 0; iq < points.size(); ++iq) {
localFiniteElem.localBasis().evaluateFunction(points[iq].position(), shapeValues);
value_type result = 0;
for (std::size_t j = 0; j < shapeValues.size(); ++j)
result += localVec[j] * shapeValues[j];
values[iq] = result;
}
}
value_type const& operator[](std::size_t iq) const
{
return values[iq];
}
int getDegree(Dune::GeometryType const&) const
{
return degree;
}
private:
typename DOFVectorType::BaseVector const& vector;
double factor;
using Basis = typename DOFVectorType::FeSpace;
typename Basis::LocalView localView;
typename Basis::LocalIndexSet localIndexSet;
int degree = getPolynomialDegree<Basis>;
std::vector<value_type> values;
};
#endif
/// \brief An expression that evalues a DOFVector at a given element, either on the
/// dofs or on the quadrature points. \see valueOf.
template <class DiscreteGlobalBasisFct>
......@@ -100,26 +25,29 @@ namespace AMDiS
public:
using value_type = typename DiscreteGlobalBasisFct::Range;
DOFVectorTerm(std::shared_ptr<DiscreteGlobalBasisFct> const& dofvector, double factor = 1.0)
explicit DOFVectorTerm(std::shared_ptr<DiscreteGlobalBasisFct> const& dofvector, double factor = 1.0)
: dofvector_(dofvector)
, localfunction_(localFunction(*dofvector_))
, factor_(factor)
{
degree_ = getPolynomialDegree(dofvector->basis().localView().tree());
}
{}
DOFVectorTerm(DiscreteGlobalBasisFct& dofvector, double factor = 1.0)
explicit DOFVectorTerm(DiscreteGlobalBasisFct& dofvector, double factor = 1.0)
: localfunction_(localFunction(dofvector))
, factor_(factor)
{}
template <class Element>
void bind(Element const& element)
{
degree_ = getPolynomialDegree(dofvector.basis().localView().tree());
auto localView = dofvector_->basis().localView();
localView.bind(element);
localfunction_.bind(element);
degree_ = getPolynomialDegree(localView.tree());
}
template <class Element, class PointList>
void init(Element const& element, PointList const& points)
{
localfunction_.bind(element);
values_.resize(points.size());
for (std::size_t iq = 0; iq < points.size(); ++iq)
values_[iq] = factor_ * localfunction_(points[iq].position());
......@@ -132,7 +60,8 @@ namespace AMDiS
return values_[iq];
}
int getDegree(Dune::GeometryType const&) const
// [expects: dofvector_->basis().localView() is bound to an element]
int getDegree() const
{
return degree_;
}
......@@ -142,7 +71,7 @@ namespace AMDiS
typename DiscreteGlobalBasisFct::LocalFunction localfunction_;
double factor_;
int degree_;
int degree_ = 0;
std::vector<value_type> values_;
};
......@@ -152,10 +81,9 @@ namespace AMDiS
* points are scaled with.
**/
template <class DiscreteGlobalBasisFct>
DOFVectorTerm<DiscreteGlobalBasisFct>
valueOf(std::shared_ptr<DiscreteGlobalBasisFct> const& vector, double factor = 1.0)
auto valueOf(std::shared_ptr<DiscreteGlobalBasisFct> const& vector, double factor = 1.0)
{
return {vector, factor};
return DOFVectorTerm<DiscreteGlobalBasisFct>{vector, factor};
}
......@@ -175,29 +103,31 @@ namespace AMDiS
DOFVectorFuncTerm(DOFVectorType const& dofvector, F_&& f, int f_deg)
: vector(dofvector.getVector())
, f(std::forward<F_>(f))
, localView(dofvector.getFeSpace().localView())
, localIndexSet(dofvector.getFeSpace().localIndexSet())
, localView_(dofvector.getFeSpace().localView())
, localIndexSet_(dofvector.getFeSpace().localIndexSet())
, f_deg(f_deg)
{}
template <class Element>
void bind(Element const& element)
{
localView_.bind(element);
localIndexSet_.bind(localView_);
degree_ = getPolynomialDegree(localView_.tree());
}
template <class Element, class PointList>
void init(Element const& element,
PointList const& points)
{
localView.bind(element);
localIndexSet.bind(localView);
const auto& localFiniteElem = localView.tree().finiteElement();
const auto& localFiniteElem = localView_.tree().finiteElement();
const std::size_t nBasisFct = localFiniteElem.size();
// update degree of basis functions
degree = localFiniteElem.localBasis().order();
std::vector<Dune::FieldVector<double,1> > shapeValues(nBasisFct);
std::vector<value_type> localVec(nBasisFct);
for (std::size_t j = 0; j < nBasisFct; ++j) {
const auto global_idx = localIndexSet.index(j);
const auto global_idx = localIndexSet_.index(j);
localVec[j] = vector[global_idx];
}
......@@ -217,9 +147,9 @@ namespace AMDiS
return values[iq];
}
int getDegree(Dune::GeometryType const&) const
int getDegree() const
{
return degree * f_deg;
return degree_ * f_deg;
}
private:
......@@ -227,11 +157,11 @@ namespace AMDiS
Func f;
using Basis = typename DOFVectorType::FeSpace;
typename Basis::LocalView localView;
typename Basis::LocalIndexSet localIndexSet;
typename Basis::LocalView localView_;
typename Basis::LocalIndexSet localIndexSet_;
int degree = getPolynomialDegree<Basis>;
int f_deg;
int degree_ = 0;
int f_deg = 1;
std::vector<value_type> values;
};
......@@ -274,33 +204,39 @@ namespace AMDiS
using value_type = GradientGlobalType;
using field_type = Value_t<DOFVectorType>;
GradientTerm(DOFVectorType const& dofvector, double factor = 1.0)
explicit GradientTerm(DOFVectorType const& dofvector, double factor = 1.0)
: vector(dofvector.getVector())
, factor(factor)
, localView(dofvector.getFeSpace().localView())
, localIndexSet(dofvector.getFeSpace().localIndexSet())
, localView_(dofvector.getFeSpace().localView())
, localIndexSet_(dofvector.getFeSpace().localIndexSet())
{}
template <class Element>
void bind(Element const& element)
{
localView_.bind(element);
localIndexSet_.bind(localView_);
int deg = getPolynomialDegree(localView_.tree());
auto t = element.type();
degree_ = t.isSimplex() ? Math::max(0, deg-1) : deg;
}