diff --git a/CMakeLists.txt b/CMakeLists.txt index ede6de2c358b43d67854d0a3d230fb97eb95c86b..7de534f765097fc7737a9b63b410f2eb95141495 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,7 +1,7 @@ cmake_minimum_required(VERSION 3.1) project(dune-amdis CXX) -set(CXX_MAX_STANDARD 14 CACHE BOOL "" FORCE) +#set(CXX_MAX_STANDARD 14 CACHE BOOL "" FORCE) if(NOT (dune-common_DIR OR dune-common_ROOT OR "${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*")) @@ -20,11 +20,34 @@ include(DuneMacros) # start a dune project with information from dune.module dune_project() -dune_enable_all_packages() +dune_enable_all_packages(MODULE_LIBRARIES amdis) + +include(AmdisMacros) add_subdirectory("src") add_subdirectory("dune") add_subdirectory("doc") add_subdirectory("cmake/modules") + +# some additional packages and flags +add_dune_alberta_flags("amdis" OBJECT USE_GENERIC) +target_compile_definitions("amdis" PUBLIC AMDIS_BACKEND_MTL=1) +target_compile_options("amdis" PUBLIC -ftemplate-backtrace-limit=0 -Wall -pedantic -Wno-unused-parameter) + +find_package(MTL REQUIRED + PATHS /usr/local/lib/mtl4 /opt/sources/mtl4 /opt/development/mtl4) +if (MTL_FOUND) + target_include_directories("amdis" PUBLIC ${MTL_INCLUDE_DIRS}) + + set (CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL") + foreach (feature ${CXX_ELEVEN_FEATURE_LIST}) + target_compile_definitions("amdis" PUBLIC MTL_WITH_${feature}) + endforeach () + + if (HAVE_UMFPACK OR ENABLE_SUITESPARSE OR SuiteSparse_FOUND) + target_compile_definitions("amdis" PUBLIC MTL_HAS_UMFPACK) + endif () +endif (MTL_FOUND) + # finalize the dune project, e.g. generating config.h etc. finalize_dune_project(GENERATE_CONFIG_H_CMAKE) diff --git a/cmake/modules/AMDiSCXXFeatures.cmake b/cmake/modules/AMDiSCXXFeatures.cmake new file mode 100644 index 0000000000000000000000000000000000000000..c6b5831f4dd093e8024b4fc8a743fd8c60c4d669 --- /dev/null +++ b/cmake/modules/AMDiSCXXFeatures.cmake @@ -0,0 +1,33 @@ +#include(CheckIncludeFileCXX) +include(CheckCXXSourceCompiles) +#include(CheckCXXSymbolExists) + +# fold expressions (a + ...) +check_cxx_source_compiles(" + template + int f(Args... args) + { + return (args + ...); + } + int main() + { + f(0,1,2,3,4,5); + } +" AMDIS_HAS_CXX_FOLD_EXPRESSION +) + +check_cxx_source_compiles(" + template + auto f() + { + if constexpr(I == 0) + return double{1.0}; + else + return int{0}; + } + int main() + { + return f<1>(); + } +" AMDIS_HAS_CXX_CONSTEXPR_IF +) \ No newline at end of file diff --git a/cmake/modules/AmdisMacros.cmake b/cmake/modules/AmdisMacros.cmake index 613dfb664b75999b008f0003a4d7c409cbf409ad..641e272763ba71437d5073fc74823335ee82240c 100644 --- a/cmake/modules/AmdisMacros.cmake +++ b/cmake/modules/AmdisMacros.cmake @@ -1 +1,3 @@ # File for module specific CMake tests. + +include(AMDiSCXXFeatures) \ No newline at end of file diff --git a/config.h.cmake b/config.h.cmake index 785daa80c98e01eeb81b556489d4f413f1facb4b..e53392bc7ccd643d244e463b52e0df06dcb2fa1b 100644 --- a/config.h.cmake +++ b/config.h.cmake @@ -40,6 +40,12 @@ /* Define to the revision of amdis */ #define AMDIS_VERSION_REVISION @AMDIS_VERSION_REVISION@ -/* end amdis + +/* some detected compiler features may be used in AMDiS */ +#cmakedefine AMDIS_HAS_CXX_FOLD_EXPRESSIONS 1 + +#cmakedefine AMDIS_HAS_CXX_CONSTEXPR_IF 1 + +/* end dune-amdis Everything below here will be overwritten */ diff --git a/dune/amdis/AdaptInfo.hpp b/dune/amdis/AdaptInfo.hpp index 5c4848e0ccc671050a4b19e138e7c7ac2c85bc65..019d94e34a5bdae145e5aea386ae684980596b3c 100644 --- a/dune/amdis/AdaptInfo.hpp +++ b/dune/amdis/AdaptInfo.hpp @@ -9,7 +9,7 @@ // AMDiS includes #include -#include +#include namespace AMDiS { diff --git a/dune/amdis/Assembler.inc.hpp b/dune/amdis/Assembler.inc.hpp index c9e19143487b065e08a00cbcd0e36efab9ce61e5..7d39929038057cc4f7c5a69da8c0f9531a2c853b 100644 --- a/dune/amdis/Assembler.inc.hpp +++ b/dune/amdis/Assembler.inc.hpp @@ -4,7 +4,7 @@ #include #include -#include +#include namespace AMDiS { diff --git a/dune/amdis/CMakeLists.txt b/dune/amdis/CMakeLists.txt index 0dda28e576d3f1bfe5f6756ab32c10dc6f3715c7..a0cfbd893df73cc693622e57b74b1837286cd9ca 100644 --- a/dune/amdis/CMakeLists.txt +++ b/dune/amdis/CMakeLists.txt @@ -1,6 +1,6 @@ #install headers -dune_add_library("duneamdis" NO_EXPORT +dune_library_add_sources(amdis SOURCES AdaptBase.cpp AdaptInfo.cpp AdaptInstationary.cpp @@ -11,33 +11,7 @@ dune_add_library("duneamdis" NO_EXPORT # ProblemInstat.cpp ProblemStat.cpp StandardProblemIteration.cpp - #linear_algebra/istl/SystemMatrix.cpp - #linear_algebra/istl/SystemVector.cpp - linear_algebra/mtl/SystemMatrix.cpp - linear_algebra/mtl/SystemVector.cpp - utility/Filesystem.cpp - ) -add_dune_alberta_flags("duneamdis" OBJECT USE_GENERIC) -target_compile_definitions("duneamdis" PUBLIC AMDIS_BACKEND_MTL=1) -target_compile_options("duneamdis" PUBLIC -ftemplate-backtrace-limit=0 -Wall -pedantic -Wno-unused-parameter) - -find_package(MTL REQUIRED - PATHS /usr/local/lib/mtl4) -if (MTL_FOUND) - target_include_directories("duneamdis" PUBLIC ${MTL_INCLUDE_DIRS}) -# target_link_libraries("duneamdis" PUBLIC ${MTL_LIBRARIES}) -# target_compile_options("duneamdis" PUBLIC ${MTL_CXX_DEFINITIONS}) - - set (CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL") - foreach (feature ${CXX_ELEVEN_FEATURE_LIST}) - target_compile_definitions("duneamdis" PUBLIC MTL_WITH_${feature}) - endforeach () - - if (HAVE_UMFPACK OR ENABLE_SUITESPARSE OR SuiteSparse_FOUND) - target_compile_definitions("duneamdis" PUBLIC MTL_HAS_UMFPACK) - endif () -endif (MTL_FOUND) - +) install(FILES AdaptBase.hpp @@ -45,73 +19,40 @@ install(FILES AdaptInstationary.hpp AdaptStationary.hpp AMDiS.hpp + Assembler.hpp + Assembler.inc.hpp + ContextGeometry.hpp CreatorInterface.hpp CreatorMap.hpp DirichletBC.hpp DirichletBC.inc.hpp FileWriter.hpp Flag.hpp + GridFunctionOperator.hpp + GridFunctions.hpp Initfile.hpp LinearAlgebra.hpp - Log.hpp - Math.hpp + LocalAssembler.hpp + LocalAssemblerBase.hpp Mesh.hpp - OperatorEvaluation.hpp - Operator.hpp - Operator.inc.hpp - OperatorTermBase.hpp - OperatorTerm.hpp - ProblemInstatBase.hpp + Operators.hpp + Output.hpp ProblemInstat.hpp ProblemInstat.inc.hpp + ProblemInstatBase.hpp ProblemInterationInterface.hpp - ProblemStatBase.hpp ProblemStat.hpp ProblemStat.inc.hpp + ProblemStatBase.hpp + ProblemStatTraits.hpp ProblemTimeInterface.hpp StandardProblemIteration.hpp - TermGenerator.hpp - common/ConceptsBase.hpp - common/ClonablePtr.hpp - common/IndexSeq.hpp - common/Literals.hpp - common/Loops.hpp - common/Mpl.hpp - common/ScalarTypes.hpp - common/Timer.hpp - common/Traits.hpp - common/TupleUtility.hpp - common/Utility.hpp - common/ValueCategory.hpp - linear_algebra/LinearAlgebraBse.hpp - linear_algebra/LinearSolverInterface.hpp - linear_algebra/PreconditionerInterface.hpp - linear_algebra/RunnerInterface.hpp - linear_algebra/SolverInfo.hpp - linear_algebra/istl/DOFMatrix.hpp - linear_algebra/istl/DOFVector.hpp - linear_algebra/istl/ISTL_Preconditioner.hpp - linear_algebra/istl/ISTLRunner.hpp - linear_algebra/istl/ISTL_Solver.hpp - linear_algebra/istl/LinearSolver.hpp - linear_algebra/istl/Preconditioner.hpp - linear_algebra/istl/SystemMatrix.hpp - linear_algebra/istl/SystemVector.hpp - linear_algebra/mtl/BITL_Preconditioner.hpp - linear_algebra/mtl/BITL_Solver.hpp - linear_algebra/mtl/BlockMTLMatrix.hpp - linear_algebra/mtl/BlockMTLVector.hpp - linear_algebra/mtl/Copy.hpp - linear_algebra/mtl/DOFMatrix.hpp - linear_algebra/mtl/DOFVector.hpp - linear_algebra/mtl/ITL_Preconditioner.hpp - linear_algebra/mtl/ITL_Solver.hpp - linear_algebra/mtl/KrylovRunner.hpp - linear_algebra/mtl/LinearSolver.hpp - linear_algebra/mtl/Mapper.hpp - linear_algebra/mtl/MTLDenseVector.hpp - linear_algebra/mtl/Preconditioner.hpp - linear_algebra/mtl/SystemMatrix.hpp - linear_algebra/mtl/SystemVector.hpp - linear_algebra/mtl/UmfpackRunner.hpp DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis) + +add_subdirectory("assembler") +add_subdirectory("common") +add_subdirectory("gridfunctions") +add_subdirectory("io") +add_subdirectory("linear_algebra") +add_subdirectory("operations") +add_subdirectory("utility") \ No newline at end of file diff --git a/dune/amdis/OperatorBase.hpp b/dune/amdis/GridFunctionOperator.hpp similarity index 58% rename from dune/amdis/OperatorBase.hpp rename to dune/amdis/GridFunctionOperator.hpp index 17d161640f6502b1b3f1ec128b9c9cc820cd32f0..5e0e704b2aedd123adddc58ed6bf9e3857bf8880 100644 --- a/dune/amdis/OperatorBase.hpp +++ b/dune/amdis/GridFunctionOperator.hpp @@ -1,15 +1,22 @@ #pragma once #include +#include -#include +#include #include namespace AMDiS { - template - class ExpressionOperatorBase + template + class GridFunctionOperatorBase { + using LocalFunction = decltype(localFunction(std::declval())); + using Element = typename GridFunction::EntitySet::Element; + using Geometry = typename Element::Geometry; + + using LocalCoordinate = typename GridFunction::EntitySet::LocalCoordinate; + public: /// \brief Constructor. Stores a copy of `expr`. /** @@ -17,8 +24,9 @@ namespace AMDiS * \ref ExpressionBase, and stores a copy. Additionally, it gets the * differentiation order, to calculate the quadrature degree in \ref getDegree. **/ - ExpressionOperatorBase(Expr const& expr, int order) - : expr_(expr) + GridFunctionOperatorBase(GridFunction const& gridFct, int order) + : gridFct_(gridFct) + , localFct_(localFunction(gridFct_)) , order_(order) {} @@ -28,13 +36,12 @@ namespace AMDiS * Since all operators need geometry information, the `Element::Geometry` * object `geometry` is created once, during grid traversal, and passed in. * - * By default, it binds the \ref expr_ to the `element`. + * By default, it binds the \ref localFct_ to the `element`. **/ - template void bind(Element const& element, Geometry const& geometry) { if (!bound_) { - expr_.bind(element); + localFct_.bind(element); isSimplex_ = geometry.type().isSimplex(); isAffine_ = geometry.affine(); bound_ = true; @@ -45,32 +52,14 @@ namespace AMDiS void unbind() { bound_ = false; - initSize_ = 0; - expr_.unbind(); - } - - /// \brief Evaluate expression at quadrature points. - /** - * Initialization at all points `quad[iq].position()` where - * `iq` in `[0, quad.size())`. - * - * \param context Container for localContext, localGeometry and geometry, see \ref LocalAssembler - * \param quad Liste of evaluations points. Follows interface of \ref Dune::QuadratureRule - **/ - template - void init(Context const& context, Quad const& quad) - { - if (initSize_ != quad.size()) { - expr_.init(context, quad); - initSize_ = quad.size(); - } + localFct_.unbind(); } /// Return expressions iq'th value. Must be initialized in \ref init before. - auto eval(std::size_t iq) const + auto coefficient(LocalCoordinate const& local) const { - assert( initSize_ > 0u ); - return expr_[iq]; + assert( bound_ ); + return localFct_(local); } @@ -86,9 +75,9 @@ namespace AMDiS assert( bound_ ); int psiDegree = getPolynomialDegree(node); - int termDegree = expr_.getDegree(); + int coeffDegree = order(localFct_); - int degree = psiDegree + termDegree; + int degree = psiDegree + coeffDegree; if (isSimplex_) degree -= order_; if (isAffine_) @@ -111,9 +100,9 @@ namespace AMDiS int psiDegree = getPolynomialDegree(rowNode); int phiDegree = getPolynomialDegree(colNode); - int termDegree = expr_.getDegree(); + int coeffDegree = order(localFct_); - int degree = psiDegree + phiDegree + termDegree; + int degree = psiDegree + phiDegree + coeffDegree; if (isSimplex_) degree -= order_; if (isAffine_) @@ -122,45 +111,31 @@ namespace AMDiS return degree; } - - protected: - - Expr& expr() - { - return expr_; - } - - Expr const& expr() const - { - return expr_; - } - - private: + GridFunction gridFct_; + LocalFunction localFct_; - Expr expr_; //< the stores expression, see \ref ExpressionBase int order_; //< the derivative order of this operator bool isSimplex_ = false; //< the bound element is a simplex bool isAffine_ = false; //< the bound geometry is affine bool bound_ = false; //< an element is bound to the operator - std::size_t initSize_ = 0; //< number of quadrature points, the expression was evaluated at. }; - /// A default implementation of an ExpressionOperator if no specialization is available. + /// A default implementation of an GridFunctionOperator if no specialization is available. /** * An operator must implement either \ref calculateElementVector, or * \ref calculateElementMatrix, if it is a vector or matrix operator, * respectively. **/ - template - class ExpressionOperator - : public ExpressionOperatorBase + template + class GridFunctionOperator + : public GridFunctionOperatorBase { public: - ExpressionOperator(Tag, Expr const& expr) - : ExpressionOperatorBase(expr, 0) + GridFunctionOperator(Tag, GridFct const& gridFct) + : GridFunctionOperatorBase(gridFct, 0) {} /// Assemble a local element vector on the element that is bound. @@ -191,21 +166,50 @@ namespace AMDiS } }; + template + struct ExpressionPreOperator + { + Tag tag; + Expr expr; + }; + + /// Store tag and expression in struct + template + auto makeOperator(Tag t, Expr const& expr) + { + return ExpressionPreOperator{t, expr}; + } + + /// Generate an \ref GridFunctionOperator from a PreOperator (tag, expr). + template + auto makeGridOperator(ExpressionPreOperator const& op, GridView const& gridView) + { + auto gridFct = makeGridFunction(op.expr, gridView, Dune::PriorityTag<42>{}); + return GridFunctionOperator{op.tag, gridFct}; + } + + template + auto makeGridOperator(std::reference_wrapper> op, GridView const& gridView) + { + ExpressionPreOperator const& op_ref = op; + auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView, Dune::PriorityTag<42>{}); + return GridFunctionOperator{op_ref.tag, gridFct}; + } - /// Generate an \ref ExpressionOperator with a given tag `t` and pre-expression `epxr`. - template - auto makeOperator(Tag t, PreExpr const& expr) + /// Generate a shared_ptr to \ref GridFunctionOperator from a PreOperator (tag, expr). + template + auto makeGridOperatorPtr(ExpressionPreOperator const& op, GridView const& gridView) { - using Expr = ToTerm_t; - return ExpressionOperator{t, toTerm(expr)}; + auto gridFct = makeGridFunction(op.expr, gridView, Dune::PriorityTag<42>{}); + return std::make_shared>(op.tag, gridFct); } - /// Generate a shared_ptr to \ref ExpressionOperator with a given tag `t` and pre-expression `epxr`. - template - auto makeOperatorPtr(Tag t, PreExpr const& expr) + template + auto makeGridOperatorPtr(std::reference_wrapper> op, GridView const& gridView) { - using Expr = ToTerm_t; - return std::make_shared>(t, toTerm(expr)); + ExpressionPreOperator const& op_ref = op; + auto gridFct = makeGridFunction(std::ref(op_ref.expr), gridView, Dune::PriorityTag<42>{}); + return std::make_shared>(op_ref.tag, gridFct); } } // end namespace AMDiS diff --git a/dune/amdis/GridFunctions.hpp b/dune/amdis/GridFunctions.hpp new file mode 100644 index 0000000000000000000000000000000000000000..acbfda0f4ea419450a1ce0593fb3a73f0a53f841 --- /dev/null +++ b/dune/amdis/GridFunctions.hpp @@ -0,0 +1,8 @@ +#pragma once + +#include +#include +#include +#include +#include +#include diff --git a/dune/amdis/Initfile.hpp b/dune/amdis/Initfile.hpp index b548e36ec5bef214861035c09c12fbe991dd434b..84480f6b77372babcc500869c4f32ed686e1c100 100644 --- a/dune/amdis/Initfile.hpp +++ b/dune/amdis/Initfile.hpp @@ -16,7 +16,7 @@ #include #include -#include +#include namespace AMDiS { diff --git a/dune/amdis/LocalAssembler.hpp b/dune/amdis/LocalAssembler.hpp index 17884fb2f330d11ad2cebfd261db17957a98c927..0519ed73f24ac52195eda2fac1f9182623bc4612 100644 --- a/dune/amdis/LocalAssembler.hpp +++ b/dune/amdis/LocalAssembler.hpp @@ -189,7 +189,7 @@ namespace AMDiS private: - std::unique_ptr storage_; //< the stored operator, implementing \ref ExpressionOperatorBase + std::unique_ptr storage_; //< the stored operator, implementing \ref GridFunctionOperatorBase Operator& op_; Element const* element_ = nullptr; diff --git a/dune/amdis/Operations.hpp b/dune/amdis/Operations.hpp new file mode 100644 index 0000000000000000000000000000000000000000..d798bcbc32046fc68386650ed46634eedf826507 --- /dev/null +++ b/dune/amdis/Operations.hpp @@ -0,0 +1,8 @@ +#pragma once + +#include +#include +#include +#include +#include +#include diff --git a/dune/amdis/Operators.hpp b/dune/amdis/Operators.hpp index 82fe1a2b12a7734097f566f9a4c924d527ec2918..6ab4d429c80aca52e8bc919eacbda5971b78c613 100644 --- a/dune/amdis/Operators.hpp +++ b/dune/amdis/Operators.hpp @@ -1,7 +1,7 @@ #pragma once #include -#include +#include /** * In the following comments we use the notation @@ -15,25 +15,25 @@ */ // zero-order operators -#include // -#include // -#include // -#include // -#include // -#include // , +#include // +#include // +#include // +#include // +#include // +#include // , // first-order operators -#include // -#include // -#include // -#include // -#include // -#include // -#include // -#include // +#include // +#include // +#include // +#include // +#include // +#include // +#include // +#include // // second-order operators -#include // -#include // -#include // +#include // +#include // +#include // diff --git a/dune/amdis/ProblemStat.hpp b/dune/amdis/ProblemStat.hpp index 2719f0ea6e7b541f95bdf1b18b7511d766caa0fe..78228307ee6279a72c1c123a6a4f85962d3f8286 100644 --- a/dune/amdis/ProblemStat.hpp +++ b/dune/amdis/ProblemStat.hpp @@ -31,6 +31,8 @@ #include #include +#include + #include #include @@ -142,21 +144,6 @@ namespace AMDiS RowTreePath row, ColTreePath col, Values const& values); - - private: // implementation detail - - template - void addDirichletBCImpl(Predicate const& predicate, - RowNode const& rowNode, ColNode const& colNode, - RangeType const& value); - - template - std::enable_if_t>::value> - addDirichletBCImpl(Predicate const& predicate, - RowNode const& rowNode, ColNode const& colNode, - Values const& values); - - public: /// Implementation of \ref ProblemStatBase::solve @@ -233,7 +220,7 @@ namespace AMDiS auto getSolution(TreePath path) const { auto tp = makeTreePath(path); - return makeDiscreteFunctionPtr(globalBasis, tp, solution); + return makeDOFVectorView(*solution, tp); } diff --git a/dune/amdis/ProblemStat.inc.hpp b/dune/amdis/ProblemStat.inc.hpp index 6cb0bfe00749a967ca66b3361ec3842e7ef6dfa6..5be3139ba10d72f3fc202e00b23a42a3839eb64b 100644 --- a/dune/amdis/ProblemStat.inc.hpp +++ b/dune/amdis/ProblemStat.inc.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include @@ -132,13 +133,14 @@ void ProblemStat::createFileWriter() template template void ProblemStat::addMatrixOperator( - Operator const& op, + Operator const& preOp, RowTreePath row, ColTreePath col, double* factor, double* estFactor) { auto i = child(globalBasis->localView().tree(), makeTreePath(row)); auto j = child(globalBasis->localView().tree(), makeTreePath(col)); + auto op = makeGridOperator(preOp, globalBasis->gridView()); auto localAssembler = makeLocalAssemblerPtr(op, i, j); matrixOperators[i][j].element.push_back({localAssembler, factor, estFactor}); @@ -150,7 +152,7 @@ template template void ProblemStat::addMatrixOperator( BoundaryType b, - Operator const& op, + Operator const& preOp, RowTreePath row, ColTreePath col, double* factor, double* estFactor) { @@ -158,6 +160,7 @@ void ProblemStat::addMatrixOperator( auto j = child(globalBasis->localView().tree(), makeTreePath(col)); using Intersection = typename GridView::Intersection; + auto op = makeGridOperator(preOp, globalBasis->gridView()); auto localAssembler = makeLocalAssemblerPtr(op, i, j); matrixOperators[i][j].boundary.push_back({localAssembler, factor, estFactor, b}); @@ -168,12 +171,13 @@ void ProblemStat::addMatrixOperator( template template void ProblemStat::addVectorOperator( - Operator const& op, + Operator const& preOp, TreePath path, double* factor, double* estFactor) { auto i = child(globalBasis->localView().tree(), makeTreePath(path)); + auto op = makeGridOperator(preOp, globalBasis->gridView()); auto localAssembler = makeLocalAssemblerPtr(op, i); rhsOperators[i].element.push_back({localAssembler, factor, estFactor}); @@ -185,13 +189,14 @@ template template void ProblemStat::addVectorOperator( BoundaryType b, - Operator const& op, + Operator const& preOp, TreePath path, double* factor, double* estFactor) { auto i = child(globalBasis->localView().tree(), makeTreePath(path)); using Intersection = typename GridView::Intersection; + auto op = makeGridOperator(preOp, globalBasis->gridView()); auto localAssembler = makeLocalAssemblerPtr(op, i); rhsOperators[i].boundary.push_back({localAssembler, factor, estFactor, b}); @@ -211,33 +216,13 @@ addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Val auto i = child(globalBasis->localView().tree(), makeTreePath(row)); auto j = child(globalBasis->localView().tree(), makeTreePath(col)); - addDirichletBCImpl(predicate, i, j, values); -} - -template - template -void ProblemStat:: -addDirichletBCImpl(Predicate const& predicate, RowNode const& row, ColNode const& col, RangeType const& value) -{ - using Range = RangeType; - using BcType = DirichletBC; - auto bc = std::make_shared(predicate, [value](auto const&) { return value; }); - constraints[row][col].push_back(bc); -} - -template - template - std::enable_if_t>::value> -ProblemStat:: -addDirichletBCImpl(Predicate const& predicate, RowNode const& row, ColNode const& col, Values const& values) -{ - using Range = RangeType; - static_assert( Concepts::Functor, - "Function passed to addDirichletBC for `values` does not model the Functor concept"); + auto valueGridFct = makeGridFunction(values, globalBasis->gridView(), Dune::PriorityTag<42>{}); + using Range = typename decltype(valueGridFct)::Range; using BcType = DirichletBC; - auto bc = std::make_shared(predicate, values); - constraints[row][col].push_back(bc); + auto bc = std::make_shared(predicate, valueGridFct); + constraints[i][j].push_back(bc); + // TODO: make DirichletBC an abstract class and add specialization with gridfunction type } diff --git a/dune/amdis/Terms.hpp b/dune/amdis/Terms.hpp deleted file mode 100644 index 79a0266f77e5f01069b4ac151ef28acad890fc01..0000000000000000000000000000000000000000 --- a/dune/amdis/Terms.hpp +++ /dev/null @@ -1,6 +0,0 @@ -#pragma once - -#include -#include -#include -#include diff --git a/dune/amdis/assembler/CMakeLists.txt b/dune/amdis/assembler/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..fd3abb09dbf748ee06c05ce4b8bfaf1750ab173e --- /dev/null +++ b/dune/amdis/assembler/CMakeLists.txt @@ -0,0 +1,21 @@ +#install headers + +install(FILES + FirstOrderDivTestvecTrial.hpp + FirstOrderGradTestTrial.hpp + FirstOrderGradTestTrialvec.hpp + FirstOrderPartialTestTrial.hpp + FirstOrderTestDivTrialvec.hpp + FirstOrderTestGradTrial.hpp + FirstOrderTestPartialTrial.hpp + FirstOrderTestvecGradTrial.hpp + SecondOrderDivTestvecDivTrialvec.hpp + SecondOrderGradTestGradTrial.hpp + SecondOrderPartialTestPartialTrial.hpp + ZeroOrderTest.hpp + ZeroOrderTestTrial.hpp + ZeroOrderTestTrialvec.hpp + ZeroOrderTestvec.hpp + ZeroOrderTestvecTrial.hpp + ZeroOrderTestvecTrialvec.hpp +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis/operators) diff --git a/dune/amdis/operators/FirstOrderDivTestvecTrial.hpp b/dune/amdis/assembler/FirstOrderDivTestvecTrial.hpp similarity index 68% rename from dune/amdis/operators/FirstOrderDivTestvecTrial.hpp rename to dune/amdis/assembler/FirstOrderDivTestvecTrial.hpp index 7d276628bc4ea388f5604864270c4da15d80aacf..50e39886b519ad510e31bf68587f4f6c157784b9 100644 --- a/dune/amdis/operators/FirstOrderDivTestvecTrial.hpp +++ b/dune/amdis/assembler/FirstOrderDivTestvecTrial.hpp @@ -2,7 +2,7 @@ #include -#include +#include namespace AMDiS { @@ -13,14 +13,14 @@ namespace AMDiS // first-order operator - template - class ExpressionOperator - : public ExpressionOperator + template + class GridFunctionOperator + : public GridFunctionOperator { - using Transposed = ExpressionOperator; + using Transposed = GridFunctionOperator; public: - ExpressionOperator(tag::divtestvec_trial, Expr const& expr) + GridFunctionOperator(tag::divtestvec_trial, GridFct const& expr) : Transposed(tag::test_divtrialvec{}, expr) {} @@ -34,6 +34,9 @@ namespace AMDiS std::integral_constant flagSameFE, std::integral_constant flagSameNode) { + static_assert(RowNode::isPower && ColNode::isLeaf, + "ColNode must be Leaf-Node and RowNode must be a Power-Node."); + auto elementMatrixTransposed = trans(elementMatrix); Transposed::calculateElementMatrix( context, quad, elementMatrixTransposed, colNode, rowNode, flagSameFE, flagSameNode); diff --git a/dune/amdis/operators/FirstOrderGradTestTrial.hpp b/dune/amdis/assembler/FirstOrderGradTestTrial.hpp similarity index 75% rename from dune/amdis/operators/FirstOrderGradTestTrial.hpp rename to dune/amdis/assembler/FirstOrderGradTestTrial.hpp index 1d818c7ebeb1371386ed51a3667e8745736c5558..a7c5de9a39b030e8248ce2fe36e48a787be52b7a 100644 --- a/dune/amdis/operators/FirstOrderGradTestTrial.hpp +++ b/dune/amdis/assembler/FirstOrderGradTestTrial.hpp @@ -2,7 +2,7 @@ #include -#include +#include namespace AMDiS { @@ -13,14 +13,14 @@ namespace AMDiS // first-order operator - template - class ExpressionOperator - : public ExpressionOperator + template + class GridFunctionOperator + : public GridFunctionOperator { - using Transposed = ExpressionOperator; + using Transposed = GridFunctionOperator; public: - ExpressionOperator(tag::gradtest_trial, Expr const& expr) + GridFunctionOperator(tag::gradtest_trial, GridFct const& expr) : Transposed(tag::test_gradtrial{}, expr) {} diff --git a/dune/amdis/operators/FirstOrderGradTestTrialvec.hpp b/dune/amdis/assembler/FirstOrderGradTestTrialvec.hpp similarity index 74% rename from dune/amdis/operators/FirstOrderGradTestTrialvec.hpp rename to dune/amdis/assembler/FirstOrderGradTestTrialvec.hpp index 98896146710907e8309858fdb23371bec8e11251..dd36c9e2dbd50a024c0e271121351e8bc9999a09 100644 --- a/dune/amdis/operators/FirstOrderGradTestTrialvec.hpp +++ b/dune/amdis/assembler/FirstOrderGradTestTrialvec.hpp @@ -2,7 +2,7 @@ #include -#include +#include namespace AMDiS { @@ -13,14 +13,14 @@ namespace AMDiS // first-order operator - template - class ExpressionOperator - : public ExpressionOperator + template + class GridFunctionOperator + : public GridFunctionOperator { - using Transposed = ExpressionOperator; + using Transposed = GridFunctionOperator; public: - ExpressionOperator(tag::gradtest_trialvec, Expr const& expr) + GridFunctionOperator(tag::gradtest_trialvec, GridFct const& expr) : Transposed(tag::testvec_gradtrial{}, expr) {} diff --git a/dune/amdis/operators/FirstOrderPartialTestTrial.hpp b/dune/amdis/assembler/FirstOrderPartialTestTrial.hpp similarity index 74% rename from dune/amdis/operators/FirstOrderPartialTestTrial.hpp rename to dune/amdis/assembler/FirstOrderPartialTestTrial.hpp index db739f39b3d8a7e10545e5f5f0a98ac04ac5dd28..e505b88bdca0f6aed75f59f28156869caac1ff4d 100644 --- a/dune/amdis/operators/FirstOrderPartialTestTrial.hpp +++ b/dune/amdis/assembler/FirstOrderPartialTestTrial.hpp @@ -2,7 +2,7 @@ #include -#include +#include namespace AMDiS { @@ -16,14 +16,14 @@ namespace AMDiS // first-order operator - template - class ExpressionOperator - : public ExpressionOperator + template + class GridFunctionOperator + : public GridFunctionOperator { - using Transposed = ExpressionOperator; + using Transposed = GridFunctionOperator; public: - ExpressionOperator(tag::partialtest_trial tag, Expr const& expr) + GridFunctionOperator(tag::partialtest_trial tag, GridFct const& expr) : Transposed(tag::test_partialtrial{tag.comp}, expr) {} diff --git a/dune/amdis/operators/FirstOrderTestDivTrialvec.hpp b/dune/amdis/assembler/FirstOrderTestDivTrialvec.hpp similarity index 78% rename from dune/amdis/operators/FirstOrderTestDivTrialvec.hpp rename to dune/amdis/assembler/FirstOrderTestDivTrialvec.hpp index f745cd82a8513784f0ca68eb0a30d5a88ef62119..2446d3aac156211ebca8b26741b5bdd399da1a2a 100644 --- a/dune/amdis/operators/FirstOrderTestDivTrialvec.hpp +++ b/dune/amdis/assembler/FirstOrderTestDivTrialvec.hpp @@ -2,7 +2,7 @@ #include -#include +#include #include #include @@ -15,16 +15,16 @@ namespace AMDiS // first-order operator - template - class ExpressionOperator - : public ExpressionOperatorBase + template + class GridFunctionOperator + : public GridFunctionOperatorBase { - using Super = ExpressionOperatorBase; + using Super = GridFunctionOperatorBase; - static_assert( Category::Scalar, "Expression must be of scalar type." ); + static_assert( Category::Scalar, "Expression must be of scalar type." ); public: - ExpressionOperator(tag::test_divtrialvec, Expr const& expr) + GridFunctionOperator(tag::test_divtrialvec, GridFct const& expr) : Super(expr, 1) {} @@ -38,8 +38,8 @@ namespace AMDiS std::integral_constant, std::integral_constant) { - static_assert( std::is_same::value && - std::is_same::value, "" ); + static_assert(RowNode::isLeaf && ColNode::isPower, + "RowNode must be Leaf-Node and ColNode must be a Power-Node."); static const std::size_t CHILDREN = ColNode::CHILDREN; static_assert( CHILDREN == Context::dow, "" ); @@ -47,8 +47,6 @@ namespace AMDiS auto const& rowLocalFE = rowNode.finiteElement(); auto const& colLocalFE = colNode.child(0).finiteElement(); - Super::init(context, quad); - std::vector > referenceGradients; for (std::size_t iq = 0; iq < quad.size(); ++iq) { @@ -59,7 +57,7 @@ namespace AMDiS const auto jacobian = context.geometry.jacobianInverseTransposed(local); // The multiplicative factor in the integral transformation formula - const double factor = Super::eval(iq) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); + double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_); diff --git a/dune/amdis/operators/FirstOrderTestGradTrial.hpp b/dune/amdis/assembler/FirstOrderTestGradTrial.hpp similarity index 77% rename from dune/amdis/operators/FirstOrderTestGradTrial.hpp rename to dune/amdis/assembler/FirstOrderTestGradTrial.hpp index 5ed708210e3d126f7cd4ec820a2fdea3dcef4a8e..2a1abf6c83326198cf0dd5bcdf952da12524ef83 100644 --- a/dune/amdis/operators/FirstOrderTestGradTrial.hpp +++ b/dune/amdis/assembler/FirstOrderTestGradTrial.hpp @@ -2,7 +2,7 @@ #include -#include +#include #include #include @@ -15,16 +15,16 @@ namespace AMDiS // first-order operator - template - class ExpressionOperator - : public ExpressionOperatorBase + template + class GridFunctionOperator + : public GridFunctionOperatorBase { - using Super = ExpressionOperatorBase; + using Super = GridFunctionOperatorBase; - static_assert( Category::Vector, "Expression must be of vector type." ); + static_assert( Category::Vector, "Expression must be of vector type." ); public: - ExpressionOperator(tag::test_gradtrial, Expr const& expr) + GridFunctionOperator(tag::test_gradtrial, GridFct const& expr) : Super(expr, 1) {} @@ -38,14 +38,12 @@ namespace AMDiS std::integral_constant, std::integral_constant) { - static_assert( std::is_same::value && - std::is_same::value, "" ); + static_assert(RowNode::isLeaf && ColNode::isLeaf, + "Operator can be applied to Leaf-Nodes only."); auto const& rowLocalFE = rowNode.finiteElement(); auto const& colLocalFE = colNode.finiteElement(); - Super::init(context, quad); - std::vector > referenceGradients; for (std::size_t iq = 0; iq < quad.size(); ++iq) { @@ -56,8 +54,8 @@ namespace AMDiS const auto jacobian = context.geometry.jacobianInverseTransposed(local); // The multiplicative factor in the integral transformation formula - const double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight(); - const auto b = Super::eval(iq); + double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight(); + const auto b = Super::coefficient(local); rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_); diff --git a/dune/amdis/operators/FirstOrderTestPartialTrial.hpp b/dune/amdis/assembler/FirstOrderTestPartialTrial.hpp similarity index 78% rename from dune/amdis/operators/FirstOrderTestPartialTrial.hpp rename to dune/amdis/assembler/FirstOrderTestPartialTrial.hpp index 4c36ce831d012a507fb027dd65f0cbbdadf88f39..8ff9a1125cc27209ddf966756935db2361d9f77c 100644 --- a/dune/amdis/operators/FirstOrderTestPartialTrial.hpp +++ b/dune/amdis/assembler/FirstOrderTestPartialTrial.hpp @@ -2,7 +2,7 @@ #include -#include +#include #include #include @@ -18,16 +18,16 @@ namespace AMDiS // first-order operator - template - class ExpressionOperator - : public ExpressionOperatorBase + template + class GridFunctionOperator + : public GridFunctionOperatorBase { - using Super = ExpressionOperatorBase; + using Super = GridFunctionOperatorBase; - static_assert( Category::Scalar, "Expression must be of scalar type." ); + static_assert( Category::Scalar, "Expression must be of scalar type." ); public: - ExpressionOperator(tag::test_partialtrial tag, Expr const& expr) + GridFunctionOperator(tag::test_partialtrial tag, GridFct const& expr) : Super(expr, 1) , comp_(tag.comp) {} @@ -42,14 +42,12 @@ namespace AMDiS std::integral_constant, std::integral_constant) { - static_assert( std::is_same::value && - std::is_same::value, "" ); + static_assert(RowNode::isLeaf && ColNode::isLeaf, + "Operator can be applied to Leaf-Nodes only."); auto const& rowLocalFE = rowNode.finiteElement(); auto const& colLocalFE = colNode.finiteElement(); - Super::init(context, quad); - for (std::size_t iq = 0; iq < quad.size(); ++iq) { // Position of the current quadrature point in the reference element decltype(auto) local = context.position(quad[iq].position()); @@ -58,8 +56,8 @@ namespace AMDiS const auto jacobian = context.geometry.jacobianInverseTransposed(local); // The multiplicative factor in the integral transformation formula - const double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight(); - const double c = Super::eval(iq); + double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight(); + double c = Super::coefficient(local); rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_); diff --git a/dune/amdis/operators/FirstOrderTestvecGradTrial.hpp b/dune/amdis/assembler/FirstOrderTestvecGradTrial.hpp similarity index 78% rename from dune/amdis/operators/FirstOrderTestvecGradTrial.hpp rename to dune/amdis/assembler/FirstOrderTestvecGradTrial.hpp index c33b00773af43ceadedee421500fa4f0a8c11e84..fe924bdd4acbd3d00580aab1b1b5c6babe196db0 100644 --- a/dune/amdis/operators/FirstOrderTestvecGradTrial.hpp +++ b/dune/amdis/assembler/FirstOrderTestvecGradTrial.hpp @@ -2,7 +2,7 @@ #include -#include +#include #include #include @@ -15,16 +15,16 @@ namespace AMDiS // first-order operator - template - class ExpressionOperator - : public ExpressionOperatorBase + template + class GridFunctionOperator + : public GridFunctionOperatorBase { - using Super = ExpressionOperatorBase; + using Super = GridFunctionOperatorBase; - static_assert( Category::Scalar, "Expression must be of scalar type." ); + static_assert( Category::Scalar, "Expression must be of scalar type." ); public: - ExpressionOperator(tag::testvec_gradtrial, Expr const& expr) + GridFunctionOperator(tag::testvec_gradtrial, GridFct const& expr) : Super(expr, 1) {} @@ -38,8 +38,8 @@ namespace AMDiS std::integral_constant, std::integral_constant) { - static_assert( std::is_same::value && - std::is_same::value, "" ); + static_assert(RowNode::isPower && ColNode::isLeaf, + "ColNode must be Leaf-Node and RowNode must be a Power-Node."); static const std::size_t CHILDREN = RowNode::CHILDREN; static_assert( CHILDREN == Context::dow, "" ); @@ -47,8 +47,6 @@ namespace AMDiS auto const& rowLocalFE = rowNode.child(0).finiteElement(); auto const& colLocalFE = colNode.finiteElement(); - Super::init(context, quad); - std::vector > referenceGradients; for (std::size_t iq = 0; iq < quad.size(); ++iq) { @@ -59,7 +57,7 @@ namespace AMDiS const auto jacobian = context.geometry.jacobianInverseTransposed(local); // The multiplicative factor in the integral transformation formula - const double factor = Super::eval(iq) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); + double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_); diff --git a/dune/amdis/operators/SecondOrderDivTestvecDivTrialvec.hpp b/dune/amdis/assembler/SecondOrderDivTestvecDivTrialvec.hpp similarity index 85% rename from dune/amdis/operators/SecondOrderDivTestvecDivTrialvec.hpp rename to dune/amdis/assembler/SecondOrderDivTestvecDivTrialvec.hpp index 556bf7fbd5c4178a368cb93553fbb8644d63b0d6..64956b5b8e907a3ae9b4c3d243b264af09e69db4 100644 --- a/dune/amdis/operators/SecondOrderDivTestvecDivTrialvec.hpp +++ b/dune/amdis/assembler/SecondOrderDivTestvecDivTrialvec.hpp @@ -2,7 +2,7 @@ #include -#include +#include #include #include @@ -15,16 +15,16 @@ namespace AMDiS // second-order operator - template - class ExpressionOperator - : public ExpressionOperatorBase + template + class GridFunctionOperator + : public GridFunctionOperatorBase { - using Super = ExpressionOperatorBase; + using Super = GridFunctionOperatorBase; - static_assert( Category::Scalar, "Expression must be of scalar type." ); + static_assert( Category::Scalar, "Expression must be of scalar type." ); public: - ExpressionOperator(tag::divtestvec_divtrialvec, Expr const& expr) + GridFunctionOperator(tag::divtestvec_divtrialvec, GridFct const& expr) : Super(expr, 2) {} @@ -38,8 +38,8 @@ namespace AMDiS std::integral_constant, std::integral_constant) { - static_assert( std::is_same::value && - std::is_same::value, "" ); + static_assert(RowNode::isPower && ColNode::isPower, + "Operator can be applied to Power-Nodes only."); static const std::size_t CHILDREN = RowNode::CHILDREN; static_assert( RowNode::CHILDREN == ColNode::CHILDREN, "" ); @@ -47,8 +47,6 @@ namespace AMDiS auto const& rowLocalFE = rowNode.child(0).finiteElement(); auto const& colLocalFE = colNode.child(0).finiteElement(); - Super::init(context, quad); - std::vector > rowReferenceGradients; std::vector > colReferenceGradients; @@ -60,7 +58,7 @@ namespace AMDiS const auto jacobian = context.geometry.jacobianInverseTransposed(local); // The multiplicative factor in the integral transformation formula - const double factor = Super::eval(iq) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); + double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); // The gradients of the shape functions on the reference element rowLocalFE.localBasis().evaluateJacobian(local, rowReferenceGradients); @@ -106,8 +104,6 @@ namespace AMDiS auto const& localFE = node.child(0).finiteElement(); - Super::init(context, quad); - std::vector > referenceGradients; for (std::size_t iq = 0; iq < quad.size(); ++iq) { @@ -118,7 +114,7 @@ namespace AMDiS const auto jacobian = context.geometry.jacobianInverseTransposed(local); // The multiplicative factor in the integral transformation formula - const double factor = Super::eval(iq) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); + double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); // The gradients of the shape functions on the reference element localFE.localBasis().evaluateJacobian(local, referenceGradients); diff --git a/dune/amdis/operators/SecondOrderGradTestGradTrial.hpp b/dune/amdis/assembler/SecondOrderGradTestGradTrial.hpp similarity index 85% rename from dune/amdis/operators/SecondOrderGradTestGradTrial.hpp rename to dune/amdis/assembler/SecondOrderGradTestGradTrial.hpp index 14803b3cc9848067c6163418fc43626ba30758f0..dc3e77cb347b1f434331877b5f2e3029cc7200e2 100644 --- a/dune/amdis/operators/SecondOrderGradTestGradTrial.hpp +++ b/dune/amdis/assembler/SecondOrderGradTestGradTrial.hpp @@ -2,7 +2,7 @@ #include -#include +#include #include #include @@ -15,18 +15,18 @@ namespace AMDiS // second-order operator , or - template - class ExpressionOperator - : public ExpressionOperatorBase + template + class GridFunctionOperator + : public GridFunctionOperatorBase { - using Super = ExpressionOperatorBase; + using Super = GridFunctionOperatorBase; - using expr_value_type = typename Expr::value_type; + using expr_value_type = typename GridFct::Range; static_assert( Category::Scalar || Category::Matrix, "Expression must be of scalar or matrix type." ); public: - ExpressionOperator(tag::gradtest_gradtrial, Expr const& expr) + GridFunctionOperator(tag::gradtest_gradtrial, GridFct const& expr) : Super(expr, 2) {} @@ -40,16 +40,14 @@ namespace AMDiS std::integral_constant, std::false_type /*sameNode*/) { - static_assert( std::is_same::value && - std::is_same::value, "" ); + static_assert(RowNode::isLeaf && ColNode::isLeaf, + "Operator can be applied to Leaf-Nodes only."); auto const& rowLocalFE = rowNode.finiteElement(); auto const& colLocalFE = colNode.finiteElement(); test_exit( sameFE, "Not implemented: currently only the implementation for equal fespaces available" ); - Super::init(context, quad); - for (std::size_t iq = 0; iq < quad.size(); ++iq) { // Position of the current quadrature point in the reference element decltype(auto) local = context.position(quad[iq].position()); @@ -58,8 +56,8 @@ namespace AMDiS const auto jacobian = context.geometry.jacobianInverseTransposed(local); // The multiplicative factor in the integral transformation formula - const double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight(); - const auto exprValue = Super::eval(iq); + double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight(); + const auto exprValue = Super::coefficient(local); // The gradients of the shape functions on the reference element std::vector > referenceGradients; @@ -92,10 +90,10 @@ namespace AMDiS std::true_type /*sameFE*/, std::true_type /*sameNode*/) { - static_assert( std::is_same::value && - std::is_same::value, "" ); + static_assert(RowNode::isLeaf && ColNode::isLeaf, + "Operator can be applied to Leaf-Nodes only."); - using Category = ValueCategory_t; + using Category = ValueCategory_t; calcElementMatrix(context, quad, elementMatrix, node, Category{}); } @@ -112,8 +110,6 @@ namespace AMDiS auto const& localFE = node.finiteElement(); std::vector > referenceGradients; - Super::init(context, quad); - for (std::size_t iq = 0; iq < quad.size(); ++iq) { // Position of the current quadrature point in the reference element decltype(auto) local = context.position(quad[iq].position()); @@ -122,7 +118,7 @@ namespace AMDiS const auto jacobian = context.geometry.jacobianInverseTransposed(local); // The multiplicative factor in the integral transformation formula - const double factor = Super::eval(iq) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); + double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); // The gradients of the shape functions on the reference element localFE.localBasis().evaluateJacobian(local, referenceGradients); @@ -161,8 +157,6 @@ namespace AMDiS auto const& localFE = node.finiteElement(); std::vector > referenceGradients; - Super::init(context, quad); - for (std::size_t iq = 0; iq < quad.size(); ++iq) { // Position of the current quadrature point in the reference element decltype(auto) local = context.position(quad[iq].position()); @@ -171,8 +165,8 @@ namespace AMDiS const auto jacobian = context.geometry.jacobianInverseTransposed(local); // The multiplicative factor in the integral transformation formula - const double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight(); - const auto exprValue = Super::eval(iq); + double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight(); + const auto exprValue = Super::coefficient(local); // The gradients of the shape functions on the reference element localFE.localBasis().evaluateJacobian(local, referenceGradients); diff --git a/dune/amdis/operators/SecondOrderPartialTestPartialTrial.hpp b/dune/amdis/assembler/SecondOrderPartialTestPartialTrial.hpp similarity index 81% rename from dune/amdis/operators/SecondOrderPartialTestPartialTrial.hpp rename to dune/amdis/assembler/SecondOrderPartialTestPartialTrial.hpp index bc5ee90274675fdd70771981e32f8ec2bced8e5c..4ef310f5058cf74d6da4e89d83f6edaa960764ab 100644 --- a/dune/amdis/operators/SecondOrderPartialTestPartialTrial.hpp +++ b/dune/amdis/assembler/SecondOrderPartialTestPartialTrial.hpp @@ -2,7 +2,7 @@ #include -#include +#include #include #include @@ -19,16 +19,16 @@ namespace AMDiS // second-order operator - template - class ExpressionOperator - : public ExpressionOperatorBase + template + class GridFunctionOperator + : public GridFunctionOperatorBase { - using Super = ExpressionOperatorBase; + using Super = GridFunctionOperatorBase; - static_assert( Category::Scalar, "Expression must be of scalar type." ); + static_assert( Category::Scalar, "Expression must be of scalar type." ); public: - ExpressionOperator(tag::partialtest_partialtrial tag, Expr const& expr) + GridFunctionOperator(tag::partialtest_partialtrial tag, GridFct const& expr) : Super(expr, 2) , compTest_(tag.comp_test) , compTrial_(tag.comp_trial) @@ -44,14 +44,12 @@ namespace AMDiS std::integral_constant, std::integral_constant) { - static_assert( std::is_same::value && - std::is_same::value, "" ); + static_assert(RowNode::isLeaf && ColNode::isLeaf, + "Operator can be applied to Leaf-Nodes only."); auto const& rowLocalFE = rowNode.finiteElement(); auto const& colLocalFE = colNode.finiteElement(); - Super::init(context, quad); - std::vector > rowReferenceGradients; std::vector > colReferenceGradients; @@ -63,8 +61,8 @@ namespace AMDiS const auto jacobian = context.geometry.jacobianInverseTransposed(local); // The multiplicative factor in the integral transformation formula - const double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight(); - const double c = Super::eval(iq); + double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight(); + double c = Super::coefficient(local); // The gradients of the shape functions on the reference element rowLocalFE.localBasis().evaluateJacobian(local, rowReferenceGradients); diff --git a/dune/amdis/operators/ZeroOrderTest.hpp b/dune/amdis/assembler/ZeroOrderTest.hpp similarity index 63% rename from dune/amdis/operators/ZeroOrderTest.hpp rename to dune/amdis/assembler/ZeroOrderTest.hpp index 2dd8dfe7347a913f5ce61425ccc58bf3740784c7..00822bf2ff76f99818959869d0920630f4726f9e 100644 --- a/dune/amdis/operators/ZeroOrderTest.hpp +++ b/dune/amdis/assembler/ZeroOrderTest.hpp @@ -2,7 +2,8 @@ #include -#include +#include +#include namespace AMDiS { @@ -13,16 +14,16 @@ namespace AMDiS // zero-order vector-operator - template - class ExpressionOperator - : public ExpressionOperatorBase + template + class GridFunctionOperator + : public GridFunctionOperatorBase { - using Super = ExpressionOperatorBase; + using Super = GridFunctionOperatorBase; - static_assert( Category::Scalar, "Expression must be of scalar type." ); + static_assert( Category::Scalar, "Expression must be of scalar type." ); public: - ExpressionOperator(tag::test, Expr const& expr) + GridFunctionOperator(tag::test, GridFct const& expr) : Super(expr, 0) {} @@ -33,19 +34,16 @@ namespace AMDiS ElementVector& elementVector, Node const& node) { - static_assert( std::is_same::value, "" ); - static_assert( std::is_same::value, "" ); + static_assert(Node::isLeaf, "Operator can be applied to Leaf-Nodes only"); auto const& localFE = node.finiteElement(); - Super::init(context, quad); - for (std::size_t iq = 0; iq < quad.size(); ++iq) { // Position of the current quadrature point in the reference element decltype(auto) local = context.position(quad[iq].position()); // The multiplicative factor in the integral transformation formula - const double factor = Super::eval(iq) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); + double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); localFE.localBasis().evaluateFunction(local, shapeValues_); diff --git a/dune/amdis/operators/ZeroOrderTestTrial.hpp b/dune/amdis/assembler/ZeroOrderTestTrial.hpp similarity index 73% rename from dune/amdis/operators/ZeroOrderTestTrial.hpp rename to dune/amdis/assembler/ZeroOrderTestTrial.hpp index 93372fc2ab49c0c8e5b0ec64e94709de9a32e055..5eb16a96dac5ff2854f2ddc8fb970e654f3645da 100644 --- a/dune/amdis/operators/ZeroOrderTestTrial.hpp +++ b/dune/amdis/assembler/ZeroOrderTestTrial.hpp @@ -2,7 +2,8 @@ #include -#include +#include +#include namespace AMDiS { @@ -13,16 +14,16 @@ namespace AMDiS // zero-order operator - template - class ExpressionOperator - : public ExpressionOperatorBase + template + class GridFunctionOperator + : public GridFunctionOperatorBase { - using Super = ExpressionOperatorBase; + using Super = GridFunctionOperatorBase; - static_assert( Category::Scalar, "Expression must be of scalar type." ); + static_assert( Category::Scalar, "Expression must be of scalar type." ); public: - ExpressionOperator(tag::test_trial, Expr const& expr) + GridFunctionOperator(tag::test_trial, GridFct const& expr) : Super(expr, 0) {} @@ -36,10 +37,8 @@ namespace AMDiS std::integral_constant, std::false_type /*sameNode*/) { - static_assert( std::is_same::value && - std::is_same::value, "" ); - - static_assert( std::is_same::value, "" ); + static_assert(RowNode::isLeaf && ColNode::isLeaf, + "Operator can be applied to Leaf-Nodes only."); auto const& rowLocalFE = rowNode.finiteElement(); auto const& colLocalFE = colNode.finiteElement(); @@ -47,14 +46,12 @@ namespace AMDiS auto* rowShapeValues = &rowShapeValues_; auto* colShapeValues = (sameFE ? &rowShapeValues_ : &colShapeValues_); - Super::init(context, quad); - for (std::size_t iq = 0; iq < quad.size(); ++iq) { // Position of the current quadrature point in the reference element decltype(auto) local = context.position(quad[iq].position()); // The multiplicative factor in the integral transformation formula - const double factor = Super::eval(iq) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); + double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_); if (!sameFE) @@ -83,22 +80,18 @@ namespace AMDiS std::true_type /*sameFE*/, std::true_type /*sameNode*/) { - static_assert( std::is_same::value && - std::is_same::value, "" ); - - static_assert( std::is_same::value, "" ); + static_assert(RowNode::isLeaf && ColNode::isLeaf, + "Operator can be applied to Leaf-Nodes only."); auto const& localFE = node.finiteElement(); auto& shapeValues = rowShapeValues_; - Super::init(context, quad); - for (std::size_t iq = 0; iq < quad.size(); ++iq) { // Position of the current quadrature point in the reference element decltype(auto) local = context.position(quad[iq].position()); // The multiplicative factor in the integral transformation formula - const double factor = Super::eval(iq) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); + const double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); localFE.localBasis().evaluateFunction(local, shapeValues); diff --git a/dune/amdis/operators/ZeroOrderTestTrialvec.hpp b/dune/amdis/assembler/ZeroOrderTestTrialvec.hpp similarity index 71% rename from dune/amdis/operators/ZeroOrderTestTrialvec.hpp rename to dune/amdis/assembler/ZeroOrderTestTrialvec.hpp index 401ac70c70e70efb591730c1b10505dbbd1a6774..2268afaf6d7f6d82b868ad66bf554c966b193939 100644 --- a/dune/amdis/operators/ZeroOrderTestTrialvec.hpp +++ b/dune/amdis/assembler/ZeroOrderTestTrialvec.hpp @@ -2,7 +2,8 @@ #include -#include +#include +#include namespace AMDiS { @@ -13,16 +14,16 @@ namespace AMDiS // zero-order operator - template - class ExpressionOperator - : public ExpressionOperatorBase + template + class GridFunctionOperator + : public GridFunctionOperatorBase { - using Super = ExpressionOperatorBase; + using Super = GridFunctionOperatorBase; - static_assert( Category::Vector, "Expression must be of vector type." ); + static_assert( Category::Vector, "Expression must be of vector type." ); public: - ExpressionOperator(tag::test_trialvec, Expr const& expr) + GridFunctionOperator(tag::test_trialvec, GridFct const& expr) : Super(expr, 0) {} @@ -36,11 +37,11 @@ namespace AMDiS std::integral_constant, std::integral_constant) { - static_assert( std::is_same::value && - std::is_same::value, "" ); + static_assert(RowNode::isLeaf && ColNode::isPower, + "RowNode must be Leaf-Node and ColNode must be a Power-Node."); static const std::size_t CHILDREN = ColNode::CHILDREN; - static_assert( std::is_same>::value, "" ); + static_assert( std::is_same>::value, "" ); auto const& rowLocalFE = rowNode.finiteElement(); auto const& colLocalFE = colNode.child(0).finiteElement(); @@ -48,15 +49,13 @@ namespace AMDiS auto* rowShapeValues = &rowShapeValues_; auto* colShapeValues = (sameFE ? &rowShapeValues_ : &colShapeValues_); - Super::init(context, quad); - for (std::size_t iq = 0; iq < quad.size(); ++iq) { // Position of the current quadrature point in the reference element decltype(auto) local = context.position(quad[iq].position()); // The multiplicative factor in the integral transformation formula - const double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight(); - const auto b = Super::eval(iq); + double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight(); + const auto b = Super::coefficient(local); rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_); if (!sameFE) diff --git a/dune/amdis/operators/ZeroOrderTestvec.hpp b/dune/amdis/assembler/ZeroOrderTestvec.hpp similarity index 65% rename from dune/amdis/operators/ZeroOrderTestvec.hpp rename to dune/amdis/assembler/ZeroOrderTestvec.hpp index e42948d8e76f681abd6ce8465ef40e71570c0b80..e7e5c3fec20786c86b8f75f7a7dddd55da75a399 100644 --- a/dune/amdis/operators/ZeroOrderTestvec.hpp +++ b/dune/amdis/assembler/ZeroOrderTestvec.hpp @@ -2,7 +2,8 @@ #include -#include +#include +#include namespace AMDiS { @@ -13,16 +14,16 @@ namespace AMDiS // zero-order vector-operator - template - class ExpressionOperator - : public ExpressionOperatorBase + template + class GridFunctionOperator + : public GridFunctionOperatorBase { - using Super = ExpressionOperatorBase; + using Super = GridFunctionOperatorBase; - static_assert( Category::Vector, "Expression must be of vector type." ); + static_assert( Category::Vector, "Expression must be of vector type." ); public: - ExpressionOperator(tag::testvec, Expr const& expr) + GridFunctionOperator(tag::testvec, GridFct const& expr) : Super(expr, 0) {} @@ -33,22 +34,20 @@ namespace AMDiS ElementVector& elementVector, Node const& node) { - static_assert( std::is_same::value, "" ); + static_assert(Node::isPower, + "Operator can be applied to Power-Nodes only."); static const std::size_t CHILDREN = Node::CHILDREN; - static_assert( std::is_same>::value, "" ); auto const& localFE = node.child(0).finiteElement(); - Super::init(context, quad); - for (std::size_t iq = 0; iq < quad.size(); ++iq) { // Position of the current quadrature point in the reference element decltype(auto) local = context.position(quad[iq].position()); // The multiplicative factor in the integral transformation formula - const double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight(); - const auto exprValue = Super::eval(iq); + double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight(); + const auto exprValue = Super::coefficient(local); localFE.localBasis().evaluateFunction(local, shapeValues_); diff --git a/dune/amdis/operators/ZeroOrderTestvecTrial.hpp b/dune/amdis/assembler/ZeroOrderTestvecTrial.hpp similarity index 75% rename from dune/amdis/operators/ZeroOrderTestvecTrial.hpp rename to dune/amdis/assembler/ZeroOrderTestvecTrial.hpp index 5b851e42a418563bd03e1754ebe59b3599e7e6e1..b356ae25a9569d86814b6c12b7df38daaa8ae397 100644 --- a/dune/amdis/operators/ZeroOrderTestvecTrial.hpp +++ b/dune/amdis/assembler/ZeroOrderTestvecTrial.hpp @@ -2,7 +2,7 @@ #include -#include +#include namespace AMDiS { @@ -13,14 +13,14 @@ namespace AMDiS // zero-order operator - template - class ExpressionOperator - : public ExpressionOperator + template + class GridFunctionOperator + : public GridFunctionOperator { - using Transposed = ExpressionOperator; + using Transposed = GridFunctionOperator; public: - ExpressionOperator(tag::testvec_trial, Expr const& expr) + GridFunctionOperator(tag::testvec_trial, GridFct const& expr) : Transposed(tag::test_trialvec{}, expr) {} diff --git a/dune/amdis/operators/ZeroOrderTestvecTrialvec.hpp b/dune/amdis/assembler/ZeroOrderTestvecTrialvec.hpp similarity index 86% rename from dune/amdis/operators/ZeroOrderTestvecTrialvec.hpp rename to dune/amdis/assembler/ZeroOrderTestvecTrialvec.hpp index 5d83e96f89f09c8f1b3b27b635994c032fdffd5a..94319155bc4febbcf1c70ab72a0f881926b9c4d3 100644 --- a/dune/amdis/operators/ZeroOrderTestvecTrialvec.hpp +++ b/dune/amdis/assembler/ZeroOrderTestvecTrialvec.hpp @@ -2,7 +2,7 @@ #include -#include +#include #include namespace AMDiS @@ -14,18 +14,18 @@ namespace AMDiS // zero-order operator , or - template - class ExpressionOperator - : public ExpressionOperatorBase + template + class GridFunctionOperator + : public GridFunctionOperatorBase { - using Super = ExpressionOperatorBase; + using Super = GridFunctionOperatorBase; - using expr_value_type = typename Expr::value_type; + using expr_value_type = typename GridFct::Range; static_assert( Category::Scalar || Category::Matrix, "Expression must be of scalar or matrix type." ); public: - ExpressionOperator(tag::testvec_trialvec, Expr const& expr) + GridFunctionOperator(tag::testvec_trialvec, GridFct const& expr) : Super(expr, 0) {} @@ -39,13 +39,13 @@ namespace AMDiS std::integral_constant, std::integral_constant) { - static_assert( std::is_same::value && - std::is_same::value, "" ); + static_assert(RowNode::isPower && ColNode::isPower, + "Operator can be applied to Power-Nodes only."); static_assert( (RowNode::CHILDREN == ColNode::CHILDREN), "" ); // theoretically the restriction of the same nr of childs would not be necessary - using Category = ValueCategory_t; + using Category = ValueCategory_t; static_assert( std::is_same::value || std::is_same::value, "" ); calcElementMatrix(context, quad, elementMatrix, rowNode, colNode, @@ -75,14 +75,12 @@ namespace AMDiS auto* rowShapeValues = &rowShapeValues_; auto* colShapeValues = (sameFE ? &rowShapeValues_ : &colShapeValues_); - Super::init(context, quad); - for (std::size_t iq = 0; iq < quad.size(); ++iq) { // Position of the current quadrature point in the reference element decltype(auto) local = context.position(quad[iq].position()); // The multiplicative factor in the integral transformation formula - const double factor = Super::eval(iq) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); + double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_); if (!sameFE) @@ -123,14 +121,12 @@ namespace AMDiS auto& shapeValues = rowShapeValues_; - Super::init(context, quad); - for (std::size_t iq = 0; iq < quad.size(); ++iq) { // Position of the current quadrature point in the reference element decltype(auto) local = context.position(quad[iq].position()); // The multiplicative factor in the integral transformation formula - const double factor = Super::eval(iq) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); + double factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); localFE.localBasis().evaluateFunction(local, shapeValues); @@ -169,7 +165,7 @@ namespace AMDiS tag::matrix) { static const std::size_t CHILDREN = RowNode::CHILDREN; - static_assert( std::is_same>::value, "" ); + static_assert( std::is_same>::value, "" ); auto const& rowLocalFE = rowNode.child(0).finiteElement(); auto const& colLocalFE = colNode.child(0).finiteElement(); @@ -177,15 +173,13 @@ namespace AMDiS auto* rowShapeValues = &rowShapeValues_; auto* colShapeValues = (sameFE ? &rowShapeValues_ : &colShapeValues_); - Super::init(context, quad); - for (std::size_t iq = 0; iq < quad.size(); ++iq) { // Position of the current quadrature point in the reference element decltype(auto) local = context.position(quad[iq].position()); // The multiplicative factor in the integral transformation formula - const double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight(); - const Dune::FieldMatrix exprValue = Super::eval(iq); + double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight(); + const Dune::FieldMatrix exprValue = Super::coefficient(local); rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_); if (!sameFE) diff --git a/dune/amdis/common/Assigner.hpp b/dune/amdis/common/Assigner.hpp new file mode 100644 index 0000000000000000000000000000000000000000..920975bfde41869eeb5de5b4b79b2ed92930a49a --- /dev/null +++ b/dune/amdis/common/Assigner.hpp @@ -0,0 +1,53 @@ +#pragma once + +namespace AMDiS +{ + namespace Assigner + { + struct assign + { + template + constexpr void operator()(T const& from, S& to) const + { + to = from; + } + }; + + struct plus_assign + { + template + constexpr void operator()(T const& from, S& to) const + { + to += from; + } + }; + + struct minus_assign + { + template + constexpr void operator()(T const& from, S& to) const + { + to -= from; + } + }; + + struct multiplies_assign + { + template + constexpr void operator()(T const& from, S& to) const + { + to *= from; + } + }; + + struct divides_assign + { + template + constexpr void operator()(T const& from, S& to) const + { + to /= from; + } + }; + + } // end namespcae Assigner +} // end namespace AMDiS diff --git a/dune/amdis/common/CMakeLists.txt b/dune/amdis/common/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..fb11b09588470fc0647c891fc4345d008c918cee --- /dev/null +++ b/dune/amdis/common/CMakeLists.txt @@ -0,0 +1,21 @@ +#install headers + +install(FILES + ClonablePtr.hpp + Concepts.hpp + ConceptsBase.hpp + FieldMatVec.hpp + IndexSeq.hpp + Literals.hpp + Loops.hpp + Math.hpp + Mpl.hpp + Operations.hpp + OptionalDelete.hpp + ScalarTypes.hpp + Timer.hpp + TupleUtility.hpp + TypeDefs.hpp + Utility.hpp + ValueCategory.hpp +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis/common) diff --git a/dune/amdis/common/Concepts.hpp b/dune/amdis/common/Concepts.hpp index 9cc20b8f448e30352bc43cde134384c2fa75c63a..3c10a0e3c96359573a690ca9f35cd4bdf6f64770 100644 --- a/dune/amdis/common/Concepts.hpp +++ b/dune/amdis/common/Concepts.hpp @@ -97,6 +97,16 @@ namespace AMDiS auto requires_(F&& f, Args&&... args) -> decltype( f(std::forward(args)...)); }; + // f(args...) + struct Evaluable + { + template + auto requires_(F&& f, Args&&... args) -> Void_t< + decltype( f(std::forward(args)...)), + std::result_of_t(std::decay_t...)> + >; + }; + } // end namespace Definition #endif // DOXYGEN @@ -138,6 +148,10 @@ namespace AMDiS template constexpr bool Callable = models; + template + constexpr bool Evaluable = models; + + /// \brief A Functor is a function `F` with signature `Signature`. /** * To be used as follows: `Concepts::Functor`. Returns true, if diff --git a/dune/amdis/common/FieldMatVec.hpp b/dune/amdis/common/FieldMatVec.hpp index ab67359aca2f74aa2061407b813ce92eb535a767..6469012201898d30b5bc9836a1dba443e39f6e6c 100644 --- a/dune/amdis/common/FieldMatVec.hpp +++ b/dune/amdis/common/FieldMatVec.hpp @@ -5,10 +5,11 @@ #include #include -#include +#include #include -#include #include +#include +#include namespace AMDiS { @@ -99,7 +100,7 @@ namespace AMDiS template T sum(Dune::FieldVector const& x) { - return Impl::accumulate(x, Operation::plus{}); + return Impl::accumulate(x, Operation::Plus{}); } @@ -115,28 +116,28 @@ namespace AMDiS template auto max(Dune::FieldVector const& x) { - return Impl::accumulate(x, Operation::maximum{}); + return Impl::accumulate(x, Operation::Max{}); } /// Minimum over all vector entries template auto min(Dune::FieldVector const& x) { - return Impl::accumulate(x, Operation::minimum{}); + return Impl::accumulate(x, Operation::Min{}); } /// Maximum of the absolute values of vector entries template auto abs_max(Dune::FieldVector const& x) { - return Impl::accumulate(x, Operation::abs_max{}); + return Impl::accumulate(x, Operation::AbsMax{}); } /// Minimum of the absolute values of vector entries template auto abs_min(Dune::FieldVector const& x) { - return Impl::accumulate(x, Operation::abs_min{}); + return Impl::accumulate(x, Operation::AbsMin{}); } // ---------------------------------------------------------------------------- diff --git a/dune/amdis/Math.hpp b/dune/amdis/common/Math.hpp similarity index 71% rename from dune/amdis/Math.hpp rename to dune/amdis/common/Math.hpp index 71d7585dd45d94a3496d3f149e230c71fa20c692..5293b0821c6f2e306c3979a5431a1b5ef3515655 100644 --- a/dune/amdis/Math.hpp +++ b/dune/amdis/common/Math.hpp @@ -73,22 +73,71 @@ namespace AMDiS /// Implementation of the minimum of two values \f$ min(a,b)\f$ of any type - /// supporting the `>` relation. + /// supporting the `<` relation. + /// @{ + template constexpr auto min(T0 a, T1 b) { - return a > b ? b : a; + return a < b ? a : b; + } + +#ifndef AMDIS_HAS_CXX_FOLD_EXPRESSIONS + template + constexpr T0 min(T0 a) + { + return a; } + template + constexpr auto min(T0 a, Ts... ts) + { + return min(a, min(ts...)); + } +#else + template + constexpr auto min(T0 a, Ts... ts) + { + std::common_type_t result = a; + return (result = min(result, ts), ...); + } +#endif + + /// @} + /// Implementation of the maximum of two values \f$ max(a,b)\f$ of any type /// supporting the `>` relation. + /// @{ template constexpr auto max(T0 a, T1 b) { - return a > b ? a : b; + return a < b ? b : a; + } + +#ifndef AMDIS_HAS_CXX_FOLD_EXPRESSIONS + template + constexpr T0 max(T0 a) + { + return a; } + template + constexpr auto max(T0 a, Ts... ts) + { + return max(a, max(ts...)); + } +#else + template + constexpr auto max(T0 a, Ts... ts) + { + std::common_type_t result = a; + return (result = max(result, ts), ...); + } +#endif + + /// @} + } // end namespace Math diff --git a/dune/amdis/common/Mpl.hpp b/dune/amdis/common/Mpl.hpp index f73c27af1b31455a01dadbe25ed29fca62b5ea05..bb3cfc4213486dcf780610853bce6fa662c9075e 100644 --- a/dune/amdis/common/Mpl.hpp +++ b/dune/amdis/common/Mpl.hpp @@ -33,6 +33,9 @@ namespace AMDiS template struct Indices {}; + template + auto get(Indices) { return std::get(std::make_tuple(index_...)); } + /// Some predefined index constants static constexpr index_t<0> _0 = {}; static constexpr index_t<1> _1 = {}; diff --git a/dune/amdis/common/Operations.hpp b/dune/amdis/common/Operations.hpp deleted file mode 100644 index ac83871d4640f70bfb8b877d673f62eaf4718f6d..0000000000000000000000000000000000000000 --- a/dune/amdis/common/Operations.hpp +++ /dev/null @@ -1,287 +0,0 @@ -#pragma once - -#include - -#include -#include -#include - -namespace AMDiS -{ - namespace Operation - { - /** \defgroup operations Operations - * \brief Functors, representing unary/binary operations used in expressions. - * @{ - **/ - - struct assign - { - template - constexpr void operator()(T const& from, S& to) const - { - to = from; - } - }; - - /// Operation that represents A+B - struct plus - { - template - constexpr auto operator()(T const& lhs, S const& rhs) const - { - return lhs + rhs; - } - - friend int order(plus, int lhs, int rhs) - { - return std::max(lhs, rhs); - } - }; - - struct plus_assign - { - template - constexpr void operator()(T const& from, S& to) const - { - to += from; - } - }; - - /// Operation that represents A-B - struct minus - { - template - constexpr auto operator()(T const& lhs, S const& rhs) const - { - return lhs - rhs; - } - - friend int order(minus, int lhs, int rhs) - { - return std::max(lhs, rhs); - } - }; - - struct minus_assign - { - template - constexpr void operator()(T const& from, S& to) const - { - to -= from; - } - }; - - /// Operation that represents A*B - struct times - { - template - constexpr auto operator()(T const& lhs, S const& rhs) const - { - return lhs * rhs; - } - - friend int order(times, int lhs, int rhs) - { - return lhs + rhs; - } - }; - - struct times_assign - { - template - constexpr void operator()(T const& from, S& to) const - { - to *= from; - } - }; - - /// Operation that represents A/B - struct divides - { - template - constexpr auto operator()(T const& lhs, S const& rhs) const - { - return lhs / rhs; - } - - // NOTE: the order is given only approximately - friend int order(divides, int lhs, int rhs) - { - return lhs + rhs + 1; - } - }; - - struct divides_assign - { - template - constexpr void operator()(T const& from, S& to) const - { - to /= from; - } - }; - - /// Operation that represents A-B - struct negate - { - template - constexpr auto operator()(T const& x) const - { - return -x; - } - - friend int order(divides, int x) - { - return x; - } - }; - - /// Operation that represents max(A,B) - struct maximum - { - template - constexpr auto operator()(T const& lhs, S const& rhs) const - { - return Math::max(lhs, rhs); - } - - // NOTE: the order is given only approximately - friend int order(maximum, int lhs, int rhs) - { - return std::max(lhs, rhs); - } - }; - - /// Operation that represents min(A,B) - struct minimum - { - template - constexpr auto operator()(T const& lhs, S const& rhs) const - { - return Math::min(lhs, rhs); - } - - // NOTE: the order is given only approximately - friend int order(minimum, int lhs, int rhs) - { - return std::max(lhs, rhs); - } - }; - - /// Operation that represents max(|A|,|B|) - struct abs_max - { - template - constexpr auto operator()(T const& lhs, S const& rhs) const - { - return Math::max(std::abs(lhs), std::abs(rhs)); - } - - // NOTE: the order is given only approximately - friend int order(abs_max, int lhs, int rhs) - { - return std::max(lhs, rhs); - } - }; - - /// Operation that represents min(|A|,|B|) - struct abs_min - { - template - constexpr auto operator()(T const& lhs, S const& rhs) const - { - return Math::min(std::abs(lhs), std::abs(rhs)); - } - - // NOTE: the order is given only approximately - friend int order(abs_min, int lhs, int rhs) - { - return std::max(lhs, rhs); - } - }; - - /// Operation that represents A*factor - template - struct times_scalar; - - template - struct times_scalar> > - { - constexpr explicit times_scalar(S const& scalar) : scalar(scalar) {} - - template ) > - constexpr auto operator()(T const& lhs) const - { - return lhs * scalar; - } - - friend int order(times_scalar, int lhs) - { - return lhs; - } - - private: - S scalar; - }; - - /// Operation that represents A/factor - template - struct divides_scalar; - - template - struct divides_scalar> > - { - constexpr explicit divides_scalar(S const& scalar) : scalar(scalar) {} - - template ) > - constexpr auto operator()(T const& lhs) const - { - return lhs / scalar; - } - - friend int order(divides_scalar, int lhs) - { - return lhs; - } - - private: - S scalar; - }; - - - /// Operation that represents conj(A) or A, if A is complex or not, respectively - template - struct hermitian - { - constexpr T operator()(T const& x) const - { - return std::conj(x); - } - - friend int order(hermitian, int x) - { - return x; - } - }; - - template - struct hermitian> > - { - constexpr T const& operator()(T const& x) const - { - return x; - } - - friend int order(hermitian, int x) - { - return x; - } - }; - - /** @} **/ - - } // end namespace Operation - -} // end namespace AMDiS diff --git a/dune/amdis/common/TupleUtility.hpp b/dune/amdis/common/TupleUtility.hpp index 33b73105d59e8ee6b038a0eec4862563ded0aafd..b03213f890b7392d333598cff6fd7da1092dce39 100644 --- a/dune/amdis/common/TupleUtility.hpp +++ b/dune/amdis/common/TupleUtility.hpp @@ -210,4 +210,5 @@ namespace AMDiS MakeSeq_t>::value>{}); } + } // end namespace AMDiS diff --git a/dune/amdis/common/ValueCategory.hpp b/dune/amdis/common/ValueCategory.hpp index 15d83a1bb41b6a231b3855788d45d561daaf3cfb..a2d78a58124edf5fb33d10bf4585895ff263c29b 100644 --- a/dune/amdis/common/ValueCategory.hpp +++ b/dune/amdis/common/ValueCategory.hpp @@ -27,7 +27,7 @@ namespace AMDiS }; template - using ValueCategory_t = typename ValueCategory::type; + using ValueCategory_t = typename ValueCategory>::type; template struct ValueCategory::value >> diff --git a/dune/amdis/gridfunctions/AnalyticGridFunction.hpp b/dune/amdis/gridfunctions/AnalyticGridFunction.hpp index eb55d0c56dd4f7a99544092e87e1803e19bb2cc2..5eef8dfb92cf4705b59c423a72c50697d3c8368d 100644 --- a/dune/amdis/gridfunctions/AnalyticGridFunction.hpp +++ b/dune/amdis/gridfunctions/AnalyticGridFunction.hpp @@ -3,11 +3,13 @@ #include #include +#include #include +#include #include -#define AMDIS_CALLABLE_DEFAULT_ORDER 1 +#define AMDIS_CALLABLE_DEFAULT_ORDER 4 namespace AMDiS { @@ -16,79 +18,116 @@ namespace AMDiS * @{ **/ + template + class AnalyticLocalFunction; - /// \brief A Gridfunction that applies a functor to the evaluated Gridfunctions - template - class AnalyticGridFunction + template + class AnalyticLocalFunction { public: - using EntitySet = Dune::Functions::GridViewEntitySet; - using Element = typename EntitySet::Element; - using Geometry = typename Element::Geometry; + using Range = R; + using Domain = D; // LocalDomain - using Domain = typename Geometry::GlobalCoordinate; - using Range = typename std::result_of::type; + using Geometry = typename LocalContext::Geometry; public: + AnalyticLocalFunction(Function const& fct) + : fct_(fct) + {} + + void bind(LocalContext const& element) + { + geometry_.emplace(element.geometry()); + } - template - class LocalFunction + void unbind() { - public: - using Range = typename std::result_of::type; - using Domain = typename Geometry::LocalCoordinate; - - public: - LocalFunction(Fct const& fct) - : fct_(fct) - {} - - void bind(Element const& element) - { - geometry_.emplace(element.geometry()); - } - - void unbind() - { - geometry_.reset(); - } - - Range operator()(Domain const& local) const - { - return fct_(geometry_.value().global(local)); - } - - template )> - friend auto derivative(LocalFunction const& glf) - { - auto df = derivative(glf.fct_); - return LocalFunction{df}; - } - - private: - Fct fct_; - Dune::Std::optional geometry_; - }; + geometry_.reset(); + } + Range operator()(Domain const& local) const + { + return fct_(geometry_.value().global(local)); + } + + Function const& fct() const + { + return fct_; + } + + private: + Function fct_; + Dune::Std::optional geometry_; + }; + + + template + auto derivative(AnalyticLocalFunction const& lf) + { + static_assert(Concepts::HasPartial, + "No partial(_0,...) defined for Functor F of AnalyticLocalFunction."); + + auto df = partial(lf.fct(), index_<0>); + + using RawSignature = typename Dune::Functions::SignatureTraits::RawSignature; + using DerivativeSignature = typename Dune::Functions::DefaultDerivativeTraits::Range(D); + return AnalyticLocalFunction{df}; + } + + /// Return the polynomial order of the function f, if a free function + /// order(f) exists, otherwise return a default value. +#ifdef AMDIS_HAS_CXX_CONSTEXPR_IF + template + int order(AnalyticLocalFunction const& lf) + { + if constexpr (Concepts::HasOrder) + return order(lf.fct(),1); + else + return AMDIS_CALLABLE_DEFAULT_ORDER; + } +#else + template + int order(AnalyticLocalFunction const& lf) + { + return Dune::Hybrid::ifElse( bool_>, + [&lf](auto id) { return id( order(lf.fct(),1) ); }, + [] (auto id) { return AMDIS_CALLABLE_DEFAULT_ORDER; } + ); + } +#endif + + + /// A Gridfunction that evaluates a function with global coordinates + template + class AnalyticGridFunction + { + public: + using EntitySet = Dune::Functions::GridViewEntitySet; + using Element = typename EntitySet::Element; + + using Domain = typename EntitySet::GlobalCoordinate; + using LocalDomain = typename EntitySet::LocalCoordinate; + using Range = std::decay_t>; public: + using LocalFunction = AnalyticLocalFunction; + public: + /// Constructor. Stores the function `fct` and creates an `EntitySet`. AnalyticGridFunction(Function const& fct, GridView const& gridView) : fct_(fct) , entitySet_(gridView) {} + /// Return the evaluated functor at global coordinates Range operator()(Domain const& x) const { return fct_(x); } - template - friend auto localFunction(AnalyticGridFunction const& gf) - -> LocalFunction + friend auto localFunction(AnalyticGridFunction const& gf) { - return LocalFunction{gf.fct_}; + return LocalFunction{gf.fct_}; } EntitySet const& entitySet() const @@ -105,46 +144,44 @@ namespace AMDiS }; + /// Return a GridFunction representing the derivative of a functor. + /// [expects: Functor f has free function derivative(f)] template - auto derivative(AnalyticGridFunction const& gf, - REQUIRES(Concepts::HasDerivative)) + auto derivative(AnalyticGridFunction const& gf) { - auto df = derivative(gf.fct()); - return AnalyticGridFunction{df, gf.entitySet().gridView()}; - } + static_assert(Concepts::HasPartial, + "No partial(_0,...) defined for Functor of AnalyticLocalFunction."); - template - int order(AnalyticGridFunction const& gf) - { - return Dune::Hybrid::ifElse( bool_>, - [&gf](auto id) { return id( order( gf.fct()) ); }, - [] (auto id) { return AMDIS_CALLABLE_DEFAULT_ORDER; } - ); + auto df = partial(gf.fct(), index_<0>); + return AnalyticGridFunction{df, gf.entitySet().gridView()}; } + /// A pre-GridFunction that just stores the function template struct AnalyticPreGridFunction { Function fct_; }; - namespace Concepts { namespace Definition + namespace Traits { template struct IsPreGridFunction> : std::true_type {}; - }} + } + + /// Generator function for \ref AnalyticPreGridFunction template )> - auto eval_(Function const& f) + auto evalAtQP(Function const& f) { return AnalyticPreGridFunction{f}; } - /// Generator function for \ref FunctorGridFunction expressions + /// Generator function for \ref AnalyticGridFunction expressions from callable functions. template )> auto makeGridFunction(Function const& f, GridView const& gridView, Dune::PriorityTag<3>) @@ -152,6 +189,7 @@ namespace AMDiS return AnalyticGridFunction{f, gridView}; } + /// Generator function for \ref AnalyticGridFunction expressions from \ref AnalyticPreGridFunction. template auto makeGridFunction(AnalyticPreGridFunction const& pre, GridView const& gridView, Dune::PriorityTag<3>) { diff --git a/dune/amdis/gridfunctions/CMakeLists.txt b/dune/amdis/gridfunctions/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..51955a7961165e0f35f8a89f70326aa83c3b2a59 --- /dev/null +++ b/dune/amdis/gridfunctions/CMakeLists.txt @@ -0,0 +1,14 @@ +#install headers + +install(FILES + AnalyticGridFunction.hpp + ConstantGridFunction.hpp + CoordsGridFunction.hpp + DerivativeGridFunction.hpp + DOFVectorView.hpp + DOFVectorView.inc.hpp + FunctorGridFunction.hpp + GridFunctionConcepts.hpp + Integrate.hpp + OperationsGridFunction.hpp +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis/gridfunctions) diff --git a/dune/amdis/gridfunctions/ConstantGridFunction.hpp b/dune/amdis/gridfunctions/ConstantGridFunction.hpp index dfced112e27270e0f5a4e6538891a6bb5e8973ce..3a32439a2ee9ec8e255bd1f7a1c5ff31e3f9e9f7 100644 --- a/dune/amdis/gridfunctions/ConstantGridFunction.hpp +++ b/dune/amdis/gridfunctions/ConstantGridFunction.hpp @@ -6,7 +6,6 @@ #include #include #include -#include #include #include @@ -20,41 +19,105 @@ namespace AMDiS * @{ **/ - template - struct ConstantFunction + /// \brief Gridfunction returning a constant value. + /** + * A stored constant is return in global and local evaluation of this functor. + * May be used with arithmetic types and vectors/matrices of arithmetic types. + * It is also allowed to pass a \ref std::reference_wrapper to allow to + * modify the value after construction. + **/ + template + class ConstantLocalFunction; + + template + class ConstantLocalFunction { + public: + using Range = R; + using Domain = D; // LocalDomain + + using Geometry = typename LocalContext::Geometry; + + public: + ConstantLocalFunction(T const& value) + : value_(value) + {} + + void bind(LocalContext const& /*element*/) {} + + void unbind() {} + + Range const& operator()(Domain const& /*local*/) const + { + return value_; + } + + friend int order(ConstantLocalFunction const& /*lf*/) + { + return 0; + } + + private: + T value_; + }; + + + template + auto derivative(ConstantLocalFunction const& lf) + { + using RawSignature = typename Dune::Functions::SignatureTraits::RawSignature; + using DerivativeRange = typename Dune::Functions::DefaultDerivativeTraits::Range; + DerivativeRange diff(0); + return ConstantLocalFunction{diff}; + } + + + /// A Gridfunction that evaluates a function with global coordinates + template + class ConstantGridFunction + { + public: + using EntitySet = Dune::Functions::GridViewEntitySet; + using Element = typename EntitySet::Element; + + using Domain = typename EntitySet::GlobalCoordinate; + using LocalDomain = typename EntitySet::LocalCoordinate; using Range = Underlying_t; - using DerivativeRange = typename Dune::Functions::DefaultDerivativeTraits::Range; - using Derivative = ConstantFunction; + public: + using LocalFunction = ConstantLocalFunction; - explicit ConstantFunction(T const& value) + public: + /// Constructor. Stores the function `fct` and creates an `EntitySet`. + ConstantGridFunction(T const& value, GridView const& gridView) : value_(value) + , entitySet_(gridView) {} - template - Range const& operator()(GlobalCoordinate const& /*x*/) const + /// Return the evaluated functor at global coordinates + Range const& operator()(Domain const& /*x*/) const { return value_; } - friend Derivative derivative(ConstantFunction const& /*f*/) + friend auto localFunction(ConstantGridFunction const& gf) { - return Derivative{DerivativeRange(0)}; + return LocalFunction{gf.value_}; } - friend int order(ConstantFunction const& /*f*/) + EntitySet const& entitySet() const { - return 0; + return entitySet_; } private: T value_; + EntitySet entitySet_; }; - /** @} **/ + namespace Concepts { /** \addtogroup Concepts @@ -79,10 +142,6 @@ namespace AMDiS struct ConstantToGridFunction> : ConstantToGridFunction {}; - template - struct ConstantToGridFunction> - : ConstantToGridFunction {}; - template struct ConstantToGridFunction> : ConstantToGridFunction {}; @@ -104,11 +163,12 @@ namespace AMDiS * @{ **/ + /// Generator for a GridFunction representing a constant value template )> auto makeGridFunction(T const& value, GridView const& gridView, Dune::PriorityTag<2>) { - return makeGridFunction(ConstantFunction{value}, gridView, Dune::PriorityTag<3>{}); + return ConstantGridFunction{value, gridView}; } /** @} **/ diff --git a/dune/amdis/gridfunctions/CoordsGridFunction.hpp b/dune/amdis/gridfunctions/CoordsGridFunction.hpp index 3afad41f9ccaa6fde40ae720181413181fb0fee5..994538419c048a933e9d28ac1f4658f0e6dd2ff5 100644 --- a/dune/amdis/gridfunctions/CoordsGridFunction.hpp +++ b/dune/amdis/gridfunctions/CoordsGridFunction.hpp @@ -2,122 +2,136 @@ #include #include -#include +#include #include namespace AMDiS { - /** - * \addtogroup GridFunctions - * @{ - **/ - - struct CoordsFunction + namespace Operation { - struct Derivative + /** \addtogroup operations + * @{ + **/ + + /// A functor that evaluates to the global coordinates + struct CoordsFunction { + struct Derivative + { + template + Dune::DiagonalMatrix const& operator()(Dune::FieldVector const& x) const + { + return Dune::DiagonalMatrix{T(1)}; + } + }; + template - Dune::IdentityMatrix const& operator()(Dune::FieldVector const& x) const + Dune::FieldVector const& operator()(Dune::FieldVector const& x) const + { + return x; + } + + friend Derivative partial(CoordsFunction const& /*f*/, index_t<0>) + { + return Derivative{}; + } + + friend int order(CoordsFunction const& /*f*/, int /*d*/) { - return Dune::IdentityMatrix{}; + return 1; } }; - template - Dune::FieldVector const& operator()(Dune::FieldVector const& x) const - { - return x; - } - friend Derivative derivative(CoordsFunction const& /*f*/) + /// A functor that evaluates to a component of the global coordinates + struct CoordsCompFunction { - return Derivative{}; - } + struct Derivative + { + explicit Derivative(int comp) + : comp_(comp) + {} - friend int order(CoordsFunction const& /*f*/) - { - return 1; - } - }; + template + Dune::FieldVector operator()(Dune::FieldVector const& x) const + { + Dune::FieldVector result(0); + result[comp_] = T(1); + return result; + } - struct CoordsCompFunction - { - struct Derivative - { - explicit Derivative(int comp) + private: + int comp_; + }; + + /// Constructor. Stores the component `comp` of the coordinates. + explicit CoordsCompFunction(int comp) : comp_(comp) {} template - Dune::FieldVector const& operator()(Dune::FieldVector const& x) const + T const& operator()(Dune::FieldVector const& x) const + { + return x[comp_]; + } + + friend Derivative partial(CoordsCompFunction const& f, index_t<0>) { - Dune::FieldVector result(0); - result[comp_] = T(1); + return Derivative{f.comp_}; + } - return result; + friend int order(CoordsCompFunction const& /*f*/, int /*d*/) + { + return 1; } private: int comp_; }; - explicit CoordsCompFunction(int comp) - : comp_(comp) - {} - - template - T const& operator()(Dune::FieldVector const& x) const - { - return x[comp_]; - } - - friend Derivative derivative(CoordsCompFunction const& f) - { - return Derivative{f.comp_}; - } - - friend int order(CoordsCompFunction const& /*f*/) - { - return 1; - } + /** @} **/ - private: - int comp_; - }; - - namespace Concepts { namespace Definition - { - template <> - struct IsPreGridFunction - : std::true_type {}; - - template <> - struct IsPreGridFunction - : std::true_type {}; - }} + } // end namespace Operation /// Generator for \ref CoordsFunction inline auto X() { - return CoordsFunction{}; + return Operation::CoordsFunction{}; } /// Generator for \ref CoordsCompFunction inline auto X(int comp) { - return CoordsCompFunction{comp}; + return Operation::CoordsCompFunction{comp}; } + namespace Traits + { + template <> + struct IsPreGridFunction + : std::true_type {}; + + template <> + struct IsPreGridFunction + : std::true_type {}; + } + + /** + * \addtogroup GridFunctions + * @{ + **/ + /// Generator for GridFunction from the pre-GridFunction \ref CoordsFunction template - auto makeGridFunction(CoordsFunction const& f, GridView const& gridView, Dune::PriorityTag<1>) + auto makeGridFunction(Operation::CoordsFunction const& f, GridView const& gridView, Dune::PriorityTag<1>) { return makeAnalyticGridFunction(f, gridView); } + /// Generator for GridFunction from the pre-GridFunction \ref CoordsCompFunction template - auto makeGridFunction(CoordsCompFunction const& f, GridView const& gridView, Dune::PriorityTag<1>) + auto makeGridFunction(Operation::CoordsCompFunction const& f, GridView const& gridView, Dune::PriorityTag<1>) { return makeAnalyticGridFunction(f, gridView); } diff --git a/dune/amdis/gridfunctions/DOFVectorView.hpp b/dune/amdis/gridfunctions/DOFVectorView.hpp index 547afdee3e4c6b3c6f8151b463c9c34b46b06bd9..7f13f72d36862e272a6a892e0bb2f275d0b0238e 100644 --- a/dune/amdis/gridfunctions/DOFVectorView.hpp +++ b/dune/amdis/gridfunctions/DOFVectorView.hpp @@ -4,11 +4,15 @@ #include #include +#include #include #include #include #include +#include +#include + namespace AMDiS { template @@ -92,6 +96,12 @@ namespace AMDiS /// Evaluate Gradient at bound element in local coordinates Range operator()(Domain const& x) const; + friend int order(GradientLocalFunction const& self) + { + assert( self.bound_ ); + return std::max(0, getPolynomialDegree(*self.subTree_)-1); + } + /// Return the bound element Element const& localContext() const { @@ -165,6 +175,12 @@ namespace AMDiS return GradientLocalFunction{*localFunction.globalFunction_}; } + friend int order(LocalFunction const& self) + { + assert( self.bound_ ); + return getPolynomialDegree(*self.subTree_); + } + /// Return the bound element Element const& localContext() const { @@ -207,12 +223,6 @@ namespace AMDiS return LocalFunction{self}; } - friend int order(DOFVectorView const& self) - { - auto node = child(self.basis().localView().tree(), self.treePath()); - return getPolynomialDegree(self.basis().gridView(), node); - } - EntitySet const& entitySet() const { return entitySet_; @@ -281,7 +291,7 @@ namespace AMDiS auto&& gridFct = makeGridFunction(std::forward(expr), basis.gridView(), Dune::PriorityTag<42>{}); DOFVector tmp(basis, "tmp"); - Dune::Functions::interpolate(basis, treePath, tmp, gridFct); + Dune::Functions::interpolate(basis, treePath, tmp, std::forward(gridFct)); // move data from temporary vector into stored DOFVector mutableDofVector_->getVector() = std::move(tmp.getVector()); diff --git a/dune/amdis/gridfunctions/DOFVectorView.inc.hpp b/dune/amdis/gridfunctions/DOFVectorView.inc.hpp index b617f8849fd408821b28e86d18fcf8b1a6ecdd79..141aac5680689f866c9f7195b40677c6ef387b81 100644 --- a/dune/amdis/gridfunctions/DOFVectorView.inc.hpp +++ b/dune/amdis/gridfunctions/DOFVectorView.inc.hpp @@ -1,5 +1,7 @@ #pragma once +#include + namespace AMDiS { template diff --git a/dune/amdis/gridfunctions/DerivativeGridFunction.hpp b/dune/amdis/gridfunctions/DerivativeGridFunction.hpp index 98505f2c6a4c6167d71bf2e55e50ba070dff38f9..48eaf4db6bcc3c8429edf94bd05021fd9c4b9b7e 100644 --- a/dune/amdis/gridfunctions/DerivativeGridFunction.hpp +++ b/dune/amdis/gridfunctions/DerivativeGridFunction.hpp @@ -13,7 +13,7 @@ namespace AMDiS * @{ **/ - /// \brief A Gridfunction that returns the derivative when calling localFunction + /// A Gridfunction that returns the derivative when calling localFunction template class DerivativeGridFunction { @@ -32,6 +32,7 @@ namespace AMDiS : gridFct_(gridFct) {} + /// no global derivative available Range operator()(Domain const& x) const { error_exit("Not implemented"); @@ -43,11 +44,6 @@ namespace AMDiS return derivative(localFunction(gf.gridFct_)); } - friend int order(DerivativeGridFunction const& gf) - { - return std::max(0, order(gf.gridFct_)-1); - } - EntitySet const& entitySet() const { return gridFct_.entitySet(); @@ -57,16 +53,43 @@ namespace AMDiS GridFct gridFct_; }; - - /// Generator function for \ref DerivativeGridFunction expressions template && - Concepts::HasLocalFunctionDerivative)> + REQUIRES(not Concepts::HasDerivative && + Concepts::HasLocalFunctionDerivative)> auto derivative(GridFct const& gridFct) { return DerivativeGridFunction{gridFct}; } + + template + struct DerivativePreGridFunction + { + Expr expr; + }; + + namespace Traits + { + template + struct IsPreGridFunction> + : std::true_type {}; + } + + + /// Generator function for \ref DerivativeGridFunction expressions + template + auto gradientAtQP(Expr const& expr) + { + return DerivativePreGridFunction{expr}; + } + + /// Generator function for \ref AnalyticGridFunction expressions from \ref AnalyticPreGridFunction. + template + auto makeGridFunction(DerivativePreGridFunction const& pre, GridView const& gridView, Dune::PriorityTag<3>) + { + return derivative(makeGridFunction(pre.expr, gridView, Dune::PriorityTag<42>{})); + } + /** @} **/ } // end namespace AMDiS diff --git a/dune/amdis/gridfunctions/FunctorGridFunction.hpp b/dune/amdis/gridfunctions/FunctorGridFunction.hpp index 2d7f579e889054b3cd07650f5e2223358dbd036d..bf7b82654cc1d927c676e4dd91e924c0a2a80c6f 100644 --- a/dune/amdis/gridfunctions/FunctorGridFunction.hpp +++ b/dune/amdis/gridfunctions/FunctorGridFunction.hpp @@ -3,12 +3,15 @@ #include #include +#include + +#include #include #include #include #include -#define AMDIS_FUNCTOR_DEFAULT_ORDER 1 +#define AMDIS_FUNCTOR_DEFAULT_ORDER 4 namespace AMDiS { @@ -35,181 +38,236 @@ namespace AMDiS "All GridFunctions must have the same EntitySet." ); }; - struct SumFunctor - { - auto operator()() const { return 0; } - template - auto operator()(T0 const& t0) const { return t0; } - template - auto operator()(T0 const& t0, T1 const& t1) const { return t0 + t1; } - template - auto operator()(T0 const& t0, Ts const&... ts) const { return t0 + operator()(ts...); } - }; - } // end namespace Impl + template + class FunctorLocalFunction; - /// \brief A Gridfunction that applies a functor to the evaluated Gridfunctions - template - class FunctorGridFunction + template + class FunctorLocalFunction { public: - using Range = typename std::result_of::type; - using Domain = typename Impl::DomainType::type; - - using EntitySet = typename Impl::EntitySetType::type; + using Range = R; + using Domain = D; public: + /// Constructor. Stores copies of the functor and localFunction(gridfunction)s. + template + FunctorLocalFunction(Functor const& fct, LocalFcts&&... localFcts) + : fct_(fct) + , localFcts_(std::forward(localFcts)...) + {} - class LocalFunction + /// Calls \ref bind for all localFunctions + template + void bind(Element const& element) { - template - using LocalFct = std::decay_t()))>; - - using Element = typename EntitySet::Element; + forEach(localFcts_, [&](auto& localFct) { + localFct.bind(element); + }); + } - public: - using Range = typename FunctorGridFunction::Range; - using Domain = typename Impl::DomainType...>::type; + /// Calls \ref unbind for all localFunctions + void unbind() + { + forEach(localFcts_, [&](auto& localFct) { + localFct.unbind(); + }); + } - public: - template - LocalFunction(Functor const& fct, GridFcts&&... gridFcts) - : fct_(fct) - , localFcts_(localFunction(std::forward(gridFcts))...) - {} + /// Applies the functor \ref fct_ to the evaluated localFunctions + Range operator()(Domain const& x) const + { + return eval([&x](auto const& localFct) { return localFct(x); }, + MakeSeq_t{}); + } + private: + template + auto eval(Op op, Indices) const + { + return fct_(op(std::get(localFcts_))...); + } - void bind(Element const& element) - { - forEach(localFcts_, [&](auto& localFct) { - localFct.bind(element); - }); - } + public: + /// Return the stored functor + Functor const& fct() const + { + return fct_; + } - void unbind() - { - forEach(localFcts_, [&](auto& localFct) { - localFct.unbind(); - }); - } + std::tuple const& localFcts() const + { + return localFcts_; + } - Range operator()(Domain const& x) const - { - return eval([&x](auto const& localFct) { return localFct(x); }, - MakeSeq_t{}); - } + private: + Functor fct_; + std::tuple localFcts_; + }; - private: - template - auto eval(Op op, Indices) const - { - return fct_(op(std::get(localFcts_))...); - } - private: - Functor fct_; - std::tuple...> localFcts_; - }; + /// \brief A Gridfunction that applies a functor to the evaluated Gridfunctions + template + class FunctorGridFunction + { + public: + using Range = typename std::result_of::type; + using Domain = typename Impl::DomainType::type; + using EntitySet = typename Impl::EntitySetType::type; + using LocalDomain = typename EntitySet::LocalCoordinate; public: + template + using LocalFct = std::decay_t()))>; + + using RawRange = std::decay_t; + using LocalFunction = FunctorLocalFunction...>; + + public: + /// Constructor. Stores copies of the functor and gridfunctions. template explicit FunctorGridFunction(Functor const& fct, GridFcts&&... gridFcts) : fct_(fct) , gridFcts_(std::forward(gridFcts)...) {} + /// Applies the functor \ref fct_ to the evaluated gridfunctions Range operator()(Domain const& x) const { return eval(fct_, [&x](auto const& gridFct) { return gridFct(x); }, MakeSeq_t{}); } + /// Creates a LocalFunction from the LocalFunctions of the GridFunctions friend LocalFunction localFunction(FunctorGridFunction const& gf) { - return gf.expand([&gf](auto const&... gridFcts) { return LocalFunction{gf.fct_, gridFcts...}; }, - MakeSeq_t{}); - } - -#if 0 - template )> - friend auto derivative(FunctorGridFunction const& gf) - { - auto index_seq = MakeSeq_t{}; - return gf.expand_idx([&](auto const... _i) + return Dune::Std::apply([&gf](auto const&... gridFcts) { - return localFunction( - makeFunctorGridFunction(Impl::SumFunctor{}, - return gf.expand([&](auto const&... gfs) { - return makeFunctorGridFunction(partial(_i, gf.fct_), gfs...); - }, index_seq)...) ); - }, index_seq); + return LocalFunction{gf.fct_, localFunction(gridFcts)...}; + }, + gf.gridFcts_); } -#endif + /// Return the stored entityset of the first GridFunction EntitySet const& entitySet() const { return std::get<0>(gridFcts_).entitySet(); } - Functor const& fct() const { return fct_; } - - + private: template auto eval(Outer outer, Inner inner, Indices) const { return outer(inner(std::get(gridFcts_))...); } - template - auto expand(Op op, Indices) const - { - return op(std::get(gridFcts_)...); - } - - template - auto expand_idx(Op op, Indices) const - { - return op(index_...); - } - private: Functor fct_; std::tuple gridFcts_; }; + /// Generator function for \ref FunctorGridFunction expressions + template + auto makeFunctorGridFunction(Functor const& f, GridFcts const&... gridFcts) + { + static_assert(all_of_v...>, + "All passed parameters must be GridFunctions."); + static_assert(Concepts::Evaluable, + "Range types of grid functions are not compatible with the functor."); + return FunctorGridFunction{f, gridFcts...}; + } + - template )> - int order(FunctorGridFunction const& gf) + /// Derivative of the LocalFunction of a \ref FunctorGridFunction, utilizing + /// the chain-rule. Only available of the functor provides partial derivatives. + /* + * d_x(f(lf(x)...)) = sum_i d_i(f)[lf(x)...] * derivative(lf[i]) + */ + template )> + auto derivative(FunctorLocalFunction const& lf) { - static_assert(all_of_v...>, - "All GridFunctions gf must have a free function order(gf)."); + auto index_seq = std::make_index_sequence{}; + + // d_i(f)[lgfs...] * lgfs_i + auto term_i = [&](auto const _i) + { + auto di_f = Dune::Std::apply([&](auto const&... lgfs) { + return makeFunctorGridFunction(partial(lf.fct(), _i), lgfs...); + }, lf.localFcts()); + + auto const& lgfs_i = std::get<_i>(lf.localFcts()); + return makeFunctorGridFunction(Operation::Multiplies{}, di_f, derivative(lgfs_i)); + }; - return gf.expand([&gf](auto const&... gfs) { return order(gf.fct(), order(gfs)...); }, - MakeSeq_t{}); + // sum_i [ d_i(f)[lgfs...] * derivative(lgfs_i) ] + auto gridFct = Dune::Std::apply([&](auto const... _i) + { + return makeFunctorGridFunction(Operation::Plus{}, term_i(_i)...); + }, index_seq); + + return localFunction(gridFct); } - template )> - int order(FunctorGridFunction const& gf) + +#ifdef AMDIS_HAS_CXX_CONSTEXPR_IF + + template + int order(FunctorLocalFunction const& lf) { - return AMDIS_FUNCTOR_DEFAULT_ORDER; + if constexpr (Concepts::HasFunctorOrder && + all_of_v...>) + { + return Dune::Std::apply([&lf](auto const&... lfs) { return order(lf.fct(), order(lfs)...); }, lf.localFcts()); + } + else if constexpr (not Concepts::HasFunctorOrder && + all_of_v...>) + { + return Dune::Std::apply([](auto const&... lfs) { return Math::max(order(lfs)...); }, lf.localFcts()); + } + else { + return AMDIS_FUNCTOR_DEFAULT_ORDER; + } } +#else + /// Calculate the polynomial order for functors f providing a free functions + /// order(f, degs...), where degs are the order of the LocalFunctions + template + && all_of_v...>)> + int order(FunctorLocalFunction const& lf) + { + return Dune::Std::apply([&lf](auto const&... gfs) { return order(lf.fct(), order(gfs)...); }, + lf.localFcts()); + } - /// Generator function for \ref FunctorGridFunction expressions - template ...>)> - auto makeFunctorGridFunction(Functor const& f, GridFcts const&... gridFcts) + /// Specialization for \ref order for functor not providing an order free function. + /// Take the maximum of all provided LocalFunction orders + template + && all_of_v...>)> + int order(FunctorLocalFunction const& lf) { - return FunctorGridFunction{f, gridFcts...}; + return Dune::Std::apply([](auto const&... gfs) { return Math::max(order(gfs)...); }, + lf.localFcts()); + } + + /// Specialization for \ref order for functor not providing an order free function + /// and some LocalFunctions do not provide an order function. Take a constant value. + template )...>)> + int order(FunctorLocalFunction const& /*lf*/) + { + return AMDIS_FUNCTOR_DEFAULT_ORDER; } +#endif + /// PreGridFunction related to \ref FunctorGridFunction, \relates FunctorGridFunction template struct FunctorPreGridFunction { @@ -221,36 +279,33 @@ namespace AMDiS Functor fct_; std::tuple gridFcts_; - - template - auto expand(Op op, Indices) const - { - return op(std::get(gridFcts_)...); - } }; - namespace Concepts { namespace Definition + namespace Traits { template struct IsPreGridFunction> : std::true_type {}; - }} + } + /// Generator function for \ref FunctorPreGridFunction, i.e., creates a + /// pre GridFunction. \relates FunctorGridFunction template - auto apply(Functor const& f, GridFcts const&... gridFcts) + auto invokeAtQP(Functor const& f, GridFcts const&... gridFcts) { return FunctorPreGridFunction{f, gridFcts...}; } + /// Generator function from a \ref FunctorPreGridFunction to a \ref FunctorGridFunction. template auto makeGridFunction(FunctorPreGridFunction const& preGridFct, GridView const& gridView, Dune::PriorityTag<5>) { - return preGridFct.expand([&](auto const&... pgf) { + return Dune::Std::apply([&](auto const&... pgf) { return makeFunctorGridFunction(preGridFct.fct_, makeGridFunction(pgf, gridView, Dune::PriorityTag<42>{})...); - }, MakeSeq_t{}); + }, preGridFct.gridFcts_); } /** @} **/ diff --git a/dune/amdis/gridfunctions/GridFunctionConcepts.hpp b/dune/amdis/gridfunctions/GridFunctionConcepts.hpp index bdd2b3ac67159c9e1fa121def86e559b41786e3f..9ecb8c9f7d3e8a812746f3e6ec0470decd515293 100644 --- a/dune/amdis/gridfunctions/GridFunctionConcepts.hpp +++ b/dune/amdis/gridfunctions/GridFunctionConcepts.hpp @@ -8,9 +8,12 @@ namespace AMDiS { - // forward declaration - template - struct FunctorPreGridFunction; + namespace Traits + { + template + struct IsPreGridFunction + : std::false_type {}; + } // end namespace Traits namespace Concepts { @@ -34,6 +37,18 @@ namespace AMDiS auto requires_(F&& f) -> decltype( derivative(localFunction(f)) ); }; + struct HasLocalFunctionOrder + { + template + auto requires_(F&& f) -> decltype( order(localFunction(f)) ); + }; + + struct HasPartial + { + template + auto requires_(F&& f) -> decltype( partial(f, index_<0>) ); + }; + struct HasOrder { template @@ -59,10 +74,6 @@ namespace AMDiS constexpr bool CallableDow = Concepts::Callable>; - template - struct IsPreGridFunction - : std::false_type {}; - } // end namespace Definition @@ -75,6 +86,12 @@ namespace AMDiS template constexpr bool HasLocalFunctionDerivative = models; + template + constexpr bool HasLocalFunctionOrder = models; + + template + constexpr bool HasPartial = models; + template constexpr bool HasOrder = models; @@ -99,7 +116,7 @@ namespace AMDiS template constexpr bool AnyGridFunction = any_of_v>...> || - any_of_v>::value...>; + any_of_v>::value...>; } // end namespace Concepts diff --git a/dune/amdis/gridfunctions/Integrate.hpp b/dune/amdis/gridfunctions/Integrate.hpp new file mode 100644 index 0000000000000000000000000000000000000000..e0d0e9f861c349a2b83374a0ff5dbb28a300c8e1 --- /dev/null +++ b/dune/amdis/gridfunctions/Integrate.hpp @@ -0,0 +1,35 @@ +#pragma once + +#include + +#include + +namespace AMDiS +{ + template + auto integrate(Expr&& expr, GridView const& gridView) + { + auto&& gridFct = makeGridFunction(std::forward(expr), gridView, Dune::PriorityTag<42>{}); + auto localFct = localFunction(gridFct); + + using QuadratureRules = Dune::QuadratureRules; + + using Range = typename std::decay_t::Range; + Range result(0); + + for (auto const& element : elements(gridView)) { + auto geometry = element.geometry(); + + localFct.bind(element); + int degree = order(localFct); + auto const& quad = QuadratureRules::rule(element.type(), degree); + + for (auto const qp : quad) + result += localFct(qp.position()) * geometry.integrationElement(qp.position()) * qp.weight(); + localFct.unbind(); + } + + return result; + } + +} // end namespace AMDiS diff --git a/dune/amdis/gridfunctions/OperationsGridFunction.hpp b/dune/amdis/gridfunctions/OperationsGridFunction.hpp index fd47b7260edc8edde6fb85f12465b4d1caa068b2..5808737454d91c922de649b3dfc42ab3e9bcfd33 100644 --- a/dune/amdis/gridfunctions/OperationsGridFunction.hpp +++ b/dune/amdis/gridfunctions/OperationsGridFunction.hpp @@ -1,7 +1,9 @@ #pragma once -#include +#include #include +#include +#include namespace AMDiS { @@ -9,35 +11,35 @@ namespace AMDiS REQUIRES(Concepts::AnyGridFunction)> auto operator-(Lhs&& lhs) { - return apply(Operation::negate{}, std::forward(lhs)); + return invokeAtQP(Operation::Negate{}, std::forward(lhs)); } template )> auto operator+(Lhs&& lhs, Rhs&& rhs) { - return apply(Operation::plus{}, std::forward(lhs), std::forward(rhs)); + return invokeAtQP(Operation::Plus{}, std::forward(lhs), std::forward(rhs)); } template )> auto operator-(Lhs&& lhs, Rhs&& rhs) { - return apply(Operation::minus{}, std::forward(lhs), std::forward(rhs)); + return invokeAtQP(Operation::Minus{}, std::forward(lhs), std::forward(rhs)); } template )> auto operator*(Lhs&& lhs, Rhs&& rhs) { - return apply(Operation::times{}, std::forward(lhs), std::forward(rhs)); + return invokeAtQP(Operation::Multiplies{}, std::forward(lhs), std::forward(rhs)); } template )> auto operator/(Lhs&& lhs, Rhs&& rhs) { - return apply(Operation::divides{}, std::forward(lhs), std::forward(rhs)); + return invokeAtQP(Operation::Divides{}, std::forward(lhs), std::forward(rhs)); } @@ -45,28 +47,106 @@ namespace AMDiS REQUIRES(Concepts::AnyGridFunction)> auto max(Lhs&& lhs, Rhs&& rhs) { - return apply(Operation::maximum{}, std::forward(lhs), std::forward(rhs)); + return invokeAtQP(Operation::Max{}, std::forward(lhs), std::forward(rhs)); } template )> auto min(Lhs&& lhs, Rhs&& rhs) { - return apply(Operation::minimum{}, std::forward(lhs), std::forward(rhs)); + return invokeAtQP(Operation::Min{}, std::forward(lhs), std::forward(rhs)); } template )> auto abs_max(Lhs&& lhs, Rhs&& rhs) { - return apply(Operation::abs_max{}, std::forward(lhs), std::forward(rhs)); + return invokeAtQP(Operation::AbsMax{}, std::forward(lhs), std::forward(rhs)); } template )> auto abs_min(Lhs&& lhs, Rhs&& rhs) { - return apply(Operation::abs_min{}, std::forward(lhs), std::forward(rhs)); + return invokeAtQP(Operation::AbsMin{}, std::forward(lhs), std::forward(rhs)); + } + + // unary vector operations + + template )> + auto sum(Vec&& vec) + { + return invokeAtQP([](auto const& v) { return sum(v); }, std::forward(vec)); + } + + template )> + auto unary_dot(Vec&& vec) + { + return invokeAtQP(Operation::UnaryDot{}, std::forward(vec)); + } + + template )> + auto one_norm(Vec&& vec) + { + return invokeAtQP([](auto const& v) { return one_norm(v); }, std::forward(vec)); + } + + template )> + auto two_norm(Vec&& vec) + { + return invokeAtQP(Operation::TwoNorm{}, std::forward(vec)); + } + + template )> + auto p_norm(Vec&& vec) + { + return invokeAtQP([](auto const& v) { return p_norm

