From c786927261d1c07223ce3f9feabd89e439f2f4bb Mon Sep 17 00:00:00 2001 From: Simon Praetorius <simon.praetorius@tu-dresden.de> Date: Tue, 1 May 2018 19:33:24 +0200 Subject: [PATCH] moved to dune-2.7 version --- examples/ellipt.cc | 8 +- examples/heat.cc | 28 +-- examples/navier_stokes.cc | 19 +- examples/stokes0.cc | 13 +- examples/stokes1.cc | 14 +- examples/stokes3.cc | 12 +- examples/vecellipt.cc | 13 +- src/amdis/Assembler.hpp | 15 -- src/amdis/Assembler.inc.hpp | 12 +- src/amdis/ProblemStat.hpp | 55 +++-- src/amdis/ProblemStat.inc.hpp | 24 +-- src/amdis/ProblemStatTraits.hpp | 203 ++++++++---------- src/amdis/common/Concepts.hpp | 2 +- src/amdis/common/ConceptsBase.hpp | 15 +- src/amdis/common/Loops.hpp | 2 + src/amdis/common/Mpl.hpp | 3 - .../gridfunctions/AnalyticGridFunction.hpp | 2 +- src/amdis/gridfunctions/DOFVectorView.hpp | 43 ++-- src/amdis/gridfunctions/DOFVectorView.inc.hpp | 55 ++--- .../gridfunctions/GridFunctionConcepts.hpp | 2 +- .../linear_algebra/HierarchicWrapper.hpp | 4 +- src/amdis/linear_algebra/istl/DOFMatrix.hpp | 14 +- src/amdis/test/kdtree.hpp | 5 +- src/amdis/test/test1.cc | 5 +- src/amdis/utility/CMakeLists.txt | 2 - src/amdis/utility/Traversal.hpp | 137 ------------ src/amdis/utility/TreeData.hpp | 2 +- test/ExpressionsTest.cpp | 13 +- 28 files changed, 250 insertions(+), 472 deletions(-) delete mode 100644 src/amdis/utility/Traversal.hpp diff --git a/examples/ellipt.cc b/examples/ellipt.cc index 0414e633..5b620582 100644 --- a/examples/ellipt.cc +++ b/examples/ellipt.cc @@ -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); diff --git a/examples/heat.cc b/examples/heat.cc index d2322641..f09f7505 100644 --- a/examples/heat.cc +++ b/examples/heat.cc @@ -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(); diff --git a/examples/navier_stokes.cc b/examples/navier_stokes.cc index 4fcf941b..4ac468a0 100644 --- a/examples/navier_stokes.cc +++ b/examples/navier_stokes.cc @@ -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; diff --git a/examples/stokes0.cc b/examples/stokes0.cc index fbcc5fb8..37588d94 100644 --- a/examples/stokes0.cc +++ b/examples/stokes0.cc @@ -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; diff --git a/examples/stokes1.cc b/examples/stokes1.cc index cd62b4cd..cee81273 100644 --- a/examples/stokes1.cc +++ b/examples/stokes1.cc @@ -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; diff --git a/examples/stokes3.cc b/examples/stokes3.cc index 7c4e3a6c..f7b29e87 100644 --- a/examples/stokes3.cc +++ b/examples/stokes3.cc @@ -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; diff --git a/examples/vecellipt.cc b/examples/vecellipt.cc index 8ea86360..78aaa0e3 100644 --- a/examples/vecellipt.cc +++ b/examples/vecellipt.cc @@ -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"); diff --git a/src/amdis/Assembler.hpp b/src/amdis/Assembler.hpp index 8473ac32..ef99fb52 100644 --- a/src/amdis/Assembler.hpp +++ b/src/amdis/Assembler.hpp @@ -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_; diff --git a/src/amdis/Assembler.inc.hpp b/src/amdis/Assembler.inc.hpp index 05501f75..7afe0dfb 100644 --- a/src/amdis/Assembler.inc.hpp +++ b/src/amdis/Assembler.inc.hpp @@ -1,8 +1,8 @@ #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(); } diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp index 47e5a636..fc5f6048 100644 --- a/src/amdis/ProblemStat.hpp +++ b/src/amdis/ProblemStat.hpp @@ -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() diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp index 4f08455f..1db46d2b 100644 --- a/src/amdis/ProblemStat.inc.hpp +++ b/src/amdis/ProblemStat.inc.hpp @@ -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; } diff --git a/src/amdis/ProblemStatTraits.hpp b/src/amdis/ProblemStatTraits.hpp index bd0c562e..6dcb9e17 100644 --- a/src/amdis/ProblemStatTraits.hpp +++ b/src/amdis/ProblemStatTraits.hpp @@ -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; - using VelocityNodeFactory = Dune::Functions::PowerNodeFactory<MultiIndex, IndexTagV, - typename LagrangeBasisSubFactory<GridView,2,MultiIndex>::type, 2>; + using VelocityPreBasis = Dune::Functions::PowerPreBasis<MultiIndex, IndexTagV, + Dune::Functions::PQkPreBasis<GridView,k+1,MultiIndex>, GridView::dimensionworld>; - using PressureNodeFactory = typename LagrangeBasisSubFactory<GridView,1,MultiIndex>::type; + using PressurePreBasis = Dune::Functions::PQkPreBasis<GridView,k,MultiIndex>; - using IndexTag = Dune::Functions::BasisBuilder::FlatLexicographic; - using type = Dune::Functions::CompositeNodeFactory<MultiIndex, IndexTag, - VelocityNodeFactory, PressureNodeFactory>; + using IndexTag = Dune::Functions::BasisFactory::FlatLexicographic; + using type = Dune::Functions::CompositePreBasis<MultiIndex, IndexTag, + VelocityPreBasis, PressurePreBasis>; }; } // end namespace Impl - - /// \brief A factory for a composite basis composed of lagrange bases of different degree. - template <class GridView, int... degrees> - struct LagrangeBasis + /// An Exemplar for ProblemStatTraits + template <class GlobalBasisType> + class DefaultProblemTraits { - using GlobalBasis - = Dune::Functions::DefaultGlobalBasis<typename Impl::LagrangeBasisBuilder<GridView, degrees...>::type>; + public: + /// A global-basis of type dune-functions DefaultGlobalBasis<NodeFactory> + using GlobalBasis = GlobalBasisType; }; + /// \brief ProblemStatTraits for a (composite) basis composed of + /// lagrange bases of different degree. + template <class GridView, int... degrees> + struct LagrangeBasis + : public DefaultProblemTraits< + Dune::Functions::DefaultGlobalBasis< + typename Impl::LagrangeBasisCreator<GridView, degrees...>::type> > + {}; - /// Specialization of \ref LagrangeBasis for Grid type \ref Dune::YaspGrid for a given dimension. + /// \brief Specialization of \ref LagrangeBasis for Grid type + /// \ref Dune::YaspGrid for a given dimension. template <int dim, int... degrees> struct YaspGridBasis - { - using GlobalBasis = typename LagrangeBasis<typename Dune::YaspGrid<dim>::LeafGridView, degrees...>::GlobalBasis; - }; + : public LagrangeBasis<typename Dune::YaspGrid<dim>::LeafGridView, degrees...> + {}; - - template <class GridView> + /// \brief ProblemStatTraits of Taylor-Hood basis of lagrange-type + /// with pressure degree k + template <class GridView, int k = 1> struct TaylorHoodBasis - { - using GlobalBasis - = Dune::Functions::DefaultGlobalBasis<typename Impl::TaylorHoodBasisBuilder<GridView>::type>; - }; - - template <class GlobalBasisType> - struct DefaultProblemTraits - { - using GlobalBasis = GlobalBasisType; - }; - -} // end namespace AMDiS - - -namespace Dune { -namespace Functions { -namespace BasisBuilder { - - namespace Imp - { - struct PQ1NodeFactoryBuilder - { - static const std::size_t requiredMultiIndexSize=1; - - template<class MultiIndex, class GridView> - auto build(const GridView& gridView) - -> PQ1NodeFactory<GridView, MultiIndex> - { - return {gridView}; - } - }; - - } // end namespace Imp + : public DefaultProblemTraits< + Dune::Functions::DefaultGlobalBasis< + typename Impl::TaylorHoodBasisCreator<GridView,k>::type> > + {}; - // define the basis-build (factory-tag) for a given node-factory - template <class NodeFactory> - struct FactoryTag; + using Dune::Functions::BasisFactory::makeBasis; + using Dune::Functions::BasisFactory::lagrange; + using Dune::Functions::BasisFactory::pq; - template <class MultiIndex, class IndexMergingStrategy, class... SubFactories> - struct FactoryTag<CompositeNodeFactory<MultiIndex, IndexMergingStrategy, SubFactories...>> + template <std::size_t k, class ChildPBF> + auto repeat(ChildPBF&& childPreBasisFactory) { - using type = Imp::CompositeNodeFactoryBuilder<IndexMergingStrategy, typename FactoryTag<SubFactories>::type...>; - }; - - template <class MultiIndex, class IndexMergingStrategy, class SubFactory, std::size_t C> - struct FactoryTag<PowerNodeFactory<MultiIndex, IndexMergingStrategy, SubFactory, C>> - { - using type = Imp::PowerNodeFactoryBuilder<C, IndexMergingStrategy, typename FactoryTag<SubFactory>::type>; - }; - - template <class GridView, int k, class MultiIndex> - struct FactoryTag<PQkNodeFactory<GridView, k, MultiIndex>> - { - using type = Imp::PQkNodeFactoryBuilder<k>; - }; - - template <class GridView, class MultiIndex> - struct FactoryTag<PQ1NodeFactory<GridView, MultiIndex>> - { - using type = Imp::PQ1NodeFactoryBuilder; - }; - -}}} // end namespace Dune::Functions::BasisBuilder - + return Dune::Functions::BasisFactory::power<k>( + std::forward<ChildPBF>(childPreBasisFactory), + Dune::Functions::BasisFactory::flatInterleaved()); + } -namespace AMDiS -{ - template <class GlobalBasis, class GridView> - GlobalBasis makeGlobalBasis(GridView const& gv) + template <class... Args> + auto compose(Args&&... args) { - using NodeFactory = typename GlobalBasis::NodeFactory; - using Builder = typename Dune::Functions::BasisBuilder::FactoryTag<NodeFactory>::type; - return makeBasis(gv, Builder{}); + return Dune::Functions::BasisFactory::composite( + std::forward<Args>(args)..., + Dune::Functions::BasisFactory::flatLexicographic()); } } // end namespace AMDiS diff --git a/src/amdis/common/Concepts.hpp b/src/amdis/common/Concepts.hpp index a04721a0..c548cc02 100644 --- a/src/amdis/common/Concepts.hpp +++ b/src/amdis/common/Concepts.hpp @@ -102,7 +102,7 @@ namespace AMDiS struct Evaluable { template <class F, class... Args> - auto requires_(F&& f, Args&&... args) -> Void_t< + auto requires_(F&& f, Args&&... args) -> void_t< decltype( f(std::forward<Args>(args)...)), std::result_of_t<std::decay_t<F>(std::decay_t<Args>...)> >; diff --git a/src/amdis/common/ConceptsBase.hpp b/src/amdis/common/ConceptsBase.hpp index 7df39159..a42cd50c 100644 --- a/src/amdis/common/ConceptsBase.hpp +++ b/src/amdis/common/ConceptsBase.hpp @@ -14,28 +14,19 @@ namespace AMDiS { - namespace Impl - { - // Workaround for MSVC (problems with alias templates in pack expansion) - template <class... Args> - struct VoidType { using type = void; }; - } - - template <class... Args> - using Void_t = typename Impl::VoidType<Args...>::type; - + using Dune::void_t; namespace Concepts { namespace Impl_ { - template <class Concept, class = Void_t<>> + template <class Concept, class = void_t<>> struct models : std::false_type {}; template <class Concept, class... Ts> - struct models<Concept(Ts...), Void_t< decltype(std::declval<Concept>().requires_(std::declval<Ts>()...)) >> + struct models<Concept(Ts...), void_t< decltype(std::declval<Concept>().requires_(std::declval<Ts>()...)) >> : std::true_type {}; diff --git a/src/amdis/common/Loops.hpp b/src/amdis/common/Loops.hpp index f58ac532..ca004320 100644 --- a/src/amdis/common/Loops.hpp +++ b/src/amdis/common/Loops.hpp @@ -1,8 +1,10 @@ #pragma once #include <dune/common/hybridutilities.hh> +#include <dune/common/rangeutilities.hh> namespace AMDiS { using Dune::Hybrid::forEach; + } // end namespace AMDiS diff --git a/src/amdis/common/Mpl.hpp b/src/amdis/common/Mpl.hpp index 716c95df..2e49f5a5 100644 --- a/src/amdis/common/Mpl.hpp +++ b/src/amdis/common/Mpl.hpp @@ -51,9 +51,6 @@ namespace AMDiS template <bool B> constexpr bool_t<B> bool_ = {}; - static constexpr bool_t<true> _true {}; - static constexpr bool_t<false> _false {}; - namespace Impl { diff --git a/src/amdis/gridfunctions/AnalyticGridFunction.hpp b/src/amdis/gridfunctions/AnalyticGridFunction.hpp index c40e5af1..00f63c5e 100644 --- a/src/amdis/gridfunctions/AnalyticGridFunction.hpp +++ b/src/amdis/gridfunctions/AnalyticGridFunction.hpp @@ -154,7 +154,7 @@ namespace AMDiS auto derivative(AnalyticGridFunction<F,GV> const& gf) { static_assert(Concepts::HasPartial<F>, - "No partial(_0,...) defined for Functor of AnalyticLocalFunction."); + "No partial(...,_0) defined for Functor of AnalyticLocalFunction."); auto df = partial(gf.fct(), index_<0>); return AnalyticGridFunction<decltype(df), GV>{df, gf.entitySet().gridView()}; diff --git a/src/amdis/gridfunctions/DOFVectorView.hpp b/src/amdis/gridfunctions/DOFVectorView.hpp index e71139e6..0cfccfc2 100644 --- a/src/amdis/gridfunctions/DOFVectorView.hpp +++ b/src/amdis/gridfunctions/DOFVectorView.hpp @@ -6,7 +6,7 @@ #include <dune/functions/common/defaultderivativetraits.hh> #include <dune/functions/common/treedata.hh> #include <dune/functions/functionspacebases/defaultnodetorangemap.hh> -#include <dune/functions/functionspacebases/flatvectorbackend.hh> +#include <dune/functions/functionspacebases/flatvectorview.hh> #include <dune/functions/gridfunctions/gridviewentityset.hh> #include <dune/typetree/childextraction.hh> @@ -46,9 +46,6 @@ namespace AMDiS using Element = typename EntitySet::Element; using Geometry = typename Element::Geometry; - template <class Block> - using Flat = Dune::Functions::FlatVectorBackend<Block>; - enum { hasDerivative = false }; public: // a local view on the gradients @@ -63,8 +60,7 @@ namespace AMDiS enum { hasDerivative = false }; private: - using LocalBasisView = typename GlobalBasis::LocalView; - using LocalIndexSet = typename GlobalBasis::LocalIndexSet; + using LocalView = typename GlobalBasis::LocalView; template <class LeafNode> using LocalBasisJacobian = typename LeafNode::FiniteElement::Traits::LocalBasisType::Traits::JacobianType; @@ -77,25 +73,22 @@ namespace AMDiS public: GradientLocalFunction(DOFVectorConstView const& globalFunction) : globalFunction_(&globalFunction) - , localBasisView_(globalFunction_->basis().localView()) - , localIndexSet_(globalFunction_->basis().localIndexSet()) - , subTree_(&Dune::TypeTree::child(localBasisView_.tree(), globalFunction_->treePath())) + , localView_(globalFunction_->basis().localView()) + , subTree_(&Dune::TypeTree::child(localView_.tree(), globalFunction_->treePath())) { referenceGradientContainer_.init(*subTree_); } void bind(Element const& element) { - localBasisView_.bind(element); - localIndexSet_.bind(localBasisView_); + localView_.bind(element); geometry_.emplace(element.geometry()); bound_ = true; } void unbind() { - localIndexSet_.unbind(); - localBasisView_.unbind(); + localView_.unbind(); geometry_.reset(); bound_ = false; } @@ -113,13 +106,12 @@ namespace AMDiS Element const& localContext() const { assert( bound_ ); - return localBasisView_.element(); + return localView_.element(); } private: DOFVectorConstView const* globalFunction_; - LocalBasisView localBasisView_; - LocalIndexSet localIndexSet_; + LocalView localView_; SubTree const* subTree_; mutable ReferenceGradientContainer referenceGradientContainer_; @@ -141,8 +133,7 @@ namespace AMDiS enum { hasDerivative = true }; private: - using LocalBasisView = typename GlobalBasis::LocalView; - using LocalIndexSet = typename GlobalBasis::LocalIndexSet; + using LocalView = typename GlobalBasis::LocalView; template <class LeafNode> using LocalBasisRange = RangeType_t<LeafNode>; @@ -156,24 +147,21 @@ namespace AMDiS public: LocalFunction(DOFVectorConstView const& globalFunction) : globalFunction_(&globalFunction) - , localBasisView_(globalFunction_->basis().localView()) - , localIndexSet_(globalFunction_->basis().localIndexSet()) - , subTree_(&Dune::TypeTree::child(localBasisView_.tree(), globalFunction_->treePath())) + , localView_(globalFunction_->basis().localView()) + , subTree_(&Dune::TypeTree::child(localView_.tree(), globalFunction_->treePath())) { shapeFunctionValueContainer_.init(*subTree_); } void bind(Element const& element) { - localBasisView_.bind(element); - localIndexSet_.bind(localBasisView_); + localView_.bind(element); bound_ = true; } void unbind() { - localIndexSet_.unbind(); - localBasisView_.unbind(); + localView_.unbind(); bound_ = false; } @@ -197,13 +185,12 @@ namespace AMDiS Element const& localContext() const { assert( bound_ ); - return localBasisView_.element(); + return localView_.element(); } private: DOFVectorConstView const* globalFunction_; - LocalBasisView localBasisView_; - LocalIndexSet localIndexSet_; + LocalView localView_; SubTree const* subTree_; mutable ShapeFunctionValueContainer shapeFunctionValueContainer_; diff --git a/src/amdis/gridfunctions/DOFVectorView.inc.hpp b/src/amdis/gridfunctions/DOFVectorView.inc.hpp index 123de812..5d122263 100644 --- a/src/amdis/gridfunctions/DOFVectorView.inc.hpp +++ b/src/amdis/gridfunctions/DOFVectorView.inc.hpp @@ -15,9 +15,6 @@ LocalFunction::operator()(LocalDomain const& x) const auto&& nodeToRangeEntry = globalFunction_->nodeToRangeEntry_; AMDiS::forEachLeafNode_(*subTree_, [&,this](auto const& node, auto) { - using Node = std::decay_t<decltype(node)>; - using LocalBasisRange = typename LocalFunction::template LocalBasisRange<Node>; - auto&& fe = node.finiteElement(); auto&& localBasis = fe.localBasis(); @@ -25,30 +22,24 @@ LocalFunction::operator()(LocalDomain const& x) const localBasis.evaluateFunction(x, shapeFunctionValues); // Get range entry associated to this node - auto&& re = nodeToRangeEntry(node, y); - using RangeBlock = std::decay_t<decltype(re)>; + auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, y)); for (std::size_t i = 0; i < localBasis.size(); ++i) { - auto&& multiIndex = localIndexSet_.index(node.localIndex(i)); + auto&& multiIndex = localView_.index(node.localIndex(i)); // Get coefficient associated to i-th shape function - auto&& c = coefficients[multiIndex]; - using CoefficientBlock = std::decay_t<decltype(c)>; + auto c = Dune::Functions::flatVectorView(coefficients[multiIndex]); // Get value of i-th shape function - auto&& v = shapeFunctionValues[i]; + auto v = Dune::Functions::flatVectorView(shapeFunctionValues[i]); - // Notice that the range entry re, the coefficient c, and the shape functions - // value v may all be scalar, vector, matrix, or general container valued. - std::size_t dimC = Flat<CoefficientBlock>::size(c); - std::size_t dimV = Flat<LocalBasisRange>::size(v); - assert(dimC*dimV == Flat<RangeBlock>::size(re)); + std::size_t dimC = c.size(); + std::size_t dimV = v.size(); + assert(dimC*dimV == std::size_t(re.size())); for(std::size_t j = 0; j < dimC; ++j) { - auto&& c_j = Flat<CoefficientBlock>::getEntry(c, j); - for(std::size_t k = 0; k < dimV; ++k) { - auto&& v_k = Flat<LocalBasisRange>::getEntry(v, k); - Flat<RangeBlock>::getEntry(re, j*dimV + k) += c_j*v_k; - } + auto&& c_j = c[j]; + for(std::size_t k = 0; k < dimV; ++k) + re[j*dimV + k] += c_j*v[k]; } } }); @@ -89,30 +80,24 @@ GradientLocalFunction::operator()(LocalDomain const& x) const multiplies_ABt(referenceGradients[i], jacobian, gradients[i]); // D[phi] * J^(-1) -> grad // Get range entry associated to this node - auto&& re = nodeToRangeEntry(node, dy); - using RangeBlock = std::decay_t<decltype(re)>; + auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, dy)); for (std::size_t i = 0; i < localBasis.size(); ++i) { - auto&& multiIndex = localIndexSet_.index(node.localIndex(i)); + auto&& multiIndex = localView_.index(node.localIndex(i)); // Get coefficient associated to i-th shape function - auto&& c = coefficients[multiIndex]; - using CoefficientBlock = std::decay_t<decltype(c)>; + auto c = Dune::Functions::flatVectorView(coefficients[multiIndex]); // Get value of i-th transformed reference gradient - auto&& grad = gradients[i]; + auto grad = Dune::Functions::flatVectorView(gradients[i]); - // Notice that the range entry re, the coefficient c, and the transformed - // reference gradient grad may all be scalar, vector, matrix, or general container valued. - std::size_t dimC = Flat<CoefficientBlock>::size(c); - std::size_t dimV = Flat<GradientBlock>::size(grad); - assert(dimC*dimV == std::size_t(Flat<RangeBlock>::size(re))); + std::size_t dimC = c.size(); + std::size_t dimV = grad.size(); + assert(dimC*dimV == std::size_t(re.size())); for(std::size_t j = 0; j < dimC; ++j) { - auto&& c_j = Flat<CoefficientBlock>::getEntry(c, j); - for(std::size_t k = 0; k < dimV; ++k) { - auto&& v_k = Flat<GradientBlock>::getEntry(grad, k); - Flat<RangeBlock>::getEntry(re, j*dimV + k) += c_j*v_k; - } + auto&& c_j = c[j]; + for(std::size_t k = 0; k < dimV; ++k) + re[j*dimV + k] += c_j*grad[k]; } } }); diff --git a/src/amdis/gridfunctions/GridFunctionConcepts.hpp b/src/amdis/gridfunctions/GridFunctionConcepts.hpp index 50d02d6a..80e277df 100644 --- a/src/amdis/gridfunctions/GridFunctionConcepts.hpp +++ b/src/amdis/gridfunctions/GridFunctionConcepts.hpp @@ -62,7 +62,7 @@ namespace AMDiS struct HasGridFunctionTypes { template <class GF> - auto requires_(GF const& /*gf*/) -> Void_t< + auto requires_(GF const& /*gf*/) -> void_t< typename GF::Range, typename GF::Domain, typename GF::EntitySet >; diff --git a/src/amdis/linear_algebra/HierarchicWrapper.hpp b/src/amdis/linear_algebra/HierarchicWrapper.hpp index 858231db..9da18972 100644 --- a/src/amdis/linear_algebra/HierarchicWrapper.hpp +++ b/src/amdis/linear_algebra/HierarchicWrapper.hpp @@ -67,14 +67,14 @@ namespace AMDiS protected: template <class Vec, class Basis> - Void_t<decltype( std::declval<Vec>().resize(std::size_t(0)) )> + void_t<decltype( std::declval<Vec>().resize(std::size_t(0)) )> resizeImpl(Vec& vec, Dune::Functions::SizeInfo<Basis> const& s) { vec.resize(std::size_t(s)); } template <class Vec, class Basis> - Void_t<decltype( std::declval<Vec>().change_dim(std::size_t(0)) )> + void_t<decltype( std::declval<Vec>().change_dim(std::size_t(0)) )> resizeImpl(Vec& vec, Dune::Functions::SizeInfo<Basis> const& s) { vec.change_dim(std::size_t(s)); diff --git a/src/amdis/linear_algebra/istl/DOFMatrix.hpp b/src/amdis/linear_algebra/istl/DOFMatrix.hpp index 808af385..e256a05e 100644 --- a/src/amdis/linear_algebra/istl/DOFMatrix.hpp +++ b/src/amdis/linear_algebra/istl/DOFMatrix.hpp @@ -110,24 +110,18 @@ namespace AMDiS // A loop over all elements of the grid auto rowLocalView = rowFeSpace.localView(); - auto rowIndexSet = rowFeSpace.localIndexSet(); - auto colLocalView = colFeSpace.localView(); - auto colIndexSet = colFeSpace.localIndexSet(); for (const auto& element : elements(meshView)) { rowLocalView.bind(element); - rowIndexSet.bind(rowLocalView); - colLocalView.bind(element); - colIndexSet.bind(colLocalView); - for (std::size_t i = 0; i < rowIndexSet.size(); ++i) { + for (std::size_t i = 0; i < rowLocalView.size(); ++i) { // The global index of the i-th vertex of the element - auto row = rowIndexSet.index(i); - for (std::size_t j = 0; j < colIndexSet.size(); ++j) { + auto row = rowLocalView.index(i); + for (std::size_t j = 0; j < colLocalView.size(); ++j) { // The global index of the j-th vertex of the element - auto col = colIndexSet.index(j); + auto col = colLocalView.index(j); occupationPattern.add(row, col); } } diff --git a/src/amdis/test/kdtree.hpp b/src/amdis/test/kdtree.hpp index d74e4669..007d74ca 100644 --- a/src/amdis/test/kdtree.hpp +++ b/src/amdis/test/kdtree.hpp @@ -46,7 +46,6 @@ namespace AMDiS const int leaf_max_size = 10) : feSpace(feSpace_) , localView(feSpace.localView()) - , localIndexSet(feSpace.localIndexSet()) { init(); @@ -168,7 +167,6 @@ namespace AMDiS LocalCoords const& lambda) { localView.bind(element); - localIndexSet.bind(localView); auto const& localFiniteElem = localView.tree().finiteElement(); auto const& localBasis = localFiniteElem.localBasis(); @@ -179,7 +177,7 @@ namespace AMDiS Value data = 0; for (std::size_t j = 0; j < nBasisFct; ++j) { - auto const global_idx = localIndexSet.index(j); + auto const global_idx = localView.index(j); data += vector[global_idx] * shapeValues[j]; } @@ -203,7 +201,6 @@ namespace AMDiS private: typename FeSpace::LocalView localView; - typename FeSpace::LocalIndexSet localIndexSet; }; // end of KDTreeMeshAdaptor diff --git a/src/amdis/test/test1.cc b/src/amdis/test/test1.cc index e7ecb748..a3a154c5 100644 --- a/src/amdis/test/test1.cc +++ b/src/amdis/test/test1.cc @@ -29,14 +29,12 @@ namespace { public: LocalEvaluation(FeSpace const& feSpace) : localView(feSpace.localView()) - , localIndexSet(feSpace.localIndexSet()) {} template <class Element> void bind(Element const& element) { localView.bind(element); - localIndexSet.bind(localView); } /// evaluate DOFVector at barycentric coordinates \p lambda in element that @@ -49,7 +47,7 @@ namespace { ValueType data = 0; for (std::size_t j = 0; j < shapeValues.size(); ++j) { - const auto global_idx = localIndexSet.index(j); + const auto global_idx = localView.index(j); data += vec[global_idx] * shapeValues[j]; } @@ -58,7 +56,6 @@ namespace { private: typename FeSpace::LocalView localView; - typename FeSpace::LocalIndexSet localIndexSet; std::vector<Dune::FieldVector<double,1> > shapeValues; }; diff --git a/src/amdis/utility/CMakeLists.txt b/src/amdis/utility/CMakeLists.txt index bef683e0..69b9b57d 100644 --- a/src/amdis/utility/CMakeLists.txt +++ b/src/amdis/utility/CMakeLists.txt @@ -10,8 +10,6 @@ install(FILES MultiIndex.hpp RangeType.hpp String.hpp - Traversal.hpp TreeData.hpp TreePath.hpp - Visitor.hpp DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/utility) diff --git a/src/amdis/utility/Traversal.hpp b/src/amdis/utility/Traversal.hpp deleted file mode 100644 index 3a02e742..00000000 --- a/src/amdis/utility/Traversal.hpp +++ /dev/null @@ -1,137 +0,0 @@ -#pragma once - -#include <dune/typetree/nodetags.hh> -#include <dune/typetree/treepath.hh> -#include <dune/typetree/visitor.hh> - -#include <amdis/common/Loops.hpp> - -namespace AMDiS -{ - // forward declaration of main engine struct - template <typename NodeTag, bool visit = true> - struct TraverseTree; - - - // Do not visit nodes the visitor is not interested in - template <typename NodeTag> - struct TraverseTree<NodeTag, false> - { - template <typename Node, typename Visitor, typename TreePath> - static void apply(const Node& node, const Visitor& visitor, TreePath const& tp) - {} - }; - -#ifndef DOXYGEN - - // some implementation details - - template <class Node, class Index> - struct HybridChildType - : HybridChildType<std::remove_const_t<Node>, std::remove_const_t<Index>> {}; - - template <class Node> - struct HybridChildType<Node, std::size_t> - { - using type = typename Node::template Child<0>::Type; - }; - - template <class Node, std::size_t K> - struct HybridChildType<Node, Dune::index_constant<K>> - { - using type = typename Node::template Child<K>::Type; - }; - - template <class NodeTag, class Node> - constexpr std::size_t hybridDegree(NodeTag, Node const& node) - { - return Dune::TypeTree::degree(node); - } - - template <class Node> - constexpr auto hybridDegree(Dune::TypeTree::CompositeNodeTag, Node const& node) - { - return Dune::index_constant<Node::CHILDREN>{}; - } - - - template <std::size_t k, std::size_t n> - constexpr bool notLastChild(Dune::index_constant<k> const&, Dune::index_constant<n> const&) - { - return k < n-1; - } - - constexpr bool notLastChild(std::size_t k, std::size_t n) - { - return k < n-1; - } - -#endif - - - template <class NodeTag> - struct TraverseTree<NodeTag, true> - { - template <typename N, typename V, typename TreePath> - static void apply(N&& n, V&& v, TreePath const& tp) - { - using Node = std::remove_reference_t<N>; - using Visitor = std::remove_reference_t<V>; - - v.pre(std::forward<N>(n),tp); - - auto const deg = hybridDegree(NodeTag{}, n); - forEach(Dune::range(deg), [&](auto const _k) - { - // always call beforeChild(), regardless of the value of visit - v.beforeChild(std::forward<N>(n),n.child(_k),tp,_k); - - // descend to child - using C = typename HybridChildType<Node, decltype(_k)>::type; - const bool visit = Visitor::template VisitChild<Node,C,TreePath>::value; - TraverseTree<Dune::TypeTree::NodeTag<C>,visit>::apply(n.child(_k),std::forward<V>(v),push_back(tp, _k)); - - // always call afterChild(), regardless of the value of visit - v.afterChild(std::forward<N>(n),n.child(_k),tp,_k); - - // if this is not the last child, call infix callback - if (notLastChild(_k, deg)) - v.in(std::forward<N>(n),tp); - }); - - v.post(std::forward<N>(n),tp); - } - }; - - // LeafNode - just call the leaf() callback - template <> - struct TraverseTree<Dune::TypeTree::LeafNodeTag, true> - { - template <typename N, typename V, typename TreePath> - static void apply(N&& n, V&& v, TreePath const& tp) - { - v.leaf(std::forward<N>(n),tp); - } - }; - - //! Apply visitor to TypeTree. - /** - * This function applies the given visitor to the given tree. Both visitor and tree may be const - * or non-const (if the compiler supports rvalue references, they may even be a non-const temporary). - * - * \note The visitor must implement the interface laid out by DefaultVisitor (most easily achieved by - * inheriting from it). - * - * \param tree The tree the visitor will be applied to. - * \param visitor The visitor to apply to the tree. - */ - template <typename Tree, typename Visitor> - void traverseTree(Tree&& tree, Visitor&& visitor) - { - using Node = std::remove_reference_t<Tree>; - using NodeTag = Dune::TypeTree::NodeTag<Node>; - using TreePath = Dune::TypeTree::HybridTreePath<>; - TraverseTree<NodeTag>::apply(std::forward<Tree>(tree), std::forward<Visitor>(visitor), TreePath{}); - } - -} // end namespace AMDiS diff --git a/src/amdis/utility/TreeData.hpp b/src/amdis/utility/TreeData.hpp index 183f5339..b913ad46 100644 --- a/src/amdis/utility/TreeData.hpp +++ b/src/amdis/utility/TreeData.hpp @@ -5,8 +5,8 @@ #include <utility> #include <vector> +#include <dune/typetree/traversal.hh> #include <dune/typetree/typetree.hh> -#include <amdis/utility/Visitor.hpp> namespace AMDiS { diff --git a/test/ExpressionsTest.cpp b/test/ExpressionsTest.cpp index bc544975..8c01c202 100644 --- a/test/ExpressionsTest.cpp +++ b/test/ExpressionsTest.cpp @@ -13,19 +13,20 @@ using namespace AMDiS; -using ElliptParam = YaspGridBasis<2, 2>; -using ElliptProblem = ProblemStat<ElliptParam>; - int main(int argc, char** argv) { AMDiS::init(argc, argv); using namespace Dune::Indices; - ElliptProblem prob("ellipt"); + using Grid = Dune::YaspGrid<2>; + auto gridPtr = MeshCreator<Grid>::create("mesh"); + auto basis = makeBasis(gridPtr->leafGridView(), lagrange<2>()); + + auto prob = makeProblemStat("ellipt", *gridPtr, basis); prob.initialize(INIT_ALL); - auto u = prob.getSolution(_0); + auto u = prob.getSolution(); // eval a functor at global coordinates (at quadrature points) auto expr1 = evalAtQP([](Dune::FieldVector<double, 2> const& x) { return x[0] + x[1]; }); @@ -42,7 +43,7 @@ int main(int argc, char** argv) auto expr8 = X(0); // Evaluation of the DOFVector (component) at quadrature points - auto expr9 = prob.getSolution(_0); + auto expr9 = prob.getSolution(); // --------------------------------- -- GitLab