Commit 384b27d4 authored by Praetorius, Simon's avatar Praetorius, Simon

compatibility with dune 2.6.0 of all dependend modules

parent c7869272
......@@ -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
......@@ -36,15 +36,11 @@ int main(int argc, char** argv)
return {-20.0 * std::exp(-10.0 * dot(x,x)) * x};
};
using Grid = Dune::YaspGrid<AMDIS_DIM>;
auto gridPtr = MeshCreator<Grid>::create("mesh");
auto basis = makeBasis(gridPtr->leafGridView(), lagrange<1>());
auto prob = makeProblemStat("ellipt", *gridPtr, basis);
ElliptProblem prob("ellipt");
prob.initialize(INIT_ALL);
auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
prob.addMatrixOperator(opL);
prob.addMatrixOperator(opL, _0, _0);
auto opForce = makeOperator(tag::test{}, f, 6);
prob.addVectorOperator(opForce, _0);
......@@ -60,7 +56,7 @@ int main(int argc, char** argv)
std::vector<double> widths; widths.reserve(numLevels);
for (int i = 0; i < numLevels; ++i) {
prob.getGrid()->globalRefine(1);
auto gridView = prob.getGlobalBasis()->gridView();
auto gridView = prob.gridView();
double h = 0;
for (auto const& e : edges(gridView))
......@@ -76,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
......
......@@ -13,44 +13,44 @@
using namespace AMDiS;
// 1 component with polynomial degree 1
//using Grid = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;
using HeatParam = YaspGridBasis<AMDIS_DIM, 2>;
using HeatProblem = ProblemStat<HeatParam>;
using HeatProblemInstat = ProblemInstat<HeatParam>;
int main(int argc, char** argv)
{
AMDiS::init(argc, argv);
using Grid = Dune::YaspGrid<AMDIS_DIM>;
auto gridPtr = MeshCreator<Grid>::create("mesh");
auto basis = makeBasis(gridPtr->leafGridView(), lagrange<2>());
auto prob = makeProblemStat("heat", *gridPtr, basis);
HeatProblem prob("heat");
prob.initialize(INIT_ALL);
auto probInstat = makeProblemInstat("heat", prob);
HeatProblemInstat probInstat("heat", prob);
probInstat.initialize(INIT_UH_OLD);
AdaptInfo adaptInfo("adapt");
double* invTau = probInstat.getInvTau();
auto _ = Dune::TypeTree::hybridTreePath();
auto opTimeLhs = makeOperator(tag::test_trial{}, std::ref(*invTau));
prob.addMatrixOperator(opTimeLhs, _, _);
prob.addMatrixOperator(opTimeLhs, 0, 0);
auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
prob.addMatrixOperator(opL, _, _);
prob.addMatrixOperator(opL, 0, 0);
auto opTimeRhs = makeOperator(tag::test{},
invokeAtQP([invTau](double u) { return u * (*invTau); }, prob.getSolution()), 2);
prob.addVectorOperator(opTimeRhs, _);
invokeAtQP([invTau](double u) { return u * (*invTau); }, prob.getSolution(0)), 2);
prob.addVectorOperator(opTimeRhs, 0);
auto opForce = makeOperator(tag::test{}, [](auto const& x) { return -1.0; }, 0);
prob.addVectorOperator(opForce, _);
prob.addVectorOperator(opForce, 0);
// set boundary condition
auto predicate = [](auto const& p){ return p[0] < 1.e-8 || p[1] < 1.e-8; };
auto dbcValues = [](auto const& p){ return 0.0; };
prob.addDirichletBC(predicate, _, _, dbcValues);
prob.addDirichletBC(predicate, 0, 0, dbcValues);
AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
adapt.adapt();
......
......@@ -12,20 +12,23 @@
using namespace AMDiS;
int main(int argc, char** argv)
struct NavierStokesBasis
{
AMDiS::init(argc, argv);
using Grid = Dune::YaspGrid<AMDIS_DIM>;
auto gridPtr = MeshCreator<Grid>::create("mesh");
using GlobalBasis = typename TaylorHoodBasis<Grid::LeafGridView>::GlobalBasis;
};
using StokesProblem = ProblemStat<NavierStokesBasis>;
using StokesProblemInstat = ProblemInstat<NavierStokesBasis>;
auto preBasis = compose(repeat<AMDIS_DOW>(lagrange<2>()), lagrange<1>());
auto basis = makeBasis(gridPtr->leafGridView(), preBasis);
int main(int argc, char** argv)
{
AMDiS::init(argc, argv);
auto prob = makeProblemStat("stokes", *gridPtr, basis);
StokesProblem prob("stokes");
prob.initialize(INIT_ALL);
auto probInstat = makeProblemInstat("stokes", prob);
StokesProblemInstat probInstat("stokes", prob);
probInstat.initialize(INIT_UH_OLD);
double viscosity = 1.0;
......
......@@ -13,16 +13,17 @@
using namespace AMDiS;
// 3 components: velocity with polynomial degree 2 and pressure with polynomial degree 1
using Grid = Dune::YaspGrid<AMDIS_DIM>;
using StokesParam = TaylorHoodBasis<Grid::LeafGridView>;
using StokesProblem = ProblemStat<StokesParam>;
int main(int argc, char** argv)
{
AMDiS::init(argc, argv);
using Grid = Dune::YaspGrid<AMDIS_DIM>;
auto gridPtr = MeshCreator<Grid>::create("mesh");
auto preBasis = compose(repeat<AMDIS_DOW>(lagrange<2>()), lagrange<1>());
auto basis = makeBasis(gridPtr->leafGridView(), preBasis);
auto prob = makeProblemStat("stokes", *gridPtr, basis);
StokesProblem prob("stokes");
prob.initialize(INIT_ALL);
double viscosity = 1.0;
......
......@@ -13,17 +13,17 @@
using namespace AMDiS;
// 3 components: velocity with polynomial degree 2 and pressure with polynomial degree 1
using Grid = Dune::YaspGrid<AMDIS_DIM>;
using StokesParam = TaylorHoodBasis<Grid::LeafGridView>;
using StokesProblem = ProblemStat<StokesParam>;
int main(int argc, char** argv)
{
AMDiS::init(argc, argv);
using Grid = Dune::YaspGrid<AMDIS_DIM>;
auto gridPtr = MeshCreator<Grid>::create("mesh");
auto preBasis = compose(repeat<AMDIS_DOW>(lagrange<2>()), lagrange<1>());
auto basis = makeBasis(gridPtr->leafGridView(), preBasis);
auto prob = makeProblemStat("stokes", *gridPtr, basis);
StokesProblem prob("stokes");
prob.initialize(INIT_ALL);
double viscosity = 1.0;
......
......@@ -9,17 +9,15 @@
using namespace AMDiS;
using Grid = Dune::YaspGrid<AMDIS_DIM>;
using StokesParam = TaylorHoodBasis<Grid::LeafGridView>;
using StokesProblem = ProblemStat<StokesParam>;
int main(int argc, char** argv)
{
AMDiS::init(argc, argv);
using Grid = Dune::YaspGrid<AMDIS_DIM>;
auto gridPtr = MeshCreator<Grid>::create("mesh");
auto preBasis = compose(repeat<AMDIS_DOW>(lagrange<2>()), lagrange<1>());
auto basis = makeBasis(gridPtr->leafGridView(), preBasis);
auto prob = makeProblemStat("stokes", *gridPtr, basis);
StokesProblem prob("stokes");
prob.initialize(INIT_ALL);
double viscosity = 1.0;
......@@ -40,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};
};
......
......@@ -6,17 +6,16 @@
using namespace AMDiS;
// 1 component with polynomial degree 1
//using Grid = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;
using ElliptParam = YaspGridBasis<AMDIS_DIM, 2,2>;
using ElliptProblem = ProblemStat<ElliptParam>;
int main(int argc, char** argv)
{
AMDiS::init(argc, argv);
using Grid = Dune::YaspGrid<AMDIS_DIM>;
auto gridPtr = MeshCreator<Grid>::create("mesh");
auto preBasis = repeat<2>(lagrange<2>());
auto basis = makeBasis(gridPtr->leafGridView(), preBasis);
auto prob = makeProblemStat("stokes", *gridPtr, basis);
ElliptProblem prob("ellipt");
prob.initialize(INIT_ALL);
AdaptInfo adaptInfo("adapt");
......
......@@ -37,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();
......@@ -55,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()) {
......@@ -95,9 +95,9 @@ void Assembler<Traits>::assemble(
}
// 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();
});
});
......@@ -151,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)
......@@ -160,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);
......@@ -183,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"));
}
......
......@@ -67,7 +67,6 @@ namespace AMDiS
using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
public:
#if 0
/**
* \brief Constructor. Takes the name of the problem that is used to
* access values corresponding to this problem in the parameter file.
......@@ -94,7 +93,6 @@ namespace AMDiS
globalBasis_ = Dune::stackobject_to_shared_ptr(globalBasis);
initGlobalBasis(*globalBasis_);
}
#endif
/// \brief Constructor taking a grid reference and a basis reference.
/// Stores pointers to both.
......@@ -211,7 +209,6 @@ namespace AMDiS
/// Return a pointer to the grid, \ref grid
std::shared_ptr<Grid> getGrid() { return grid_; }
#if 0
/// Set the mesh. Stores pointer to passed reference and initializes feSpaces
/// matrices and vectors, as well as the file-writer.
void setGrid(Grid& grid)
......@@ -222,13 +219,12 @@ namespace AMDiS
createMatricesAndVectors();
createFileWriter();
}
#endif
/// Return the gridView of the leaf-level
auto const& gridView() { return globalBasis_->gridView(); }
/// Return the \ref feSpaces
std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis_; }
std::shared_ptr<GlobalBasis> const& globalBasis() { return globalBasis_; }
/// Implementation of \ref ProblemStatBase::getName
......@@ -239,6 +235,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("Can not create GlobalBasis from type. Pass a BasisCreator instead!");
}
void createGrid()
{
gridName_ = "";
......
......@@ -57,12 +57,10 @@ void ProblemStat<Traits>::initialize(
warning("globalBasis already created");
}
else {
#if 0
if (initFlag.isSet(INIT_FE_SPACE) ||
(initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
createGlobalBasis();
}
#endif
if (adoptProblem &&
(adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
......
#pragma once
#include <tuple>
#include <dune/functions/functionspacebases/basistags.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/flatmultiindex.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/pqknodalbasis.hh>
#include <dune/grid/yaspgrid.hh>
namespace Dune {
namespace Functions {
namespace BasisFactory {
// define the basis-build (factory-tag) for a given pre-basis
template <class PreBasis>
struct FactoryTag;
template <class MultiIndex, class IndexMergingStrategy, class... ChildPreBasisFactory>
struct FactoryTag<CompositePreBasis<MultiIndex, IndexMergingStrategy, ChildPreBasisFactory...>>
{
using type = Imp::CompositePreBasisFactory<IndexMergingStrategy, typename FactoryTag<ChildPreBasisFactory>::type...>;
};
template <class MultiIndex, class IndexMergingStrategy, class ChildPreBasisFactory, std::size_t C>
struct FactoryTag<PowerPreBasis<MultiIndex, IndexMergingStrategy, ChildPreBasisFactory, C>>
{
using type = Imp::PowerPreBasisFactory<C, IndexMergingStrategy, typename FactoryTag<ChildPreBasisFactory>::type>;
};
template <class GridView, int k, class MultiIndex>
struct FactoryTag<PQkPreBasis<GridView, k, MultiIndex>>
{
using type = Imp::PQkPreBasisFactory<k>;
};
}}} // end namespace Dune::Functions::BasisFactory
#include <amdis/common/Mpl.hpp>
namespace AMDiS
{
namespace Impl
{
template <class GridView, bool same, int... deg>
struct LagrangeBasisCreatorImpl;
template <bool same, int... degs>
struct LagrangePreBasisCreatorImpl;
// specialization for composite basis
template <class GridView, int... deg>
struct LagrangeBasisCreatorImpl<GridView, false, deg...>
template <int... degs>
struct LagrangePreBasisCreatorImpl<false, degs...>
{
using MultiIndex = Dune::ReservedVector<std::size_t,1>;
using IndexTag = Dune::Functions::BasisFactory::FlatLexicographic;
using type = Dune::Functions::CompositePreBasis<MultiIndex, IndexTag,
Dune::Functions::PQkPreBasis<GridView,deg,MultiIndex>...>;
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 LagrangeBasisCreatorImpl<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::BasisFactory::FlatLexicographic;
using type = Dune::Functions::PowerPreBasis<MultiIndex, IndexTag,
Dune::Functions::PQkPreBasis<GridView,deg,MultiIndex>, 1 + sizeof...(degs)>;
static auto create()
{
using namespace Dune::Functions::BasisFactory;
return power<1+sizeof...(degs)>(lagrange<deg>(), flatLexicographic());
}
};
// factory to construct a global basis of several lagrange bases, with flat indexing.
template <class GridView, int deg, int... degs>
struct LagrangeBasisCreator
{
using type = typename LagrangeBasisCreatorImpl<GridView, and_t<(deg == degs)...>::value, deg, degs...>::type;
};
template <int deg, int... degs>
struct LagrangePreBasisCreator
: public LagrangePreBasisCreatorImpl<all_of_v<(deg == degs)...>, deg, degs...>
{};
template <class GridView, int k>
struct TaylorHoodBasisCreator
template <int dow, int k = 1>
struct TaylorHoodPreBasisCreator
{
using MultiIndex = Dune::ReservedVector<std::size_t,1>;
using IndexTagV = Dune::Functions::BasisFactory::FlatInterleaved;