PythonBoundaryConditionLocalAssembler.h

Go to the documentation of this file.

 
#pragma once
 
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
#include "NumLib/DOF/DOFTableUtil.h"
#include "ProcessLib/BoundaryConditionAndSourceTerm/GenericNaturalBoundaryConditionLocalAssembler.h"
#include "PythonBoundaryCondition.h"
#include "PythonBoundaryConditionPythonSideInterface.h"
 
namespace ProcessLib
{
struct MethodNotOverriddenInDerivedClassException
{
};
 
template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
class PythonBoundaryConditionLocalAssembler final
    : public GenericNaturalBoundaryConditionLocalAssembler<
          ShapeFunction, IntegrationMethod, GlobalDim>
{
    using Base = GenericNaturalBoundaryConditionLocalAssembler<
        ShapeFunction, IntegrationMethod, GlobalDim>;
 
public:
    PythonBoundaryConditionLocalAssembler(
        MeshLib::Element const& e,
        std::size_t const /*local_matrix_size*/,
        bool is_axially_symmetric,
        unsigned const integration_order,
        PythonBoundaryConditionData const& data)
        : Base(e, is_axially_symmetric, integration_order), _data(data)
    {
    }
 
    void assemble(std::size_t const boundary_element_id,
                  NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
                  double const t, std::vector<GlobalVector*> const& x,
                  int const process_id, GlobalMatrix& /*K*/, GlobalVector& b,
                  GlobalMatrix* Jac) override
    {
        using ShapeMatricesType =
            ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
        auto const fe = NumLib::createIsoparametricFiniteElement<
            ShapeFunction, ShapeMatricesType>(Base::_element);
 
        unsigned const num_integration_points =
            Base::_integration_method.getNumberOfPoints();
        auto const num_var = _data.dof_table_bulk.getNumberOfVariables();
        auto const num_nodes = Base::_element.getNumberOfNodes();
        auto const num_comp_total =
            _data.dof_table_bulk.getNumberOfGlobalComponents();
 
        auto const& bulk_node_ids_map =
            *_data.boundary_mesh.getProperties()
                 .template getPropertyVector<std::size_t>(
                     "bulk_node_ids", MeshLib::MeshItemType::Node, 1);
 
        // gather primary variables
        Eigen::MatrixXd primary_variables_mat(num_nodes, num_comp_total);
        for (int var = 0; var < num_var; ++var)
        {
            auto const num_comp =
                _data.dof_table_bulk.getNumberOfVariableComponents(var);
            for (int comp = 0; comp < num_comp; ++comp)
            {
                auto const global_component =
                    dof_table_boundary.getGlobalComponent(var, comp);
 
                for (unsigned element_node_id = 0; element_node_id < num_nodes;
                     ++element_node_id)
                {
                    auto const* const node =
                        Base::_element.getNode(element_node_id);
                    auto const boundary_node_id = node->getID();
                    auto const bulk_node_id =
                        bulk_node_ids_map[boundary_node_id];
                    MeshLib::Location loc{_data.bulk_mesh_id,
                                          MeshLib::MeshItemType::Node,
                                          bulk_node_id};
                    auto const dof_idx =
                        _data.dof_table_bulk.getGlobalIndex(loc, var, comp);
                    if (dof_idx == NumLib::MeshComponentMap::nop)
                    {
                        // TODO extend Python BC to mixed FEM ansatz functions
                        OGS_FATAL(
                            "No d.o.f. found for (node={:d}, var={:d}, "
                            "comp={:d}).  "
                            "That might be due to the use of mixed FEM ansatz "
                            "functions, which is currently not supported by "
                            "the implementation of Python BCs. That excludes, "
                            "e.g., the HM process.",
                            bulk_node_id, var, comp);
                    }
                    primary_variables_mat(element_node_id, global_component) =
                        (*x[process_id])[dof_idx];
                }
            }
        }
 
        Eigen::VectorXd local_rhs = Eigen::VectorXd::Zero(num_nodes);
        Eigen::MatrixXd local_Jac =
            Eigen::MatrixXd::Zero(num_nodes, num_nodes * num_comp_total);
        std::vector<double> prim_vars_data(num_comp_total);
        auto prim_vars = MathLib::toVector(prim_vars_data);
 
        for (unsigned ip = 0; ip < num_integration_points; ip++)
        {
            auto const& N = Base::_ns_and_weights[ip].N;
            auto const& w = Base::_ns_and_weights[ip].weight;
            auto const coords = fe.interpolateCoordinates(N);
            prim_vars =
                N * primary_variables_mat;  // Assumption: all primary variables
                                            // have same shape functions.
            auto const flag_flux_dFlux =
                _data.bc_object->getFlux(t, coords, prim_vars_data);
            if (!_data.bc_object->isOverriddenNatural())
            {
                // getFlux() is not overridden in Python, so we can skip the
                // whole BC assembly (i.e., for all boundary elements).
                throw MethodNotOverriddenInDerivedClassException{};
            }
 
            if (!std::get<0>(flag_flux_dFlux))
            {
                // No flux value for this integration point. Skip assembly of
                // the entire element.
                return;
            }
            auto const flux = std::get<1>(flag_flux_dFlux);
            auto const& dFlux = std::get<2>(flag_flux_dFlux);
 
            local_rhs.noalias() += N * (flux * w);
 
            if (static_cast<int>(dFlux.size()) != num_comp_total)
            {
                // This strict check is technically mandatory only if a Jacobian
                // is assembled. However, it is done as a consistency check also
                // for cases without Jacobian assembly.
                OGS_FATAL(
                    "The Python BC must return the derivative of the flux "
                    "w.r.t. each primary variable. {:d} components expected. "
                    "{:d} components returned from Python.",
                    num_comp_total, dFlux.size());
            }
 
            if (Jac)
            {
                for (int comp = 0; comp < num_comp_total; ++comp)
                {
                    auto const top = 0;
                    auto const left = comp * num_nodes;
                    auto const width = num_nodes;
                    auto const height = num_nodes;
                    // The assignment -= takes into account the sign convention
                    // of 1st-order in time ODE systems in OpenGeoSys.
                    local_Jac.block(top, left, width, height).noalias() -=
                        N.transpose() * (dFlux[comp] * w) * N;
                }
            }
        }
 
        auto const& indices_specific_component =
            dof_table_boundary(boundary_element_id, _data.global_component_id)
                .rows;
        b.add(indices_specific_component, local_rhs);
 
        if (Jac)
        {
            // only assemble a block of the Jacobian, not the whole local matrix
            auto const indices_all_components =
                NumLib::getIndices(boundary_element_id, dof_table_boundary);
            MathLib::RowColumnIndices<GlobalIndexType> rci{
                indices_specific_component, indices_all_components};
            Jac->add(rci, local_Jac);
        }
    }
 
private:
    PythonBoundaryConditionData const& _data;
};
 
}  // namespace ProcessLib