diff --git a/doc/Mainpage.md b/doc/Mainpage.md index 3797679e04930f9431853114358495e70ddecc13..74909596b3dfe606455fc62f2efbacffe350fe6c 100644 --- a/doc/Mainpage.md +++ b/doc/Mainpage.md @@ -29,7 +29,7 @@ The corresponding weak form of the equation reads: with \f$ u\in V_g \f$. Thus, we need to define a grid (discretization of \f$ \Omega \f$) and a finite -element space (discretization of \f$ V_0 \f$ and \f$ V_g \f$). This needs to be +element space (discretization of \f$ V_0 \f$ and \f$ V_g \f$). This cab be provided as `Traits` parameter in the ProblemStat class: ~~~~~~~~~~~~~~~{.cpp} @@ -37,7 +37,7 @@ using Grid = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>; using Traits = LagrangeBasis<Grid::LeafGridView, 1>; ~~~~~~~~~~~~~~~ -This defines and `AlbertaGrid` as grid and a Lagrange Finite-Element space with +where `AlbertaGrid` defines a grid and `Traits` a Lagrange Finite-Element space with local polynomial degree 1 on the elements. All AMDiS programs start with initialization of the library, using `AMDiS::init` @@ -62,7 +62,7 @@ prob.initialize(INIT_ALL); Operators specify the (bi-)linear-form and the coefficient function in the term, see \ref operators for a list of possible types. The bilinear-form in the Poisson -equation consists of second-order term with coefficitn = 1 and the linear-form +equation consists of a second-order term with coefficient = 1 and the linear-form includes the function \f$ f(x)=-1 \f$: ~~~~~~~~~~~~~~~{.cpp} @@ -74,8 +74,8 @@ prob.addVectorOperator(opF, 0); ~~~~~~~~~~~~~~~ Boundary conditions, in the example above a Dirichlet condition, is specified by -defining a predicate for the boundary \f$ \partial\Omega \f$ and the values on the -boundary \f$ g(x) = 0 \f$: +defining a predicate for the boundary \f$ \Gamma\subset\partial\Omega \f$ and the +values on the boundary \f$ g(x) = 0 \f$: ~~~~~~~~~~~~~~~{.cpp} auto predicate = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8; }; @@ -93,7 +93,7 @@ adapt.adapt(); // assemble and solve Finally, finish the AMDiS program with `AMDiS::finish()`. -In complete program then reads: +The complete program then reads: ~~~~~~~~~~~~~~~{.cpp} #include <dune/amdis/AMDiS.hpp> #include <dune/amdis/AdaptInfo.hpp> diff --git a/dune/amdis/FileWriter.hpp b/dune/amdis/FileWriter.hpp index eedf48ff2fdcd40e0a2c8d66bdc60c29ac6bffed..2ef389ed12aeda33f7bd5da5453ca3ca42afc0f8 100644 --- a/dune/amdis/FileWriter.hpp +++ b/dune/amdis/FileWriter.hpp @@ -10,20 +10,48 @@ #include <dune/typetree/childextraction.hh> #include <dune/amdis/Initfile.hpp> +#include <dune/amdis/common/Size.hpp> +#include <dune/amdis/common/ValueCategory.hpp> +#include <dune/amdis/gridfunctions/DOFVectorView.hpp> #include <dune/amdis/io/FileWriterInterface.hpp> #include <dune/amdis/utility/Filesystem.hpp> namespace AMDiS { + namespace Impl + { + template <class Tag> struct VTKFieldTypeImpl; + template <> + struct VTKFieldTypeImpl<tag::scalar> { + static const Dune::VTK::FieldInfo::Type value = Dune::VTK::FieldInfo::Type::scalar; + }; + template <> + struct VTKFieldTypeImpl<tag::vector> { + static const Dune::VTK::FieldInfo::Type value = Dune::VTK::FieldInfo::Type::vector; + }; + template <> + struct VTKFieldTypeImpl<tag::matrix> { + static const Dune::VTK::FieldInfo::Type value = Dune::VTK::FieldInfo::Type::tensor; + }; + + } // end namespace Impl + + template <class Range> + constexpr Dune::VTK::FieldInfo::Type VTKFieldType = Impl::VTKFieldTypeImpl<ValueCategory_t<Range>>::value; + + template <class Range> + constexpr std::size_t VTKFieldSize = Size<Range>; + + template <class Traits, class TreePath> class FileWriter : public FileWriterInterface { private: // typedefs and static constants - using GlobalBasis = typename Traits::GlobalBasis; - using GridView = typename GlobalBasis::GridView; + using Vector = DOFVectorConstView<GlobalBasis,TreePath>; + using Range = typename Vector::Range; /// Dimension of the mesh static constexpr int dim = GridView::dimension; @@ -31,117 +59,57 @@ namespace AMDiS /// Dimension of the world static constexpr int dow = GridView::dimensionworld; - using Vector = DOFVector<GlobalBasis>; - - - public: /// Constructor. FileWriter(std::string baseName, - std::shared_ptr<GlobalBasis> const& basis, - TreePath const& tp, - std::shared_ptr<Vector> const& vector) + Vector const& dofvector) : FileWriterInterface(baseName) - , basis_(basis) - , treePath_(tp) - , vector_(vector) + , dofvector_(dofvector) { - //int subSampling = Parameters::get<int>(baseName + "->subsampling").value_or(0); + int m = Parameters::get<int>(baseName + "->ParaView mode").value_or(0); + mode_ = (m == 0 ? Dune::VTK::ascii : Dune::VTK::appendedraw); - //if (subSampling > 0) - // vtkWriter_ = std::make_shared<Dune::SubsamplingVTKWriter<GridView>>(basis_->gridView(), subSampling); - //else - vtkWriter_ = std::make_shared<Dune::VTKWriter<GridView>>(basis_->gridView()); + int subSampling = Parameters::get<int>(baseName + "->subsampling").value_or(0); - vtkSeqWriter_ = std::make_shared<Dune::VTKSequenceWriter<GridView>>(vtkWriter_, filename_, dir_, ""); + if (subSampling > 0) + vtkWriter_ = std::make_shared<Dune::SubsamplingVTKWriter<GridView>>(gridView(), subSampling); + else + vtkWriter_ = std::make_shared<Dune::VTKWriter<GridView>>(gridView()); - mode_ = Parameters::get<int>(baseName + "->ParaView mode").value_or(0); + vtkSeqWriter_ = std::make_shared<Dune::VTKSequenceWriter<GridView>>(vtkWriter_, filename_, dir_, ""); + vtkWriter_->addVertexData(dofvector_, Dune::VTK::FieldInfo(name_, VTKFieldType<Range>, VTKFieldSize<Range>)); } /// Implements \ref FileWriterInterface::writeFiles virtual void writeFiles(AdaptInfo& adaptInfo, bool force) override { - using Tree = typename GlobalBasis::LocalView::Tree; - using Node = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>; - - writeVertexData(typename Node::NodeTag{}, index_<Node::CHILDREN>, [&,this]() - { - if (mode_ == 0) - vtkSeqWriter_->write(adaptInfo.getTime(), Dune::VTK::ascii); - else if (mode_ == 1) - vtkSeqWriter_->write(adaptInfo.getTime(), Dune::VTK::appendedraw); - }); + vtkSeqWriter_->write(adaptInfo.getTime(), mode_); } protected: - template <class W> - void writeVertexData(Dune::TypeTree::LeafNodeTag, index_t<0>, W write) + GridView const& gridView() const { - using Dune::Functions::BasisBuilder::makeBasis; - using Dune::Functions::BasisBuilder::lagrange; - auto fct = makeDiscreteFunctionPtr(basis_,treePath_,vector_); - - auto p1basis = makeBasis(basis_->gridView(), lagrange<1>()); - auto data = makeDOFVector(p1basis, name_); - interpolate(p1basis, data, *fct); - - auto dataFct = Dune::Functions::makeDiscreteGlobalBasisFunction<double>(p1basis,data); - - vtkWriter_->addVertexData(dataFct, Dune::VTK::FieldInfo(name_, Dune::VTK::FieldInfo::Type::scalar, 1)); - write(); - vtkWriter_->clear(); - } - - template <std::size_t C, class W> - void writeVertexData(Dune::TypeTree::PowerNodeTag, index_t<C>, W write) - { - using Dune::Functions::BasisBuilder::makeBasis; - using Dune::Functions::BasisBuilder::lagrange; - using Dune::Functions::BasisBuilder::power; - using Dune::Functions::BasisBuilder::flatLexicographic; - - assert( C == dow ); - auto fct = makeDiscreteFunctionPtr(basis_,treePath_,vector_); - - auto p1basis = makeBasis(basis_->gridView(), power<C>(lagrange<1>(), flatLexicographic())); - auto data = makeDOFVector(p1basis, name_); - interpolate(p1basis, data, *fct); - - using Range = Dune::FieldVector<double,C>; - auto dataFct = Dune::Functions::makeDiscreteGlobalBasisFunction<Range>(p1basis,data); - - vtkWriter_->addVertexData(dataFct, Dune::VTK::FieldInfo(name_, Dune::VTK::FieldInfo::Type::vector, C)); - write(); - vtkWriter_->clear(); + return dofvector_.basis().gridView(); } - template <class NodeTag, std::size_t C, class W> - void writeVertexData(NodeTag, index_t<C>, W) {} - private: - std::shared_ptr<GlobalBasis> basis_; - TreePath treePath_; - std::shared_ptr<Vector> vector_; + Vector dofvector_; std::shared_ptr<Dune::VTKWriter<GridView>> vtkWriter_; std::shared_ptr<Dune::VTKSequenceWriter<GridView>> vtkSeqWriter_; - // std::vector<double> data_vector; - - // represents VTK::OutputType: 0...ascii, 1...appendedraw - int mode_; + // represents VTK::OutputType: ascii, appendedraw + Dune::VTK::OutputType mode_; }; - template <class Traits, class GlobalBasis, class TreePath, class Vector> + template <class Traits, class GlobalBasis, class TreePath> std::shared_ptr<FileWriter<Traits,TreePath>> makeFileWriterPtr( std::string baseName, - std::shared_ptr<GlobalBasis> const& basis, - TreePath const& tp, - std::shared_ptr<Vector> const& vector) + DOFVectorConstView<GlobalBasis,TreePath> const& dofvector) { - return std::make_shared<FileWriter<Traits,TreePath>>(baseName, basis, tp, vector); + return std::make_shared<FileWriter<Traits,TreePath>>(baseName, dofvector); } } // end namespace AMDiS diff --git a/dune/amdis/GridFunctionOperator.hpp b/dune/amdis/GridFunctionOperator.hpp index 2a2197f214eec52328c0c4eeb99ee7c9e563b7db..c60074feccf67ce20bc32b96f4395affd4da5bff 100644 --- a/dune/amdis/GridFunctionOperator.hpp +++ b/dune/amdis/GridFunctionOperator.hpp @@ -173,29 +173,42 @@ namespace AMDiS : GridFunctionOperatorBase<GridFct, QuadratureCreator>(gridFct, 0, quadCreator) {} - /// Assemble a local element vector on the element that is bound. + /// \brief Assemble a local element vector on the element that is bound. + /** + * \param contextGeometry Container for context related data + * \param quad A quadrature formula + * \param elementVector The output element vector + * \param node The node of the test-basis + **/ template <class Context, class QuadratureRule, class ElementVector, class Node> - void calculateElementVector(Context const&, //< container for context related data - QuadratureRule const&, //< a quadrature formula - ElementVector&, //< the output element vector - Node const& //< the node of the test-basis - ) + void calculateElementVector(Context const& contextGeometry, + QuadratureRule const& quad, + ElementVector& elementVector, + Node const& node) { error_exit("Needs to be implemented!"); } - /// Assemble a local element matrix on the element that is bound. + /// \brief Assemble a local element matrix on the element that is bound. + /** + * \param contextGeometry Container for context related data + * \param quad A quadrature formula + * \param elementMatrix The output element matrix + * \param rowNode The node of the test-basis + * \param colNode The node of the trial-basis + * \param flag1 finiteelement is the same for test and trial + * \param flag2 nodes are the same in the basis-tree + **/ template <class Context, class QuadratureRule, class ElementMatrix, class RowNode, class ColNode, bool sameFE, bool sameNode> - void calculateElementMatrix(Context const&, //< container for context related data - QuadratureRule const&, //< a quadrature formula - ElementMatrix&, //< the output element matrix - RowNode const&, //< the node of the test-basis - ColNode const&, //< the node of the trial-basis - std::integral_constant<bool, sameFE>, //< finiteelement is the same for test and trial - std::integral_constant<bool, sameNode> //< nodes are the same in the basis-tree - ) + void calculateElementMatrix(Context const& contextGeometry, + QuadratureRule const& quad, + ElementMatrix& elementMatrix, + RowNode const& rowNode, + ColNode const& colNode, + std::integral_constant<bool, sameFE> flag1, + std::integral_constant<bool, sameNode> flag2) { error_exit("Needs to be implemented!"); } @@ -261,28 +274,28 @@ namespace AMDiS /// Store tag and expression to create a \ref GridFunctionOperator template <class Tag, class Expr> - auto makeOperator(Tag t, Expr const& expr) + auto makeOperator(Tag tag, Expr const& expr) { using RawExpr = Underlying_t<Expr>; static_assert(Concepts::HasGridFunctionOrder<RawExpr>, "Polynomial degree of expression can not be deduced. You need to provide an explicit value for polynomial order or a quadrature rule in `makeOperator()`."); - return ExpressionPreOperator<Tag, Expr, tag::deduce>{t, expr}; + return ExpressionPreOperator<Tag, Expr, tag::deduce>{tag, expr}; } /// Store tag and expression and polynomial order of expression to create a \ref GridFunctionOperator template <class Tag, class Expr> - auto makeOperator(Tag t, Expr const& expr, int order, + auto makeOperator(Tag tag, Expr const& expr, int order, Dune::QuadratureType::Enum qt = Dune::QuadratureType::GaussLegendre) { - return ExpressionPreOperator<Tag, Expr, tag::order>{t, expr, order, qt}; + return ExpressionPreOperator<Tag, Expr, tag::order>{tag, expr, order, qt}; } /// Store tag and expression and a quadrature rule to create a \ref GridFunctionOperator template <class Tag, class Expr, class ctype, int dim> - auto makeOperator(Tag t, Expr const& expr, Dune::QuadratureRule<ctype,dim> const& rule) + auto makeOperator(Tag tag, Expr const& expr, Dune::QuadratureRule<ctype,dim> const& rule) { - return ExpressionPreOperator<Tag, Expr, tag::rule, Dune::QuadratureRule<ctype,dim>>{t, expr, rule}; + return ExpressionPreOperator<Tag, Expr, tag::rule, Dune::QuadratureRule<ctype,dim>>{tag, expr, rule}; } /** @} **/ diff --git a/dune/amdis/Operators.hpp b/dune/amdis/Operators.hpp index 01092cfb4c9ff951b733b0bde601396e15eb70a4..19f930d2e74a620199e7fe0ef7e055aef5358994 100644 --- a/dune/amdis/Operators.hpp +++ b/dune/amdis/Operators.hpp @@ -20,17 +20,17 @@ * \brief Defines operators to be assembled in the matrix/vector * * An `Operator` is a class providing methods necessary for assembling: - * - `void bind(Element, Geometry)` and `void unbind()` for binding an unbinding the + * - `bind(Element, Geometry)` and `unbind()` for binding an unbinding the * element to (from) an GridView entity of codim 0. Additionally the Geometry * object of the element is provided. - * - `Dune::QuadratureRule<ctype,dim> getQuadratureRule(Nodes...)` factory for the + * - `getQuadratureRule(Nodes...)` factory for the * quadrature rules used in assembling the operator on the element. `Nodes...` * is either `{RowNode, ColNode}` for Matrix-Operators or `{Node}` for a * Vector-Operator. - * - `void calculateElementVector(ContextGeometry, QuadratureRule, ElementVector, Node)` + * - `calculateElementVector(ContextGeometry, QuadratureRule, ElementVector, Node)` * where the `ContextGeometry` provides a reference to the ElementGeometry and * geometry of the LocalContext (that can be different), *or* - * - `void calculateElementMatrix(ContextGeometry, QuadratureRule, ElementMatrix, RowNode, ColNode, Flags...)` + * - `calculateElementMatrix(ContextGeometry, QuadratureRule, ElementMatrix, RowNode, ColNode, Flags...)` * Same as for `calculateElementVector` but additionally two optimization flags * are provided as `bool_t<...>` type: * + `sameFE`: the FiniteElementSpace of `RowNode` and `ColNode` are the same. diff --git a/dune/amdis/ProblemStat.inc.hpp b/dune/amdis/ProblemStat.inc.hpp index d8fb2227bddd768b68d22fa617d6ded17936a311..2ab868c47c34f2cb04fdb0a1a55404ed14a6bdcb 100644 --- a/dune/amdis/ProblemStat.inc.hpp +++ b/dune/amdis/ProblemStat.inc.hpp @@ -122,8 +122,8 @@ void ProblemStat<Traits>::createFileWriter() if (!Parameters::get<std::string>(componentName + "->filename")) return; - auto writer = makeFileWriterPtr<Traits>(componentName, globalBasis, treePath, solution); - filewriter.push_back(writer); + auto writer = makeFileWriterPtr<Traits>(componentName, this->getSolution(treePath)); + filewriter.push_back(std::move(writer)); }); } diff --git a/dune/amdis/common/FieldMatVec.hpp b/dune/amdis/common/FieldMatVec.hpp index 6469012201898d30b5bc9836a1dba443e39f6e6c..36703b412c4ec987041c9ebb2671ef624d62a02c 100644 --- a/dune/amdis/common/FieldMatVec.hpp +++ b/dune/amdis/common/FieldMatVec.hpp @@ -13,45 +13,48 @@ namespace AMDiS { - // some arithmetic operations with Dune::FieldVector + using Dune::FieldVector; + using Dune::FieldMatrix; + + // some arithmetic operations with FieldVector template <class T, int N, class S, REQUIRES(Concepts::Arithmetic<S>) > - Dune::FieldVector<T,N> operator*(Dune::FieldVector<T,N> v, S factor) + FieldVector<T,N> operator*(FieldVector<T,N> v, S factor) { return v *= factor; } template <class S, class T, int N, REQUIRES(Concepts::Arithmetic<S>) > - Dune::FieldVector<T,N> operator*(S factor, Dune::FieldVector<T,N> v) + FieldVector<T,N> operator*(S factor, FieldVector<T,N> v) { return v *= factor; } template <class T, int N, class S, REQUIRES(Concepts::Arithmetic<S>) > - Dune::FieldVector<T,N> operator/(Dune::FieldVector<T,N> v, S factor) + FieldVector<T,N> operator/(FieldVector<T,N> v, S factor) { return v /= factor; } - // some arithmetic operations with Dune::FieldMatrix + // some arithmetic operations with FieldMatrix template <class T, int N, int M> - Dune::FieldMatrix<T,N,M> operator+(Dune::FieldMatrix<T,N,M> A, Dune::FieldMatrix<T,N,M> const& B) + FieldMatrix<T,N,M> operator+(FieldMatrix<T,N,M> A, FieldMatrix<T,N,M> const& B) { return A += B; } template <class T, int N, int M> - Dune::FieldMatrix<T,N,M> operator-(Dune::FieldMatrix<T,N,M> A, Dune::FieldMatrix<T,N,M> const& B) + FieldMatrix<T,N,M> operator-(FieldMatrix<T,N,M> A, FieldMatrix<T,N,M> const& B) { return A -= B; } template <class T, int N, int M> - Dune::FieldVector<T,N> operator*(Dune::FieldMatrix<T,N,M> const& mat, Dune::FieldVector<T,M> const& vec) + FieldVector<T,N> operator*(FieldMatrix<T,N,M> const& mat, FieldVector<T,M> const& vec) { return Dune::FMatrixHelp::mult(mat, vec); } @@ -60,14 +63,14 @@ namespace AMDiS /// Cross-product a 2d-vector = orthogonal vector template <class T> - Dune::FieldVector<T, 2> cross(Dune::FieldVector<T, 2> const& a) + FieldVector<T, 2> cross(FieldVector<T, 2> const& a) { return {{ a[1], -a[0] }}; } /// Cross-product of two vectors (in 3d only) template <class T> - Dune::FieldVector<T, 3> cross(Dune::FieldVector<T, 3> const& a, Dune::FieldVector<T, 3> const& b) + FieldVector<T, 3> cross(FieldVector<T, 3> const& a, FieldVector<T, 3> const& b) { return {{ a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], @@ -76,7 +79,7 @@ namespace AMDiS /// Dot product (vec1^T * vec2) template <class T, class S, int N> - auto dot(Dune::FieldVector<T,N> const& vec1, Dune::FieldVector<S,N> const& vec2) + auto dot(FieldVector<T,N> const& vec1, FieldVector<S,N> const& vec2) { return vec1.dot(vec2); } @@ -86,7 +89,7 @@ namespace AMDiS namespace Impl { template <class T, int N, class Operation> - T accumulate(Dune::FieldVector<T, N> const& x, Operation op) + T accumulate(FieldVector<T, N> const& x, Operation op) { T result = 0; for (int i = 0; i < N; ++i) @@ -98,7 +101,7 @@ namespace AMDiS /// Sum of vector entires. template <class T, int N> - T sum(Dune::FieldVector<T, N> const& x) + T sum(FieldVector<T, N> const& x) { return Impl::accumulate(x, Operation::Plus{}); } @@ -106,7 +109,7 @@ namespace AMDiS /// Dot-product with the vector itself template <class T, int N> - auto unary_dot(Dune::FieldVector<T, N> const& x) + auto unary_dot(FieldVector<T, N> const& x) { auto op = [](auto const& a, auto const& b) { return a + Math::sqr(std::abs(b)); }; return Impl::accumulate(x, op); @@ -114,28 +117,28 @@ namespace AMDiS /// Maximum over all vector entries template <class T, int N> - auto max(Dune::FieldVector<T, N> const& x) + auto max(FieldVector<T, N> const& x) { return Impl::accumulate(x, Operation::Max{}); } /// Minimum over all vector entries template <class T, int N> - auto min(Dune::FieldVector<T, N> const& x) + auto min(FieldVector<T, N> const& x) { return Impl::accumulate(x, Operation::Min{}); } /// Maximum of the absolute values of vector entries template <class T, int N> - auto abs_max(Dune::FieldVector<T, N> const& x) + auto abs_max(FieldVector<T, N> const& x) { return Impl::accumulate(x, Operation::AbsMax{}); } /// Minimum of the absolute values of vector entries template <class T, int N> - auto abs_min(Dune::FieldVector<T, N> const& x) + auto abs_min(FieldVector<T, N> const& x) { return Impl::accumulate(x, Operation::AbsMin{}); } @@ -146,7 +149,7 @@ namespace AMDiS * \brief The 1-norm of a vector = sum_i |x_i| **/ template <class T, int N> - auto one_norm(Dune::FieldVector<T, N> const& x) + auto one_norm(FieldVector<T, N> const& x) { auto op = [](auto const& a, auto const& b) { return a + std::abs(b); }; return Impl::accumulate(x, op); @@ -156,7 +159,7 @@ namespace AMDiS * \brief The euklidean 2-norm of a vector = sqrt( sum_i |x_i|^2 ) **/ template <class T, int N> - auto two_norm(Dune::FieldVector<T, N> const& x) + auto two_norm(FieldVector<T, N> const& x) { return std::sqrt(unary_dot(x)); } @@ -165,7 +168,7 @@ namespace AMDiS * \brief The p-norm of a vector = ( sum_i |x_i|^p )^(1/p) **/ template <int p, class T, int N> - auto p_norm(Dune::FieldVector<T, N> const& x) + auto p_norm(FieldVector<T, N> const& x) { auto op = [](auto const& a, auto const& b) { return a + Math::pow<p>(std::abs(b)); }; return std::pow( Impl::accumulate(x, op), 1.0/p ); @@ -175,7 +178,7 @@ namespace AMDiS * \brief The infty-norm of a vector = max_i |x_i| = alias for \ref abs_max **/ template <class T, int N> - auto infty_norm(Dune::FieldVector<T, N> const& x) + auto infty_norm(FieldVector<T, N> const& x) { return abs_max(x); } @@ -184,7 +187,7 @@ namespace AMDiS /// The euklidean distance between two vectors = |lhs-rhs|_2 template <class T, int N> - T distance(Dune::FieldVector<T, N> const& lhs, Dune::FieldVector<T, N> const& rhs) + T distance(FieldVector<T, N> const& lhs, FieldVector<T, N> const& rhs) { T result = 0; for (int i = 0; i < N; ++i) @@ -196,9 +199,9 @@ namespace AMDiS /// Outer product (vec1 * vec2^T) template <class T, class S, int N, int M, int K> - auto outer(Dune::FieldMatrix<T,N,K> const& vec1, Dune::FieldMatrix<S,M,K> const& vec2) + auto outer(FieldMatrix<T,N,K> const& vec1, FieldMatrix<S,M,K> const& vec2) { - using result_type = Dune::FieldMatrix<decltype( std::declval<T>() * std::declval<S>() ), N, M>; + using result_type = FieldMatrix<decltype( std::declval<T>() * std::declval<S>() ), N, M>; result_type mat; for (int i = 0; i < N; ++i) for (int j = 0; j < M; ++j) @@ -209,28 +212,28 @@ namespace AMDiS // ---------------------------------------------------------------------------- template <class T> - T det(Dune::FieldMatrix<T, 0, 0> const& /*mat*/) + T det(FieldMatrix<T, 0, 0> const& /*mat*/) { return 0; } /// Determinant of a 1x1 matrix template <class T> - T det(Dune::FieldMatrix<T, 1, 1> const& mat) + T det(FieldMatrix<T, 1, 1> const& mat) { return mat[0][0]; } /// Determinant of a 2x2 matrix template <class T> - T det(Dune::FieldMatrix<T, 2, 2> const& mat) + T det(FieldMatrix<T, 2, 2> const& mat) { return mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0]; } /// Determinant of a 3x3 matrix template <class T> - T det(Dune::FieldMatrix<T, 3, 3> const& mat) + T det(FieldMatrix<T, 3, 3> const& mat) { return mat[0][0]*mat[1][1]*mat[2][2] + mat[0][1]*mat[1][2]*mat[2][0] + mat[0][2]*mat[1][0]*mat[2][1] - (mat[0][2]*mat[1][1]*mat[2][0] + mat[0][1]*mat[1][0]*mat[2][2] + mat[0][0]*mat[1][2]*mat[2][1]); @@ -238,14 +241,14 @@ namespace AMDiS /// Determinant of a NxN matrix template <class T, int N> - T det(Dune::FieldMatrix<T, N, N> const& mat) + T det(FieldMatrix<T, N, N> const& mat) { return mat.determinant(); } /// Return the inverse of the matrix `mat` template <class T, int N> - auto inv(Dune::FieldMatrix<T, N, N> mat) + auto inv(FieldMatrix<T, N, N> mat) { mat.invert(); return mat; @@ -253,7 +256,7 @@ namespace AMDiS /// Solve the linear system A*x = b template <class T, int N> - void solve(Dune::FieldMatrix<T, N, N> const& A, Dune::FieldVector<T, N>& x, Dune::FieldVector<T, N> const& b) + void solve(FieldMatrix<T, N, N> const& A, FieldVector<T, N>& x, FieldVector<T, N> const& b) { A.solve(x, b); } @@ -261,7 +264,7 @@ namespace AMDiS /// Gramian determinant = sqrt( det( DT^T * DF ) ) template <class T, int N, int M> - T gramian(Dune::FieldMatrix<T,N,M> const& DF) + T gramian(FieldMatrix<T,N,M> const& DF) { using std::sqrt; return sqrt( det(outer(DF, DF)) ); @@ -269,16 +272,16 @@ namespace AMDiS /// Gramian determinant, specialization for 1 column matrices template <class T, int M> - T gramian(Dune::FieldMatrix<T, 1, M> const& DF) + T gramian(FieldMatrix<T, 1, M> const& DF) { using std::sqrt; return sqrt(dot(DF[0], DF[0])); } template <class T, int M, int N> - Dune::FieldMatrix<T,N,M> trans(Dune::FieldMatrix<T, M, N> const& A) + FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A) { - Dune::FieldMatrix<T,N,M> At; + FieldMatrix<T,N,M> At; for (int i = 0; i < M; ++i) for (int j = 0; j < N; ++j) At[j][i] = A[i][j]; @@ -288,7 +291,7 @@ namespace AMDiS template <class T, int M, int N, int L> - Dune::FieldMatrix<T,M,N> multiplies(Dune::FieldMatrix<T, M, L> const& A, Dune::FieldMatrix<T, L, N> const& B) + FieldMatrix<T,M,N> multiplies(FieldMatrix<T, M, L> const& A, FieldMatrix<T, L, N> const& B) { return A.rightmultiplyany(B); } @@ -296,9 +299,9 @@ namespace AMDiS template <class T, int M, int N, int L> - Dune::FieldMatrix<T,M,N> multiplies_AtB(Dune::FieldMatrix<T, L, M> const& A, Dune::FieldMatrix<T, N, L> const& B) + FieldMatrix<T,M,N> multiplies_AtB(FieldMatrix<T, L, M> const& A, FieldMatrix<T, N, L> const& B) { - Dune::FieldMatrix<T,M,N> C; + FieldMatrix<T,M,N> C; for (int m = 0; m < M; ++m) { for (int n = 0; n < N; ++n) { @@ -311,9 +314,9 @@ namespace AMDiS } template <class T, int M, int N, int L> - Dune::FieldMatrix<T,M,N> multiplies_ABt(Dune::FieldMatrix<T, M, L> const& A, Dune::FieldMatrix<T, N, L> const& B) + FieldMatrix<T,M,N> multiplies_ABt(FieldMatrix<T, M, L> const& A, FieldMatrix<T, N, L> const& B) { - Dune::FieldMatrix<T,M,N> C; + FieldMatrix<T,M,N> C; for (int m = 0; m < M; ++m) { for (int n = 0; n < N; ++n) { @@ -326,7 +329,7 @@ namespace AMDiS } template <class T, int M, int N, int L> - Dune::FieldMatrix<T,M,N>& multiplies_ABt(Dune::FieldMatrix<T, M, L> const& A, Dune::FieldMatrix<T, N, L> const& B, Dune::FieldMatrix<T,M,N>& C) + FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, L> const& A, FieldMatrix<T, N, L> const& B, FieldMatrix<T,M,N>& C) { for (int m = 0; m < M; ++m) { for (int n = 0; n < N; ++n) { diff --git a/dune/amdis/common/Overloaded.hpp b/dune/amdis/common/Overloaded.hpp new file mode 100644 index 0000000000000000000000000000000000000000..f60381815e68c238b9ae34ec366e888aa33e643d --- /dev/null +++ b/dune/amdis/common/Overloaded.hpp @@ -0,0 +1,31 @@ +#pragma once + +#include <typet_traits> + +namespace AMDiS +{ +#ifdef AMDIS_HAS_CXX_VARIADIC_USING + template <class... Fs> + struct overloaded : Fs... + { + template <class... Fs_> + overloaded(Fs_&&... fs) + : Fs(std::forward<Fs_>(fs))... + {} + + using Fs::operator()...; + }; + +#ifdef DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION + template <class... Fs_> + overloaded(Fs_&&... fs) -> overloaded<std::decay_t<Fs_>...> +#endif + + template <class... Fs> + overloaded<std::decay_t<Fs>...> overload(Fs&&... fs) + { + return {std::forward<Fs>(fs)...}; + } +#endif + +} // end namespace AMDiS diff --git a/dune/amdis/common/Size.hpp b/dune/amdis/common/Size.hpp new file mode 100644 index 0000000000000000000000000000000000000000..fca48b7fa3d5c77f9b217af3839dd8f13fc74bf5 --- /dev/null +++ b/dune/amdis/common/Size.hpp @@ -0,0 +1,63 @@ +#pragma once + +#include <array> +#include <tuple> +#include <type_traits> + +#include <dune/amdis/common/ConceptsBase.hpp> +#include <dune/amdis/common/Mpl.hpp> +#include <dune/amdis/common/Utility.hpp> + +namespace Dune +{ + template <class T, int N> + class FieldVector; + + template <class T, int N, int M> + class FieldMatrix; +} + +namespace AMDiS +{ + // Get the number of elements in a tuple / pair / array / ... + //@{ + namespace Impl + { + template <class Tuple, class = Void_t<>> + struct SizeImpl : index_t<0> {}; + + template <class... Args> + struct SizeImpl<std::tuple<Args...>> + : index_t< sizeof...(Args) > {}; + + template <class Arg0, class Arg1> + struct SizeImpl<std::pair<Arg0, Arg1>> + : index_t<2> {}; + + template <class T, std::size_t N> + struct SizeImpl<std::array<T,N>> + : index_t<N> {}; + + template <class T, int N> + struct SizeImpl<Dune::FieldVector<T,N>, std::enable_if_t<(N!=1)>> + : index_t<N> {}; + + template <class T, int N, int M> + struct SizeImpl<Dune::FieldMatrix<T,N,M>, std::enable_if_t<(N!=1 || M!=1)>> + : index_t<N*M> {}; + + + // Specialization for arithmetic types + template <class T> + struct SizeImpl<T, std::enable_if_t<std::is_arithmetic<T>::value> > + : index_t<1> {}; + + } // end namespace Impl + //@} + + + template <class T> + constexpr std::size_t Size = Impl::SizeImpl<std::decay_t<T>>::value; + + +} // end namespace AMDiS diff --git a/dune/amdis/common/Utility.hpp b/dune/amdis/common/Utility.hpp index e2f012656c524259b2488e01d40eec1d6ddda60a..be6a05b5f690c738bd94ee7bc4a97c4af329da93 100644 --- a/dune/amdis/common/Utility.hpp +++ b/dune/amdis/common/Utility.hpp @@ -10,37 +10,6 @@ namespace AMDiS { - namespace Impl - { - // workaround for MSVC (problems with alias templates in pack expansion) - template <class, class T> - struct InvokeType { using type = T; }; - - template <class, class, class T> - struct InvokeType2 { using type = T; }; - - template <class T> - struct ValueType - { - using type = typename InvokeType<T, typename T::value_type>::type; - }; - } - - template <class T> - using Value_t = typename Impl::ValueType<T>::type; - - template <class T> - using Size_t = typename Impl::InvokeType<T, typename T::size_type>::type; - - template <class T> - using Result_t = typename Impl::InvokeType<T, typename T::result_type>::type; - - template <class T> - using Decay_t = typename Impl::InvokeType<T, typename std::decay<T>::type>::type; - - template <class T1, class T2> - using Common_t = typename Impl::InvokeType2<T1, T2, typename std::common_type<T1,T2>::type>::type; - namespace Impl { template <class T> @@ -68,53 +37,13 @@ namespace AMDiS struct Types {}; template <class... Ts> - using Types_t = Types<Decay_t<Ts>...>; - - - /// An identity type wrapper, represents the type itself - template <class T> - struct Id - { - using type = T; - }; - - template <class T> - using Id_t = typename Id<T>::type; + using Types_t = Types<std::decay_t<Ts>...>; /// Alias that indicates ownership of resources template <class T> using owner = T; - /// Dummy type for unknown return type - struct no_valid_type {}; - - - template <class F1, class F2> - struct overloaded : F1, F2 - { - overloaded(F1 x1, F2 x2) - : F1(x1), F2(x2) - {} - - using F1::operator(); - using F2::operator(); - }; - - template <class F> - F overload(F f) { return f; } - - template <class F1, class F2> - overloaded<F1, F2> overload(F1 f1, F2 f2) - { - return {f1, f2}; - } - - template <class F1, class... Fs> - auto overload(F1 f1, Fs... fs) - { - return overload(f1, overload(fs...)); - } // --------------------------------------------------------------------------- diff --git a/dune/amdis/common/ValueCategory.hpp b/dune/amdis/common/ValueCategory.hpp index a2d78a58124edf5fb33d10bf4585895ff263c29b..108d2f38edc7310ad54c0c08d93e05ef66d9501e 100644 --- a/dune/amdis/common/ValueCategory.hpp +++ b/dune/amdis/common/ValueCategory.hpp @@ -19,6 +19,7 @@ namespace AMDiS } // end namespace tag + /// Category of type T, e.g. scalar, vector matrix, specified by a tag template <class T, class = void> struct ValueCategory @@ -60,11 +61,12 @@ namespace AMDiS }; template <class T> - struct ValueCategory<std::reference_wrapper<T>> + struct ValueCategory< std::reference_wrapper<T> > { using type = typename ValueCategory<T>::type; }; + namespace Category { template <class T> diff --git a/dune/amdis/gridfunctions/DOFVectorView.hpp b/dune/amdis/gridfunctions/DOFVectorView.hpp index a43158dca671c7cba1efb5f6d36135a4f6de9171..06ecfab48e3a1fdcecbc4a8abc95bd230d4dcd83 100644 --- a/dune/amdis/gridfunctions/DOFVectorView.hpp +++ b/dune/amdis/gridfunctions/DOFVectorView.hpp @@ -20,13 +20,8 @@ namespace AMDiS * @{ **/ -#ifndef DOXYGEN - template <class GlobalBasisType, class TreePathType, bool isConst = true> - class DOFVectorView; -#endif - template <class GlobalBasisType, class TreePathType> - class DOFVectorView<GlobalBasisType, TreePathType, true> + class DOFVectorConstView { public: using GlobalBasis = GlobalBasisType; @@ -41,7 +36,7 @@ namespace AMDiS using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>; using Domain = typename EntitySet::GlobalCoordinate; - using Range = typename Vector::value_type; + using Range = RangeType<SubTree>; using DerivativeTraits = Dune::Functions::DefaultDerivativeTraits<Range(Domain)>; using DerivativeRange = typename DerivativeTraits::Range; @@ -76,7 +71,7 @@ namespace AMDiS using ReferenceGradientContainer = Dune::Functions::TreeData<SubTree, NodeData, true>; public: - GradientLocalFunction(DOFVectorView const& globalFunction) + GradientLocalFunction(DOFVectorConstView const& globalFunction) : globalFunction_(&globalFunction) , localBasisView_(globalFunction_->basis().localView()) , localIndexSet_(globalFunction_->basis().localIndexSet()) @@ -118,7 +113,7 @@ namespace AMDiS } private: - DOFVectorView const* globalFunction_; + DOFVectorConstView const* globalFunction_; LocalBasisView localBasisView_; LocalIndexSet localIndexSet_; @@ -136,8 +131,8 @@ namespace AMDiS class LocalFunction { public: - using Domain = typename DOFVectorView::LocalDomain; - using Range = typename DOFVectorView::Range; + using Domain = typename DOFVectorConstView::LocalDomain; + using Range = typename DOFVectorConstView::Range; private: using LocalBasisView = typename GlobalBasis::LocalView; @@ -152,7 +147,7 @@ namespace AMDiS using ShapeFunctionValueContainer = Dune::Functions::TreeData<SubTree, NodeData, true>; public: - LocalFunction(DOFVectorView const& globalFunction) + LocalFunction(DOFVectorConstView const& globalFunction) : globalFunction_(&globalFunction) , localBasisView_(globalFunction_->basis().localView()) , localIndexSet_(globalFunction_->basis().localIndexSet()) @@ -198,7 +193,7 @@ namespace AMDiS } private: - DOFVectorView const* globalFunction_; + DOFVectorConstView const* globalFunction_; LocalBasisView localBasisView_; LocalIndexSet localIndexSet_; SubTree const* subTree_; @@ -211,7 +206,7 @@ namespace AMDiS public: /// Constructor. Stores a pointer to the dofVector and a copy of the treePath. - DOFVectorView(DOFVector<GlobalBasis> const& dofVector, TreePath const& treePath) + DOFVectorConstView(DOFVector<GlobalBasis> const& dofVector, TreePath const& treePath) : dofVector_(&dofVector) , treePath_(treePath) , entitySet_(dofVector.getFeSpace().gridView()) @@ -226,7 +221,7 @@ namespace AMDiS } /// \brief Create a local function for this view on the DOFVector. \relates LocalFunction - friend LocalFunction localFunction(DOFVectorView const& self) + friend LocalFunction localFunction(DOFVectorConstView const& self) { return LocalFunction{self}; } @@ -266,28 +261,28 @@ namespace AMDiS // A mutable version of DOFVectorView template <class GlobalBasisType, class TreePathType> - class DOFVectorView<GlobalBasisType, TreePathType, false> - : public DOFVectorView<GlobalBasisType, TreePathType, true> + class DOFVectorMutableView + : public DOFVectorConstView<GlobalBasisType, TreePathType> { - using DOFVectorConstView = DOFVectorView<GlobalBasisType, TreePathType, true>; + using Super = DOFVectorConstView<GlobalBasisType, TreePathType>; using GlobalBasis = GlobalBasisType; using TreePath = TreePathType; public: /// Constructor. Stores a pointer to the mutable `dofvector`. - DOFVectorView(DOFVector<GlobalBasis>& dofVector, TreePath const& treePath) - : DOFVectorConstView(dofVector, treePath) + DOFVectorMutableView(DOFVector<GlobalBasis>& dofVector, TreePath const& treePath) + : Super(dofVector, treePath) , mutableDofVector_(&dofVector) {} public: /// Interpolation of GridFunction to DOFVector template <class Expr> - DOFVectorView& interpolate(Expr&& expr) + DOFVectorMutableView& interpolate(Expr&& expr) { - auto const& basis = DOFVectorConstView::basis(); - auto const& treePath = DOFVectorConstView::treePath(); + auto const& basis = Super::basis(); + auto const& treePath = Super::treePath(); auto&& gridFct = makeGridFunction(std::forward<Expr>(expr), basis.gridView()); @@ -303,7 +298,7 @@ namespace AMDiS DOFVector<GlobalBasis>& coefficients() { return *mutableDofVector_; } /// Return the const DOFVector - using DOFVectorConstView::coefficients; + using Super::coefficients; protected: DOFVector<GlobalBasis>* mutableDofVector_; @@ -317,14 +312,14 @@ namespace AMDiS template <class GlobalBasis, class TreePath> auto makeDOFVectorView(DOFVector<GlobalBasis> const& dofVector, TreePath const& treePath) { - return DOFVectorView<GlobalBasis, TreePath, true>{dofVector, treePath}; + return DOFVectorConstView<GlobalBasis, TreePath>{dofVector, treePath}; } // A Generator for a mutable \ref DOFVectorView. template <class GlobalBasis, class TreePath> auto makeDOFVectorView(DOFVector<GlobalBasis>& dofVector, TreePath const& treePath) { - return DOFVectorView<GlobalBasis, TreePath, false>{dofVector, treePath}; + return DOFVectorMutableView<GlobalBasis, TreePath>{dofVector, treePath}; } @@ -333,7 +328,7 @@ namespace AMDiS auto makeDOFVectorView(DOFVector<GlobalBasis> const& dofVector) { auto treePath = Dune::TypeTree::hybridTreePath(); - return DOFVectorView<GlobalBasis, decltype(treePath), true>{dofVector, treePath}; + return DOFVectorConstView<GlobalBasis, decltype(treePath)>{dofVector, treePath}; } // A Generator for a mutable \ref DOFVectorView. @@ -341,7 +336,7 @@ namespace AMDiS auto makeDOFVectorView(DOFVector<GlobalBasis>& dofVector) { auto treePath = Dune::TypeTree::hybridTreePath(); - return DOFVectorView<GlobalBasis, decltype(treePath), false>{dofVector, treePath}; + return DOFVectorMutableView<GlobalBasis, decltype(treePath)>{dofVector, treePath}; } #endif diff --git a/dune/amdis/gridfunctions/DOFVectorView.inc.hpp b/dune/amdis/gridfunctions/DOFVectorView.inc.hpp index 141aac5680689f866c9f7195b40677c6ef387b81..4b2c9fed92de728ca6c8db23cc3309813e6fd5ba 100644 --- a/dune/amdis/gridfunctions/DOFVectorView.inc.hpp +++ b/dune/amdis/gridfunctions/DOFVectorView.inc.hpp @@ -5,7 +5,7 @@ namespace AMDiS { template <class GlobalBasis, class TreePath> -typename DOFVectorView<GlobalBasis, TreePath, true>::Range DOFVectorView<GlobalBasis, TreePath, true>:: +typename DOFVectorConstView<GlobalBasis, TreePath>::Range DOFVectorConstView<GlobalBasis, TreePath>:: LocalFunction::operator()(LocalDomain const& x) const { assert( bound_ ); @@ -59,7 +59,7 @@ LocalFunction::operator()(LocalDomain const& x) const template <class GlobalBasis, class TreePath> -typename DOFVectorView<GlobalBasis, TreePath, true>::DerivativeRange DOFVectorView<GlobalBasis, TreePath, true>:: +typename DOFVectorConstView<GlobalBasis, TreePath>::DerivativeRange DOFVectorConstView<GlobalBasis, TreePath>:: GradientLocalFunction::operator()(LocalDomain const& x) const { assert( bound_ ); diff --git a/dune/amdis/operations/Arithmetic.hpp b/dune/amdis/operations/Arithmetic.hpp index 5f7a45684296cbc71ff7c5781aaf20e1abb6e197..95f1159d9964771d730716d6b5825116ad88167f 100644 --- a/dune/amdis/operations/Arithmetic.hpp +++ b/dune/amdis/operations/Arithmetic.hpp @@ -194,15 +194,17 @@ namespace AMDiS using Sqr = Pow<2>; + /// \see Pow template <> struct Pow<1> : public Id {}; + /// \see Pow template <> struct Pow<0> : public Zero {}; - /// Functor that represents x^p + /// Functor that represents x^p, \see \ref Pow struct Pow_ { Pow_(int p) diff --git a/dune/amdis/utility/TreePath.hpp b/dune/amdis/utility/TreePath.hpp index 264513808a939677b3601ed6053585d710b0b292..2fb4ce335493988bd520a57a1b971bcddca81e58 100644 --- a/dune/amdis/utility/TreePath.hpp +++ b/dune/amdis/utility/TreePath.hpp @@ -3,9 +3,11 @@ #include <sstream> #include <string> +#include <dune/common/std/apply.hh> #include <dune/typetree/treepath.hh> #include <dune/typetree/typetraits.hh> + namespace AMDiS { auto makeTreePath(int i) { return Dune::TypeTree::hybridTreePath(std::size_t(i)); } @@ -33,6 +35,16 @@ namespace AMDiS return Dune::TypeTree::hybridTreePath(); } + template <class... T> + auto makeDynamicTreePath(Dune::TypeTree::HybridTreePath<T...> const& tp) + { + auto tmp = Dune::TypeTree::MakeableDynamicTreePath<sizeof...(T)>().mutablePath(); + forEach(range_<0,sizeof...(T)>, [&tp,&tmp](auto _i) mutable { + tmp.push_back(std::size_t(std::get<_i>(tp._data))); + }); + return tmp.view(); + } + namespace Impl {