(v); }, std::forward(vec)); + } + + template )> + auto infty_norm(Vec&& vec) + { + return invokeAtQP([](auto const& v) { return infty_norm(v); }, std::forward(vec)); + } + + template )> + auto trans(Mat&& mat) + { + return invokeAtQP([](auto const& m) { return trans(m); }, std::forward(mat)); + } + + // binary vector operations + + template )> + auto dot(Lhs&& lhs, Rhs&& rhs) + { + return invokeAtQP(Operation::Dot{}, + std::forward(lhs), std::forward(rhs)); + } + + template )> + auto distance(Lhs&& lhs, Rhs&& rhs) + { + using namespace Operation; + return invokeAtQP(compose(TwoNorm{}, Minus{}), + std::forward(lhs), std::forward(rhs)); + } + + template )> + auto outer(Lhs&& lhs, Rhs&& rhs) + { + return invokeAtQP([](auto const& v, auto const& w) { return outer(v,w); }, + std::forward(lhs), std::forward(rhs)); } } // end namespace AMDiS diff --git a/dune/amdis/io/CMakeLists.txt b/dune/amdis/io/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..2f050333cbbe94f837fa94021c2f673cba323864 --- /dev/null +++ b/dune/amdis/io/CMakeLists.txt @@ -0,0 +1,5 @@ +#install headers + +install(FILES + FileWriterInterface.hpp +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis/io) diff --git a/dune/amdis/linear_algebra/CMakeLists.txt b/dune/amdis/linear_algebra/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..1fdefd5df83cf04b611e520c25964c35fc847d0a --- /dev/null +++ b/dune/amdis/linear_algebra/CMakeLists.txt @@ -0,0 +1,29 @@ +#install headers + +install(FILES + HierarchicWrapper.hpp + LinearAlgebraBse.hpp + LinearSolverInterface.hpp + PreconditionerInterface.hpp + RunnerInterface.hpp + SolverInfo.hpp + mtl/BITL_Preconditioner.hpp + mtl/BITL_Solver.hpp + mtl/BlockMTLMatrix.hpp + mtl/BlockMTLVector.hpp + mtl/Copy.hpp + mtl/DOFMatrix.hpp + mtl/DOFVector.hpp + mtl/ITL_Preconditioner.hpp + mtl/ITL_Solver.hpp + mtl/KrylovRunner.hpp + mtl/LinearSolver.hpp + mtl/Mapper.hpp + mtl/MTLDenseVector.hpp + mtl/Preconditioner.hpp + mtl/SystemMatrix.hpp + mtl/SystemVector.hpp + mtl/UmfpackRunner.hpp +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis/linear_algebra) + +add_subdirectory("mtl") \ No newline at end of file diff --git a/dune/amdis/linear_algebra/istl/CMakeLists.txt b/dune/amdis/linear_algebra/istl/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..3213d6b309e5cd2379ea3bf408e6285f9d485ba1 --- /dev/null +++ b/dune/amdis/linear_algebra/istl/CMakeLists.txt @@ -0,0 +1,17 @@ +#dune_library_add_sources(amdis SOURCES +# SystemMatrix.cpp +# SystemVector.cpp +#) + +install(FILES + DOFMatrix.hpp + DOFVector.hpp + ISTL_Preconditioner.hpp + ISTL_Solver.hpp + ISTLRunner.hpp + LinearSolver.hpp + Preconditioner.hpp + SystemMatrix.hpp + SystemVector.hpp + UmfpackRunner.hpp +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis/linear_algebra/istl) diff --git a/dune/amdis/linear_algebra/mtl/CMakeLists.txt b/dune/amdis/linear_algebra/mtl/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..490d448f2e1d9ccb42215f55c430a2d7d566a88e --- /dev/null +++ b/dune/amdis/linear_algebra/mtl/CMakeLists.txt @@ -0,0 +1,24 @@ +#dune_library_add_sources(amdis SOURCES +# SystemMatrix.cpp +# SystemVector.cpp +#) + +install(FILES + BITL_Preconditioner.hpp + BITL_Solver.hpp + BlockMTLMatrix.hpp + BlockMTLVector.hpp + Copy.hpp + DOFMatrix.hpp + DOFVector.hpp + ITL_Preconditioner.hpp + ITL_Solver.hpp + KrylovRunner.hpp + LinearSolver.hpp + Mapper.hpp + MTLDenseVector.hpp + Preconditioner.hpp + SystemMatrix.hpp + SystemVector.hpp + UmfpackRunner.hpp +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis/linear_algebra/mtl) diff --git a/dune/amdis/operations/Arithmetic.hpp b/dune/amdis/operations/Arithmetic.hpp new file mode 100644 index 0000000000000000000000000000000000000000..d7d1a2816eae0caedcf8cc35d83a7229659f297b --- /dev/null +++ b/dune/amdis/operations/Arithmetic.hpp @@ -0,0 +1,256 @@ +#pragma once + +#include + +#include +#include +#include +#include + +namespace AMDiS +{ + namespace Operation + { + /** \addtogroup operations + * @{ + **/ + + /// Operation that represents A+B + struct Plus + { +#ifdef AMDIS_HAS_CXX_FOLD_EXPRESSIONS + template + constexpr auto operator()(Ts const&... ts) const + { + return (ts + ...); + } +#else + template + constexpr auto operator()(T0 const& t0, T1 const& t1) const + { + return t0 + t1; + } + + template + constexpr auto operator()(T0 const& t0, Ts const&... ts) const + { + return t0 + operator()(ts...); + } +#endif + }; + + // [0] + g + ... => g + ... + template + struct ComposerBuilder + : ComposerBuilder {}; + + // ...+ g + [0] => ... + g + template + struct ComposerBuilder + : ComposerBuilder {}; + + // [0] + [0] => [0] + template <> + struct ComposerBuilder + : ComposerBuilder {}; + + + template + int order(Plus, Int... orders) + { + return Math::max(int(orders)...); + } + + template + auto partial(Plus, index_t) + { + static_assert((I < 2), "Derivatives of `Plus` only defined for the binary case."); + return One{}; + } + + // ------------------------------------------------------------------------- + + /// Operation that represents A-B + struct Minus + { + template + constexpr auto operator()(T const& lhs, S const& rhs) const + { + return lhs - rhs; + } + + friend int order(Minus, int lhs, int rhs) + { + return Math::max(lhs, rhs); + } + + friend auto partial(Minus, index_t<0>) + { + return One{}; + } + + friend auto partial(Minus, index_t<1>) + { + return StaticConstant{}; + } + }; + + // [0] - g => -g + // template + // struct ComposerBuilder + // : ComposerBuilder {}; + + // g - [0] => g + template + struct ComposerBuilder + : ComposerBuilder {}; + + // [0] - [0] => [0] + template <> + struct ComposerBuilder + : ComposerBuilder {}; + + + // ------------------------------------------------------------------------- + + /// Operation that represents A*B + struct Multiplies + { +#ifdef AMDIS_HAS_CXX_FOLD_EXPRESSIONS + template + constexpr auto operator()(Ts const&... ts) const + { + return (ts * ...); + } +#else + template + constexpr auto operator()(T const& lhs, S const& rhs) const + { + return lhs * rhs; + } +#endif + }; + + // [0] * g * ... => [0] + template + struct ComposerBuilder + : ComposerBuilder {}; + + // ...* g * [0] => [0] + template + struct ComposerBuilder + : ComposerBuilder {}; + + // [0] * [0] => [0] + template <> + struct ComposerBuilder + : ComposerBuilder {}; + + + template + int order(Multiplies, Int... orders) + { + return Plus{}(int(orders)...); + } + + // only for binary * + // d_0 (x * y) = y, d_1 (x * y) = x + template + auto partial(Multiplies, index_t) + { + static_assert((I < 2), "Derivatives of `Multiplies` only defined for the binary case."); + return Component<1-I>{}; + } + + // ------------------------------------------------------------------------- + + /// Operation that represents x^p + template + struct Pow + { + static_assert(p > 1, "Exponent in power must be non-negative!"); + + template + constexpr auto operator()(T const& x) const + { + return Math::pow

