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

moved to dune-2.7 version

parent 3d0d7673
......@@ -36,11 +36,15 @@ int main(int argc, char** argv)
return {-20.0 * std::exp(-10.0 * dot(x,x)) * x};
};
ElliptProblem prob("ellipt");
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);
prob.initialize(INIT_ALL);
auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
prob.addMatrixOperator(opL, _0, _0);
prob.addMatrixOperator(opL);
auto opForce = makeOperator(tag::test{}, f, 6);
prob.addVectorOperator(opForce, _0);
......
......@@ -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);
HeatProblem prob("heat");
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);
prob.initialize(INIT_ALL);
HeatProblemInstat probInstat("heat", prob);
auto probInstat = makeProblemInstat("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, 0, 0);
prob.addMatrixOperator(opTimeLhs, _, _);
auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
prob.addMatrixOperator(opL, 0, 0);
prob.addMatrixOperator(opL, _, _);
auto opTimeRhs = makeOperator(tag::test{},
invokeAtQP([invTau](double u) { return u * (*invTau); }, prob.getSolution(0)), 2);
prob.addVectorOperator(opTimeRhs, 0);
invokeAtQP([invTau](double u) { return u * (*invTau); }, prob.getSolution()), 2);
prob.addVectorOperator(opTimeRhs, _);
auto opForce = makeOperator(tag::test{}, [](auto const& x) { return -1.0; }, 0);
prob.addVectorOperator(opForce, 0);
prob.addVectorOperator(opForce, _);
// 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, 0, 0, dbcValues);
prob.addDirichletBC(predicate, _, _, dbcValues);
AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
adapt.adapt();
......
......@@ -12,23 +12,20 @@
using namespace AMDiS;
struct NavierStokesBasis
{
using Grid = Dune::YaspGrid<AMDIS_DIM>;
using GlobalBasis = typename TaylorHoodBasis<Grid::LeafGridView>::GlobalBasis;
};
using StokesProblem = ProblemStat<NavierStokesBasis>;
using StokesProblemInstat = ProblemInstat<NavierStokesBasis>;
int main(int argc, char** argv)
{
AMDiS::init(argc, argv);
StokesProblem prob("stokes");
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);
prob.initialize(INIT_ALL);
StokesProblemInstat probInstat("stokes", prob);
auto probInstat = makeProblemInstat("stokes", prob);
probInstat.initialize(INIT_UH_OLD);
double viscosity = 1.0;
......
......@@ -13,17 +13,16 @@
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);
StokesProblem prob("stokes");
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);
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);
StokesProblem prob("stokes");
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);
prob.initialize(INIT_ALL);
double viscosity = 1.0;
......
......@@ -9,15 +9,17 @@
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);
StokesProblem prob("stokes");
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);
prob.initialize(INIT_ALL);
double viscosity = 1.0;
......
......@@ -6,16 +6,17 @@
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);
ElliptProblem prob("ellipt");
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);
prob.initialize(INIT_ALL);
AdaptInfo adaptInfo("adapt");
......
......@@ -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_;
......
#pragma once
#include <dune/functions/functionspacebases/subspacebasis.hh>
#include <dune/typetree/traversal.hh>
#include <amdis/utility/TreePath.hpp>
#include <amdis/utility/Visitor.hpp>
#include <amdis/common/Math.hpp>
......@@ -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();
......@@ -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,7 +89,7 @@ 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];
}
}
......@@ -105,7 +102,6 @@ void Assembler<Traits>::assemble(
});
});
localIndexSet.unbind();
localView.unbind();
}
......
......@@ -67,13 +67,14 @@ 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.
**/
explicit ProblemStat(std::string name)
: StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
, name_(std::move(name))
, name(std::move(name))
{}
/// Constructor taking additionally a reference to a grid that is used
......@@ -93,6 +94,20 @@ 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.
ProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
: StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
, name(std::move(name))
{
gridName = "";
Parameters::get(name + "->mesh", gridName);
globalBasis_ = Dune::stackobject_to_shared_ptr(globalBasis);
initGlobalBasis(*globalBasis_);
}
/**
......@@ -108,27 +123,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 = {});
/** @} */
......@@ -200,7 +211,8 @@ namespace AMDiS
/// Return a pointer to the grid, \ref grid
std::shared_ptr<Grid> getGrid() { return grid_; }
/// Set the grid. Stores pointer to passed reference and initializes feSpaces
#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)
{
......@@ -210,7 +222,10 @@ 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_; }
......@@ -239,18 +254,12 @@ 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());
matrixOperators_.init(localView_->tree(), tag::store{});
rhsOperators_.init(localView_->tree(), tag::store{});
constraints_.init(localView_->tree(), tag::store{});
matrixOperators.init(localView_->tree(), tag::store{});
rhsOperators.init(localView_->tree(), tag::store{});
constraints.init(localView_->tree(), tag::store{});
}
void createMatricesAndVectors()
......
......@@ -57,10 +57,12 @@ 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))) {
......@@ -146,7 +148,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 +167,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 +180,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 +190,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 +204,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 +213,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 +223,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 +233,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 +244,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;
}
......
......@@ -3,173 +3,148 @@
#include <tuple>
#include <dune/functions/functionspacebases/basistags.hh>
#include <dune/functions/functionspacebases/flatmultiindex.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/functions/functionspacebases/pq1nodalbasis.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
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;
struct LagrangeBasisCreatorImpl;
// specialization for composite basis
template <class GridView, int... deg>
struct LagrangeBasisBuilderImpl<GridView, false, deg...>
struct LagrangeBasisCreatorImpl<GridView, false, deg...>
{
using MultiIndex = Dune::ReservedVector<std::size_t,1>;
using IndexTag = Dune::Functions::BasisBuilder::FlatLexicographic;
using IndexTag = Dune::Functions::BasisFactory::FlatLexicographic;
using type = Dune::Functions::CompositeNodeFactory<MultiIndex, IndexTag,
typename LagrangeBasisSubFactory<GridView,deg,MultiIndex>::type...>;
using type = Dune::Functions::CompositePreBasis<MultiIndex, IndexTag,
Dune::Functions::PQkPreBasis<GridView,deg,MultiIndex>...>;
};
// specialization for power basis
template <class GridView, int deg, int... degs>
struct LagrangeBasisBuilderImpl<GridView, true, deg, degs...>
struct LagrangeBasisCreatorImpl<GridView, true, deg, degs...>
{
using MultiIndex = Dune::ReservedVector<std::size_t,1>;
using IndexTag = Dune::Functions::BasisBuilder::FlatLexicographic;
using IndexTag = Dune::Functions::BasisFactory::FlatLexicographic;
using type = Dune::Functions::PowerNodeFactory<MultiIndex, IndexTag,
typename LagrangeBasisSubFactory<GridView,deg,MultiIndex>::type, 1 + sizeof...(degs)>;
using type = Dune::Functions::PowerPreBasis<MultiIndex, IndexTag,
Dune::Functions::PQkPreBasis<GridView,deg,MultiIndex>, 1 + sizeof...(degs)>;
};
// factory to construct a global basis of several lagrange bases, with flat indexing.
template <class GridView, int deg, int... degs>
struct LagrangeBasisBuilder
struct LagrangeBasisCreator
{
using type = typename LagrangeBasisBuilderImpl<GridView, and_t<(deg == degs)...>::value, deg, degs...>::type;
using type = typename LagrangeBasisCreatorImpl<GridView, and_t<(deg == degs)...>::value, deg, degs...>::type;
};
template <class GridView>
struct TaylorHoodBasisBuilder
template <class GridView, int k>
struct TaylorHoodBasisCreator
{
using MultiIndex = Dune::ReservedVector<std::size_t,1>;
using IndexTagV = Dune::Functions::BasisBuilder::FlatInterleaved;
using IndexTagV = Dune::Functions::BasisFactory::FlatInterleaved;