#pragma once #include #include #include #include namespace AMDiS { /** * \addtogroup operators * @{ **/ namespace tag { struct test_partialtrial { std::size_t comp; }; } /// first-order operator \f$\langle\psi, c\,\partial_i\phi\rangle \f$ template class GridFunctionOperator : public GridFunctionOperatorBase { using Super = GridFunctionOperatorBase; static_assert( Category::Scalar, "Expression must be of scalar type." ); public: GridFunctionOperator(tag::test_partialtrial tag, GridFct const& expr, QuadCreator const& quadCreator) : Super(expr, quadCreator, 1) , comp_(tag.comp) {} 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(); 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(); double c = Super::coefficient(local); rowLocalFE.localBasis().evaluateFunction(local, rowShapeValues_); // The gradients of the shape functions on the reference element std::vector > referenceGradients; colLocalFE.localBasis().evaluateJacobian(local, referenceGradients); // Compute the shape function gradients on the real element std::vector colPartial(referenceGradients.size()); for (std::size_t i = 0; i < colPartial.size(); ++i) { colPartial[i] = jacobian[comp_][0] * referenceGradients[i][0][0]; for (std::size_t j = 1; j < jacobian.M(); ++j) colPartial[i] += jacobian[comp_][j] * referenceGradients[i][0][j]; } for (std::size_t j = 0; j < colLocalFE.size(); ++j) { const int local_j = colNode.localIndex(j); double value = factor * (c * colPartial[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::size_t comp_; std::vector> rowShapeValues_; }; /** @} **/ } // end namespace AMDiS