Commit a6593546 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Merge branch 'issue/update_dune' into 'develop'

Issue/update dune

See merge request !27
parents 3d0d7673 b6a20add
Pipeline #1119 passed with stage
in 3 minutes and 4 seconds
......@@ -7,5 +7,5 @@ Module: amdis
Version: 0.1
Maintainer: simon.praetorius@tu-dresden.de
#depending on
Depends: dune-common dune-geometry dune-localfunctions dune-typetree dune-grid dune-functions
Depends: dune-common (>= 2.6) dune-geometry (>= 2.6) dune-localfunctions (>= 2.6) dune-typetree (>= 2.6) dune-grid (>= 2.6) dune-functions (>= 2.6)
#Suggests: dune-uggrid dune-alugrid dune-foamgrid
......@@ -55,8 +55,8 @@ int main(int argc, char** argv)
std::vector<double> errH1; errH1.reserve(numLevels);
std::vector<double> widths; widths.reserve(numLevels);
for (int i = 0; i < numLevels; ++i) {
prob.getGrid()->globalRefine(1);
auto gridView = prob.getGlobalBasis()->gridView();
prob.grid().globalRefine(1);
auto gridView = prob.gridView();
double h = 0;
for (auto const& e : edges(gridView))
......@@ -72,7 +72,7 @@ int main(int argc, char** argv)
errH1.push_back(std::sqrt(errorH1));
#if WRITE_FILES
Dune::VTKWriter<typename ElliptProblem::GridView> vtkWriter(prob.getGlobalBasis()->gridView());
Dune::VTKWriter<typename ElliptProblem::GridView> vtkWriter(gridView);
vtkWriter.addVertexData(prob.getSolution(_0), Dune::VTK::FieldInfo("u", Dune::VTK::FieldInfo::Type::scalar, 1));
vtkWriter.write("u_" + std::to_string(i));
#endif
......
......@@ -38,12 +38,12 @@ int main(int argc, char** argv)
auto not_left = [](auto const& x) { return x[0] > 1.0 - 1.e-8 || x[1] < 1.e-8 || x[1] > 1.0 - 1.e-8; };
// define boundary values
auto parabolic_y = [](auto const& x) -> Dune::FieldVector<double,2>
auto parabolic_y = [](auto const& x) -> FieldVector<double,2>
{
return {0.0, x[1]*(1.0 - x[1])};
};
auto zero = [](auto const& x) -> Dune::FieldVector<double,2>
auto zero = [](auto const& x) -> FieldVector<double,2>
{
return {0.0, 0.0};
};
......
......@@ -72,21 +72,6 @@ namespace AMDiS
bool asmMatrix, bool asmVector) const;
/// Return the element the LocalViews are bound to
template <class LocalView0, class... LovalViews>
auto const& getElement(LocalView0 const& localView, LovalViews const&...) const
{
return localView.element();
}
/// Return the gridView the localViews are bound to
template <class LocalView0, class... LovalViews>
auto const& getGridView(LocalView0 const& localView, LovalViews const&...) const
{
return globalBasis_.gridView();
}
private:
GlobalBasis& globalBasis_;
......
......@@ -20,7 +20,6 @@ void Assembler<Traits>::assemble(
initMatrixVector(matrix, solution, rhs, asmMatrix, asmVector);
auto localView = globalBasis_.localView();
auto localIndexSet = globalBasis_.localIndexSet();
// 2. create a local matrix and vector
std::size_t localSize = localView.maxSize();
......@@ -38,7 +37,7 @@ void Assembler<Traits>::assemble(
auto geometry = element.geometry();
// traverse type-tree of global-basis
AMDiS::forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
{
auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
auto rowLocalView = rowBasis.localView();
......@@ -56,7 +55,7 @@ void Assembler<Traits>::assemble(
this->assembleElementOperators(element, rhsOp, vecAssembler);
}
AMDiS::forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
{
auto& matOp = matrixOperators_[rowNode][colNode];
if (matOp.doAssemble(asmMatrix) && !matOp.empty()) {
......@@ -76,14 +75,12 @@ void Assembler<Traits>::assemble(
});
});
localIndexSet.bind(localView);
// add element-matrix to system-matrix
for (std::size_t i = 0; i < localView.size(); ++i) {
auto const row = localIndexSet.index(i);
auto const row = localView.index(i);
for (std::size_t j = 0; j < localView.size(); ++j) {
if (std::abs(elementMatrix(i,j)) > threshold<double>) {
auto const col = localIndexSet.index(j);
auto const col = localView.index(j);
matrix(row,col) += elementMatrix(i,j);
}
}
......@@ -92,20 +89,19 @@ void Assembler<Traits>::assemble(
// add element-vector to system-vector
for (std::size_t i = 0; i < localView.size(); ++i) {
if (std::abs(elementVector[i]) > threshold<double>) {
auto const idx = localIndexSet.index(i);
auto const idx = localView.index(i);
rhs[idx] += elementVector[i];
}
}
// unbind all operators
AMDiS::forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto&&) {
forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto&&) {
rhsOperators_[rowNode].unbind();
AMDiS::forEachNode_(localView.tree(), [&,this](auto const& colNode, auto&&) {
forEachNode_(localView.tree(), [&,this](auto const& colNode, auto&&) {
matrixOperators_[rowNode][colNode].unbind();
});
});
localIndexSet.unbind();
localView.unbind();
}
......@@ -155,7 +151,7 @@ void Assembler<Traits>::initMatrixVector(
rhs = 0;
auto localView = globalBasis_.localView();
AMDiS::forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
{
#ifdef HAVE_EXTENDED_DUNE_FUNCTIONS
if (rowNode.isLeaf)
......@@ -164,7 +160,7 @@ void Assembler<Traits>::initMatrixVector(
auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
AMDiS::forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
{
auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
......@@ -187,14 +183,14 @@ std::size_t Assembler<Traits>::finishMatrixVector(
matrix.finish();
auto localView = globalBasis_.localView();
AMDiS::forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
{
auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
auto& rhsOp = rhsOperators_[rowNode];
if (rhsOp.doAssemble(asmVector))
rhsOp.assembled = true;
AMDiS::forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTreePath)
{
auto colBasis = Dune::Functions::subspaceBasis(globalBasis_, colTreePath);
auto& matOp = matrixOperators_[rowNode][colNode];
......
......@@ -5,9 +5,13 @@
#include <type_traits>
#include <vector>
#include <dune/common/hybridutilities.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <amdis/Output.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/common/ValueCategory.hpp>
#include <amdis/linear_algebra/mtl/MTLDenseVector.hpp>
#include <amdis/utility/RangeType.hpp>
#include <amdis/utility/TreeData.hpp>
......@@ -75,7 +79,18 @@ namespace AMDiS
test_exit_dbg(initialized_, "Boundary condition not initialized!");
auto columns = matrix.applyDirichletBC(dirichletNodes_);
finishImpl(matrix, rhs, solution, rowBasis, colBasis, ValueCategory_t<Range>{});
using Dune::Functions::interpolate;
Dune::Hybrid::ifElse(std::is_same<RangeType_t<typename RowBasis::LocalView::Tree>, Range>{},
[&](auto id) {
auto rhsWrapper = wrapper(rhs.getVector());
interpolate(id(rowBasis), rhsWrapper, values_, dirichletNodes_);
});
Dune::Hybrid::ifElse(std::is_same<RangeType_t<typename ColBasis::LocalView::Tree>, Range>{},
[&](auto id) {
auto solutionWrapper = wrapper(solution.getVector());
interpolate(id(colBasis), solutionWrapper, values_, dirichletNodes_);
});
// subtract columns of dirichlet nodes from rhs
// for (auto const& triplet : columns)
......@@ -96,16 +111,6 @@ namespace AMDiS
template <class RB, class RowNodeTag>
void initImpl(RB const&, RowNodeTag) {}
template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis, class ValueCat>
void finishImpl(Matrix& matrix, VectorX& rhs, VectorB& solution,
RowBasis const& rowBasis, ColBasis const& colBasis,
ValueCat);
template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis>
void finishImpl(Matrix& matrix, VectorX& rhs, VectorB& solution,
RowBasis const& rowBasis, ColBasis const& colBasis,
tag::unknown) {}
private:
std::function<bool(Domain)> predicate_;
std::function<Range(Domain)> values_;
......
......@@ -4,6 +4,7 @@
#include <dune/functions/functionspacebases/subspacebasis.hh>
#include <amdis/LinearAlgebra.hpp>
#include <amdis/linear_algebra/HierarchicWrapper.hpp>
#include <amdis/linear_algebra/mtl/MTLDenseVector.hpp>
namespace AMDiS
{
......@@ -13,7 +14,7 @@ namespace AMDiS
RowBasis const& rowBasis, Dune::TypeTree::LeafNodeTag)
{
using Dune::Functions::interpolate;
interpolate(rowBasis, hierarchicVectorWrapper(dirichletNodes_), predicate_);
interpolate(rowBasis, dirichletNodes_, predicate_);
}
template <class WorldVector, class Range>
......@@ -46,19 +47,4 @@ namespace AMDiS
});
}
template <class WorldVector, class Range>
template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis, class ValueCat>
void DirichletBC<WorldVector, Range>::finishImpl(
Matrix& matrix, VectorX& solution, VectorB& rhs,
RowBasis const& rowBasis, ColBasis const& colBasis, ValueCat)
{
using Dune::Functions::interpolate;
interpolate(rowBasis, hierarchicVectorWrapper(rhs.getVector()), values_,
hierarchicVectorWrapper(dirichletNodes_));
interpolate(colBasis, hierarchicVectorWrapper(solution.getVector()), values_,
hierarchicVectorWrapper(dirichletNodes_));
}
} // end namespace AMDiS
......@@ -45,7 +45,7 @@ void ProblemInstat<Traits>::createUhOld()
if (oldSolution)
warning("oldSolution already created\n");
else // create oldSolution
oldSolution.reset(new SystemVector(*problemStat.getGlobalBasis(), name + "_uOld"));
oldSolution.reset(new SystemVector(problemStat.globalBasis(), name + "_uOld"));
}
......
......@@ -40,6 +40,10 @@
namespace AMDiS
{
// forward declaration
template <class Traits>
class ProblemInstat;
template <class Traits>
class ProblemStat
: public ProblemStatBase
......@@ -47,6 +51,8 @@ namespace AMDiS
{
using Self = ProblemStat;
friend class ProblemInstat<Traits>;
public: // typedefs and static constants
using GlobalBasis = typename Traits::GlobalBasis;
......@@ -108,27 +114,23 @@ namespace AMDiS
/// Adds an operator to \ref A.
/** @{ */
template <class Operator, class RowTreePath, class ColTreePath>
void addMatrixOperator(Operator const& op, RowTreePath, ColTreePath,
double* factor = nullptr, double* estFactor = nullptr);
template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
void addMatrixOperator(Operator const& op, RowTreePath = {}, ColTreePath = {});
// operator evaluated on the boundary
template <class Operator, class RowTreePath, class ColTreePath>
void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath, ColTreePath,
double* factor = nullptr, double* estFactor = nullptr);
template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath = {}, ColTreePath = {});
/** @} */
/// Adds an operator to \ref rhs.
/** @{ */
template <class Operator, class TreePath>
void addVectorOperator(Operator const& op, TreePath,
double* factor = nullptr, double* estFactor = nullptr);
template <class Operator, class TreePath = RootTreePath>
void addVectorOperator(Operator const& op, TreePath = {});
// operator evaluated on the boundary
template <class Operator, class TreePath>
void addVectorOperator(BoundaryType b, Operator const& op, TreePath,
double* factor = nullptr, double* estFactor = nullptr);
template <class Operator, class TreePath = RootTreePath>
void addVectorOperator(BoundaryType b, Operator const& op, TreePath = {});
/** @} */
......@@ -157,15 +159,18 @@ namespace AMDiS
public: // get-methods
/// Returns a pointer to system-matrix, \ref systemMatrix_
std::shared_ptr<SystemMatrix> getSystemMatrix() { return systemMatrix_; }
std::shared_ptr<SystemMatrix> getSystemMatrix() const { return systemMatrix_; }
/// Returns a reference to system-matrix, \ref systemMatrix_
SystemMatrix& getSystemMatrix() { return *systemMatrix_; }
SystemMatrix const& getSystemMatrix() const { return *systemMatrix_; }
/// Returns a reference to the solution vector, \ref solution_
SystemVector& getSolutionVector() { return *solution_; }
SystemVector const& getSolutionVector() const { return *solution_; }
/// Return a reference to the rhs system-vector, \ref rhs
SystemVector& getRhsVector() { return *rhs_; }
SystemVector const& getRhsVector() const { return *rhs_; }
/// Returns a pointer to the solution vector, \ref solution_
std::shared_ptr<SystemVector> getSolutionVector() const
{
return solution_;
}
/// Return a mutable view to a solution component
template <class TreePath = RootTreePath>
......@@ -184,36 +189,36 @@ namespace AMDiS
}
/// Return a point to the rhs system-vector, \ref rhs
std::shared_ptr<SystemVector> getRhsVector() { return rhs_; }
std::shared_ptr<SystemVector> getRhsVector() const { return rhs_; }
/// Return a pointer to the linear solver, \ref linearSolver
std::shared_ptr<LinearSolverType> getSolver() { return linearSolver_; }
/// Return a reference to the linear solver, \ref linearSolver
LinearSolverType& getSolver() { return *linearSolver_; }
LinearSolverType const& getSolver() const { return *linearSolver_; }
/// Set a new linear solver for the problem
void setSolver(std::shared_ptr<LinearSolverType> const& solver)
{
linearSolver_ = solver;
}
/// Return a pointer to the grid, \ref grid
std::shared_ptr<Grid> getGrid() { return grid_; }
/// Return a reference to the grid, \ref grid
Grid& grid() { return *grid_; }
Grid const& grid() const { return *grid_; }
/// Set the grid. Stores pointer to passed reference and initializes feSpaces
/// Set the grid. Stores pointer and initializes feSpaces
/// matrices and vectors, as well as the file-writer.
void setGrid(Grid& grid)
void setGrid(std::shared_ptr<Grid> const& grid)
{
grid_ = Dune::stackobject_to_shared_ptr(grid);
grid_ = grid;
createGlobalBasis();
createMatricesAndVectors();
createFileWriter();
}
/// Return the gridView of the leaf-level
GridView const& gridView() { return globalBasis_->gridView(); }
/// Return the \ref feSpaces
std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis_; }
GlobalBasis const& globalBasis() { return *globalBasis_; }
/// Implementation of \ref ProblemStatBase::getName
......@@ -224,6 +229,27 @@ namespace AMDiS
protected: // initialization methods
template <class T, class GV>
using HasCreate = decltype(T::create(std::declval<GV>()));
void createGlobalBasis()
{
createGlobalBasis(Dune::Std::is_detected<HasCreate,Traits,GridView>{});
initGlobalBasis(*globalBasis_);
}
void createGlobalBasis(std::true_type)
{
assert( bool(grid_) );
static_assert(std::is_same<GridView, typename Grid::LeafGridView>::value, "");
globalBasis_ = std::make_shared<GlobalBasis>(Traits::create(grid_->leafGridView()));
}
void createGlobalBasis(std::false_type)
{
error_exit("Cannot create GlobalBasis from type. Pass a BasisCreator instead!");
}
void createGrid()
{
gridName_ = "";
......@@ -239,15 +265,9 @@ namespace AMDiS
msg("");
}
void createGlobalBasis()
{
globalBasis_ = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(grid_->leafGridView()));
initGlobalBasis(*globalBasis_);
}
void initGlobalBasis(GlobalBasis const& globalBasis)
{
localView_ = std::make_shared<typename GlobalBasis::LocalView>(globalBasis.localView());
localView_ = std::make_shared<typename GlobalBasis::LocalView>(globalBasis_->localView());
matrixOperators_.init(localView_->tree(), tag::store{});
rhsOperators_.init(localView_->tree(), tag::store{});
constraints_.init(localView_->tree(), tag::store{});
......
......@@ -146,7 +146,7 @@ void ProblemStat<Traits>::createMarker()
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
{
AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
{
std::string componentName = name_ + "->output[" + to_string(treePath) + "]";
......@@ -165,8 +165,7 @@ template <class Traits>
template <class Operator, class RowTreePath, class ColTreePath>
void ProblemStat<Traits>::addMatrixOperator(
Operator const& preOp,
RowTreePath row, ColTreePath col,
double* factor, double* estFactor)
RowTreePath row, ColTreePath col)
{
static_assert( Concepts::PreTreePath<RowTreePath>,
"row must be a valid treepath, or an integer/index-constant");
......@@ -179,7 +178,7 @@ void ProblemStat<Traits>::addMatrixOperator(
auto op = makeLocalOperator<Element>(preOp, globalBasis_->gridView());
auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i, j);
matrixOperators_[i][j].element.push_back({localAssembler, factor, estFactor});
matrixOperators_[i][j].element.push_back({localAssembler, nullptr, nullptr});
matrixOperators_[i][j].changing = true;
}
......@@ -189,8 +188,7 @@ template <class Traits>
void ProblemStat<Traits>::addMatrixOperator(
BoundaryType b,
Operator const& preOp,
RowTreePath row, ColTreePath col,
double* factor, double* estFactor)
RowTreePath row, ColTreePath col)
{
static_assert( Concepts::PreTreePath<RowTreePath>,
"row must be a valid treepath, or an integer/index-constant");
......@@ -204,7 +202,7 @@ void ProblemStat<Traits>::addMatrixOperator(
auto op = makeLocalOperator<Intersection>(preOp, globalBasis_->gridView());
auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i, j);
matrixOperators_[i][j].boundary.push_back({localAssembler, factor, estFactor, b});
matrixOperators_[i][j].boundary.push_back({localAssembler, nullptr, nullptr, b});
matrixOperators_[i][j].changing = true;
}
......@@ -213,8 +211,7 @@ template <class Traits>
template <class Operator, class TreePath>
void ProblemStat<Traits>::addVectorOperator(
Operator const& preOp,
TreePath path,
double* factor, double* estFactor)
TreePath path)
{
static_assert( Concepts::PreTreePath<TreePath>,
"path must be a valid treepath, or an integer/index-constant");
......@@ -224,7 +221,7 @@ void ProblemStat<Traits>::addVectorOperator(
auto op = makeLocalOperator<Element>(preOp, globalBasis_->gridView());
auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i);
rhsOperators_[i].element.push_back({localAssembler, factor, estFactor});
rhsOperators_[i].element.push_back({localAssembler, nullptr, nullptr});
rhsOperators_[i].changing = true;
}
......@@ -234,8 +231,7 @@ template <class Traits>
void ProblemStat<Traits>::addVectorOperator(
BoundaryType b,
Operator const& preOp,
TreePath path,
double* factor, double* estFactor)
TreePath path)
{
static_assert( Concepts::PreTreePath<TreePath>,
"path must be a valid treepath, or an integer/index-constant");
......@@ -246,7 +242,7 @@ void ProblemStat<Traits>::addVectorOperator(
auto op = makeLocalOperator<Intersection>(preOp, globalBasis_->gridView());
auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i);
rhsOperators_[i].boundary.push_back({localAssembler, factor, estFactor, b});
rhsOperators_[i].boundary.push_back({localAssembler, nullptr, nullptr, b});
rhsOperators_[i].changing = true;
}
......
#pragma once
#include <tuple>
#include <dune/functions/functionspacebases/basistags.hh>
#include <dune/functions/functionspacebases/flatmultiindex.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/pqknodalbasis.hh>
#include <dune/functions/functionspacebases/pq1nodalbasis.hh>
#include <dune/grid/yaspgrid.hh>
#include <amdis/common/Mpl.hpp>
namespace AMDiS
{
namespace Impl
{
template <class GridView, int k, class MultiIndex>
struct LagrangeBasisSubFactory
{
using type = Dune::Functions::PQkNodeFactory<GridView, k, MultiIndex>;
};
template <class GridView, class MultiIndex>
struct LagrangeBasisSubFactory<GridView, 1, MultiIndex>
{
using type = Dune::Functions::PQ1NodeFactory<GridView, MultiIndex>;
};
template <class GridView, bool same, int... deg>
struct LagrangeBasisBuilderImpl;
template <bool same, int... degs>
struct LagrangePreBasisCreatorImpl;
// specialization for composite basis
template <class GridView, int... deg>
struct LagrangeBasisBuilderImpl<GridView, false, deg...>
template <int... degs>
struct LagrangePreBasisCreatorImpl<false, degs...>
{
using MultiIndex = Dune::ReservedVector<std::size_t,1>;
using IndexTag = Dune::Functions::BasisBuilder::FlatLexicographic;
using type = Dune::Functions::CompositeNodeFactory<MultiIndex, IndexTag,
typename LagrangeBasisSubFactory<GridView,deg,MultiIndex>::type...>;
static auto create()
{
using namespace Dune::Functions::BasisFactory;
return composite(lagrange<degs>()..., flatLexicographic());
}
};
// specialization for power basis
template <class GridView, int deg, int... degs>
struct LagrangeBasisBuilderImpl<GridView, true, deg, degs...>
template <int deg, int... degs>
struct LagrangePreBasisCreatorImpl<true, deg, degs...>
{
using MultiIndex = Dune::ReservedVector<std::size_t,1>;
using IndexTag = Dune::Functions::BasisBuilder::FlatLexicographic;
using type = Dune::Functions::PowerNodeFactory<MultiIndex, IndexTag,