#pragma once #include #include #include #include namespace AMDiS { /** * \addtogroup operators * @{ **/ namespace tag { struct test_gradtrial {}; } /// first-order operator \f$\langle\psi, \mathbf{b}\cdot\nabla\phi\rangle \f$ template class GridFunctionOperator : public GridFunctionOperatorBase { using Super = GridFunctionOperatorBase; static_assert( Category::Vector, "Expression must be of vector type." ); public: GridFunctionOperator(tag::test_gradtrial, GridFct const& expr, QuadCreator const& quadCreator) : Super(expr, quadCreator, 1) {} template void calculateElementMatrix(Context const& context, QuadratureRule const& quad, ElementMatrix& elementMatrix, RowNode const& rowNode, ColNode const& colNode, std::integral_constant, std::integral_constant) { static_assert(RowNode::isLeaf && ColNode::isLeaf, "Operator can be applied to Leaf-Nodes only."); auto const& rowLocalFE = rowNode.finiteElement(); auto const& colLocalFE = colNode.finiteElement(); std::vector > referenceGradients; 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 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 double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight(); const auto b = Super::coefficient(local); rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_); // The gradients of the shape functions on the reference element colLocalFE.localBasis().evaluateJacobian(local, referenceGradients); // Compute the shape function gradients on the real element std::vector > colGradients(referenceGradients.size()); for (std::size_t i = 0; i < colGradients.size(); ++i) jacobian.mv(referenceGradients[i][0], colGradients[i]); for (std::size_t j = 0; j < colLocalFE.size(); ++j) { const int local_j = colNode.localIndex(j); double value = factor * (b * colGradients[j]); for (std::size_t i = 0; i < rowLocalFE.size(); ++i) { const int local_i = rowNode.localIndex(i); elementMatrix[local_i][local_j] += value * rowShapeValues_[i]; } } } } private: std::vector> rowShapeValues_; }; /** @} **/ } // end namespace AMDiS