Commit a1c2c4f1 authored by Praetorius, Simon's avatar Praetorius, Simon

new global assembler class finished, needs some more corrections

parent 5de54da3
......@@ -6,9 +6,10 @@
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/amdis/FiniteElementSpaces.hpp>
#include <dune/amdis/LinearAlgebra.hpp>
#include <dune/amdis/LocalFiniteElemenet.hpp>
#include <dune/amdis/common/Mpl.hpp>
#include <dune/amdis/common/TypeDefs.hpp>
namespace AMDiS
{
......@@ -27,6 +28,9 @@ namespace AMDiS
template <class T>
using VectorEntries = Dune::FieldVector<T, nComponents>;
using ElementVector = Impl::ElementVector;
using ElementMatrix = Impl::ElementMatrix;
public:
/// Constructor, stores a shared-pointer to the feSpaces
Assembler(std::shared_ptr<FeSpaces> const& feSpaces_)
......@@ -42,92 +46,103 @@ namespace AMDiS
/// Assemble the linear system
template <class Operators>
void assemble(SystemMatrixType& matrix,
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators,
bool asmMatrix_, bool asmVector_);
void assemble(
SystemMatrixType& matrix,
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators,
bool asmMatrix_, bool asmVector_);
private:
// Return whether the matrix-block needs to be assembled
template <int R, int C, class Operators>
bool assembleMatrix(const index_t<R>, const index_t<C>,
MatrixEntries<Operators> const& matrix_operators, bool asmMatrix_) const
{
return asmMatrix_ && (!matrix_operators[R][C].assembled || matrix_operators[R][C].changing);
}
// Return whether the vector-block needs to be assembled
template <int R, class Operators>
bool assembleVector(const index_t<R>, VectorEntries<Operators> const& rhs_operators, bool asmVector_) const
{
return asmVector_ && (!rhs_operators[R].assembled || rhs_operators[R].changing);
}
// Sets the system to zero and initializes all operators and boundary conditions
/// Sets the system to zero and initializes all operators and boundary conditions
template <class Operators>
void initMatrixVector(SystemMatrixType& matrix,
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators,
bool asmMatrix_, bool asmVector_);
void initMatrixVector(
SystemMatrixType& matrix,
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators,
bool asmMatrix_, bool asmVector_);
/// Assemble a block-matrix of element-matrices and return a matrix of flags, whether
/// a given block has received some entries.
template <class Operators>
MatrixEntries<bool> assembleElementMatrices(MatrixEntries<Operators>& operators,
MatrixEntries<ElementMatrix>& elementMatrix,
LocalFiniteElement<FeSpaces> const& localFiniteElem,
bool asmMatrix_);
MatrixEntries<bool> assembleElementMatrices(
MatrixEntries<Operators>& operators,
MatrixEntries<ElementMatrix>& elementMatrix,
FiniteElementSpaces<FeSpaces> const& elementBases,
bool asmMatrix_);
/// Assemble a block-vector of element-vectors and return a vector of flags, whether
/// a given block has received some entries.
template <class Operators>
VectorEntries<bool> assembleElementVectors(VectorEntries<Operators>& operators,
VectorEntries<ElementVector>& elementVector,
LocalFiniteElement<FeSpaces> const& localFiniteElem,
bool asmVector_);
VectorEntries<bool> assembleElementVectors(
VectorEntries<Operators>& operators,
VectorEntries<ElementVector>& elementVector,
FiniteElementSpaces<FeSpaces> const& elementBases,
bool asmVector_);
/// Assemble one block of the block-element-matrix
// The MatrixData argument stores all matrix-operators
template <class Operators, class RowView, class ColView>
bool assembleElementMatrix(Operators& operators,
ElementMatrix& elementMatrix,
RowView const& rowLocalView, ColView const& colLocalView);
bool assembleElementMatrix(
Operators& operators,
ElementMatrix& elementMatrix,
RowView const& rowLocalView,
ColView const& colLocalView);
/// Assemble one block of the block-element-vector
// The VectorData argument stores all vector-operators
template <class Operators, class RowView>
bool assembleElementVector(Operators& operators,
ElementVector& elementVector,
RowView const& rowLocalView);
bool assembleElementVector(
Operators& operators,
ElementVector& elementVector,
RowView const& rowLocalView);
/// Add the block-element-matrix to the system-matrix
void addElementMatrices(SystemMatrixType& dofmatrix,
LocalFiniteElement<FeSpaces> const& localFiniteElem,
MatrixEntries<bool> const& addMat,
MatrixEntries<ElementMatrix> const& elementMatrix);
void addElementMatrices(
SystemMatrixType& dofmatrix,
FiniteElementSpaces<FeSpaces> const& elementBases,
MatrixEntries<bool> const& addMat,
MatrixEntries<ElementMatrix> const& elementMatrix);
/// Add the block-element-vector to the system-vector
void addElementVectors(SystemVectorType& dofvector,
LocalFiniteElement<FeSpaces> const& localFiniteElem,
VectorEntries<bool> const& addVec,
VectorEntries<ElementVector> const& elementVector);
void addElementVectors(
SystemVectorType& dofvector,
FiniteElementSpaces<FeSpaces> const& elementBases,
VectorEntries<bool> const& addVec,
VectorEntries<ElementVector> const& elementVector);
/// Finish insertion into the matrix and assembles boundary conditions
/// Return the number of nonzeros assembled into the matrix
template <class Operators>
std::size_t finishMatrixVector(SystemMatrixType& matrix,
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators,
bool asmMatrix_, bool asmVector_);
std::size_t finishMatrixVector(
SystemMatrixType& matrix,
SystemVectorType& solution,
SystemVectorType& rhs,
MatrixEntries<Operators>& matrix_operators,
VectorEntries<Operators>& rhs_operators,
bool asmMatrix_, bool asmVector_);
/// Return whether the matrix-block needs to be assembled
template <std::size_t R, std::size_t C, class Operators>
bool assembleMatrix(const index_t<R>, const index_t<C>,
MatrixEntries<Operators> const& matrix_operators, bool asmMatrix_) const
{
return asmMatrix_ && (!matrix_operators[R][C].assembled || matrix_operators[R][C].changing);
}
/// Return whether the vector-block needs to be assembled
template <std::size_t R, class Operators>
bool assembleVector(const index_t<R>, VectorEntries<Operators> const& rhs_operators, bool asmVector_) const
{
return asmVector_ && (!rhs_operators[R].assembled || rhs_operators[R].changing);
}
private:
std::shared_ptr<FeSpaces> feSpaces;
......
This diff is collapsed.
......@@ -21,25 +21,6 @@ add_dune_alberta_flags("duneamdis" OBJECT USE_GENERIC)
target_compile_definitions("duneamdis" PUBLIC AMDIS_BACKEND_MTL=1)
target_compile_options("duneamdis" PUBLIC -ftemplate-backtrace-limit=0)
# set(BOOST_VERSION "1.54")
# set(BOOST_LIBS_REQUIRED system program_options)
# if (NOT BUILD_SHARED_LIBS)
# set(Boost_USE_STATIC_LIBS ON)
# endif (NOT BUILD_SHARED_LIBS)
# find_package(Boost ${BOOST_VERSION} REQUIRED ${BOOST_LIBS_REQUIRED})
# if (Boost_FOUND)
# add_library(boost INTERFACE)
# target_include_directories(boost INTERFACE ${Boost_INCLUDE_DIR})
# target_link_libraries(boost INTERFACE ${Boost_LIBRARIES})
# target_link_libraries(duneamdis INTERFACE boost)
# if (MSVC_SHARED_LIBS)
# link_directories(${Boost_LIBRARY_DIRS})
# target_compile_definitions("duneamdis" INTERFACE ${Boost_LIB_DIAGNOSTIC_DEFINITIONS})
# endif (MSVC_SHARED_LIBS)
# endif (Boost_FOUND)
find_package(MTL REQUIRED
PATHS /usr/local/lib/mtl4)
if (MTL_FOUND)
......
#pragma once
#include <tuple>
#include <dune/amdis/common/TupleUtility.hpp>
#include <dune/amdis/common/IndexSeq.hpp>
namespace AMDiS
{
namespace Impl
{
template <template<class> class Base, class Tuple, class Indices> struct MakeTupleType;
template <template<class> class Base, class Tuple, std::size_t... I>
struct MakeTupleType<Base, Tuple, Indices<I...>>
{
using type = std::tuple<typename Base<std::tuple_element_t<I, Tuple>>::type...>;
};
/// Constructs tuple type, by wrapping Base around the tuple elements.
/// Returns tuple type tuple<Base<tuple_element_t<0>>::type, Base<tuple_element_t<1>>::type, ...>
template <template<class> class Base, class Tuple>
using MakeTupleType_t =
typename MakeTupleType<Base, Tuple, MakeSeq_t<std::tuple_size<Tuple>::value>>::type;
template <class Functor, class Tuple, std::size_t... I>
auto makeTupleType(Functor f, Tuple&& tuple, Indices<I...>)
{
return std::make_tuple(f(std::get<I>(std::forward<Tuple>(tuple>)))...);
}
// construction method to construct a tuple of DOFVectors
template <class Functor, class Tuple>
auto makeTupleType(Functor f, Tuple&& tuple)
{
return makeTupleType(f, std::forward<Tuple>(tuple), MakeSeq_t<std::tuple_size<Tuple>::value>{});
}
} // end namespace Impl
template <class FeSpace>
struct LocalView
{
......@@ -52,15 +19,17 @@ namespace AMDiS
template <class FeSpaces>
class LocalFiniteElement
class FiniteElementSpaces
{
template <std::size_t I>
using FeSpace = std::tuple_element_t<I, FeSpaces>;
static constexpr int nComponents = std::tuple_size<FeSpaces>::value;
using LocalViews = Impl::MakeTupleType_t<LocalView, FeSpaces>;
using LocalIndexSets = Impl::MakeTupleType_t<LocalIndexSet, FeSpaces>;
static_assert( nComponents > 0, "" );
using LocalViews = MapTuple_t<LocalView, FeSpaces>;
using LocalIndexSets = MapTuple_t<LocalIndexSet, FeSpaces>;
/// The grid view the global FE basis lives on
using GridView = typename FeSpace<0>::GridView;
......@@ -69,9 +38,10 @@ namespace AMDiS
using Element = typename GridView::template Codim<0>::Entity;
public:
LocalFiniteElement(FeSpaces const& feSpaces)
: localViews(Impl::makeTupleType([](auto const& basis) { return basis.localView(); }, feSpaces))
, localIndexSets(Impl::makeTupleType([](auto const& basis) { return basis.localIndexSet(); }, feSpaces))
FiniteElementSpaces(FeSpaces const& feSpaces)
: localViews_(mapTuple([](auto const& basis) { return basis.localView(); }, feSpaces))
, localIndexSets_(mapTuple([](auto const& basis) { return basis.localIndexSet(); }, feSpaces))
, gridView_(std::get<0>(feSpaces).gridView())
{}
/// Bind the \ref localViews and \ref localIndexSets to a grid element
......@@ -80,8 +50,8 @@ namespace AMDiS
forEach(range_<0, nComponents>, [&,this](auto const _i)
{
static const int I = decltype(_i)::value;
std::get<I>(localViews).bind(element);
std::get<I>(localIndexSets).bind(localView);
std::get<I>(localViews_).bind(element);
std::get<I>(localIndexSets_).bind(std::get<I>(localViews_));
});
// NOTE: maybe create element-geometry here
}
......@@ -92,29 +62,36 @@ namespace AMDiS
forEach(range_<0, nComponents>, [&,this](auto const _i)
{
static const int I = decltype(_i)::value;
std::get<I>(localIndexSets).unbind();
std::get<I>(localViews).unbind();
std::get<I>(localIndexSets_).unbind();
std::get<I>(localViews_).unbind();
});
}
template <std::size_t I>
auto const& localView(const index_t<I> _i = {}) const
{
return std::get<I>(localViews);
return std::get<I>(localViews_);
}
template <std::size_t I>
auto const& localIndexSet(const index_t<I> _i = {}) const
{
return std::get<I>(localIndexSets);
return std::get<I>(localIndexSets_);
}
auto const& gridView() const
{
return gridView_;
}
private:
/// tuple of localView objects, obtained from the tuple of global bases
LocalViews localViews;
LocalViews localViews_;
// tuple of localIndexSet objects, obtained from the tuple of global bases
LocalIndexSets localIndexSets;
LocalIndexSets localIndexSets_;
GridView gridView_;
};
} // end namespace AMDiS
......@@ -119,6 +119,9 @@ namespace AMDiS
void init(RowFeSpace const& rowFeSpace,
ColFeSpace const& colFeSpace);
template <class RowFeSpace>
void init(RowFeSpace const& rowFeSpace);
/// Calculates the needed quadrature degree for the given term-order \p order.
/// This takes into account the polynomial degree of trial and test functions
......
......@@ -24,6 +24,18 @@ namespace AMDiS
}
template <class MeshView, class Element>
template <class RowFeSpace>
void Operator<MeshView, Element>::
init(RowFeSpace const& rowFeSpace)
{
using IdType = typename Operator<MeshView, Element>::IdType;
lastVectorId = std::numeric_limits<IdType>::max();
psiDegree = getPolynomialDegree<RowFeSpace>;
}
template <class MeshView, class Element>
template <class Geometry>
int Operator<MeshView, Element>::
......
......@@ -14,6 +14,7 @@
#include <dune/grid/common/grid.hh>
#include <dune/amdis/AdaptInfo.hpp>
#include <dune/amdis/Assembler.hpp>
#include <dune/amdis/CreatorInterface.hpp>
#include <dune/amdis/DirichletBC.hpp>
#include <dune/amdis/FileWriter.hpp>
......@@ -52,9 +53,6 @@ namespace AMDiS
using Codim0 = typename MeshView::template Codim<0>;
using Element = typename Codim0::Entity;
using ElementVector = Impl::ElementVector;
using ElementMatrix = Impl::ElementMatrix;
/// Dimension of the mesh
static constexpr int dim = Traits::dim;
......@@ -80,17 +78,6 @@ namespace AMDiS
using LinearSolverType = LinearSolverInterface<typename SystemMatrixType::MultiMatrix,
typename SystemVectorType::MultiVector>;
/// Helper class to store an operator with optionally factors and
/// a boundary indicator
template <class OperatorType>
struct Scaled
{
std::shared_ptr<OperatorType> op;
double* factor = nullptr;
double* estFactor = nullptr;
BoundaryType b = {0};
};
public:
/**
* \brief Constructor. Takes the name of the problem that is used to
......@@ -134,8 +121,10 @@ namespace AMDiS
/// Adds an operator to \ref A.
/** @{ */
template <class OperatorType>
void addMatrixOperator(OperatorType const& op, int i, int j,
void addMatrixOperator(ElementOperator const& op, int i, int j,
double* factor = nullptr, double* estFactor = nullptr);
void addMatrixOperator(IntersectionOperator const& op, int i, int j,
double* factor = nullptr, double* estFactor = nullptr);
// operator evaluated on the boundary
......@@ -149,8 +138,10 @@ namespace AMDiS
/// Adds an operator to \ref rhs.
/** @{ */
template <class OperatorType>
void addVectorOperator(OperatorType const& op, int i,
void addVectorOperator(ElementOperator const& op, int i,
double* factor = nullptr, double* estFactor = nullptr);
void addVectorOperator(IntersectionOperator const& op, int i,
double* factor = nullptr, double* estFactor = nullptr);
// operator evaluated on the boundary
......@@ -161,6 +152,7 @@ namespace AMDiS
void addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops);
/** @} */
/// Adds a Dirichlet boundary condition
template <class Predicate, class Values>
void addDirichletBC(Predicate const& predicate,
......@@ -364,39 +356,6 @@ namespace AMDiS
}
protected: // sub-methods to assemble DOFMatrix
template <class LhsData, class RhsData, class GV>
void assemble(LhsData lhs, RhsData rhs, GV const& gridView);
template <class RowView, class ColView>
bool getElementMatrix(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
std::list<Scaled<ElementOperator>>& operators,
std::list<Scaled<IntersectionOperator>>& boundary_operators,
std::list<Scaled<IntersectionOperator>>& intersection_operators);
template <class RowView>
bool getElementVector(RowView const& rowView,
ElementVector& elementVector,
std::list<Scaled<ElementOperator>>& operators,
std::list<Scaled<IntersectionOperator>>& boundary_operators,
std::list<Scaled<IntersectionOperator>>& intersection_operators);
template <class Matrix, class RowIndexSet, class ColIndexSet>
void addElementMatrix(Matrix& matrix,
RowIndexSet const& rowIndexSet,
ColIndexSet const& colIndexSet,
ElementMatrix const& elementMatrix);
template <class Vector, class IndexSet>
void addElementVector(Vector& vector,
IndexSet const& indexSet,
ElementVector const& elementvector);
private:
/// Name of this problem.
std::string name;
......
......@@ -7,9 +7,10 @@
namespace AMDiS {
template <class Traits>
void ProblemStat<Traits>::initialize(Flag initFlag,
Self* adoptProblem,
Flag adoptFlag)
void ProblemStat<Traits>::initialize(
Flag initFlag,
Self* adoptProblem,
Flag adoptFlag)
{
// create meshes
if (mesh) {
......@@ -155,24 +156,24 @@ template <class Traits>
void ProblemStat<Traits>::
addVectorOperator(ElementOperator const& op, int i, double* factor, double* estFactor)
{
vectorOperators[i].element.push_back({std::make_shared<ElementOperator>(op), factor, estFactor});
vectorOperators[i].changing = true;
rhsOperators[i].element.push_back({std::make_shared<ElementOperator>(op), factor, estFactor});
rhsOperators[i].changing = true;
}
template <class Traits>
void ProblemStat<Traits>::
addVectorOperator(IntersectionOperator const& op, int i, double* factor, double* estFactor)
{
vectorOperators[i].intersection.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor});
vectorOperators[i].changing = true;
rhsOperators[i].intersection.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor});
rhsOperators[i].changing = true;
}
template <class Traits>
void ProblemStat<Traits>::
addVectorOperator(BoundaryType b, IntersectionOperator const& op, int i, double* factor, double* estFactor)
{
vectorOperators[i].boundary.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor, b});
vectorOperators[i].changing = true;
rhsOperators[i].boundary.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor, b});
rhsOperators[i].changing = true;
}
template <class Traits>
......@@ -180,8 +181,8 @@ void ProblemStat<Traits>::
addVectorOperator(std::map< int, ElementOperator> ops)
{
for (auto op : ops) {
vectorOperators[op.first].element.push_back({std::make_shared<ElementOperator>(op.second)});
vectorOperators[op.first].changing = true;
rhsOperators[op.first].element.push_back({std::make_shared<ElementOperator>(op.second)});
rhsOperators[op.first].changing = true;
}
}
......@@ -190,8 +191,8 @@ void ProblemStat<Traits>::
addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops)
{
for (auto op : ops) {
vectorOperators[op.first].element.push_back({std::make_shared<ElementOperator>(op.second.first)});
vectorOperators[op.first].changing = vectorOperators[op.first].changing || op.second.second;
rhsOperators[op.first].element.push_back({std::make_shared<ElementOperator>(op.second.first)});
rhsOperators[op.first].changing = rhsOperators[op.first].changing || op.second.second;
}
}
......
......@@ -166,4 +166,48 @@ namespace AMDiS
using Repeat_t =
typename Impl::RepeatedTuple<T, MakeSeq_t<N>>::type;
namespace Impl
{
template <template<class> class Base, class Tuple, class Indices> struct MapTuple;
template <template<class> class Base, class Tuple, std::size_t... I>
struct MapTuple<Base, Tuple, Indices<I...>>
{
using type = std::tuple<typename Base<std::tuple_element_t<I, Tuple>>::type...>;
};
} // end namespace Impl
/// Constructs tuple type, by wrapping Base around the tuple elements.
/// Returns tuple type tuple<Base<tuple_element_t<0>>::type, Base<tuple_element_t<1>>::type, ...>
template <template<class> class Base, class Tuple>
using MapTuple_t =
typename Impl::MapTuple<Base, Tuple, MakeSeq_t<std::tuple_size<Tuple>::value>>::type;
template <class Functor, class Tuples, std::size_t I, std::size_t... J>
auto mapTuple(Functor f, Tuples&& tuples, index_t<I>, Indices<J...>)
{
return f(std::get<I>(std::get<J>(std::forward<Tuples>(tuples)))...);
}
template <class Functor, class Tuples, std::size_t... I>
auto mapTuple(Functor f, Tuples&& tuples, Indices<I...>)
{
return std::make_tuple(mapTuple(f, std::forward<Tuples>(tuples), index_<I>,
MakeSeq_t<std::tuple_size<std::decay_t<Tuples>>::value>{})...);
}
// construction method to construct a new tuple by applying a functor to the elements
// of a given tuple
template <class Functor, class Tuple0, class... Tuples>
auto mapTuple(Functor f, Tuple0&& tuple0, Tuples&&... tuples)
{
return mapTuple(f, std::make_tuple(std::forward<Tuple0>(tuple0), std::forward<Tuples>(tuples)...),
MakeSeq_t<std::tuple_size<std::decay_t<Tuple0>>::value>{});
}
} // end namespace AMDiS
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