PythonBoundaryCondition.cpp

Go to the documentation of this file.

 
#include "PythonBoundaryCondition.h"
 
#include <pybind11/pybind11.h>
 
#include <iostream>
 
#include "BaseLib/ConfigTree.h"
#include "FlushStdoutGuard.h"
#include "MeshLib/MeshSearch/NodeSearch.h"
#include "ProcessLib/BoundaryConditionAndSourceTerm/Utils/CreateLocalAssemblers.h"
#include "PythonBoundaryConditionLocalAssembler.h"
 
namespace ProcessLib
{
PythonBoundaryCondition::PythonBoundaryCondition(
    PythonBoundaryConditionData&& bc_data,
    unsigned const integration_order,
    unsigned const shapefunction_order,
    unsigned const global_dim,
    bool const flush_stdout)
    : _bc_data(std::move(bc_data)), _flush_stdout(flush_stdout)
{
    std::vector<MeshLib::Node*> const& bc_nodes =
        _bc_data.boundary_mesh.getNodes();
    MeshLib::MeshSubset bc_mesh_subset(_bc_data.boundary_mesh, bc_nodes);
 
    // Create local DOF table from the bc mesh subset for the given variable and
    // component id.
    _dof_table_boundary = _bc_data.dof_table_bulk.deriveBoundaryConstrainedMap(
        std::move(bc_mesh_subset));
 
    BoundaryConditionAndSourceTerm::createLocalAssemblers<
        PythonBoundaryConditionLocalAssembler>(
        global_dim, _bc_data.boundary_mesh.getElements(), *_dof_table_boundary,
        shapefunction_order, _local_assemblers,
        _bc_data.boundary_mesh.isAxiallySymmetric(), integration_order,
        _bc_data);
}
 
void PythonBoundaryCondition::getEssentialBCValues(
    const double t, GlobalVector const& x,
    NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
{
    FlushStdoutGuard guard(_flush_stdout);
    (void)guard;
 
    auto const nodes = _bc_data.boundary_mesh.getNodes();
 
    auto const& bulk_node_ids_map =
        *_bc_data.boundary_mesh.getProperties().getPropertyVector<std::size_t>(
            "bulk_node_ids", MeshLib::MeshItemType::Node, 1);
 
    bc_values.ids.clear();
    bc_values.values.clear();
 
    bc_values.ids.reserve(_bc_data.boundary_mesh.getNumberOfNodes());
    bc_values.values.reserve(_bc_data.boundary_mesh.getNumberOfNodes());
 
    std::vector<double> primary_variables;
 
    for (auto const* node : _bc_data.boundary_mesh.getNodes())
    {
        auto const boundary_node_id = node->getID();
        auto const bulk_node_id = bulk_node_ids_map[boundary_node_id];
 
        // gather primary variable values
        primary_variables.clear();
        auto const num_var = _dof_table_boundary->getNumberOfVariables();
        for (int var = 0; var < num_var; ++var)
        {
            auto const num_comp =
                _dof_table_boundary->getNumberOfVariableComponents(var);
            for (int comp = 0; comp < num_comp; ++comp)
            {
                MeshLib::Location loc{_bc_data.bulk_mesh_id,
                                      MeshLib::MeshItemType::Node,
                                      bulk_node_id};
                auto const dof_idx =
                    _bc_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.push_back(x[dof_idx]);
            }
        }
 
        auto* xs = nodes[boundary_node_id]->getCoords();  // TODO DDC problems?
        auto pair_flag_value = _bc_data.bc_object->getDirichletBCValue(
            t, {xs[0], xs[1], xs[2]}, boundary_node_id, primary_variables);
        if (!_bc_data.bc_object->isOverriddenEssential())
        {
            DBUG(
                "Method `getDirichletBCValue' not overridden in Python "
                "script.");
            return;
        }
 
        if (!pair_flag_value.first)
        {
            continue;
        }
 
        MeshLib::Location l(_bc_data.bulk_mesh_id, MeshLib::MeshItemType::Node,
                            bulk_node_id);
        const auto dof_idx = _bc_data.dof_table_bulk.getGlobalIndex(
            l, _bc_data.global_component_id);
        if (dof_idx == NumLib::MeshComponentMap::nop)
        {
            OGS_FATAL(
                "Logic error. This error should already have occurred while "
                "gathering primary variables. Something nasty is going on!");
        }
 
        // For the DDC approach (e.g. with PETSc option), the negative
        // index of g_idx means that the entry by that index is a ghost
        // one, which should be dropped. Especially for PETSc routines
        // MatZeroRows and MatZeroRowsColumns, which are called to apply
        // the Dirichlet BC, the negative index is not accepted like
        // other matrix or vector PETSc routines. Therefore, the
        // following if-condition is applied.
        if (dof_idx >= 0)
        {
            bc_values.ids.emplace_back(dof_idx);
            bc_values.values.emplace_back(pair_flag_value.second);
        }
    }
}
 
void PythonBoundaryCondition::applyNaturalBC(
    const double t, std::vector<GlobalVector*> const& x, int const process_id,
    GlobalMatrix& K, GlobalVector& b, GlobalMatrix* Jac)
{
    FlushStdoutGuard guard(_flush_stdout);
 
    try
    {
        GlobalExecutor::executeMemberOnDereferenced(
            &GenericNaturalBoundaryConditionLocalAssemblerInterface::assemble,
            _local_assemblers, *_dof_table_boundary, t, x, process_id, K, b,
            Jac);
    }
    catch (MethodNotOverriddenInDerivedClassException const& /*e*/)
    {
        DBUG("Method `getFlux' not overridden in Python script.");
    }
}
 
std::unique_ptr<PythonBoundaryCondition> createPythonBoundaryCondition(
    BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh,
    NumLib::LocalToGlobalIndexMap const& dof_table, std::size_t bulk_mesh_id,
    int const variable_id, int const component_id,
    unsigned const integration_order, unsigned const shapefunction_order,
    unsigned const global_dim)
{
    config.checkConfigParameter("type", "Python");
 
    auto const bc_object = config.getConfigParameter<std::string>("bc_object");
    auto const flush_stdout = config.getConfigParameter("flush_stdout", false);
 
    // Evaluate Python code in scope of main module
    pybind11::object scope =
        pybind11::module::import("__main__").attr("__dict__");
 
    if (!scope.contains(bc_object))
    {
        OGS_FATAL(
            "Function `{:s}' is not defined in the python script file, or "
            "there was no python script file specified.",
            bc_object);
    }
 
    auto* bc = scope[bc_object.c_str()]
                   .cast<PythonBoundaryConditionPythonSideInterface*>();
 
    if (variable_id >= static_cast<int>(dof_table.getNumberOfVariables()) ||
        component_id >= dof_table.getNumberOfVariableComponents(variable_id))
    {
        OGS_FATAL(
            "Variable id or component id too high. Actual values: ({:d}, "
            "{:d}), maximum values: ({:d}, {:d}).",
            variable_id, component_id, dof_table.getNumberOfVariables(),
            dof_table.getNumberOfVariableComponents(variable_id));
    }
 
    // In case of partitioned mesh the boundary could be empty, i.e. there is no
    // boundary condition.
#ifdef USE_PETSC
    // This can be extracted to createBoundaryCondition() but then the config
    // parameters are not read and will cause an error.
    // TODO (naumov): Add a function to ConfigTree for skipping the tags of the
    // subtree and move the code up in createBoundaryCondition().
    if (boundary_mesh.getDimension() == 0 &&
        boundary_mesh.getNumberOfNodes() == 0 &&
        boundary_mesh.getNumberOfElements() == 0)
    {
        return nullptr;
    }
#endif  // USE_PETSC
 
    return std::make_unique<PythonBoundaryCondition>(
        PythonBoundaryConditionData{
            bc, dof_table, bulk_mesh_id,
            dof_table.getGlobalComponent(variable_id, component_id),
            boundary_mesh},
        integration_order, shapefunction_order, global_dim, flush_stdout);
}
 
}  // namespace ProcessLib