diff --git a/src/amdis/Operators.hpp b/src/amdis/Operators.hpp
index 9bfde78d9f97c0338aa836a4eace6961800ac41a..5fba1925cbfaa813aad9e9d7e395ecb8618bd7e9 100644
--- a/src/amdis/Operators.hpp
+++ b/src/amdis/Operators.hpp
@@ -46,6 +46,10 @@
#include <amdis/assembler/ZeroOrderTestvecTrialvec.hpp> // <Psi, A * Phi>, <Psi, c * Phi>
// first-order operators
+#include <amdis/assembler/FirstOrderPartialTest.hpp> // <d_i(psi), c>
+#include <amdis/assembler/FirstOrderGradTest.hpp> // <grad(psi), b>
+#include <amdis/assembler/FirstOrderDivTestvec.hpp> // <div(Psi), c>
+
#include <amdis/assembler/FirstOrderDivTestvecTrial.hpp> // <div(Psi), c * phi>
#include <amdis/assembler/FirstOrderGradTestTrial.hpp> // <grad(psi), b * phi>
#include <amdis/assembler/FirstOrderGradTestTrialvec.hpp> // <grad(psi), c * Phi>
diff --git a/src/amdis/assembler/ConvectionDiffusionOperator.hpp b/src/amdis/assembler/ConvectionDiffusionOperator.hpp
index ca49b35f03048d45555aeae9c4a1f99d57255a27..57d465d32b7f76cf9c2cd8f9fe919b5b027d4e60 100644
--- a/src/amdis/assembler/ConvectionDiffusionOperator.hpp
+++ b/src/amdis/assembler/ConvectionDiffusionOperator.hpp
@@ -108,11 +108,11 @@ namespace AMDiS
const auto c = localFctC(local);
IF_CONSTEXPR(conserving) {
- WorldVector gradAi, gradBi;
+ WorldVector gradAi;
for (std::size_t i = 0; i < numLocalFe; ++i) {
const int local_i = rowNode.localIndex(i);
gradAi = A * gradients[i];
- gradBi = b * gradients[i];
+ auto gradBi = b * gradients[i];
for (std::size_t j = 0; j < numLocalFe; ++j) {
const int local_j = colNode.localIndex(j);
elementMatrix[local_i][local_j] += (dot(gradAi, gradients[j]) + (c * shapeValues[i] - gradBi) * shapeValues[j]) * factor;
diff --git a/src/amdis/assembler/FirstOrderDivTestvec.hpp b/src/amdis/assembler/FirstOrderDivTestvec.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b92ab4e54fb69a5575ff99ec6f8068147fe9bfb8
--- /dev/null
+++ b/src/amdis/assembler/FirstOrderDivTestvec.hpp
@@ -0,0 +1,93 @@
+#pragma once
+
+#include <type_traits>
+
+#include <amdis/GridFunctionOperator.hpp>
+#include <amdis/LocalBasisCache.hpp>
+#include <amdis/Output.hpp>
+#include <amdis/common/Size.hpp>
+
+namespace AMDiS
+{
+ /**
+ * \addtogroup operators
+ * @{
+ **/
+
+ namespace tag
+ {
+ struct divtestvec {};
+ }
+
+
+ /// first-order operator \f$ \langle\nabla\cdot\Psi, c\rangle \f$
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::divtestvec, LocalContext, GridFct>
+ : public GridFunctionOperatorBase<GridFunctionOperator<tag::divtestvec, LocalContext, GridFct>,
+ LocalContext, GridFct>
+ {
+ using Self = GridFunctionOperator;
+ using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+
+ static_assert( Size_v<typename GridFct::Range> == 1, "Expression must be of scalar type." );
+
+ public:
+ GridFunctionOperator(tag::divtestvec, GridFct const& expr)
+ : Super(expr, 1)
+ {}
+
+ template <class Context, class Node, class ElementVector>
+ void getElementVector(Context const& context,
+ Node const& node,
+ ElementVector& elementVector)
+ {
+ static_assert(Node::isPower, "Node must be a Power-Node.");
+
+ static const std::size_t CHILDREN = Node::CHILDREN;
+ static_assert( CHILDREN == Context::dow, "" );
+
+ auto const& quad = this->getQuadratureRule(context.type(), node);
+ auto const& localFE = node.child(0).finiteElement();
+ std::size_t feSize = localFE.size();
+
+ using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
+ using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
+
+ LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
+
+ auto const& shapeGradientsCache = cache.evaluateJacobianAtQp(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.local(quad[iq].position());
+
+ // The transposed inverse Jacobian of the map from the reference element to the element
+ const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+
+ // The multiplicative factor in the integral transformation formula
+ const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
+
+ // The gradients of the shape functions on the reference element
+ auto const& shapeGradients = shapeGradientsCache[iq];
+
+ // Compute the shape function gradients on the real element
+ using WorldVector = FieldVector<RangeFieldType,Context::dow>;
+ thread_local std::vector<WorldVector> gradients;
+ gradients.resize(shapeGradients.size());
+
+ for (std::size_t i = 0; i < gradients.size(); ++i)
+ jacobian.mv(shapeGradients[i][0], gradients[i]);
+
+ for (std::size_t j = 0; j < feSize; ++j) {
+ for (std::size_t k = 0; k < CHILDREN; ++k) {
+ const auto local_kj = node.child(k).localIndex(j);
+ elementVector[local_kj] += factor * gradients[j][k];
+ }
+ }
+ }
+
+ }
+ };
+
+ /** @} **/
+
+} // end namespace AMDiS
diff --git a/src/amdis/assembler/FirstOrderGradTest.hpp b/src/amdis/assembler/FirstOrderGradTest.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cfb1e7adac531aa4a1e01fb0759246f2698870e0
--- /dev/null
+++ b/src/amdis/assembler/FirstOrderGradTest.hpp
@@ -0,0 +1,90 @@
+#pragma once
+
+#include <type_traits>
+
+#include <amdis/GridFunctionOperator.hpp>
+#include <amdis/LocalBasisCache.hpp>
+#include <amdis/Output.hpp>
+#include <amdis/common/Size.hpp>
+
+namespace AMDiS
+{
+ /**
+ * \addtogroup operators
+ * @{
+ **/
+
+ namespace tag
+ {
+ struct gradtest {};
+ }
+
+
+ /// first-order operator \f$ \langle\Psi, c\,\nabla\phi\rangle \f$
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::gradtest, LocalContext, GridFct>
+ : public GridFunctionOperatorBase<GridFunctionOperator<tag::gradtest, LocalContext, GridFct>,
+ LocalContext, GridFct>
+ {
+ static const int dow = ContextGeometry<LocalContext>::dow;
+ using Self = GridFunctionOperator;
+ using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+
+ static_assert( Size_v<typename GridFct::Range> == dow, "Expression must be of vector type." );
+
+ public:
+ GridFunctionOperator(tag::gradtest, GridFct const& expr)
+ : Super(expr, 1)
+ {}
+
+ template <class Context, class Node, class ElementVector>
+ void getElementVector(Context const& context,
+ Node const& node,
+ ElementVector& elementVector)
+ {
+ static_assert(Node::isLeaf, "Node must be Leaf-Node.");
+
+ auto const& quad = this->getQuadratureRule(context.type(), node);
+ auto const& localFE = node.finiteElement();
+ std::size_t feSize = localFE.size();
+
+ using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
+ using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
+
+ LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
+
+ auto const& shapeGradientsCache = cache.evaluateJacobianAtQp(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.local(quad[iq].position());
+
+ // The transposed inverse Jacobian of the map from the reference element to the element
+ const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+
+ // The multiplicative factor in the integral transformation formula
+ const auto factor = Super::coefficient(local);
+ const auto dx = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+
+ // The gradients of the shape functions on the reference element
+ auto const& shapeGradients = shapeGradientsCache[iq];
+
+ // Compute the shape function gradients on the real element
+ using WorldVector = FieldVector<RangeFieldType,Context::dow>;
+ thread_local std::vector<WorldVector> gradients;
+ gradients.resize(shapeGradients.size());
+
+ for (std::size_t i = 0; i < gradients.size(); ++i)
+ jacobian.mv(shapeGradients[i][0], gradients[i]);
+
+ for (std::size_t i = 0; i < feSize; ++i) {
+ const auto local_i = node.localIndex(i);
+ elementVector[local_i] += dx * (factor * gradients[i]);
+ }
+ }
+
+ }
+ };
+
+ /** @} **/
+
+} // end namespace AMDiS
diff --git a/src/amdis/assembler/FirstOrderPartialTest.hpp b/src/amdis/assembler/FirstOrderPartialTest.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5d7e02bf71b51e4d6f01f3dd7c5be17f5bd93f06
--- /dev/null
+++ b/src/amdis/assembler/FirstOrderPartialTest.hpp
@@ -0,0 +1,95 @@
+#pragma once
+
+#include <type_traits>
+
+#include <amdis/assembler/FirstOrderTestPartialTrial.hpp>
+
+namespace AMDiS
+{
+ /**
+ * \addtogroup operators
+ * @{
+ **/
+
+ namespace tag
+ {
+ struct partialtest
+ {
+ std::size_t comp;
+ };
+ }
+
+
+ /// first-order operator \f$ \langle\partial_i\psi, c\rangle \f$
+
+ template <class LocalContext, class GridFct>
+ class GridFunctionOperator<tag::partialtest, LocalContext, GridFct>
+ : public GridFunctionOperatorBase<GridFunctionOperator<tag::partialtest, LocalContext, GridFct>,
+ LocalContext, GridFct>
+ {
+ using Self = GridFunctionOperator;
+ using Super = GridFunctionOperatorBase<Self, LocalContext, GridFct>;
+
+ static_assert( Size_v<typename GridFct::Range> == 1, "Expression must be of scalar type." );
+
+ public:
+ GridFunctionOperator(tag::partialtest tag, GridFct const& expr)
+ : Super(expr, 1)
+ , comp_(tag.comp)
+ {}
+
+ template <class Context, class Node, class ElementVector>
+ void getElementVector(Context const& context,
+ Node const& node,
+ ElementVector& elementVector)
+ {
+ static_assert(Node::isLeaf, "Operator can be applied to Leaf-Nodes only.");
+
+ auto const& quad = this->getQuadratureRule(context.type(), node);
+ auto const& localFE = node.finiteElement();
+ std::size_t feSize = localFE.size();
+
+ using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
+ using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
+
+ LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
+
+ auto const& shapeGradientsCache = cache.evaluateJacobianAtQp(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.local(quad[iq].position());
+
+ // The transposed inverse Jacobian of the map from the reference element to the element
+ const auto jacobian = context.geometry().jacobianInverseTransposed(local);
+ assert(jacobian.M() == Context::dim);
+
+ // The multiplicative factor in the integral transformation formula
+ const auto dx = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+ const auto c = Super::coefficient(local);
+
+ // The gradients of the shape functions on the reference element
+ auto const& shapeGradients = shapeGradientsCache[iq];
+
+ // Compute the shape function gradients on the real element
+ thread_local std::vector<RangeFieldType> partial;
+ partial.resize(shapeGradients.size());
+ for (std::size_t i = 0; i < partial.size(); ++i) {
+ partial[i] = jacobian[comp_][0] * shapeGradients[i][0][0];
+ for (std::size_t j = 1; j < jacobian.M(); ++j)
+ partial[i] += jacobian[comp_][j] * shapeGradients[i][0][j];
+ }
+
+ for (std::size_t j = 0; j < feSize; ++j) {
+ const auto local_j = node.localIndex(j);
+ elementVector[local_j] += dx * (c * partial[j]);
+ }
+ }
+ }
+
+ private:
+ std::size_t comp_;
+ };
+
+ /** @} **/
+
+} // end namespace AMDiS
diff --git a/src/amdis/assembler/SecondOrderPartialTestPartialTrial.hpp b/src/amdis/assembler/SecondOrderPartialTestPartialTrial.hpp
index 434a9b539cf5096f5fda71d7413f81a0a6be85f0..bb558164985986bc3863bff493c2d61298615dd3 100644
--- a/src/amdis/assembler/SecondOrderPartialTestPartialTrial.hpp
+++ b/src/amdis/assembler/SecondOrderPartialTestPartialTrial.hpp
@@ -88,7 +88,7 @@ namespace AMDiS
rowPartial.resize(rowShapeGradients.size());
for (std::size_t i = 0; i < rowPartial.size(); ++i) {
rowPartial[i] = jacobian[compTest_][0] * rowShapeGradients[i][0][0];
- for (std::size_t j = 1; j < jacobian.cols(); ++j)
+ for (std::size_t j = 1; j < jacobian.cols; ++j)
rowPartial[i] += jacobian[compTest_][j] * rowShapeGradients[i][0][j];
}
@@ -96,7 +96,7 @@ namespace AMDiS
colPartial.resize(colShapeGradients.size());
for (std::size_t i = 0; i < colPartial.size(); ++i) {
colPartial[i] = jacobian[compTrial_][0] * colShapeGradients[i][0][0];
- for (std::size_t j = 1; j < jacobian.cols(); ++j)
+ for (std::size_t j = 1; j < jacobian.cols; ++j)
colPartial[i] += jacobian[compTrial_][j] * colShapeGradients[i][0][j];
}
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 57754d5dc8a495c3869ad293c8482b29928803de..f7819bc79529290f40abf4a17a9b798a99d0f011 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -40,6 +40,9 @@ dune_add_test(SOURCES MultiTypeVectorTest.cpp
dune_add_test(SOURCES MultiTypeMatrixTest.cpp
LINK_LIBRARIES amdis)
+dune_add_test(SOURCES OperatorsTest.cpp
+ LINK_LIBRARIES amdis)
+
dune_add_test(SOURCES RangeTypeTest.cpp
LINK_LIBRARIES amdis)
diff --git a/test/OperatorsTest.cpp b/test/OperatorsTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..259eaddd1870b2a470c14d61aa1ea6986f695867
--- /dev/null
+++ b/test/OperatorsTest.cpp
@@ -0,0 +1,116 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include "config.h"
+#include <amdis/AMDiS.hpp>
+#include <amdis/ProblemStat.hpp>
+#include <amdis/Operators.hpp>
+#include <amdis/assembler/ConvectionDiffusionOperator.hpp>
+#include <amdis/assembler/StokesOperator.hpp>
+
+
+using namespace AMDiS;
+
+
+using Grid = Dune::YaspGrid<2>;
+using Param = TaylorHoodBasis<typename Grid::LeafGridView>;
+using Problem = ProblemStat<Param>;
+
+int main(int argc, char** argv)
+{
+ AMDiS::init(argc, argv);
+
+ using namespace Dune::Indices;
+
+ Problem prob("ellipt");
+ prob.initialize(INIT_ALL);
+
+ auto _v = _0;
+ auto _p = _1;
+
+ FieldVector<double,1> scalar = 1.0;
+ FieldVector<double,2> vec; vec = 1.0;
+ FieldMatrix<double,2,2> mat; mat = 1.0;
+
+ auto op1 = makeOperator(tag::divtestvec{}, scalar);
+ prob.addVectorOperator(op1, _v);
+
+ auto op2 = makeOperator(tag::divtestvec_trial{}, scalar);
+ prob.addMatrixOperator(op2, _v, _p);
+
+ auto op3 = makeOperator(tag::gradtest{}, vec);
+ prob.addVectorOperator(op3, _p);
+
+ auto op4 = makeOperator(tag::gradtest_trial{}, vec);
+ prob.addMatrixOperator(op4, _p, _p);
+
+ auto op5a = makeOperator(tag::gradtest_trialvec{}, scalar);
+ // auto op5b = makeOperator(tag::gradtest_trialvec{}, mat); // TODO: needs to be implemented
+ prob.addMatrixOperator(op5a, _p, _v);
+
+ auto op6 = makeOperator(tag::partialtest{0}, scalar);
+ prob.addVectorOperator(op6, _p);
+
+ auto op7 = makeOperator(tag::partialtest_trial{0}, scalar);
+ prob.addMatrixOperator(op7, _p, _p);
+
+ auto op8 = makeOperator(tag::test_divtrialvec{}, scalar);
+ prob.addMatrixOperator(op8, _p, _v);
+
+ auto op9 = makeOperator(tag::test_gradtrial{}, vec);
+ prob.addMatrixOperator(op9, _p, _p);
+
+ auto op10 = makeOperator(tag::test_partialtrial{0}, scalar);
+ prob.addMatrixOperator(op10, _p, _p);
+
+ auto op11a = makeOperator(tag::testvec_gradtrial{}, scalar);
+ // auto op11b = makeOperator(tag::testvec_gradtrial{}, mat); // TODO: needs to be implemented
+ prob.addMatrixOperator(op11a, _v, _p);
+
+ auto op12 = makeOperator(tag::divtestvec_divtrialvec{}, scalar);
+ prob.addMatrixOperator(op12, _v, _v);
+
+ auto op13a = makeOperator(tag::gradtest_gradtrial{}, scalar);
+ auto op13b = makeOperator(tag::gradtest_gradtrial{}, mat);
+ prob.addMatrixOperator(op13a, _p, _p);
+ prob.addMatrixOperator(op13b, _p, _p);
+
+ auto op14 = makeOperator(tag::partialtest_partialtrial{0,0}, scalar);
+ prob.addMatrixOperator(op14, _p, _p);
+
+ auto op15 = makeOperator(tag::test{}, scalar);
+ prob.addVectorOperator(op15, _p);
+
+ auto op16 = makeOperator(tag::test_trial{}, scalar);
+ prob.addMatrixOperator(op16, _p, _p);
+
+ auto op17 = makeOperator(tag::test_trialvec{}, vec);
+ prob.addMatrixOperator(op17, _p, _v);
+
+ auto op18 = makeOperator(tag::testvec{}, vec);
+ prob.addVectorOperator(op18, _v);
+
+ auto op19 = makeOperator(tag::testvec_trial{}, vec);
+ prob.addMatrixOperator(op19, _v, _p);
+
+ auto op20a = makeOperator(tag::testvec_trialvec{}, scalar);
+ auto op20b = makeOperator(tag::testvec_trialvec{}, mat);
+ prob.addMatrixOperator(op20a, _v, _v);
+ prob.addMatrixOperator(op20b, _v, _v);
+
+ // some special operators
+
+ auto opStokes = makeOperator(tag::stokes{}, scalar);
+ prob.addMatrixOperator(opStokes);
+ prob.addMatrixOperator(opStokes, {}, {});
+
+ auto opCDa = convectionDiffusion(scalar, scalar, scalar, scalar);
+ prob.addMatrixOperator(opCDa, _p, _p);
+ prob.addVectorOperator(opCDa, _p);
+
+ auto opCDb = convectionDiffusion(mat, vec, scalar, scalar, std::false_type{});
+ prob.addMatrixOperator(opCDb, _p, _p);
+ prob.addVectorOperator(opCDb, _p);
+
+ AMDiS::finalize();
+ return 0;
+}