diff --git a/Applications/CLI/CMakeLists.txt b/Applications/CLI/CMakeLists.txt
index f4e4a94213a15deeb51a0f63eb7f6b1a620a471d..74c05b7a45e9a407333d00bd035eefbc70145034 100644
--- a/Applications/CLI/CMakeLists.txt
+++ b/Applications/CLI/CMakeLists.txt
@@ -39,7 +39,7 @@ if(OGS_USE_PYTHON)
target_link_libraries(ogs_embedded_python PUBLIC pybind11::embed)
target_compile_definitions(ogs_embedded_python PUBLIC OGS_USE_PYTHON)
target_link_libraries(ogs_embedded_python PRIVATE
- ProcessLibBoundaryConditionPythonModule)
+ ProcessLibBoundaryConditionPythonModule ProcessLibSourceTermPythonModule)
target_link_libraries(ogs PRIVATE ogs_embedded_python)
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index 61f2ee6746a5e216943913aa368515cdf45e6308..162ba4aacc5185ee79e8b4ad4f3c706984316932 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/ProcessLib/CMakeLists.txt
@@ -32,6 +32,10 @@ if(OGS_USE_PYTHON)
add_subdirectory(BoundaryCondition/Python)
target_link_libraries(ProcessLib
PUBLIC ProcessLibBoundaryConditionPython)
+
+ add_subdirectory(SourceTerms/Python)
+ target_link_libraries(ProcessLib
+ PUBLIC ProcessLibSourceTermPython)
endif()
if(OGS_INSITU)
diff --git a/ProcessLib/SourceTerms/CreateSourceTerm.cpp b/ProcessLib/SourceTerms/CreateSourceTerm.cpp
index 4e50e55f62343ea6eed0a69cb418af8ae9893279..820dff5831b76540cf2ac1b54d0d8a35001fcb9e 100644
--- a/ProcessLib/SourceTerms/CreateSourceTerm.cpp
+++ b/ProcessLib/SourceTerms/CreateSourceTerm.cpp
@@ -11,6 +11,9 @@
#include "CreateNodalSourceTerm.h"
#include "CreateVolumetricSourceTerm.h"
+#ifdef OGS_USE_PYTHON
+#include "Python/CreatePythonSourceTerm.h"
+#endif
#include "SourceTerm.h"
#include "SourceTermConfig.h"
@@ -41,6 +44,18 @@ std::unique_ptr<SourceTerm> createSourceTerm(
*config.component_id);
}
+ if (type == "Python")
+ {
+#ifdef OGS_USE_PYTHON
+ return ProcessLib::createPythonSourceTerm(
+ config.config, config.mesh, dof_table, mesh.getID(), variable_id,
+ *config.component_id, integration_order, shapefunction_order,
+ mesh.getDimension());
+#else
+ OGS_FATAL("OpenGeoSys has not been built with Python support.");
+#endif
+ }
+
OGS_FATAL("Unknown source term type: `%s'.", type.c_str());
}
} // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/Python/CMakeLists.txt b/ProcessLib/SourceTerms/Python/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8ecddaea2e59a5715ef26188d8e27f4795415512
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/CMakeLists.txt
@@ -0,0 +1,32 @@
+add_library(ProcessLibSourceTermPython
+ CreatePythonSourceTerm.cpp
+ CreatePythonSourceTerm.h
+ PythonSourceTerm.cpp
+ PythonSourceTerm.h
+ PythonSourceTermLocalAssembler.h
+ PythonSourceTermPythonSideInterface.h)
+if(BUILD_SHARED_LIBS)
+ install(TARGETS ProcessLibSourceTermPython
+ LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR})
+endif()
+
+target_compile_definitions(ProcessLibSourceTermPython
+ PUBLIC OGS_USE_PYTHON)
+
+target_link_libraries(ProcessLibSourceTermPython
+ PUBLIC BaseLib MathLib MeshLib NumLib logog
+ PRIVATE pybind11::pybind11)
+
+# For the embedded Python module
+add_library(ProcessLibSourceTermPythonModule
+ PythonSourceTermModule.cpp
+ PythonSourceTermModule.h)
+if(BUILD_SHARED_LIBS)
+ install(TARGETS ProcessLibSourceTermPythonModule
+ LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR})
+endif()
+
+target_link_libraries(ProcessLibSourceTermPythonModule
+ PUBLIC
+ ProcessLibSourceTermPython
+ pybind11::pybind11)
diff --git a/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.cpp b/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..491107abba3fc8a1e1e49f7dcd3ef9e62724b896
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.cpp
@@ -0,0 +1,87 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CreatePythonSourceTerm.h"
+
+#include <pybind11/pybind11.h>
+
+#include "BaseLib/ConfigTree.h"
+#include "MeshLib/Mesh.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "ProcessLib/SourceTerms/SourceTerm.h"
+#include "PythonSourceTerm.h"
+
+namespace ProcessLib
+{
+std::unique_ptr<SourceTerm> createPythonSourceTerm(
+ BaseLib::ConfigTree const& config, MeshLib::Mesh const& source_term_mesh,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::size_t const bulk_mesh_id, int const variable_id,
+ int const component_id, unsigned const integration_order,
+ unsigned const shapefunction_order, unsigned const global_dim)
+{
+ DBUG("Constructing PythonSourceTerm from config.");
+ //! \ogs_file_param{prj__process_variables__process_variable__source_term__source_term__type}
+ config.checkConfigParameter("type", "Python");
+
+ auto const source_term_object =
+ //! \ogs_file_param{prj__process_variables__process_variable__source_term__source_term__Python__source_term_object}
+ config.getConfigParameter<std::string>("source_term_object");
+ //! \ogs_file_param{prj__process_variables__process_variable__source_term__source_term__Python__flush_stdout}
+ 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(source_term_object))
+ OGS_FATAL(
+ "Function `%s' is not defined in the python script file, or there "
+ "was no python script file specified.",
+ source_term_object.c_str());
+
+ auto* source_term = scope[source_term_object.c_str()]
+ .cast<ProcessLib::SourceTerms::Python::
+ PythonSourceTermPythonSideInterface*>();
+
+ 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 source_term could be empty, i.e. there is
+ // no source_term condition.
+#ifdef USE_PETSC
+ // This can be extracted to createSourceTerm() 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 createSourceTerm().
+ if (source_term_mesh.getDimension() == 0 &&
+ source_term_mesh.getNumberOfNodes() == 0 &&
+ source_term_mesh.getNumberOfElements() == 0)
+ {
+ return nullptr;
+ }
+#endif // USE_PETSC
+
+ return std::make_unique<ProcessLib::SourceTerms::Python::PythonSourceTerm>(
+ dof_table,
+ ProcessLib::SourceTerms::Python::PythonSourceTermData{
+ source_term, dof_table, bulk_mesh_id,
+ dof_table.getGlobalComponent(variable_id, component_id),
+ source_term_mesh},
+ integration_order, shapefunction_order, global_dim, flush_stdout);
+}
+
+} // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.h b/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.h
new file mode 100644
index 0000000000000000000000000000000000000000..b7b2d6a4b7b328ba235a38ef7f6e12714b428f1a
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.h
@@ -0,0 +1,38 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MeshLib
+{
+class Mesh;
+}
+
+namespace NumLib
+{
+class LocalToGlobalIndexMap;
+}
+
+namespace ProcessLib
+{
+class SourceTerm;
+
+std::unique_ptr<SourceTerm> createPythonSourceTerm(
+ BaseLib::ConfigTree const& config, MeshLib::Mesh const& source_term_mesh,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::size_t const bulk_mesh_id, int const variable_id,
+ int const component_id, unsigned const integration_order,
+ unsigned const shapefunction_order, unsigned const global_dim);
+
+} // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTerm.cpp b/ProcessLib/SourceTerms/Python/PythonSourceTerm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1497c70e53117f1f5a4595d4334b3bd5654ec730
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTerm.cpp
@@ -0,0 +1,98 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "PythonSourceTerm.h"
+
+#include <pybind11/pybind11.h>
+#include <iostream>
+
+#include "MeshLib/MeshSearch/NodeSearch.h"
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+#include "PythonSourceTermLocalAssembler.h"
+
+namespace
+{
+//! Optionally flushes the standard output upon creation and destruction.
+//! Can be used to improve the debug output readability when printing debug
+//! messages both from OGS and from Python.
+class FlushStdoutGuard final
+{
+public:
+ //! Optionally flushes C++ stdout before running Python code.
+ explicit FlushStdoutGuard(bool const flush) : _flush(flush)
+ {
+ if (!flush)
+ return;
+
+ LOGOG_COUT << std::flush;
+ }
+
+ //! Optionally flushes Python's stdout after running Python code.
+ ~FlushStdoutGuard()
+ {
+ if (!_flush)
+ return;
+
+ using namespace pybind11::literals;
+ pybind11::print("end"_a = "", "flush"_a = true);
+ }
+
+private:
+ //! To flush or not to flush.
+ const bool _flush;
+};
+} // anonymous namespace
+
+namespace ProcessLib
+{
+namespace SourceTerms
+{
+namespace Python
+{
+PythonSourceTerm::PythonSourceTerm(
+ NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+ PythonSourceTermData&& source_term_data, unsigned const integration_order,
+ unsigned const shapefunction_order, unsigned const global_dim,
+ bool const flush_stdout)
+ : SourceTerm(source_term_dof_table),
+ _source_term_data(std::move(source_term_data)),
+ _flush_stdout(flush_stdout)
+{
+ std::vector<MeshLib::Node*> const& source_term_nodes =
+ _source_term_data.source_term_mesh.getNodes();
+ MeshLib::MeshSubset source_term_mesh_subset(
+ _source_term_data.source_term_mesh, source_term_nodes);
+
+ // Create local DOF table from the source term mesh subset for the given
+ // variable and component id.
+ _dof_table_source_term =
+ _source_term_data.dof_table_bulk.deriveBoundaryConstrainedMap(
+ std::move(source_term_mesh_subset));
+
+ createLocalAssemblers<PythonSourceTermLocalAssembler>(
+ global_dim, _source_term_data.source_term_mesh.getElements(),
+ *_dof_table_source_term, shapefunction_order, _local_assemblers,
+ _source_term_data.source_term_mesh.isAxiallySymmetric(),
+ integration_order, _source_term_data);
+}
+
+void PythonSourceTerm::integrate(const double t, const GlobalVector& x,
+ GlobalVector& b, GlobalMatrix* Jac) const
+{
+ FlushStdoutGuard guard(_flush_stdout);
+
+ GlobalExecutor::executeMemberOnDereferenced(
+ &PythonSourceTermLocalAssemblerInterface::assemble, _local_assemblers,
+ *_dof_table_source_term, t, x, b, Jac);
+}
+
+} // namespace Python
+} // namespace SourceTerms
+} // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTerm.h b/ProcessLib/SourceTerms/Python/PythonSourceTerm.h
new file mode 100644
index 0000000000000000000000000000000000000000..512c260d93f6ad7f2cbc75bd7a8723436787c110
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTerm.h
@@ -0,0 +1,77 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/IndexValueVector.h"
+#include "ProcessLib/SourceTerms/SourceTerm.h"
+
+#include "PythonSourceTermLocalAssemblerInterface.h"
+#include "PythonSourceTermPythonSideInterface.h"
+
+namespace ProcessLib
+{
+namespace SourceTerms
+{
+namespace Python
+{
+//! Groups data used by source terms, in particular by the local assemblers.
+struct PythonSourceTermData final
+{
+ //! Python object computing source term values.
+ PythonSourceTermPythonSideInterface* source_term_object;
+
+ //! DOF table of the entire domain.
+ NumLib::LocalToGlobalIndexMap const& dof_table_bulk;
+
+ //! Mesh ID of the entire domain.
+ std::size_t const bulk_mesh_id;
+
+ //! Global component ID of the (variable, component) to which this source
+ //! term is applied.
+ int const global_component_id;
+
+ //! The source term mesh, i.e., the (sub-) domain of this source term.
+ const MeshLib::Mesh& source_term_mesh;
+};
+
+//! A source term whose values are computed by a Python script.
+class PythonSourceTerm final : public ProcessLib::SourceTerm
+{
+public:
+ explicit PythonSourceTerm(
+ NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+ PythonSourceTermData&& source_term_data,
+ unsigned const integration_order, unsigned const shapefunction_order,
+ unsigned const global_dim, bool const flush_stdout);
+
+ void integrate(const double t, GlobalVector const& x, GlobalVector& b,
+ GlobalMatrix* jac) const override;
+
+private:
+ //! Auxiliary data.
+ PythonSourceTermData _source_term_data;
+
+ //! Local dof table for the boundary mesh.
+ std::unique_ptr<NumLib::LocalToGlobalIndexMap> _dof_table_source_term;
+
+ //! Local assemblers for all elements of the source term mesh.
+ std::vector<std::unique_ptr<PythonSourceTermLocalAssemblerInterface>>
+ _local_assemblers;
+
+ //! Whether or not to flush standard output before and after each call to
+ //! Python code. Ensures right order of output messages and therefore
+ //! simplifies debugging.
+ bool const _flush_stdout;
+};
+
+} // namespace Python
+} // namespace SourceTerms
+} // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h b/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h
new file mode 100644
index 0000000000000000000000000000000000000000..9a8039ad6ce99ca501a78ba53069e4eb54b05da3
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h
@@ -0,0 +1,234 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "PythonSourceTerm.h"
+
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+#include "PythonSourceTermLocalAssemblerInterface.h"
+#include "PythonSourceTermPythonSideInterface.h"
+
+namespace ProcessLib
+{
+namespace SourceTerms
+{
+namespace Python
+{
+//! Can be thrown to indicate that a member function is not overridden in a
+//! derived class (in particular, if a Python class inherits from a C++ class).
+struct MethodNotOverriddenInDerivedClassException
+{
+};
+
+template <typename NodalRowVectorType>
+struct IntegrationPointData final
+{
+ IntegrationPointData(NodalRowVectorType const& N_,
+ double const& integration_weight_)
+ : N(N_), integration_weight(integration_weight_)
+ {
+ }
+
+ NodalRowVectorType const N;
+ double const integration_weight;
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+class PythonSourceTermLocalAssembler final
+ : public PythonSourceTermLocalAssemblerInterface
+{
+public:
+ using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+ using NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType;
+ using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
+
+ explicit PythonSourceTermLocalAssembler(
+ MeshLib::Element const& e,
+ std::size_t const /*local_matrix_size*/,
+ bool is_axially_symmetric,
+ unsigned const integration_order,
+ PythonSourceTermData const& data)
+ : _data(data),
+ _element(e),
+ _is_axially_symmetric(is_axially_symmetric),
+ _integration_method(integration_order)
+ {
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ auto const shape_matrices =
+ initShapeMatrices<ShapeFunction, ShapeMatricesType,
+ IntegrationMethod, GlobalDim>(
+ _element, _is_axially_symmetric, _integration_method);
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ _ip_data.emplace_back(
+ shape_matrices[ip].N,
+ _integration_method.getWeightedPoint(ip).getWeight() *
+ shape_matrices[ip].integralMeasure *
+ shape_matrices[ip].detJ);
+ }
+ }
+
+ void assemble(std::size_t const source_term_element_id,
+ NumLib::LocalToGlobalIndexMap const& dof_table_source_term,
+ double const t, const GlobalVector& x, GlobalVector& b,
+ GlobalMatrix* Jac) override
+ {
+ using ShapeMatricesType =
+ ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+ using FemType =
+ NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
+ FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
+ &_element));
+
+ unsigned const num_integration_points =
+ _integration_method.getNumberOfPoints();
+ auto const num_var = _data.dof_table_bulk.getNumberOfVariables();
+ auto const num_nodes = ShapeFunction::NPOINTS;
+ auto const num_comp_total =
+ _data.dof_table_bulk.getNumberOfComponents();
+
+ auto const& bulk_node_ids_map =
+ *_data.source_term_mesh.getProperties()
+ .template getPropertyVector<std::size_t>("bulk_node_ids");
+
+ // gather primary variables
+ typename ShapeMatricesType::template MatrixType<ShapeFunction::NPOINTS,
+ Eigen::Dynamic>
+ 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_source_term.getGlobalComponent(var, comp);
+
+ for (unsigned element_node_id = 0; element_node_id < num_nodes;
+ ++element_node_id)
+ {
+ auto const* const node = _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[dof_idx];
+ }
+ }
+ }
+
+ NodalRowVectorType local_rhs = Eigen::VectorXd::Zero(num_nodes);
+ NodalMatrixType 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& ip_data = _ip_data[ip];
+ auto const& N = ip_data.N;
+ auto const& coords = fe.interpolateCoordinates(N);
+ // Assumption: all primary variables have same shape functions.
+ prim_vars = N * primary_variables_mat;
+
+ auto const& w = ip_data.integration_weight;
+ auto const flux_dflux =
+ _data.source_term_object->getFlux(t, coords, prim_vars_data);
+ auto const& flux = flux_dflux.first;
+ auto const& dflux = flux_dflux.second;
+ 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 source term 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 assignement -= takes into account the sign convention
+ // of 1st-order in time ODE systems in OpenGeoSys.
+ local_Jac.block(top, left, width, height).noalias() -=
+ ip_data.N.transpose() * (dflux[comp] * w) * N;
+ }
+ }
+ }
+
+ auto const& indices_specific_component =
+ dof_table_source_term(source_term_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(
+ source_term_element_id, dof_table_source_term);
+ MathLib::RowColumnIndices<GlobalIndexType> rci{
+ indices_specific_component, indices_all_components};
+ Jac->add(rci, local_Jac);
+ }
+ }
+
+private:
+ PythonSourceTermData const& _data;
+ MeshLib::Element const& _element;
+ bool _is_axially_symmetric;
+
+ IntegrationMethod const _integration_method;
+ std::vector<
+ IntegrationPointData<NodalRowVectorType>,
+ Eigen::aligned_allocator<IntegrationPointData<NodalRowVectorType>>>
+ _ip_data;
+};
+
+} // namespace Python
+} // namespace SourceTerms
+} // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssemblerInterface.h b/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssemblerInterface.h
new file mode 100644
index 0000000000000000000000000000000000000000..c006156afab1df39c8bf05df39a5d3ebed8b9151
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssemblerInterface.h
@@ -0,0 +1,33 @@
+
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+namespace ProcessLib
+{
+namespace SourceTerms
+{
+namespace Python
+{
+class PythonSourceTermLocalAssemblerInterface
+{
+public:
+ virtual void assemble(
+ std::size_t const source_term_element_id,
+ NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+ double const t, const GlobalVector& x, GlobalVector& b,
+ GlobalMatrix* Jac) = 0;
+
+ virtual ~PythonSourceTermLocalAssemblerInterface() = default;
+};
+
+} // namespace Python
+} // namespace SourceTerms
+} // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTermModule.cpp b/ProcessLib/SourceTerms/Python/PythonSourceTermModule.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d018837ad6978644cbec25fe0b3e0971b3479236
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTermModule.cpp
@@ -0,0 +1,57 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "PythonSourceTermModule.h"
+
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include "PythonSourceTermPythonSideInterface.h"
+
+namespace ProcessLib
+{
+namespace SourceTerms
+{
+namespace Python
+{
+//! Trampoline class allowing methods of class
+//! PythonSourceTermPythonSideInterface to be overridden on the Python
+//! side. Cf. https://pybind11.readthedocs.io/en/stable/advanced/classes.html
+class PythonSourceTermPythonSideInterfaceTrampoline
+ : public PythonSourceTermPythonSideInterface
+{
+public:
+ using PythonSourceTermPythonSideInterface::
+ PythonSourceTermPythonSideInterface;
+
+ std::pair<double, std::array<double, 3>> getFlux(
+ double t, std::array<double, 3> x,
+ std::vector<double> const& primary_variables) const override
+ {
+ using Ret = std::pair<double, std::array<double, 3>>;
+ PYBIND11_OVERLOAD(Ret, PythonSourceTermPythonSideInterface, getFlux, t,
+ x, primary_variables);
+ }
+};
+
+void pythonBindSourceTerm(pybind11::module& m)
+{
+ namespace py = pybind11;
+
+ py::class_<PythonSourceTermPythonSideInterface,
+ PythonSourceTermPythonSideInterfaceTrampoline>
+ pybc(m, "SourceTerm");
+
+ pybc.def(py::init());
+
+ pybc.def("getFlux", &PythonSourceTermPythonSideInterface::getFlux);
+}
+
+} // namespace Python
+} // namespace SourceTerms
+} // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTermModule.h b/ProcessLib/SourceTerms/Python/PythonSourceTermModule.h
new file mode 100644
index 0000000000000000000000000000000000000000..47b279fc0a34f28c3a400d44c488a6372d36b932
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTermModule.h
@@ -0,0 +1,25 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <pybind11/pybind11.h>
+#include "BaseLib/ExportSymbol.h"
+
+namespace ProcessLib
+{
+namespace SourceTerms
+{
+namespace Python
+{
+//! Creates Python bindings for the Python source term class.
+OGS_EXPORT_SYMBOL void pythonBindSourceTerm(pybind11::module& m);
+} // namespace Python
+} // namespace SourceTerms
+} // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTermPythonSideInterface.h b/ProcessLib/SourceTerms/Python/PythonSourceTermPythonSideInterface.h
new file mode 100644
index 0000000000000000000000000000000000000000..e2bf933785b2091b84eb4406b0d42c64bcb8b769
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTermPythonSideInterface.h
@@ -0,0 +1,51 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+namespace ProcessLib
+{
+namespace SourceTerms
+{
+namespace Python
+{
+//! Base class for source terms.
+//! This class will get Python bindings and is intended to be to be derived in
+//! Python.
+class PythonSourceTermPythonSideInterface
+{
+public:
+ /*!
+ * Computes the flux for the provided arguments (time, position of the node,
+ * primary variables at the node).
+ *
+ * \return flux Flux of the source term at that node and derivative of the
+ * flux w.r.t. all primary variables.
+ */
+ virtual std::pair<double, std::array<double, 3>> getFlux(
+ double /*t*/, std::array<double, 3> /*x*/,
+ std::vector<double> const& /*primary_variables*/) const
+ {
+ return {std::numeric_limits<double>::quiet_NaN(), {}};
+ }
+
+ //! Tells if getFlux() has been overridden in the derived class in Python.
+ //!
+ //! \pre getFlux() must already have been called once.
+ bool isOverriddenGetFlux() const { return _overridden_get_flux; }
+
+ virtual ~PythonSourceTermPythonSideInterface() = default;
+
+private:
+ //! Tells if getFlux() has been overridden in the derived class in Python.
+ mutable bool _overridden_get_flux = true;
+};
+} // namespace Python
+} // namespace SourceTerms
+} // namespace ProcessLib