(x); + } + + friend int order(Pow, int d) + { + return p*d; + } + + friend auto partial(Pow, index_t<0>) + { + return compose(Multiplies{}, StaticConstant{}, Pow{}); + } + }; + + using Sqr = Pow<2>; + + template <> + struct Pow<1> + : public Id {}; + + template <> + struct Pow<0> + : public Zero {}; + + // ------------------------------------------------------------------------- + + /// Operation that represents A/B + struct Divides + { + template + constexpr auto operator()(T const& lhs, S const& rhs) const + { + return lhs / rhs; + } + + friend int order(Divides, int lhs, int rhs) + { + // [[warning::order_not_exact]] + return lhs + rhs + 1; + } + + // d_0 f(x,y) = 1 / y + friend auto partial(Divides, index_t<0>) + { + return compose(Divides{}, One{}, Component<1>{}); + } + + // d_1 f(x,y) = (y - x)/y^2 + friend auto partial(Divides, index_t<1>) + { + return compose(Divides{}, compose(Minus{}, Component<1>{}, Component<0>{}), + compose(Pow<2>{}, Component<1>{})); + } + }; + + // ------------------------------------------------------------------------- + + /// Operation that represents A-B + struct Negate + { + template + constexpr auto operator()(T const& x) const + { + return -x; + } + + friend int order(Negate, int d) + { + return d; + } + + friend auto partial(Negate, index_t<0>) + { + return StaticConstant{}; + } + }; + + /** @} **/ + + } // end namespace Operation +} // end namespace AMDiS diff --git a/dune/amdis/operations/Basic.hpp b/dune/amdis/operations/Basic.hpp new file mode 100644 index 0000000000000000000000000000000000000000..0de39d3384f90639af54acf9495f6b5379536af1 --- /dev/null +++ b/dune/amdis/operations/Basic.hpp @@ -0,0 +1,127 @@ +#pragma once + +#include + +#include +#include +#include + +namespace AMDiS +{ + namespace Operation + { + /** \defgroup operations Operations + * \brief Functors, representing unary/binary operations used in expressions. + * @{ + **/ + + /// Functor representing a static constant value + template + struct StaticConstant + { + template + constexpr T operator()(Ts const&... /*args*/) const + { + return value; + } + }; + + using Zero = StaticConstant; + using One = StaticConstant; + + template + constexpr int order(StaticConstant const&, Int...) + { + return 0; + } + + template + constexpr auto partial(StaticConstant, index_t) + { + return Zero{}; + } + + // ------------------------------------------------------------------------- + + /// (Unary-)Functor representing the identity + struct Id + { + template + constexpr T const& operator()(T const& x) const + { + return x; + } + + friend constexpr int order(Id, int d) + { + return d; + } + + friend auto partial(Id, index_t<0>) + { + return Zero{}; + } + }; + + // ------------------------------------------------------------------------- + + /// Functor representing a constant value + template + struct Constant + { + constexpr Constant(T value) + : value_(value) + {} + + template + constexpr T const& operator()(Ts const&... /*args*/) const + { + return value_; + } + + private: + T value_; + }; + + template + constexpr int order(Constant const&, Int...) + { + return 0; + } + + template + constexpr auto partial(Constant const&, index_t) + { + return Zero{}; + } + + // ------------------------------------------------------------------------- + + template + struct Component + { + template + constexpr auto&& operator()(Ts&&... args) const + { + return std::get(std::forward_as_tuple(std::forward(args)...)); + } + }; + + template + constexpr int order(Component const&, Int... degrees) + { + return std::get(std::tie(degrees...)); + } + + template + constexpr auto partial(Component, index_t) + { + return StaticConstant{}; + } + + // ------------------------------------------------------------------------- + + /** @} **/ + + } // end namespace Operation +} // end namespace AMDiS diff --git a/dune/amdis/operations/CMakeLists.txt b/dune/amdis/operations/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..624d7666b5f147e73576fb6ca1572d38d0df2c4e --- /dev/null +++ b/dune/amdis/operations/CMakeLists.txt @@ -0,0 +1,10 @@ +#install headers + +install(FILES + Arithmetic.hpp + Basic.hpp + CMath.hpp + Composer.hpp + FieldMatVec.hpp + MaxMin.hpp +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis/operations) diff --git a/dune/amdis/operations/CMath.hpp b/dune/amdis/operations/CMath.hpp new file mode 100644 index 0000000000000000000000000000000000000000..5fbd8a83e3746b03d698e5c30ae96ab12068cdfc --- /dev/null +++ b/dune/amdis/operations/CMath.hpp @@ -0,0 +1,136 @@ +#pragma once + +#include + +#include + +/// Macro that generates a unary functor with explicitly given function and derivative. +/** + * \p NAME Name of the class. + * \p DEGREE Expression in 'd0' that gives the polynomial degree, with + * 'd' the degree of the argument passed to the functor. + * \p FCT Name of a unary c++-function that represents the functor. + * \p DERIV Name of a unary c++-function that represents the derivative of this functor. + */ +#ifndef AMDIS_MAKE_UNARY_FUNCTOR +#define AMDIS_MAKE_UNARY_FUNCTOR( NAME, DEGREE, FCT, DERIV ) \ + struct NAME \ + { \ + template \ + constexpr auto operator() (T const& x) const \ + { \ + return FCT ; \ + } \ + friend constexpr int order(NAME, int d) \ + { \ + return DEGREE ; \ + } \ + struct Derivative \ + { \ + template \ + constexpr auto operator() (T const& x) const \ + { \ + return DERIV ; \ + } \ + }; \ + friend constexpr auto partial(NAME, index_t<0>) \ + { \ + return Derivative{} ; \ + } \ + }; +#endif + + +namespace AMDiS +{ + namespace Operation + { + /** \addtogroup operations + * @{ + **/ + + /// Functor that represents the signum function + struct Signum + { + template + constexpr auto operator()(T const& x) const + { + return (x > T{0} ? T{1} : T{-1}); + } + + constexpr friend int order(Signum, int /*d*/) + { + return 0; + } + + constexpr friend auto partial(Signum, index_t<0>) + { + return Zero{}; + } + }; + + // generated unary functors using a macro... + // approximate polynomial order + + AMDIS_MAKE_UNARY_FUNCTOR( Ceil, d, std::ceil(x), 0.0 ) + AMDIS_MAKE_UNARY_FUNCTOR( Floor, d, std::floor(x), 0.0 ) + AMDIS_MAKE_UNARY_FUNCTOR( Sqrt, 3*d, std::sqrt(x), 1.0/(2.0 * std::sqrt(x)) ) + AMDIS_MAKE_UNARY_FUNCTOR( Exp, 2*d, std::exp(x), std::exp(x) ) + AMDIS_MAKE_UNARY_FUNCTOR( Log, 2*d, std::log(x), 1.0/x ) + AMDIS_MAKE_UNARY_FUNCTOR( Sin, 2*d, std::sin(x), std::cos(x) ) + AMDIS_MAKE_UNARY_FUNCTOR( Cos, 2*d, std::cos(x), -std::sin(x) ) + AMDIS_MAKE_UNARY_FUNCTOR( Tan, 2*d, std::tan(x), 1.0 + Math::sqr(std::tan(x)) ) + AMDIS_MAKE_UNARY_FUNCTOR( Asin, 2*d, std::asin(x), 1.0/std::sqrt(1.0 - Math::sqr(x)) ) + AMDIS_MAKE_UNARY_FUNCTOR( Acos, 2*d, std::acos(x), -1.0/std::sqrt(1.0 - Math::sqr(x)) ) + AMDIS_MAKE_UNARY_FUNCTOR( Atan, 2*d, std::atan(x), 1.0/(1.0 + Math::sqr(x)) ) + AMDIS_MAKE_UNARY_FUNCTOR( Sinh, 2*d, std::sinh(x), std::cosh(x) ) + AMDIS_MAKE_UNARY_FUNCTOR( Cosh, 2*d, std::cosh(x), std::sinh(x) ) + AMDIS_MAKE_UNARY_FUNCTOR( Tanh, 2*d, std::tanh(x), 1.0 - Math::sqr(std::tanh(x)) ) + AMDIS_MAKE_UNARY_FUNCTOR( Asinh, 2*d, std::asinh(x), 1.0/std::sqrt(1.0 + Math::sqr(x)) ) + AMDIS_MAKE_UNARY_FUNCTOR( Acosh, 2*d, std::acosh(x), 1.0/std::sqrt(-1.0 + Math::sqr(x)) ) + AMDIS_MAKE_UNARY_FUNCTOR( Atanh, 2*d, std::atanh(x), 1.0/(1.0 - Math::sqr(x)) ) + + + /// Binary functor representing the cmath function std::atan2(v) + struct Atan2 + { + template + constexpr auto operator() (T1 const& x, T2 const& y) const + { + return std::atan2(x, y); + } + + constexpr friend int order(int d0, int d1) + { + return 2*(d0+d1); + } + + struct D0 { + template + constexpr auto operator() (T1 const& x, T2 const& y) const + { + return -y/(Math::sqr(x) + Math::sqr(y)); + } + }; + constexpr friend auto partial(Atan2, index_t<0>) + { + return D0{}; + } + + struct D1 { + template + constexpr auto operator() (T1 const& x, T2 const& y) const + { + return x/(Math::sqr(x) + Math::sqr(y)); + } + }; + constexpr friend auto partial(Atan2, index_t<1>) + { + return D1{}; + } + }; + + /** @} **/ + + } // end namespace Operation +} // end namespace AMDiS diff --git a/dune/amdis/operations/Composer.hpp b/dune/amdis/operations/Composer.hpp new file mode 100644 index 0000000000000000000000000000000000000000..1de7a784117fcca9750d306102e0c9bf2695e577 --- /dev/null +++ b/dune/amdis/operations/Composer.hpp @@ -0,0 +1,116 @@ +#pragma once + +#include +#include + +#include + +#include + +namespace AMDiS +{ + namespace Operation + { + template + struct Composer + { + template + constexpr Composer(F_&& f, Gs_&&... gs) + : f_(std::forward(f)) + , gs_(std::forward(gs)...) + {} + + template + constexpr auto operator()(Ts const&... args) const + { + auto eval = [&](auto const& g) { return g(args...); }; + return Dune::Std::apply([&,this](auto const&... gs) { return f_(eval(gs)...); }, gs_); + } + + F f_; + std::tuple gs_; + }; + +#ifdef DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION + // deduction guide for the class \ref composer + template + Composer(F_&&, Gs_&&...) + -> Composer, std::decay_t...>; +#endif + + + template + struct ComposerBuilder + { + using type = Composer; + + template + static constexpr type build(F_&& f, Gs_&&... gs) + { + return type{std::forward(f), std::forward(gs)...}; + } + }; + + /// Generator function for \ref composer + template + auto compose(F&& f, Gs&&... gs) + { + return ComposerBuilder, std::decay_t...>::build( + std::forward(f), std::forward(gs)...); + } + + // Polynomial order or composed function combines the orders of the sub-functions + template + int order(Composer const& c, Int... degrees) + { + auto deg = [&](auto const& g) { return order(g, int(degrees)...); }; + return Dune::Std::apply([&](auto const&... gs) { return order(c.f_, deg(gs)...); }, c.gs_); + } + + // forward declaration + struct Multiplies; + struct Plus; + + /// Partial derivative of composed function: + /// Implements: // sum_i [ d_i(f)[g...] * d_j(g_i) ] + template + auto partial(Composer const& c, index_t _j) + { + auto index_seq = MakeSeq_t{}; + + // d_i(f)[g...] * d_j(g_i) + auto term_i = [&](auto const _i) + { + auto di_f = Dune::Std::apply([&](auto const&... gs) { + return compose(partial(c.f_, _i), gs...); + }, c.gs_); + + auto const& g_i = std::get<_i>(c.gs_); + return compose(Multiplies{}, di_f, partial(g_i, _j)); + }; + + // sum_i [ d_i(f)[g...] * d_j(g_i) ] + return Dune::Std::apply([&](auto const... _i) + { + return compose(Plus{}, term_i(_i)...); + }, index_seq); + } + + + // some specialization for the composer + + // id(g) => g + template + struct ComposerBuilder + { + using type = F; + + template + static constexpr F build(F_&& f) + { + return F{std::forward(f)}; + } + }; + + } // end namespace Operation +} // end namespace AMDiS diff --git a/dune/amdis/operations/FieldMatVec.hpp b/dune/amdis/operations/FieldMatVec.hpp new file mode 100644 index 0000000000000000000000000000000000000000..64e14c165faba94c17c771f3602d511c831f8e47 --- /dev/null +++ b/dune/amdis/operations/FieldMatVec.hpp @@ -0,0 +1,85 @@ +#pragma once + +#include +#include +#include +#include +#include + +namespace AMDiS +{ + namespace Operation + { + /** \addtogroup operations + * @{ + **/ + + /// (Binary-)Functor representing the euclidean dot-product + struct Dot + { + template + constexpr auto operator()(Dune::FieldVector const& lhs, Dune::FieldVector const& rhs) const + { + return dot(lhs, rhs); + } + + friend constexpr int order(Dot, int d1, int d2) + { + return d1 + d2; + } + + friend constexpr auto partial(Dot, index_t<0>) + { + return Component<1>{}; + } + + friend constexpr auto partial(Dot, index_t<1>) + { + return Component<0>{}; + } + }; + + // ------------------------------------------------------------------------- + + /// (Unary-)Functor representing the euclidean dot-product + struct UnaryDot + { + template + constexpr auto operator()(Dune::FieldVector const& vec) const + { + return unary_dot(vec); + } + + friend constexpr int order(UnaryDot, int d) + { + return 2*d; + } + + friend auto partial(UnaryDot, index_t<0>) + { + return compose(Multiplies{}, StaticConstant{}, Id{}); + } + }; + + // ------------------------------------------------------------------------- + + /// (Unary-)Functor representing the euclidean 2-norm + struct TwoNorm + : public Composer + { + TwoNorm() + : Composer(Sqrt{}, UnaryDot{}) + {} + + template + constexpr auto operator()(Dune::FieldVector const& vec) const + { + return two_norm(vec); + } + }; + + /** @} **/ + + } // end namespace Operation + +} // end namespace AMDiS diff --git a/dune/amdis/operations/MaxMin.hpp b/dune/amdis/operations/MaxMin.hpp new file mode 100644 index 0000000000000000000000000000000000000000..10e71a43db5d31e6459599e78248767a562c5f85 --- /dev/null +++ b/dune/amdis/operations/MaxMin.hpp @@ -0,0 +1,80 @@ +#pragma once + +namespace AMDiS +{ + namespace Operation + { + /// Operation that represents max(A,B) + struct Max + { + template + constexpr auto operator()(T const& lhs, S const& rhs) const + { + static_assert(Concepts::LessThanComparable, "Argument types are not comparable by <."); + return Math::max(lhs, rhs); + } + + friend int order(Max, int lhs, int rhs) + { + // [[warning::order_not_exact]] + return Math::max(lhs, rhs); + } + }; + + // ------------------------------------------------------------------------- + + /// Operation that represents min(A,B) + struct Min + { + template + constexpr auto operator()(T const& lhs, S const& rhs) const + { + static_assert(Concepts::LessThanComparable, "Argument types are not comparable by <."); + return Math::min(lhs, rhs); + } + + friend int order(Min, int lhs, int rhs) + { + // [[warning::order_not_exact]] + return Math::max(lhs, rhs); + } + }; + + // ------------------------------------------------------------------------- + + /// Operation that represents max(|A|,|B|) + struct AbsMax + { + template + constexpr auto operator()(T const& lhs, S const& rhs) const + { + return Math::max(std::abs(lhs), std::abs(rhs)); + } + + friend int order(AbsMax, int lhs, int rhs) + { + // [[warning::order_not_exact]] + return Math::max(lhs, rhs); + } + }; + + // ------------------------------------------------------------------------- + + /// Operation that represents min(|A|,|B|) + struct AbsMin + { + template + constexpr auto operator()(T const& lhs, S const& rhs) const + { + return Math::min(std::abs(lhs), std::abs(rhs)); + } + + friend int order(AbsMin, int lhs, int rhs) + { + // [[warning::order_not_exact]] + return Math::max(lhs, rhs); + } + }; + + } // end namespace Operation +} // end namespace AMDiS diff --git a/dune/amdis/terms/ConstantTerm.hpp b/dune/amdis/terms/ConstantTerm.hpp deleted file mode 100644 index 7b213af09035badb368bf56abc31999294fdd811..0000000000000000000000000000000000000000 --- a/dune/amdis/terms/ConstantTerm.hpp +++ /dev/null @@ -1,58 +0,0 @@ -#pragma once - -#include - -#include -#include -#include - -namespace AMDiS -{ - /** - * \defgroup OperatorTerms Expressions and terms to be used in GenericOperatorTerm. - * @{ - **/ - - /// An expression representing a constant (arithmetic) value, \see constant - template - class ConstantTerm - : public ExpressionBase - { - public: - // TODO: possibly convert to plain type, in case of expression templates. - using value_type = Underlying_t; - - explicit ConstantTerm(value_type const& value) - : ExpressionBase(0) - , value_(value) - {} - - template - value_type const& local(Context const& /*context*/, LocalCoordinate const& /*lambda*/) const - { - return value_; - } - - template - value_type const& global(GlobalCoordinate const& /*x*/) const - { - return value_; - } - - value_type const& operator[](std::size_t /*iq*/) const - { - return value_; - } - - private: - value_type value_; - }; - - - /// generator function for \ref ConstantTerm expressions - template - auto constant(T const& value) { return ConstantTerm{value}; } - - /** @} **/ - -} // end namespace AMDiS diff --git a/dune/amdis/terms/CoordsTerm.hpp b/dune/amdis/terms/CoordsTerm.hpp deleted file mode 100644 index 503dc8ee515d30703d576dc3bb53e2f64af206d2..0000000000000000000000000000000000000000 --- a/dune/amdis/terms/CoordsTerm.hpp +++ /dev/null @@ -1,153 +0,0 @@ -#pragma once - -#include -#include - -#include - -#include -#include - -namespace AMDiS -{ - /** - * \addtogroup OperatorTerms - * @{ - **/ - - namespace Impl - { - template - using CoordsFunctorResult - = std::result_of)>; - - template - constexpr bool CallableDow = - Concepts::Callable>; - -#ifndef AMDIS_DOW - template - using CoordsFunctorResult_t = typename std::conditional_t< - CallableDow, CoordsFunctorResult, std::conditional_t< - CallableDow, CoordsFunctorResult, std::conditional_t< - CallableDow, CoordsFunctorResult, Id> > >::type; -#else - template - using CoordsFunctorResult_t = typename std::conditional_t< - CallableDow, CoordsFunctorResult, Id>::type; -#endif - } - - - /// \brief An expression that evaluates to the current coordinate of a dof or - /// quadrature point with given index. \see eval. - template - class CoordsFunctorTerm - : public ExpressionBase - { - public: - using value_type = Impl::CoordsFunctorResult_t; - static_assert( !std::is_same::value, - "Template `Functor` can not be used as a functor."); - - template > > - explicit CoordsFunctorTerm(F&& f, int degree = 1) - : ExpressionBase(degree) - , fct_(std::forward(f)) - {} - - /// \brief Evaluates the functor \ref fct at all quadrature points and - /// stores the result in a vector \ref values. - template - void init(Context const& context, PointList const& points) - { - values_.resize(points.size()); - for (std::size_t iq = 0; iq < points.size(); ++iq) - values_[iq] = local(context, points[iq].position()); - } - - template - value_type local(Context const& context, LocalCoordinate const& lambda) const - { - return fct_(context.global(lambda)); - } - - template - value_type global(GlobalCoordinate const& x) const - { - return fct_(x); - } - - value_type const& operator[](std::size_t iq) const - { - return values_[iq]; - } - - private: - Functor fct_; - - std::vector values_; - }; - - - /// Generator function for \ref CoordsFunctorTerm expressions - template - auto eval(F&& f) { return CoordsFunctorTerm>{std::forward(f)}; } - - - /// \brief An expression that evaluates to the current coordinate of a dof or - /// quadrature point with given index. \see eval. - template - class CoordsTerm - : public ExpressionBase - { - public: - using value_type = Dune::FieldVector; - - CoordsTerm() - : ExpressionBase(1) - {} - - /// \brief Evaluates the functor \ref fct at all quadrature points and - /// stores the result in a vector \ref values. - template - void init(Context const& context, - PointList const& points) - { - values_.resize(points.size()); - for (std::size_t iq = 0; iq < points.size(); ++iq) - values_[iq] = context.global( points[iq].position() ); - } - - template - value_type const& local(Context const& context, LocalCoordinate const& lambda) const - { - return context.global(lambda); - } - - template - value_type const& global(GlobalCoordinate const& x) const - { - return x; - } - - value_type const& operator[](std::size_t iq) const - { - return values_[iq]; - } - - private: - std::vector values_; - }; - - - /// Generator function for \ref CoordsTerm expressions - template - CoordsTerm X() { return {}; } - - - - /** @} **/ - -} // end namespace AMDiS diff --git a/dune/amdis/terms/DOFVectorTerm.hpp b/dune/amdis/terms/DOFVectorTerm.hpp deleted file mode 100644 index 7d357f00a8371d9efc3b6e35dfb2b9652cc70aea..0000000000000000000000000000000000000000 --- a/dune/amdis/terms/DOFVectorTerm.hpp +++ /dev/null @@ -1,331 +0,0 @@ -#pragma once - -#include -#include -#include - -#include -#include - -#include -#include - -namespace AMDiS -{ - /** - * \addtogroup OperatorTerms - * @{ - **/ - - /// \brief An expression that evalues a DOFVector at a given element, either on the - /// dofs or on the quadrature points. \see valueOf. - template - class DOFVectorTerm - : ExpressionBase - { - public: - using value_type = typename DiscreteGlobalBasisFct::Range; - - explicit DOFVectorTerm(std::shared_ptr const& dofvector, double factor = 1.0) - : dofvector_(dofvector) - , localFunction_(localFunction(*dofvector_)) - , factor_(factor) - {} - - explicit DOFVectorTerm(DiscreteGlobalBasisFct& dofvector, double factor = 1.0) - : localFunction_(localFunction(dofvector)) - , factor_(factor) - {} - - template - void bind(Element const& element) - { - localFunction_.bind(element); - - // TODO: get localView from localFunction - auto localView = dofvector_->basis().localView(); - localView.bind(element); - degree_ = getPolynomialDegree(localView.tree()); - } - - void unbind() - { - localFunction_.unbind(); - } - - template - void init(Context const& context, PointList const& points) - { - values_.resize(points.size()); - for (std::size_t iq = 0; iq < points.size(); ++iq) - values_[iq] = factor_ * local(context, points[iq].position()); - } - - template - value_type local(Context const& /*context*/, LocalCoordinate const& lambda) const - { - return localFunction_(lambda); - } - - template - value_type global(GlobalCoordinate const& x) const - { - error_exit("Not implemented."); - return value_type(0); - } - - value_type const& operator[](std::size_t iq) const - { - return values_[iq]; - } - - // [expects: dofvector_->basis().localView() is bound to an element] - int getDegree() const - { - return degree_; - } - - private: - std::shared_ptr dofvector_ = nullptr; - typename DiscreteGlobalBasisFct::LocalFunction localFunction_; - double factor_; - - int degree_ = 0; - std::vector values_; - }; - - /// Generator function for \ref DOFVectorTerm expressions. - /** Generates a \ref DOFVectorTerm, that takes a \ref DOFVector \p vector, - * an optional factor \p factor the value of the DOFVector at the quadrature - * points are scaled with. - **/ - template - auto valueOf(std::shared_ptr> const& dofvector, - double factor = 1.0) - { - using DiscreteGlobalBasisFct = Dune::Functions::DiscreteGlobalBasisFunction; - return DOFVectorTerm{dofvector, factor}; - } - - -// ----------------------------------------------------------------------------- - -#if 0 - - /// \brief An expression that evalues the gradient of a DOFVector at a given - /// element, either on the DOFs or on the quadrature points. \see gradientOf. - template - class GradientTerm - { - using LocalFE = typename DOFVectorType::FeSpace::LocalView::Tree::FiniteElement; - using LocalBasis = typename LocalFE::Traits::LocalBasisType; - using Jacobian = typename LocalBasis::Traits::JacobianType; - - static constexpr int dim = DOFVectorType::FeSpace::GridView::dimension; - static constexpr int dow = DOFVectorType::FeSpace::GridView::dimensionworld; - - using R = typename LocalBasis::Traits::RangeFieldType; - using GradientLocalType = Dune::FieldVector; - using GradientGlobalType = Dune::FieldVector; - - public: - using value_type = GradientGlobalType; - using field_type = Value_t; - - explicit GradientTerm(DOFVectorType const& dofvector, double factor = 1.0) - : vector(dofvector.getVector()) - , factor(factor) - , localView_(dofvector.getFeSpace().localView()) - , localIndexSet_(dofvector.getFeSpace().localIndexSet()) - {} - - template - void bind(Element const& element) - { - localView_.bind(element); - localIndexSet_.bind(localView_); - int deg = getPolynomialDegree(localView_.tree()); - - auto t = element.type(); - degree_ = t.isSimplex() ? Math::max(0, deg-1) : deg; - } - - template - void init(Element const& element, - PointList const& points) - { - AMDIS_FUNCNAME("DOFVectorTerm::init()"); - - auto geometry = element.geometry(); - const auto& localFiniteElem = localView_.tree().finiteElement(); - const std::size_t nBasisFct = localFiniteElem.size(); - - std::vector shapeGradients(nBasisFct); - std::vector localVec(nBasisFct); - - for (std::size_t j = 0; j < nBasisFct; ++j) { - const auto global_idx = localIndexSet_.index(j); - localVec[j] = factor * vector[global_idx]; - } - - values.resize(points.size()); - for (std::size_t iq = 0; iq < points.size(); ++iq) { - const auto jacobian = geometry.jacobianInverseTransposed(points[iq].position()); - - localFiniteElem.localBasis().evaluateJacobian(points[iq].position(), shapeGradients); - GradientLocalType grdLocal(R(0)); - for (std::size_t j = 0; j < nBasisFct; ++j) - grdLocal.axpy(localVec[j], shapeGradients[j][0]); - - jacobian.mv(grdLocal, values[iq]); - } - } - - value_type const& operator[](std::size_t iq) const - { - return values[iq]; - } - - int getDegree() const - { - return degree_; - } - - private: - typename DOFVectorType::BaseVector const& vector; - double factor; - - using FeSpace = typename DOFVectorType::FeSpace; - typename FeSpace::LocalView localView_; - typename FeSpace::LocalIndexSet localIndexSet_; - - int degree_ = 0; - - std::vector values; - }; - - - /// \brief Generator function for an expression of gradients of \ref DOFVectors - /** Generates a \ref GradientTerm, that takes a \ref DOFVector \p vector, - * and an optional factor \p factor the value of the gradient of the - * DOFVector at the quadrature points are scaled with. - **/ - template - auto gradientOf(DOFVectorType&& vector, double factor = 1.0) - { - return GradientTerm>{std::forward(vector), factor}; - } - - - /// \brief An expression that evalues a partial derivative of a DOFVector - /// at a given element, either on the DOFs or on the quadrature points. \see derivativeOf. - template - class DerivativeTerm - { - using LocalBasis = typename DOFVectorType::FeSpace::LocalView::Tree::FiniteElement::Traits::LocalBasisType; - using Jacobian = typename LocalBasis::Traits::JacobianType; - - static constexpr int dim = DOFVectorType::FeSpace::GridView::dimension; - static constexpr int dow = DOFVectorType::FeSpace::GridView::dimensionworld; - - using R = typename LocalBasis::Traits::RangeFieldType; - using GradientLocalType = Dune::FieldVector; - using Derivative = R; - - public: - using value_type = Derivative; - using field_type = Value_t; - - DerivativeTerm(DOFVectorType const& dofvector, int component, double factor = 1.0) - : vector(dofvector.getVector()) - , component(component) - , factor(factor) - , localView_(dofvector.getFeSpace().localView()) - , localIndexSet_(dofvector.getFeSpace().localIndexSet()) - {} - - template - void bind(Element const& element) - { - localView_.bind(element); - localIndexSet_.bind(localView_); - int deg = getPolynomialDegree(localView_.tree()); - - auto t = element.type(); - degree_ = t.isSimplex() ? Math::max(0, deg-1) : deg; - } - - template - void init(Element const& element, - PointList const& points) - { - auto geometry = element.geometry(); - const auto& localFiniteElem = localView_.tree().finiteElement(); - const std::size_t nBasisFct = localFiniteElem.size(); - - std::vector shapeGradients(nBasisFct); - std::vector localVec(nBasisFct); - - for (std::size_t j = 0; j < nBasisFct; ++j) { - const auto global_idx = localIndexSet_.index(j); - localVec[j] = factor * vector[global_idx]; - } - - values.resize(points.size()); - for (std::size_t iq = 0; iq < points.size(); ++iq) { - const auto jacobian = geometry.jacobianInverseTransposed(points[iq].position()); - localFiniteElem.localBasis().evaluateJacobian(points[iq].position(), shapeGradients); - - GradientLocalType grdLocal(R(0)); - for (std::size_t j = 0; j < nBasisFct; ++j) - grdLocal.axpy(localVec[j], shapeGradients[j][0]); - - values[iq] = 0; - for (std::size_t j = 0; j < dim; ++j) - values[iq] += jacobian[component][j] * grdLocal[j]; - } - } - - value_type const& operator[](std::size_t iq) const - { - return values[iq]; - } - - int getDegree() const - { - return degree_; - } - - private: - typename DOFVectorType::BaseVector const& vector; - int component; - double factor; - - using FeSpace = typename DOFVectorType::FeSpace; - typename FeSpace::LocalView localView_; - typename FeSpace::LocalIndexSet localIndexSet_; - - int degree_ = 0; - - std::vector values; - }; - - - /// \brief Generator function for an expression of partial derivatives of \ref DOFVectors - /** Generates a \ref DerivativeTerm, that takes a \ref DOFVector \p vector, - * a index representing the component of the gradient that should be returned, - * and an optional factor \p factor the value of the partial derivative of the - * DOFVector at the quadrature points are scaled with. - **/ - template - DerivativeTerm> - derivativeOf(DOFVectorType&& vector, int component, double factor = 1.0) - { - return {std::forward(vector), component, factor}; - } - - #endif - - /** @} **/ - -} // end namespace AMDiS diff --git a/dune/amdis/terms/ExpressionBase.hpp b/dune/amdis/terms/ExpressionBase.hpp deleted file mode 100644 index 08a8a84bbb4b42cef63a359c1fea6a30a5082e90..0000000000000000000000000000000000000000 --- a/dune/amdis/terms/ExpressionBase.hpp +++ /dev/null @@ -1,29 +0,0 @@ -#pragma once - -namespace AMDiS -{ - class ExpressionBase - { - public: - ExpressionBase(int degree = 0) - : degree_(degree) - {} - - template - void bind(Element const&) { /* do nothing */ } - - void unbind() { /* do nothing */ } - - template - void init(Context const& context, PointList const& points) { /* do nothing */ } - - int getDegree() const - { - return degree_; - } - - private: - int degree_ = 0; - }; - -} // end namespace AMDiS diff --git a/dune/amdis/terms/FunctorTerm.hpp b/dune/amdis/terms/FunctorTerm.hpp deleted file mode 100644 index 24b92c8a5d2ab6a159522abb7ff38fe975ef827e..0000000000000000000000000000000000000000 --- a/dune/amdis/terms/FunctorTerm.hpp +++ /dev/null @@ -1,105 +0,0 @@ -#pragma once - -#include -#include -#include - -#include - -#include -#include -#include -#include - -namespace AMDiS -{ - /** - * \addtogroup OperatorTerms - * @{ - **/ - - - /// \brief An expression that evaluates to the current coordinate of a dof or - /// quadrature point with given index. \see eval. - template - class FunctorTerm - : public ExpressionBase - { - public: - using value_type = typename std::result_of...)>::type; - - template - explicit FunctorTerm(Functor const& fct, Ts&&... ts) - : ExpressionBase(2) // TODO: provide functor degree - , fct_(fct) - , terms_(std::forward(ts)...) - {} - - template - void bind(Element const& element) - { - forEach(terms_, [&,this](auto& term) { - term.bind(element); - }); - } - - void unbind() - { - forEach(terms_, [&,this](auto& term) { - term.unbind(); - }); - } - - /// \brief Evaluates the functor \ref fct at all quadrature points and - /// stores the result in a vector \ref values. - template - void init(Context const& context, PointList const& points) - { - forEach(terms_, [&,this](auto& term) { - term.init(context, points); - }); - } - - template - value_type local(Context const& context, LocalCoordinate const& lambda) const - { - return eval([&context,&lambda](auto const& term) { return term.local(context, lambda); }, - MakeSeq_t{}); - } - - template - value_type global(GlobalCoordinate const& x) const - { - return eval([&x](auto const& term) { return term.global(x); }, - MakeSeq_t{}); - } - - auto operator[](std::size_t iq) const - { - return eval([iq](auto const& term) { return term[iq]; }, - MakeSeq_t{}); - } - - private: - template - auto eval(Op op, Indices) const - { - return fct_(op(std::get(terms_))...); - } - - private: - Functor fct_; - std::tuple terms_; - }; - - - /// Generator function for \ref FunctorTerm expressions - template - auto func(Functor const& f, Ts&&... terms) - { - return FunctorTerm...>{f, toTerm(std::forward(terms))...}; - } - - /** @} **/ - -} // end namespace AMDiS diff --git a/dune/amdis/terms/GridViewExpression.hpp b/dune/amdis/terms/GridViewExpression.hpp deleted file mode 100644 index 82041a09d803fea307f6b117069b067aea0d92b1..0000000000000000000000000000000000000000 --- a/dune/amdis/terms/GridViewExpression.hpp +++ /dev/null @@ -1,56 +0,0 @@ -#pragma once - -#include - -#include -#include - -namespace AMDiS -{ - template - class GridViewExpression - { - using Element = typename GridView::template Codim<0>::Entity; - using Domain = typename Element::Geometry::GlobalCoordinate; - using Range = typename Expr::value_type; - - using EntitySet = Dune::Functions::GridViewEntitySet; - - using LocalExpr = LocalExpression; - - public: - GridViewExpression(Expr const& expr, GridView const& gridView) - : expr_(expr) - , entitySet_(gridView) - {} - - Range operator()(Domain const& x) const - { - return expr_.global(x); - } - - friend LocalExpr localFunction(GridViewExpression const& self) - { - return LocalExpr{self.expr_}; - } - - EntitySet const& entitySet() const - { - return entitySet_; - } - - - private: - - Expr expr_; - EntitySet entitySet_; - }; - - template ::value, int> = 0> - auto makeGridViewExpression(Expr const& expr, GridView const& gridView) - { - return GridViewExpression{expr, gridView}; - } - -} // end namespace AMDiS diff --git a/dune/amdis/terms/Interpolation.hpp b/dune/amdis/terms/Interpolation.hpp deleted file mode 100644 index ac4348514633efe167803db053f65448e32db559..0000000000000000000000000000000000000000 --- a/dune/amdis/terms/Interpolation.hpp +++ /dev/null @@ -1,47 +0,0 @@ -#pragma once - -#include - -#include -#include - -#include -#include - -namespace AMDiS -{ - template - auto operator<<(std::shared_ptr> const& dofvector, - PreExpr&& preExpr) - { - using DiscreteGlobalBasisFct = Dune::Functions::DiscreteGlobalBasisFunction; - typename ToTerm::type expr = toTerm(std::forward(preExpr)); - - using Vector = typename DiscreteGlobalBasisFct::Vector; - Vector vec(dofvector->basis(), "tmp"); - - auto gridViewExpr = makeGridViewExpression(expr, dofvector->basis().gridView()); - Dune::Functions::interpolate(dofvector->basis(), dofvector->treePath(), vec, gridViewExpr); - - const_cast(dofvector->dofs()).getVector() = vec.getVector(); - return dofvector; - } - - template - auto& operator<<(Dune::Functions::DiscreteGlobalBasisFunction& dofvector, - PreExpr& preExpr) - { - using DiscreteGlobalBasisFct = Dune::Functions::DiscreteGlobalBasisFunction; - typename ToTerm::type expr = toTerm(std::forward(preExpr)); - - using Vector = typename DiscreteGlobalBasisFct::Vector; - Vector vec(dofvector.basis(), "tmp"); - - auto gridViewExpr = makeGridViewExpression(expr, dofvector.basis().gridView()); - Dune::Functions::interpolate(dofvector.basis(), dofvector.treePath(), vec, gridViewExpr); - - const_cast(dofvector.dofs()).getVector() = vec.getVector(); - return dofvector; - } - -} // end namespace AMDiS diff --git a/dune/amdis/terms/LocalExpression.hpp b/dune/amdis/terms/LocalExpression.hpp deleted file mode 100644 index 61582397ae10fe6f824458bc72cb413e08573a2b..0000000000000000000000000000000000000000 --- a/dune/amdis/terms/LocalExpression.hpp +++ /dev/null @@ -1,78 +0,0 @@ -#pragma once - -#include - -#include -#include - -#include - -namespace AMDiS -{ - template - class LocalExpression - { - using EntitySet = Dune::Functions::GridViewEntitySet; - - using Element = typename EntitySet::Element; - using Geometry = typename Element::Geometry; - - using Domain = typename EntitySet::LocalCoordinate; - using Range = typename Expr::value_type; - - using LocalContext = Element; - using LocalGeometry = Geometry; - using GeometryWrapper = ContextGeometry; // TODO: allow intersections! - - public: - - /// Constructor. Stores a copy to the expression. - explicit LocalExpression(Expr const& expr) - : expr_(expr) - {} - - - /// Evaluate expression `expr_` at local coordinate `lambda` - /// [expects: LocalExpression is bound to element] - Range operator()(Domain const& lambda) const - { - assert( bound_ ); - return expr_.local(context_.value(), lambda); - } - - /// Bind expression to a local context. Must be called before you can evaluate it. - /// [expects: local function is not bound to other element] - void bind(Element const& element) - { - assert( not bound_ ); // NOTE: interpolateTreeSubset(...) does not unbind the element! - geometry_.emplace(element.geometry()); - context_.emplace(element, *geometry_, *geometry_); - - expr_.bind(element); - bound_ = true; - } - - /// Unbind from local context - void unbind() - { - context_.reset(); - geometry_.reset(); - bound_ = false; - } - - /// Obtain local contex this LocalFunction is bound to - LocalContext const& localContext() const - { - assert( bound_ ); - return context_.value().localContext; - } - - private: - Expr expr_; - - Dune::Std::optional context_; - Dune::Std::optional geometry_; - bool bound_ = false; - }; - -} // end namespace AMDiS diff --git a/dune/amdis/terms/TermGenerator.hpp b/dune/amdis/terms/TermGenerator.hpp deleted file mode 100644 index b568b9a9201fb5895584cb1e9e27884136dffffb..0000000000000000000000000000000000000000 --- a/dune/amdis/terms/TermGenerator.hpp +++ /dev/null @@ -1,117 +0,0 @@ -#pragma once - -// std c++ includes -#include - -// AMDiS includes -#include -#include - -#include -#include - -namespace AMDiS -{ - namespace Impl - { - // annonymous namespace to hide tags - namespace - { - struct _constant {}; - struct _functor {}; - struct _default {}; - } - - - template - constexpr bool ConstantCategory = - !std::is_same, tag::unknown>::value; - - -#ifndef AMDIS_DOW - template - constexpr bool FunctorCategory = - Concepts::Callable> || - Concepts::Callable> || - Concepts::Callable>; -#else - template - constexpr bool FunctorCategory = - Concepts::Callable>; -#endif - - - // TODO: add concept for terms - template - using TermCategory = - std::conditional_t< ConstantCategory, _constant, - std::conditional_t< FunctorCategory, _functor, - /* OperatorTerms */ _default > >; - - - // forward declaration - template struct ToTerm; - - template - struct ToTerm - { - using type = ConstantTerm; - - static type get(T const& t) - { - return type{t}; - } - }; - - // if the argument is a functor, that takes a WorldVector - // and returns a double, see concept \ref CoordsFunctor. - template - struct ToTerm - { - using type = CoordsFunctorTerm; - - template - static type get(F_&& f) - { - return type{std::forward(f)}; - } - }; - - // if the argument is already a term return it directly. - template - struct ToTerm - { - using type = T const&; - - static constexpr type get(T const& t) - { - return t; - } - }; - - } // end namespace Impl - - - /// Generator class that takes a type as argument and returns the corresponding term-type. - /** The arguments are transformed into terms in the following way: - * - If type `T` is arithmetic or any small vector/matrix type, generate a - * constant term, see \ref ConstantTerm - * - If the type `T` models a Functor concepts, that maps a WorldVector to a - * double type, construct a \ref CoordsTerm - * - In all other cases, the argument is returned without modification. - **/ - template - using ToTerm = Impl::ToTerm, Impl::TermCategory> >; - - template - using ToTerm_t = std::decay_t::type>; - - - /// Generator function for term transformation, \see ToTerm - template - inline typename ToTerm::type toTerm(T&& t) - { - return ToTerm::get(std::forward(t)); - } - -} // end namespace AMDiS diff --git a/dune/amdis/utility/CMakeLists.txt b/dune/amdis/utility/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..0ed1a7a2e069298344dd8a85eb6bc831edf5fb50 --- /dev/null +++ b/dune/amdis/utility/CMakeLists.txt @@ -0,0 +1,21 @@ +#install headers + + +dune_library_add_sources(amdis SOURCES + Filesystem.cpp +) + +install(FILES + DiscreteFunction.hpp + Filesystem.hpp + FiniteElementType.hpp + GetDegree.hpp + GetEntity.hpp + MultiIndex.hpp + RangeType.hpp + String.hpp + Traversal.hpp + TreeData.hpp + TreePath.hpp + Visitor.hpp +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/amdis/utility) diff --git a/dune/amdis/utility/GetDegree.hpp b/dune/amdis/utility/GetDegree.hpp index e6f05852b4f1df458a916b7f6e6a8fcc160ea635..2e56076d0cf8be2d56754af312cd7facd9a2d6ed 100644 --- a/dune/amdis/utility/GetDegree.hpp +++ b/dune/amdis/utility/GetDegree.hpp @@ -52,7 +52,7 @@ namespace AMDiS } template - int getPolynomialDegree(GridView const& gridView, Node const& node) + int getPolynomialDegree(GridView const& gridView, Node& node) { node.bind(*gridView.template begin<0>()); return Impl::getPolynomialDegreeImpl(node, typename Node::NodeTag{}); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f3a01d98285bd409e61c6ab4dd0be947af8d1f7f..1c3410e644d49916a83015cb942b2f829ab87320 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -5,7 +5,7 @@ foreach(project ${projects2d}) add_executable(${project}2d ${project}.cc) add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 ${project}2d) target_link_dune_default_libraries(${project}2d) - target_link_libraries(${project}2d "duneamdis") + target_link_libraries(${project}2d "amdis") target_compile_definitions(${project}2d PRIVATE AMDIS_DIM=2 AMDIS_DOW=2) endforeach() @@ -15,6 +15,6 @@ foreach(project ${projects3d}) add_executable(${project}3d ${project}.cc) add_dune_alberta_flags(GRIDDIM 3 WORLDDIM 3 ${project}3d) target_link_dune_default_libraries(${project}3d) - target_link_libraries(${project}3d "duneamdis") + target_link_libraries(${project}3d "amdis") target_compile_definitions(${project}3d PRIVATE AMDIS_DIM=3 AMDIS_DOW=3) endforeach() \ No newline at end of file diff --git a/src/ellipt.cc b/src/ellipt.cc index 90480e21aa643bbc4573f17bf1f5915c81d295be..608c0101e36253152ddcfe4a21c087d1dd7353ee 100644 --- a/src/ellipt.cc +++ b/src/ellipt.cc @@ -6,7 +6,6 @@ #include #include #include -#include #include using namespace AMDiS; diff --git a/src/expressions.cc b/src/expressions.cc index 58120101f3f0a384505920adae1358a0862b15aa..856ed1cf90d34f53f70a56079f16b7c591f509d9 100644 --- a/src/expressions.cc +++ b/src/expressions.cc @@ -6,18 +6,10 @@ #include #include #include -#include #include #include -//#include +#include -#include -#include -#include -#include -#include -#include -#include using namespace AMDiS; @@ -33,19 +25,85 @@ int main(int argc, char** argv) ElliptProblem prob("ellipt"); prob.initialize(INIT_ALL); - AdaptInfo adaptInfo("adapt"); + auto u = prob.getSolution(_0); - // prob.getSolution(_0) << [](auto const& x) { return x[0] + x[1]; }; - // prob.getSolution(_0) << func([](double u) { return 2.0 * u; }, valueOf(prob.getSolution(_0))); + // eval a functor at global coordinates (at quadrature points) + auto expr1 = evalAtQP([](Dune::FieldVector const& x) { return x[0] + x[1]; }); + auto expr2 = [](Dune::FieldVector const& x) { return x[0] + x[1]; }; + auto expr3 = [](auto const& x) { return x[0] + x[1]; }; - auto u0 = makeDOFVectorView(*prob.getSolution(), treepath(_0)); - u0.interpolate([](auto const& x) { return std::sin(5*x[0]) + x[1]; }); + // constant values at quadrature points + auto expr4 = 1.0; + auto expr5 = Dune::FieldVector{1.0, 2.0}; + auto expr6 = std::ref(expr4); - u0.interpolate( - eval_([](auto const& x) { return two_norm(x); }) + - apply([](auto const& x) { return two_norm(x); }, X()) + - apply([](auto const& grd) { return two_norm(grd); }, derivative(u0))); + // Coordinate vector and component + auto expr7 = X(); + auto expr8 = X(0); + // Evaluation of the DOFVector (component) at quadrature points + auto expr9 = prob.getSolution(_0); + + // --------------------------------- + + // derivative of expressions + + auto diff4 = gradientAtQP(expr4); + auto diff5 = gradientAtQP(expr5); + auto diff6 = gradientAtQP(expr6); + + auto diff7 = gradientAtQP(expr7); + auto diff8 = gradientAtQP(expr8); + + auto diff9 = gradientAtQP(expr9); + + // --------------------------------- + + u.interpolate(expr1); + u.interpolate(expr2); + u.interpolate(expr3); + + u.interpolate(expr4); + u.interpolate(expr6); + + u.interpolate(expr8); + + u.interpolate(expr9); + + // --------------------------------- + + // operations with expressions + + auto op1 = expr1 + expr2; + auto op2 = expr1 * expr4; + auto op3 = two_norm(expr5); + auto op4 = min(expr6, expr8); + auto op5 = one_norm(expr7); + + auto op6 = invokeAtQP([](double v) { return 2*v; }, op5); + + + u.interpolate(two_norm(diff4)); + u.interpolate(two_norm(diff6)); + + u.interpolate(two_norm(diff8)); + + u.interpolate(two_norm(diff9)); + + // --------------------------------- + + // integration of expressions + + auto gv = u.basis().gridView(); + + auto int1 = integrate(op1, gv); + auto int2 = integrate(op2, gv); + auto int3 = integrate(op3, gv); + auto int4 = integrate(op4, gv); + auto int5 = integrate(op5, gv); + auto int6 = integrate(op6, gv); + + AdaptInfo adaptInfo("adapt", 1); prob.writeFiles(adaptInfo, true); AMDiS::finalize(); return 0; diff --git a/src/heat.cc b/src/heat.cc index e79825d4cb91d9369b28e73c294639cc9a03d996..647ae60dce9cacb1d66f0c9b4fca22e4a27c57a6 100644 --- a/src/heat.cc +++ b/src/heat.cc @@ -8,7 +8,7 @@ #include #include #include -#include +#include #include using namespace AMDiS; @@ -40,7 +40,7 @@ int main(int argc, char** argv) prob.addMatrixOperator(opL, 0, 0); auto opTimeRhs = makeOperator(tag::test{}, - func([invTau](double u) { return u * (*invTau); }, valueOf(prob.getSolution(0))) ); + invokeAtQP([invTau](double u) { return u * (*invTau); }, prob.getSolution(0)) ); prob.addVectorOperator(opTimeRhs, 0); auto opForce = makeOperator(tag::test{}, [](auto const& x) { return -1.0; }); diff --git a/src/stokes1.cc b/src/stokes1.cc index 5dc4e8f5075cf23a8e1be583e5249216e44b3d5d..b65d540972e9f80d5ea171eb64f75166cc18e7dc 100644 --- a/src/stokes1.cc +++ b/src/stokes1.cc @@ -5,7 +5,6 @@ #include #include #include -#include #ifdef DOW #undef DOW diff --git a/src/vecellipt.cc b/src/vecellipt.cc index c3b8041e1f5f480b7a41a2b5db7af8df5cca3c20..7488631a9c6a5c8624b9a81d5a814bfed9c19f79 100644 --- a/src/vecellipt.cc +++ b/src/vecellipt.cc @@ -3,7 +3,6 @@ #include #include #include -#include using namespace AMDiS;