#pragma once #include #include #include #include #include #include namespace AMDiS { /** * \addtogroup operators * @{ **/ /// convection-diffusion operator, see \ref convectionDiffusion /// - + = (conserving) or /// + + = (non conserving) template class ConvectionDiffusionOperator : public LocalOperator, LocalContext> { static const int dow = ContextGeometry::dow; using A_range = typename GridFctA::Range; static_assert( Size_v == 1 || (Rows_v == dow && Cols_v == dow), "Expression A(x) must be of scalar or matrix type." ); using b_range = typename GridFctB::Range; static_assert( Size_v == 1 || Size_v == dow, "Expression b(x) must be of scalar or vector type." ); using c_range = typename GridFctC::Range; static_assert( Size_v == 1, "Expression c(x) must be of scalar type." ); using f_range = typename GridFctF::Range; static_assert( Size_v == 1, "Expression f(x) must be of scalar type." ); public: ConvectionDiffusionOperator(GridFctA const& gridFctA, GridFctB const& gridFctB, GridFctC const& gridFctC, GridFctF const& gridFctF) : gridFctA_(gridFctA) , gridFctB_(gridFctB) , gridFctC_(gridFctC) , gridFctF_(gridFctF) {} template void getElementMatrix(Context const& context, RowNode const& rowNode, ColNode const& colNode, ElementMatrix& elementMatrix) { static_assert(RowNode::isLeaf && ColNode::isLeaf, "Operator can be applied to Leaf-Nodes only." ); static_assert(std::is_same, FiniteElementType_t>{}, "Galerkin operator requires equal finite elements for test and trial space." ); using LocalBasisType = typename FiniteElementType_t::Traits::LocalBasisType; using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType; auto localFctA = localFunction(gridFctA_); localFctA.bind(context.element()); auto localFctB = localFunction(gridFctB_); localFctB.bind(context.element()); auto localFctC = localFunction(gridFctC_); localFctC.bind(context.element()); auto const& localFE = rowNode.finiteElement(); std::size_t numLocalFe = localFE.size(); int quadDeg = std::max({this->getDegree(2,coeffOrder(localFctA),rowNode,colNode), this->getDegree(1,coeffOrder(localFctB),rowNode,colNode), this->getDegree(0,coeffOrder(localFctC),rowNode,colNode)}); using QuadratureRules = Dune::QuadratureRules; auto const& quad = QuadratureRules::rule(context.type(), quadDeg); LocalBasisCache cache(localFE.localBasis()); auto const& shapeGradientsCache = cache.evaluateJacobianAtQp(context, quad); auto const& shapeValuesCache = cache.evaluateFunctionAtQp(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 = context.integrationElement(quad[iq].position()) * quad[iq].weight(); // the values of the shape functions on the reference element at the quadrature point auto const& shapeValues = shapeValuesCache[iq]; // 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; thread_local std::vector gradients; gradients.resize(shapeGradients.size()); for (std::size_t i = 0; i < gradients.size(); ++i) jacobian.mv(shapeGradients[i][0], gradients[i]); const auto A = localFctA(local); WorldVector b; b = localFctB(local); const auto c = localFctC(local); IF_CONSTEXPR(conserving) { WorldVector gradAi; for (std::size_t i = 0; i < numLocalFe; ++i) { const int local_i = rowNode.localIndex(i); gradAi = A * 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; } } } else { WorldVector grad_i; for (std::size_t i = 0; i < numLocalFe; ++i) { const int local_i = rowNode.localIndex(i); grad_i = A * gradients[i]; grad_i+= b * shapeValues[i]; for (std::size_t j = 0; j < numLocalFe; ++j) { const int local_j = colNode.localIndex(j); elementMatrix[local_i][local_j] += (dot(grad_i, gradients[j]) + c * shapeValues[i] * shapeValues[j]) * factor; } } } } localFctA.unbind(); localFctB.unbind(); localFctC.unbind(); } template void getElementVector(Context const& context, Node const& node, ElementVector& elementVector) { static_assert(Node::isLeaf, "Operator can be applied to Leaf-Nodes only." ); using LocalBasisType = typename FiniteElementType_t::Traits::LocalBasisType; auto localFctF = localFunction(gridFctF_); localFctF.bind(context.element()); auto const& localFE = node.finiteElement(); std::size_t numLocalFe = localFE.size(); int quad_order = this->getDegree(0,coeffOrder(localFctF),node); using QuadratureRules = Dune::QuadratureRules; auto const& quad = QuadratureRules::rule(context.type(), quad_order); LocalBasisCache cache(localFE.localBasis()); auto const& shapeValuesCache = cache.evaluateFunctionAtQp(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 values of the shape functions on the reference element at the quadrature point auto const& shapeValues = shapeValuesCache[iq]; // The multiplicative factor in the integral transformation formula const auto factor = localFctF(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight(); for (std::size_t i = 0; i < numLocalFe; ++i) { const int local_i = node.localIndex(i); elementVector[local_i] += shapeValues[i] * factor; } } localFctF.unbind(); } private: template using HasLocalFunctionOrder = decltype( order(std::declval()) ); template int coeffOrder(LocalFct const& localFct) { using Concept = Dune::Std::is_detected; return Dune::Hybrid::ifElse(Concept{}, [&](auto id) { return order(id(localFct)); }, [] (auto) { return 0; }); } private: GridFctA gridFctA_; GridFctB gridFctB_; GridFctC gridFctC_; GridFctF gridFctF_; }; template struct PreConvectionDiffusionOperator { static constexpr bool conserving = c::value; PreGridFctA gridFctA; PreGridFctB gridFctB; PreGridFctC gridFctC; PreGridFctF gridFctF; }; template auto convectionDiffusion(PreGridFctA const& gridFctA, PreGridFctB const& gridFctB, PreGridFctC const& gridFctC, PreGridFctF const& gridFctF, bool_t = {}) { using Pre = PreConvectionDiffusionOperator>; return Pre{gridFctA, gridFctB, gridFctC, gridFctF}; } template auto makeLocalOperator(PreConvectionDiffusionOperator const& pre, GridView const& gridView) { using Pre = PreConvectionDiffusionOperator; auto gridFctA = makeGridFunction(pre.gridFctA, gridView); auto gridFctB = makeGridFunction(pre.gridFctB, gridView); auto gridFctC = makeGridFunction(pre.gridFctC, gridView); auto gridFctF = makeGridFunction(pre.gridFctF, gridView); using GridFctOp = ConvectionDiffusionOperator; GridFctOp localOperator{gridFctA, gridFctB, gridFctC, gridFctF}; return localOperator; } /** @} **/ } // end namespace AMDiS