diff --git a/.gitmodules b/.gitmodules
index c51df9e2766afb892e5f8fd513617c44eafa9704..e22e17ff10b1e6b1d8b9880c78b2d9b2086644b7 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -22,3 +22,6 @@
[submodule "ThirdParty/cmake-modules"]
path = ThirdParty/cmake-modules
url = https://github.com/bilke/cmake-modules.git
+[submodule "ThirdParty/pybind11"]
+ path = ThirdParty/pybind11
+ url = https://github.com/pybind/pybind11.git
diff --git a/Applications/ApplicationsLib/CMakeLists.txt b/Applications/ApplicationsLib/CMakeLists.txt
index 701c66a4c58918e3921934ed248a42f9ba381ba7..9d734bcb839430e8e13bc73f572b84b1cf756a79 100644
--- a/Applications/ApplicationsLib/CMakeLists.txt
+++ b/Applications/ApplicationsLib/CMakeLists.txt
@@ -22,3 +22,7 @@ endforeach()
if(OGS_USE_PCH)
cotire(ApplicationsLib)
endif()
+
+if(OGS_USE_PYTHON)
+ target_link_libraries(ApplicationsLib PRIVATE pybind11::pybind11)
+endif()
diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 58b077d4a0af2460de993c2f215b890aaed47997..d3d07d07bd6c47cdf07a5820f03a169a120152a5 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -18,6 +18,10 @@
#include <logog/include/logog.hpp>
+#ifdef OGS_USE_PYTHON
+#include <pybind11/eval.h>
+#endif
+
#include "BaseLib/Algorithm.h"
#include "BaseLib/FileTools.h"
@@ -162,6 +166,23 @@ ProjectData::ProjectData(BaseLib::ConfigTree const& project_config,
{
_mesh_vec = readMeshes(project_config, project_directory);
+ if (auto const python_script =
+ //! \ogs_file_param{prj__python_script}
+ project_config.getConfigParameterOptional<std::string>("python_script"))
+ {
+#ifdef OGS_USE_PYTHON
+ namespace py = pybind11;
+ auto const script_path =
+ BaseLib::copyPathToFileName(*python_script, project_directory);
+
+ // Evaluate in scope of main module
+ py::object scope = py::module::import("__main__").attr("__dict__");
+ py::eval_file(script_path, scope);
+#else
+ OGS_FATAL("OpenGeoSys has not been built with Python support.");
+#endif // OGS_USE_PYTHON
+ }
+
//! \ogs_file_param{prj__curves}
parseCurves(project_config.getConfigSubtreeOptional("curves"));
diff --git a/Applications/CLI/CMakeLists.txt b/Applications/CLI/CMakeLists.txt
index c2410af0547a690a67d35040746b99cf947c35b4..f15fcb737427d8013770e431caf2a3db3901a992 100644
--- a/Applications/CLI/CMakeLists.txt
+++ b/Applications/CLI/CMakeLists.txt
@@ -5,6 +5,51 @@ target_link_libraries(ogs
PRIVATE BaseLib ApplicationsLib
)
+if(OGS_USE_PYTHON)
+ # Troubleshooting:
+ # If you get linker errors, such as ogs.cpp:(.text+0xb4): undefined reference to `_Py_ZeroStruct'
+ # it could be that OGS is compiled with the wrong Python version.
+ # I (Ch. Leh.) observed the following: The symbol _Py_ZeroStruct could not be found in /usr/lib/libpython3.6m.so (I intended to compile OGS with Python3).
+ # It's apparently a Python2 symbol (present in /usr/lib/libpython2.7.so)
+ # The compiler command-line was the following:
+ # /usr/bin/g++ ... -DvtkRenderingVolume_AUTOINIT="1(vtkRenderingVolumeOpenGL2)" \
+ # -I/usr/include/vtk -I/usr/include/python2.7 -I/usr/include/freetype2 \
+ # -I/usr/include/libxml2 ... -I/.../BaseLib ... \
+ # -isystem /usr/include/python3.6m ... -o CMakeFiles/ogs.dir/ogs.cpp.o \
+ # -c /.../Applications/CLI/ogs.cpp
+ # In particular, the Python2 include path comes before the Python3 include path.
+ # Compiling OGS with Python2 solved the issue.
+ # I assume (this is only a guess!) that VTK pulls in Python2 dependencies (on my system).
+ # I assume that this is related to https://github.com/ufz/ogs/pull/2158.
+ # Workaround: Always make sure that OGS is compiled with the same Python version as VTK.
+ # The error described above should be detected automatically by cmake and an
+ # appropriate message should be presented. The note is kept for the case
+ # that the automatic detection does not work due to whatever reason.
+
+ add_library(ogs_embedded_python STATIC ogs_embedded_python.cpp)
+
+ # Performance warning from
+ # https://github.com/pybind/pybind11/blob/master/docs/compiling.rst:
+ # Since pybind11 is a metatemplate library, it is crucial that certain compiler
+ # flags are provided to ensure high quality code generation. In contrast to the
+ # pybind11_add_module() command, the CMake interface library only provides the
+ # minimal set of parameters to ensure that the code using pybind11 compiles, but
+ # it does not pass these extra compiler flags (i.e. this is up to you).
+ # TODO: Enable further compiler/linker flags.
+ 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)
+
+ target_link_libraries(ogs PRIVATE ogs_embedded_python)
+
+ if(BUILD_SHARED_LIBS)
+ # Add macro definition, because static libs make special handling necessary
+ # s.t. the embedded OpenGeoSys Python module won't be removed by the linker.
+ target_compile_definitions(ogs_embedded_python OGS_BUILD_SHARED_LIBS)
+ endif()
+endif()
+
if(OGS_USE_PETSC)
target_link_libraries(ogs PRIVATE ${PETSC_LIBRARIES})
endif()
diff --git a/Applications/CLI/ogs.cpp b/Applications/CLI/ogs.cpp
index 1a7dbc096a267bedfe71c3135d74e77dfddc1e36..87694bcc50aa67da5077c77b5e3b087df8213583 100644
--- a/Applications/CLI/ogs.cpp
+++ b/Applications/CLI/ogs.cpp
@@ -10,8 +10,8 @@
*
*/
-#include <chrono>
#include <tclap/CmdLine.h>
+#include <chrono>
#ifndef _WIN32
#ifdef __APPLE__
@@ -42,7 +42,11 @@
#include "NumLib/NumericsConfig.h"
-int main(int argc, char *argv[])
+#ifdef OGS_USE_PYTHON
+#include "ogs_embedded_python.h"
+#endif
+
+int main(int argc, char* argv[])
{
// Parse CLI arguments.
TCLAP::CmdLine cmd(
@@ -65,34 +69,32 @@ int main(int argc, char *argv[])
"PROJECT FILE");
cmd.add(project_arg);
- TCLAP::ValueArg<std::string> outdir_arg(
- "o", "output-directory",
- "the output directory to write to",
- false,
- "",
- "output directory");
+ TCLAP::ValueArg<std::string> outdir_arg("o", "output-directory",
+ "the output directory to write to",
+ false, "", "output directory");
cmd.add(outdir_arg);
- TCLAP::ValueArg<std::string> log_level_arg(
- "l", "log-level",
- "the verbosity of logging messages: none, error, warn, info, debug, all",
- false,
+ TCLAP::ValueArg<std::string> log_level_arg("l", "log-level",
+ "the verbosity of logging "
+ "messages: none, error, warn, "
+ "info, debug, all",
+ false,
#ifdef NDEBUG
- "info",
+ "info",
#else
- "all",
+ "all",
#endif
- "log level");
+ "log level");
cmd.add(log_level_arg);
TCLAP::SwitchArg nonfatal_arg("",
- "config-warnings-nonfatal",
- "warnings from parsing the configuration file will not trigger program abortion");
+ "config-warnings-nonfatal",
+ "warnings from parsing the configuration "
+ "file will not trigger program abortion");
cmd.add(nonfatal_arg);
- TCLAP::SwitchArg unbuffered_cout_arg("",
- "unbuffered-std-out",
- "use unbuffered standard output");
+ TCLAP::SwitchArg unbuffered_cout_arg("", "unbuffered-std-out",
+ "use unbuffered standard output");
cmd.add(unbuffered_cout_arg);
#ifndef _WIN32 // TODO: On windows floating point exceptions are not handled
@@ -124,6 +126,11 @@ int main(int argc, char *argv[])
#endif // __APPLE__
#endif // _WIN32
+#ifdef OGS_USE_PYTHON
+ pybind11::scoped_interpreter guard = ApplicationsLib::setupEmbeddedPython();
+ (void)guard;
+#endif
+
BaseLib::RunTime run_time;
{
@@ -167,7 +174,9 @@ int main(int argc, char *argv[])
if (auto t = project_config->getConfigSubtreeOptional("insitu"))
{
//! \ogs_file_param{prj__insitu__scripts}
- InSituLib::Initialize(t->getConfigSubtree("scripts"), BaseLib::extractPath(project_arg.getValue()));
+ InSituLib::Initialize(
+ t->getConfigSubtree("scripts"),
+ BaseLib::extractPath(project_arg.getValue()));
isInsituConfigured = true;
}
#else
@@ -196,11 +205,11 @@ int main(int argc, char *argv[])
#ifdef USE_INSITU
if (isInsituConfigured)
InSituLib::Finalize();
- #endif
+#endif
INFO("[time] Execution took %g s.", run_time.elapsed());
#if defined(USE_PETSC)
- controller->Finalize(1) ;
+ controller->Finalize(1);
#endif
} // This nested scope ensures that everything that could possibly
// possess a ConfigTree is destructed before the final check below is
@@ -209,7 +218,9 @@ int main(int argc, char *argv[])
BaseLib::ConfigTree::assertNoSwallowedErrors();
ogs_status = solver_succeeded ? EXIT_SUCCESS : EXIT_FAILURE;
- } catch (std::exception& e) {
+ }
+ catch (std::exception& e)
+ {
ERR(e.what());
ogs_status = EXIT_FAILURE;
}
diff --git a/Applications/CLI/ogs_embedded_python.cpp b/Applications/CLI/ogs_embedded_python.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2da7a23d1a518b6df3296b914b195bdc6169b530
--- /dev/null
+++ b/Applications/CLI/ogs_embedded_python.cpp
@@ -0,0 +1,52 @@
+/**
+ * \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 <pybind11/embed.h>
+#include <logog/include/logog.hpp>
+
+#include "ogs_embedded_python.h"
+
+#include "ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionModule.h"
+
+PYBIND11_EMBEDDED_MODULE(OpenGeoSys, m)
+{
+ DBUG("Binding Python module OpenGeoSys.");
+
+ ProcessLib::pythonBindBoundaryCondition(m);
+}
+
+#ifndef OGS_BUILD_SHARED_LIBS
+
+// Hackish trick that hopefully ensures that the linker won't strip the symbol
+// pointed to by p from the library being built.
+template <typename T>
+void mark_used(T p)
+{
+ volatile T vp = (volatile T)p;
+ vp = vp;
+}
+
+#endif // OGS_BUILD_SHARED_LIBS
+
+namespace ApplicationsLib
+{
+pybind11::scoped_interpreter setupEmbeddedPython()
+{
+#ifndef OGS_BUILD_SHARED_LIBS
+ // pybind11_init_impl_OpenGeoSys is the function initializing the embedded
+ // OpenGeoSys Python module. The name is generated by pybind11. If it is not
+ // obvious that this symbol is actually used, the linker might remove it
+ // under certain circumstances.
+ mark_used(&pybind11_init_impl_OpenGeoSys);
+#endif
+
+ return pybind11::scoped_interpreter{};
+}
+
+} // namespace ApplicationsLib
diff --git a/Applications/CLI/ogs_embedded_python.h b/Applications/CLI/ogs_embedded_python.h
new file mode 100644
index 0000000000000000000000000000000000000000..758718298ba1cb3e00d47338f58901d01ea82074
--- /dev/null
+++ b/Applications/CLI/ogs_embedded_python.h
@@ -0,0 +1,20 @@
+/**
+ * \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/embed.h>
+
+namespace ApplicationsLib
+{
+/// Sets up an embedded Python interpreter and makes sure that the OpenGeoSys
+/// Python module is not removed by the linker.
+pybind11::scoped_interpreter setupEmbeddedPython();
+
+} // namespace ApplicationsLib
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3d87e8ba12738e631d4dfbb34c53744ff2ab4dc4..add4d4668f48f6f2453db13e5984dc560592ff0e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -210,6 +210,8 @@ option(OGS_ENABLE_ELEMENT_PYRAMID "Build FEM elements for pyramids." ON)
option(OGS_CHECK_HEADER_COMPILATION "Check header for standalone compilation." OFF)
+option(OGS_USE_PYTHON "Interface with Python" OFF)
+
###################
### Definitions ###
###################
@@ -259,7 +261,6 @@ if(OGS_USE_EIGEN)
add_definitions(-DEIGEN_INITIALIZE_MATRICES_BY_NAN)
endif()
endif()
-
if(MSVC AND OGS_32_BIT)
add_definitions(-DEIGEN_MAX_ALIGN_BYTES=0 -DEIGEN_DONT_ALIGN)
endif()
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/Python/c_Python.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/Python/c_Python.md
new file mode 100644
index 0000000000000000000000000000000000000000..f0a0824aed56cce00616ae2bb55ffe11123102d2
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/Python/c_Python.md
@@ -0,0 +1,4 @@
+Boundary condition using a Python script to compute its values.
+
+See ProcessLib::PythonBoundaryConditionPythonSideInterface for the Python-side
+BC interface.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/Python/t_bc_object.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/Python/t_bc_object.md
new file mode 100644
index 0000000000000000000000000000000000000000..5933158c0315078f23036aff899e4e947ff483a3
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/Python/t_bc_object.md
@@ -0,0 +1,2 @@
+The object (variable name) in the provided Python script that will be used for
+this BC.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/Python/t_flush_stdout.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/Python/t_flush_stdout.md
new file mode 100644
index 0000000000000000000000000000000000000000..8edea230fc2f64e19d3923da1df77da3136ca5c4
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/Python/t_flush_stdout.md
@@ -0,0 +1,3 @@
+Determines if stdout should be flushed each time before and after passing control to the Python script.
+This might slow down operation, but ensures the right order of messages printed
+to stdout from Python and C++.
diff --git a/Documentation/ProjectFile/prj/t_python_script.md b/Documentation/ProjectFile/prj/t_python_script.md
new file mode 100644
index 0000000000000000000000000000000000000000..95a9399605d58a98b6796d9d18c6cea5b9dd03bc
--- /dev/null
+++ b/Documentation/ProjectFile/prj/t_python_script.md
@@ -0,0 +1,5 @@
+Path to a Python script file (absolute or relative to the prj file) that is
+executed by OpenGeoSys if OpenGeoSys has been built with Python support.
+
+This script can be used, e.g., to define boundary conditions that can be
+referenced from the prj file and used by OpenGeoSys.
diff --git a/Jenkinsfile b/Jenkinsfile
index dcbbc0b0be74e0743adde3f0cf9f9f68604ec01d..cb9fe4e1f1a7126602833f8752319cedc02c31c1 100644
--- a/Jenkinsfile
+++ b/Jenkinsfile
@@ -85,7 +85,8 @@ pipeline {
cmakeOptions =
'-DOGS_USE_CONAN=ON ' +
'-DOGS_CPU_ARCHITECTURE=generic ' +
- '-DDOCS_GENERATE_LOGFILE=ON ' // redirects to build/DoxygenWarnings.log
+ '-DDOCS_GENERATE_LOGFILE=ON ' + // redirects to build/DoxygenWarnings.log
+ '-DOGS_USE_PYTHON=ON '
}
build { }
build { target="tests" }
@@ -240,7 +241,8 @@ pipeline {
configure {
cmakeOptions =
'-DOGS_USE_CONAN=ON ' +
- '-DOGS_DOWNLOAD_ADDITIONAL_CONTENT=ON '
+ '-DOGS_DOWNLOAD_ADDITIONAL_CONTENT=ON ' +
+ '-DOGS_USE_PYTHON=ON '
}
build { }
build { target="tests" }
diff --git a/NumLib/DOF/LocalToGlobalIndexMap.cpp b/NumLib/DOF/LocalToGlobalIndexMap.cpp
index 40034c2eb29109184704af316912eb01d828b773..d1a91bacf81e8b908e388d2a40d9daf14b5185f0 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.cpp
+++ b/NumLib/DOF/LocalToGlobalIndexMap.cpp
@@ -191,7 +191,8 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
std::vector<int> const& global_component_ids,
std::vector<int> const& variable_component_offsets,
std::vector<MeshLib::Element*> const& elements,
- NumLib::MeshComponentMap&& mesh_component_map)
+ NumLib::MeshComponentMap&& mesh_component_map,
+ LocalToGlobalIndexMap::ConstructorTag)
: _mesh_subsets(std::move(mesh_subsets)),
_mesh_component_map(std::move(mesh_component_map)),
_variable_component_offsets(variable_component_offsets)
@@ -249,7 +250,42 @@ LocalToGlobalIndexMap* LocalToGlobalIndexMap::deriveBoundaryConstrainedMap(
return new LocalToGlobalIndexMap(
std::move(all_mesh_subsets), global_component_ids,
- _variable_component_offsets, elements, std::move(mesh_component_map));
+ _variable_component_offsets, elements, std::move(mesh_component_map),
+ ConstructorTag{});
+}
+
+std::unique_ptr<LocalToGlobalIndexMap>
+LocalToGlobalIndexMap::deriveBoundaryConstrainedMap(
+ MeshLib::MeshSubset&& new_mesh_subset) const
+{
+ DBUG("Construct reduced local to global index map.");
+
+ // Create a subset of the current mesh component map.
+ std::vector<int> global_component_ids;
+
+ for (int i = 0; i < getNumberOfComponents(); ++i)
+ {
+ global_component_ids.push_back(i);
+ }
+
+ // Elements of the new_mesh_subset's mesh.
+ std::vector<MeshLib::Element*> const& elements =
+ new_mesh_subset.getMesh().getElements();
+
+ auto mesh_component_map = _mesh_component_map.getSubset(
+ _mesh_subsets, new_mesh_subset, global_component_ids);
+
+ // Create copies of the new_mesh_subset for each of the global components.
+ // The last component is moved after the for-loop.
+ std::vector<MeshLib::MeshSubset> all_mesh_subsets;
+ for (int i = 0; i < static_cast<int>(global_component_ids.size()) - 1; ++i)
+ all_mesh_subsets.emplace_back(new_mesh_subset);
+ all_mesh_subsets.emplace_back(std::move(new_mesh_subset));
+
+ return std::make_unique<LocalToGlobalIndexMap>(
+ std::move(all_mesh_subsets), global_component_ids,
+ _variable_component_offsets, elements, std::move(mesh_component_map),
+ ConstructorTag{});
}
std::size_t LocalToGlobalIndexMap::dofSizeWithGhosts() const
diff --git a/NumLib/DOF/LocalToGlobalIndexMap.h b/NumLib/DOF/LocalToGlobalIndexMap.h
index 41c1c02456af72fb3db23199b875237a7eb01f3c..8100ac00c155e200f072c91b1115e18e3051eede 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.h
+++ b/NumLib/DOF/LocalToGlobalIndexMap.h
@@ -41,6 +41,14 @@ namespace NumLib
/// mesh item.
class LocalToGlobalIndexMap final
{
+ // Enables using std::make_unique with private constructors from within
+ // member functions of LocalToGlobalIndexMap. Cf.
+ // http://seanmiddleditch.com/enabling-make_unique-with-private-constructors/
+ struct ConstructorTag
+ {
+ explicit ConstructorTag() = default;
+ };
+
public:
using RowColumnIndices = MathLib::RowColumnIndices<GlobalIndexType>;
using LineIndex = RowColumnIndices::LineIndex;
@@ -90,6 +98,13 @@ public:
std::vector<int> const& component_ids,
MeshLib::MeshSubset&& mesh_subset) const;
+ /// Derive a LocalToGlobalIndexMap constrained to the mesh subset and mesh
+ /// subset's elements. A new mesh component map will be constructed using
+ /// the passed mesh_subset for all variables and components of the current
+ /// LocalToGlobalIndexMap.
+ std::unique_ptr<LocalToGlobalIndexMap> deriveBoundaryConstrainedMap(
+ MeshLib::MeshSubset&& new_mesh_subset) const;
+
/// Returns total number of degrees of freedom including those located in
/// the ghost nodes.
std::size_t dofSizeWithGhosts() const;
@@ -144,10 +159,14 @@ public:
MeshLib::MeshSubset const& getMeshSubset(
int const global_component_id) const;
-private:
- /// Private constructor used by internally created local-to-global index
- /// maps. The mesh_component_map is passed as argument instead of being
- /// created by the constructor.
+ /// The global component id for the specific variable (like velocity) and a
+ /// component (like x, or y, or z).
+ int getGlobalComponent(int const variable_id, int const component_id) const;
+
+ /// Private constructor (ensured by ConstructorTag) used by internally
+ /// created local-to-global index maps. The mesh_component_map is passed as
+ /// argument instead of being created by the constructor.
+ ///
/// \attention The passed mesh_component_map is in undefined state after
/// this construtor.
explicit LocalToGlobalIndexMap(
@@ -155,8 +174,9 @@ private:
std::vector<int> const& global_component_ids,
std::vector<int> const& variable_component_offsets,
std::vector<MeshLib::Element*> const& elements,
- NumLib::MeshComponentMap&& mesh_component_map);
+ NumLib::MeshComponentMap&& mesh_component_map, ConstructorTag);
+private:
template <typename ElementIterator>
void findGlobalIndices(ElementIterator first, ElementIterator last,
std::vector<MeshLib::Node*> const& nodes,
@@ -169,11 +189,6 @@ private:
std::vector<MeshLib::Node*> const& nodes, std::size_t const mesh_id,
const int component_id, const int comp_id_write);
- /// The global component id for the specific variable (like velocity) and a
- /// component (like x, or y, or z).
- int getGlobalComponent(int const variable_id, int const component_id) const;
-
-private:
/// A vector of mesh subsets for each process variables' components.
std::vector<MeshLib::MeshSubset> _mesh_subsets;
NumLib::MeshComponentMap _mesh_component_map;
diff --git a/NumLib/Fem/FiniteElement/TemplateIsoparametric.h b/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
index 8484ea1ddf26b2aa3fe4af5997cd96514d71e1d5..c54a99aec4226b4f51f340ce779bfe80d45194f7 100644
--- a/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
+++ b/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
@@ -12,25 +12,21 @@
#pragma once
-#include <cassert>
#include <boost/math/constants/constants.hpp>
+#include <cassert>
-#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
#include "NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.h"
+#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
namespace NumLib
{
-
/**
* \brief Template class for isoparametric elements
*
* \tparam ShapeFunctionType_ The shape function type.
* \tparam ShapeMatrixTypes_ An aggregate of shape matrix types.
*/
-template <
- class ShapeFunctionType_,
- class ShapeMatrixTypes_
- >
+template <class ShapeFunctionType_, class ShapeMatrixTypes_>
class TemplateIsoparametric
{
public:
@@ -44,36 +40,29 @@ public:
/// Natural coordinates mapping tools specialization for specific
/// MeshElement and ShapeFunction types.
- using NaturalCoordsMappingType = NaturalCoordinatesMapping<
- MeshElementType, ShapeFunctionType, ShapeMatrices>;
+ using NaturalCoordsMappingType =
+ NaturalCoordinatesMapping<MeshElementType, ShapeFunctionType,
+ ShapeMatrices>;
/**
- * Constructor without specifying a mesh element. setMeshElement() must be called afterwards.
+ * Constructor without specifying a mesh element. setMeshElement() must be
+ * called afterwards.
*/
- TemplateIsoparametric()
- : _ele(nullptr)
- {
- }
+ TemplateIsoparametric() : _ele(nullptr) {}
/**
* Construct this object for the given mesh element.
*
* @param e Mesh element object
*/
- TemplateIsoparametric(const MeshElementType &e)
- : _ele(&e)
- {
- }
+ TemplateIsoparametric(const MeshElementType& e) : _ele(&e) {}
~TemplateIsoparametric() = default;
/// return current mesh element
- const MeshElementType* getMeshElement() const {return _ele;}
+ const MeshElementType* getMeshElement() const { return _ele; }
/// Sets the mesh element
- void setMeshElement(const MeshElementType &e)
- {
- this->_ele = &e;
- }
+ void setMeshElement(const MeshElementType& e) { this->_ele = &e; }
/**
* compute shape functions
*
@@ -112,23 +101,45 @@ public:
computeIntegralMeasure(is_axially_symmetric, shape);
}
+ /// Interpolates the x coordinate of the element with the given shape
+ /// matrix.
double interpolateZerothCoordinate(
typename ShapeMatrices::ShapeType const& N) const
{
auto* nodes = _ele->getNodes();
typename ShapeMatrices::ShapeType rs(N.size());
- for (int i=0; i<rs.size(); ++i) {
+ for (int i = 0; i < rs.size(); ++i)
+ {
rs[i] = (*nodes[i])[0];
}
auto const r = N.dot(rs);
return r;
}
+ /// Interpolates the coordinates of the element with the given shape matrix.
+ std::array<double, 3> interpolateCoordinates(
+ typename ShapeMatrices::ShapeType const& N) const
+ {
+ auto* nodes = _ele->getNodes();
+ typename ShapeMatrices::ShapeType rs(N.size());
+ std::array<double, 3> interpolated_coords;
+ for (int d = 0; d < 3; ++d)
+ {
+ for (int i = 0; i < rs.size(); ++i)
+ {
+ rs[i] = (*nodes[i])[d];
+ }
+ interpolated_coords[d] = N.dot(rs);
+ }
+ return interpolated_coords;
+ }
+
private:
void computeIntegralMeasure(bool is_axially_symmetric,
ShapeMatrices& shape) const
{
- if (!is_axially_symmetric) {
+ if (!is_axially_symmetric)
+ {
shape.integralMeasure = 1.0;
return;
}
@@ -140,11 +151,10 @@ private:
// located at edge of the unit triangle, it is possible that
// r becomes zero.
auto const r = interpolateZerothCoordinate(shape.N);
- shape.integralMeasure =
- boost::math::constants::two_pi<double>() * r;
+ shape.integralMeasure = boost::math::constants::two_pi<double>() * r;
}
const MeshElementType* _ele;
};
-} // NumLib
+} // namespace NumLib
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index f1e4c472c962cfe950d9d76d941f3ae5da8517a7..ba32dd836d7c5b19c54abe5f20dd98f2686609af 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -209,7 +209,7 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
BaseLib::RunTime time_dirichlet;
time_dirichlet.start();
- sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
+ sys.applyKnownSolutionsNewton(J, res, minus_delta_x, x);
INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet.elapsed());
if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
@@ -230,8 +230,7 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
// cf.
// http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecWAXPY.html
auto& x_new =
- NumLib::GlobalVectorProvider::provider.getVector(
- x, _x_new_id);
+ NumLib::GlobalVectorProvider::provider.getVector(x, _x_new_id);
LinAlg::axpy(x_new, -_damping, minus_delta_x);
if (postIterationCallback)
@@ -252,16 +251,14 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
"iteration"
" has to be repeated.");
// TODO introduce some onDestroy hook.
- NumLib::GlobalVectorProvider::provider
- .releaseVector(x_new);
+ NumLib::GlobalVectorProvider::provider.releaseVector(x_new);
continue; // That throws the iteration result away.
}
// TODO could be done via swap. Note: that also requires swapping
// the ids. Same for the Picard scheme.
LinAlg::copy(x_new, x); // copy new solution to x
- NumLib::GlobalVectorProvider::provider.releaseVector(
- x_new);
+ NumLib::GlobalVectorProvider::provider.releaseVector(x_new);
}
if (!iteration_succeeded)
diff --git a/NumLib/ODESolver/NonlinearSystem.h b/NumLib/ODESolver/NonlinearSystem.h
index 91269cb7151948a6d1006ce9919977e880651880..b3bf35c59a6a99be5ab17ef1d0f69226b62f9efc 100644
--- a/NumLib/ODESolver/NonlinearSystem.h
+++ b/NumLib/ODESolver/NonlinearSystem.h
@@ -59,7 +59,8 @@ public:
//! Apply known solutions to the linearized equation system
//! \f$ \mathit{Jac} \cdot (-\Delta x) = \mathit{res} \f$.
virtual void applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
- GlobalVector& minus_delta_x) = 0;
+ GlobalVector& minus_delta_x,
+ GlobalVector& x) = 0;
};
/*! A System of nonlinear equations to be solved with the Picard fixpoint
diff --git a/NumLib/ODESolver/ODESystem.h b/NumLib/ODESolver/ODESystem.h
index 8d16dd2b235dc9848eb252feddbef614267eaada..64005e7206fbdf952cf05678ab4030cd3974b0d7 100644
--- a/NumLib/ODESolver/ODESystem.h
+++ b/NumLib/ODESolver/ODESystem.h
@@ -48,18 +48,19 @@ public:
virtual void preAssemble(const double t, GlobalVector const& x) = 0;
//! Assemble \c M, \c K and \c b at the provided state (\c t, \c x).
- virtual void assemble(
- const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
- GlobalVector& b) = 0;
+ virtual void assemble(const double t, GlobalVector const& x,
+ GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b) = 0;
using Index = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
//! Provides known solutions (Dirichlet boundary conditions) vector for
//! the ode system at the given time \c t.
virtual std::vector<NumLib::IndexValueVector<Index>> const*
- getKnownSolutions(double const t) const
+ getKnownSolutions(double const t, GlobalVector const& x) const
{
(void)t;
+ (void)x;
return nullptr; // by default there are no known solutions
}
};
@@ -76,7 +77,6 @@ class ODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Picard>
{
public:
-
//! Calls process' pre-assembly with the provided state (\c t, \c x).
virtual void preAssemble(const double t, GlobalVector const& x) = 0;
@@ -124,11 +124,12 @@ public:
* \mathtt{Jac} \f$.
* \endparblock
*/
- virtual void assembleWithJacobian(
- const double t, GlobalVector const& x, GlobalVector const& xdot,
- const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
- GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) = 0;
+ virtual void assembleWithJacobian(const double t, GlobalVector const& x,
+ GlobalVector const& xdot,
+ const double dxdot_dx, const double dx_dx,
+ GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b, GlobalMatrix& Jac) = 0;
};
//! @}
-}
+} // namespace NumLib
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index 5c597a88e5f55e557feb2cf37291a4d7d642ae02..617c722a65fe4f41124d7b7712fcf457b35bb1e7 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -34,14 +34,14 @@ void applyKnownSolutions(std::vector<Solutions> const* const known_solutions,
}
MathLib::LinAlg::finalizeAssembly(x);
}
-}
+} // namespace detail
namespace NumLib
{
TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Newton>::
- TimeDiscretizedODESystem(const int process_id,
- ODE& ode, TimeDisc& time_discretization)
+ TimeDiscretizedODESystem(const int process_id, ODE& ode,
+ TimeDisc& time_discretization)
: _ode(ode),
_time_disc(time_discretization),
_mat_trans(createMatrixTranslator<ODETag>(time_discretization))
@@ -66,9 +66,9 @@ TimeDiscretizedODESystem<
NumLib::GlobalVectorProvider::provider.releaseVector(*_b);
}
-void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
- NonlinearSolverTag::Newton>::
- assemble(const GlobalVector& x_new_timestep)
+void TimeDiscretizedODESystem<
+ ODESystemTag::FirstOrderImplicitQuasilinear,
+ NonlinearSolverTag::Newton>::assemble(const GlobalVector& x_new_timestep)
{
namespace LinAlg = MathLib::LinAlg;
@@ -125,16 +125,16 @@ void TimeDiscretizedODESystem<
NonlinearSolverTag::Newton>::applyKnownSolutions(GlobalVector& x) const
{
::detail::applyKnownSolutions(
- _ode.getKnownSolutions(_time_disc.getCurrentTime()), x);
+ _ode.getKnownSolutions(_time_disc.getCurrentTime(), x), x);
}
void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Newton>::
applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
- GlobalVector& minus_delta_x)
+ GlobalVector& minus_delta_x, GlobalVector& x)
{
auto const* known_solutions =
- _ode.getKnownSolutions(_time_disc.getCurrentTime());
+ _ode.getKnownSolutions(_time_disc.getCurrentTime(), x);
if (!known_solutions || known_solutions->empty())
return;
@@ -149,6 +149,8 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
// For the Newton method the values must be zero
std::vector<double> values(ids.size(), 0);
MathLib::applyKnownSolution(Jac, res, minus_delta_x, ids, values);
+
+ ::detail::applyKnownSolutions(known_solutions, x);
}
TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
@@ -176,9 +178,9 @@ TimeDiscretizedODESystem<
NumLib::GlobalVectorProvider::provider.releaseVector(*_b);
}
-void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
- NonlinearSolverTag::Picard>::
- assemble(const GlobalVector& x_new_timestep)
+void TimeDiscretizedODESystem<
+ ODESystemTag::FirstOrderImplicitQuasilinear,
+ NonlinearSolverTag::Picard>::assemble(const GlobalVector& x_new_timestep)
{
namespace LinAlg = MathLib::LinAlg;
@@ -202,7 +204,7 @@ void TimeDiscretizedODESystem<
NonlinearSolverTag::Picard>::applyKnownSolutions(GlobalVector& x) const
{
::detail::applyKnownSolutions(
- _ode.getKnownSolutions(_time_disc.getCurrentTime()), x);
+ _ode.getKnownSolutions(_time_disc.getCurrentTime(), x), x);
}
void TimeDiscretizedODESystem<
@@ -212,7 +214,7 @@ void TimeDiscretizedODESystem<
GlobalVector& x)
{
auto const* known_solutions =
- _ode.getKnownSolutions(_time_disc.getCurrentTime());
+ _ode.getKnownSolutions(_time_disc.getCurrentTime(), x);
if (known_solutions)
{
@@ -229,4 +231,4 @@ void TimeDiscretizedODESystem<
}
}
-} // NumLib
+} // namespace NumLib
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
index 30e8ab1948be14bf453a3963ead8c43a36a792f1..59ead507782483e9247ee37c7fd2d91f13e20287 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.h
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -92,7 +92,8 @@ public:
void applyKnownSolutions(GlobalVector& x) const override;
void applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
- GlobalVector& minus_delta_x) override;
+ GlobalVector& minus_delta_x,
+ GlobalVector& x) override;
bool isLinear() const override
{
diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.h b/ProcessLib/BoundaryCondition/BoundaryCondition.h
index f20c2ced8db05f9f4f5fd2565e2ba9fd8d6308fd..298609592a53a09b3e72a9911de394f81f7fbd20 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.h
@@ -35,7 +35,8 @@ public:
//! Applies natural BCs (i.e. non-Dirichlet BCs) to the stiffness matrix
//! \c K and the vector \c b.
virtual void applyNaturalBC(const double /*t*/, GlobalVector const& /*x*/,
- GlobalMatrix& /*K*/, GlobalVector& /*b*/)
+ GlobalMatrix& /*K*/, GlobalVector& /*b*/,
+ GlobalMatrix* /*Jac*/)
{
// By default it is assumed that the BC is not a natural BC. Therefore
// there is nothing to do here.
@@ -43,7 +44,7 @@ public:
//! Writes the values of essential BCs to \c bc_values.
virtual void getEssentialBCValues(
- const double /*t*/,
+ const double /*t*/, GlobalVector const& /*x*/,
NumLib::IndexValueVector<GlobalIndexType>& /*bc_values*/) const
{
// By default it is assumed that the BC is not an essential BC.
diff --git a/ProcessLib/BoundaryCondition/BoundaryConditionCollection.cpp b/ProcessLib/BoundaryCondition/BoundaryConditionCollection.cpp
index 63293a99ab7efe5da2acb545b45f7fd6c6549afc..f902b7b35ef11acf9ca40287a2ce6e4414de5ef0 100644
--- a/ProcessLib/BoundaryCondition/BoundaryConditionCollection.cpp
+++ b/ProcessLib/BoundaryCondition/BoundaryConditionCollection.cpp
@@ -14,10 +14,11 @@ namespace ProcessLib
void BoundaryConditionCollection::applyNaturalBC(const double t,
GlobalVector const& x,
GlobalMatrix& K,
- GlobalVector& b)
+ GlobalVector& b,
+ GlobalMatrix* Jac)
{
for (auto const& bc : _boundary_conditions)
- bc->applyNaturalBC(t, x, K, b);
+ bc->applyNaturalBC(t, x, K, b, Jac);
}
void BoundaryConditionCollection::addBCsForProcessVariables(
@@ -43,4 +44,4 @@ void BoundaryConditionCollection::addBCsForProcessVariables(
// object if needed.
_dirichlet_bcs.resize(_boundary_conditions.size());
}
-}
+} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h b/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h
index 26adc6ae3ee1f781a0047f06ae10aea9f5bc2330..3f4af9810ebe01ae24fa564999af34619e7e183e 100644
--- a/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h
+++ b/ProcessLib/BoundaryCondition/BoundaryConditionCollection.h
@@ -26,16 +26,16 @@ public:
}
void applyNaturalBC(const double t, GlobalVector const& x, GlobalMatrix& K,
- GlobalVector& b);
+ GlobalVector& b, GlobalMatrix* Jac);
std::vector<NumLib::IndexValueVector<GlobalIndexType>> const*
- getKnownSolutions(double const t) const
+ getKnownSolutions(double const t, GlobalVector const& x) const
{
auto const n_bcs = _boundary_conditions.size();
for (std::size_t i=0; i<n_bcs; ++i) {
auto const& bc = *_boundary_conditions[i];
auto& dirichlet_storage = _dirichlet_bcs[i];
- bc.getEssentialBCValues(t, dirichlet_storage);
+ bc.getEssentialBCValues(t, x, dirichlet_storage);
}
return &_dirichlet_bcs;
}
diff --git a/ProcessLib/BoundaryCondition/BoundaryConditionConfig.h b/ProcessLib/BoundaryCondition/BoundaryConditionConfig.h
index f1f7091ea0d7e07a6e5b97de0ad3b1ce5d7e0e3a..4a537eca95272c791cd86df6ec972babba4ee576 100644
--- a/ProcessLib/BoundaryCondition/BoundaryConditionConfig.h
+++ b/ProcessLib/BoundaryCondition/BoundaryConditionConfig.h
@@ -19,19 +19,21 @@ struct BoundaryConditionConfig final
BoundaryConditionConfig(BaseLib::ConfigTree&& config_,
MeshLib::Mesh const& mesh_,
boost::optional<int> const component_id_)
- : config(std::move(config_)), mesh(mesh_), component_id(component_id_)
+ : config(std::move(config_)),
+ boundary_mesh(mesh_),
+ component_id(component_id_)
{
}
BoundaryConditionConfig(BoundaryConditionConfig&& other)
: config(std::move(other.config)),
- mesh(other.mesh),
+ boundary_mesh(other.boundary_mesh),
component_id(other.component_id)
{
}
BaseLib::ConfigTree config;
- MeshLib::Mesh const& mesh;
+ MeshLib::Mesh const& boundary_mesh;
boost::optional<int> const component_id;
};
diff --git a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
index acce3726620c45b75644b929f8c681063e7d7717..33567c6e25e44812d3e77c65358d1685b02d91c2 100644
--- a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.cpp
@@ -139,7 +139,8 @@ void ConstraintDirichletBoundaryCondition::preTimestep(double t,
}
void ConstraintDirichletBoundaryCondition::getEssentialBCValues(
- const double t, NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
+ const double t, const GlobalVector& /*x*/,
+ NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
{
SpatialPosition pos;
diff --git a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.h
index 253abfceb7fddd23c46f41ad64de4d8e63e7b1f8..04d3b4095d6906f3ecb0b7d30441e6df668ad92a 100644
--- a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryCondition.h
@@ -63,7 +63,7 @@ public:
void preTimestep(double t, GlobalVector const& x) override;
void getEssentialBCValues(
- const double t,
+ const double t, const GlobalVector& x,
NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
private:
diff --git a/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp
index 35948a821937ebded3d48920295f9fbc5355a08d..fb9d8e9f091c66ae932fe6f7d2cb82623507bd01 100644
--- a/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp
@@ -19,14 +19,17 @@
#include "NormalTractionBoundaryCondition.h"
#include "PhaseFieldIrreversibleDamageOracleBoundaryCondition.h"
#include "RobinBoundaryCondition.h"
+#ifdef OGS_USE_PYTHON
+#include "Python/PythonBoundaryCondition.h"
+#endif
namespace ProcessLib
{
std::unique_ptr<BoundaryCondition> createBoundaryCondition(
const BoundaryConditionConfig& config,
- const NumLib::LocalToGlobalIndexMap& dof_table, const MeshLib::Mesh& mesh,
- const int variable_id, const unsigned integration_order,
- const unsigned shapefunction_order,
+ const NumLib::LocalToGlobalIndexMap& dof_table,
+ const MeshLib::Mesh& bulk_mesh, const int variable_id,
+ const unsigned integration_order, const unsigned shapefunction_order,
const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>& parameters,
const Process& process)
{
@@ -36,33 +39,47 @@ std::unique_ptr<BoundaryCondition> createBoundaryCondition(
if (type == "Dirichlet")
{
return ProcessLib::createDirichletBoundaryCondition(
- config.config, config.mesh, dof_table, mesh.getID(), variable_id,
- *config.component_id, parameters);
+ config.config, config.boundary_mesh, dof_table, bulk_mesh.getID(),
+ variable_id, *config.component_id, parameters);
}
if (type == "Neumann")
{
return ProcessLib::createNeumannBoundaryCondition(
- config.config, config.mesh, dof_table, variable_id,
- *config.component_id, mesh.isAxiallySymmetric(), integration_order,
- shapefunction_order, mesh.getDimension(), parameters);
+ config.config, config.boundary_mesh, dof_table, variable_id,
+ *config.component_id, bulk_mesh.isAxiallySymmetric(),
+ integration_order, shapefunction_order, bulk_mesh.getDimension(),
+ parameters);
}
if (type == "Robin")
{
return ProcessLib::createRobinBoundaryCondition(
- config.config, config.mesh, dof_table, variable_id,
- *config.component_id, mesh.isAxiallySymmetric(), integration_order,
- shapefunction_order, mesh.getDimension(), parameters);
+ config.config, config.boundary_mesh, dof_table, variable_id,
+ *config.component_id, bulk_mesh.isAxiallySymmetric(),
+ integration_order, shapefunction_order, bulk_mesh.getDimension(),
+ parameters);
}
if (type == "NonuniformDirichlet")
{
return ProcessLib::createNonuniformDirichletBoundaryCondition(
- config.config, dof_table, variable_id, *config.component_id, mesh);
+ config.config, dof_table, variable_id, *config.component_id,
+ bulk_mesh);
}
if (type == "NonuniformNeumann")
{
return ProcessLib::createNonuniformNeumannBoundaryCondition(
config.config, dof_table, variable_id, *config.component_id,
- integration_order, shapefunction_order, mesh);
+ integration_order, shapefunction_order, bulk_mesh);
+ }
+ if (type == "Python")
+ {
+#ifdef OGS_USE_PYTHON
+ return ProcessLib::createPythonBoundaryCondition(
+ config.config, config.boundary_mesh, dof_table, bulk_mesh.getID(),
+ variable_id, *config.component_id, integration_order,
+ shapefunction_order, bulk_mesh.getDimension());
+#else
+ OGS_FATAL("OpenGeoSys has not been built with Python support.");
+#endif
}
//
@@ -71,22 +88,22 @@ std::unique_ptr<BoundaryCondition> createBoundaryCondition(
if (type == "ConstraintDirichlet")
{
return createConstraintDirichletBoundaryCondition(
- config.config, config.mesh, dof_table, variable_id,
+ config.config, config.boundary_mesh, dof_table, variable_id,
integration_order, *config.component_id, parameters, process);
}
if (type == "NormalTraction")
{
return ProcessLib::NormalTractionBoundaryCondition::
createNormalTractionBoundaryCondition(
- config.config, config.mesh, dof_table, variable_id,
- mesh.isAxiallySymmetric(), integration_order,
- shapefunction_order, mesh.getDimension(), parameters);
+ config.config, config.boundary_mesh, dof_table, variable_id,
+ bulk_mesh.isAxiallySymmetric(), integration_order,
+ shapefunction_order, bulk_mesh.getDimension(), parameters);
}
if (type == "PhaseFieldIrreversibleDamageOracleBoundaryCondition")
{
return ProcessLib::
createPhaseFieldIrreversibleDamageOracleBoundaryCondition(
- config.config, dof_table, mesh, variable_id,
+ config.config, dof_table, bulk_mesh, variable_id,
*config.component_id);
}
OGS_FATAL("Unknown boundary condition type: `%s'.", type.c_str());
diff --git a/ProcessLib/BoundaryCondition/CreateBoundaryCondition.h b/ProcessLib/BoundaryCondition/CreateBoundaryCondition.h
index 97478a17d2ccef2d057d9778ffbe13827df352af..683951c8b89d300d6ea31526fcf336043fc17e3d 100644
--- a/ProcessLib/BoundaryCondition/CreateBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/CreateBoundaryCondition.h
@@ -31,9 +31,9 @@ class Process;
std::unique_ptr<BoundaryCondition> createBoundaryCondition(
const BoundaryConditionConfig& config,
- const NumLib::LocalToGlobalIndexMap& dof_table, const MeshLib::Mesh& mesh,
- const int variable_id, const unsigned integration_order,
- const unsigned shapefunction_order,
+ const NumLib::LocalToGlobalIndexMap& dof_table,
+ const MeshLib::Mesh& bulk_mesh, const int variable_id,
+ const unsigned integration_order, const unsigned shapefunction_order,
const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>& parameters,
const Process& process);
diff --git a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp
index a81977cdc17c99b7ca976249843cd510cdf8b6f3..38a7e6d08c23e77390792e9855f37e2348c72aac 100644
--- a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp
@@ -10,14 +10,15 @@
#include "DirichletBoundaryCondition.h"
#include <algorithm>
-#include <vector>
#include <logog/include/logog.hpp>
+#include <vector>
#include "ProcessLib/Utils/ProcessUtils.h"
namespace ProcessLib
{
void DirichletBoundaryCondition::getEssentialBCValues(
- const double t, NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
+ const double t, GlobalVector const& /*x*/,
+ NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
{
SpatialPosition pos;
@@ -48,7 +49,8 @@ void DirichletBoundaryCondition::getEssentialBCValues(
// 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 (g_idx >= 0) {
+ if (g_idx >= 0)
+ {
bc_values.ids.emplace_back(g_idx);
bc_values.values.emplace_back(_parameter(t, pos).front());
}
diff --git a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.h
index 721ea5d3ad1d3d4daa67782a1386c2e8720de25a..f00681518cfe05cd57d18b8f6ed1f3e1c7c6ec41 100644
--- a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.h
@@ -58,7 +58,7 @@ public:
}
void getEssentialBCValues(
- const double t,
+ const double t, GlobalVector const& x,
NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
private:
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
index b164ed2b8c0f9bb12760060be1aa709f9322b2b2..4da13051714b8404c9380b8f25c46baf3e280c4c 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
@@ -70,11 +70,12 @@ void GenericNaturalBoundaryCondition<
LocalAssemblerImplementation>::applyNaturalBC(const double t,
const GlobalVector& x,
GlobalMatrix& K,
- GlobalVector& b)
+ GlobalVector& b,
+ GlobalMatrix* Jac)
{
GlobalExecutor::executeMemberOnDereferenced(
&GenericNaturalBoundaryConditionLocalAssemblerInterface::assemble,
- _local_assemblers, *_dof_table_boundary, t, x, K, b);
+ _local_assemblers, *_dof_table_boundary, t, x, K, b, Jac);
}
} // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
index f68304dc2f85a46f3d314637f961cf812d25d9f8..90b0d164f9ade9ccea7a830c2359d6803c5dbf62 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition.h
@@ -9,8 +9,8 @@
#pragma once
-#include "MeshLib/MeshSubset.h"
#include "BoundaryCondition.h"
+#include "MeshLib/MeshSubset.h"
namespace ProcessLib
{
@@ -38,10 +38,8 @@ public:
/// Calls local assemblers which calculate their contributions to the global
/// matrix and the right-hand-side.
- void applyNaturalBC(const double t,
- GlobalVector const& x,
- GlobalMatrix& K,
- GlobalVector& b) override;
+ void applyNaturalBC(const double t, GlobalVector const& x, GlobalMatrix& K,
+ GlobalVector& b, GlobalMatrix* Jac) override;
private:
/// Data used in the assembly of the specific boundary condition.
@@ -63,6 +61,6 @@ private:
_local_assemblers;
};
-} // ProcessLib
+} // namespace ProcessLib
#include "GenericNaturalBoundaryCondition-impl.h"
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h
index 219085cfc062f1d2a36ed66822919d96fc5f4128..1721ef2589d6f2632c401fce3fda976076932bc0 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h
@@ -23,7 +23,8 @@ public:
virtual void assemble(
std::size_t const id,
NumLib::LocalToGlobalIndexMap const& dof_table_boundary, double const t,
- const GlobalVector& x, GlobalMatrix& K, GlobalVector& b) = 0;
+ const GlobalVector& x, GlobalMatrix& K, GlobalVector& b,
+ GlobalMatrix* Jac) = 0;
};
template <typename ShapeFunction, typename IntegrationMethod,
@@ -55,4 +56,4 @@ protected:
_shape_matrices;
};
-} // ProcessLib
+} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h
index a7b50d8957aa117cd9f5483c814f4947151d8775..bf3e122f78bd1fa56384a69dc6b8e4ddf8d6a45a 100644
--- a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h
@@ -89,12 +89,13 @@ void GenericNonuniformNaturalBoundaryCondition<
LocalAssemblerImplementation>::applyNaturalBC(const double t,
const GlobalVector& x,
GlobalMatrix& K,
- GlobalVector& b)
+ GlobalVector& b,
+ GlobalMatrix* Jac)
{
GlobalExecutor::executeMemberOnDereferenced(
&GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface::
assemble,
- _local_assemblers, *_dof_table_boundary, t, x, K, b);
+ _local_assemblers, *_dof_table_boundary, t, x, K, b, Jac);
}
-} // ProcessLib
+} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h
index d6753d967192bde3de7360f2cb11c9db3b74e56c..37249a52c5a904032091b7f21a64d75602997445 100644
--- a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h
@@ -33,10 +33,8 @@ public:
/// Calls local assemblers which calculate their contributions to the global
/// matrix and the right-hand-side.
- void applyNaturalBC(const double t,
- GlobalVector const& x,
- GlobalMatrix& K,
- GlobalVector& b) override;
+ void applyNaturalBC(const double t, GlobalVector const& x, GlobalMatrix& K,
+ GlobalVector& b, GlobalMatrix* Jac) override;
private:
void constructDofTable();
diff --git a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryConditionLocalAssembler.h
index 72929d11a5e3fb97a71905bb7c08464627d09230..3cf49814d088bfa59c0da7b0b6c4665d72293ddf 100644
--- a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryConditionLocalAssembler.h
@@ -24,7 +24,8 @@ public:
virtual void assemble(
std::size_t const id,
NumLib::LocalToGlobalIndexMap const& dof_table_boundary, double const t,
- const GlobalVector& x, GlobalMatrix& K, GlobalVector& b) = 0;
+ const GlobalVector& x, GlobalMatrix& K, GlobalVector& b,
+ GlobalMatrix* Jac) = 0;
};
template <typename ShapeFunction, typename IntegrationMethod,
@@ -88,4 +89,4 @@ protected:
_ns_and_weights;
};
-} // ProcessLib
+} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
index cface4f4cdb2dfa365e8caccaa5e57ea4e10d080..a88661d70cd8166e41f262d197d3fa9d788e203e 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
@@ -9,9 +9,9 @@
#pragma once
+#include "GenericNaturalBoundaryConditionLocalAssembler.h"
#include "NumLib/DOF/DOFTableUtil.h"
#include "ProcessLib/Parameter/Parameter.h"
-#include "GenericNaturalBoundaryConditionLocalAssembler.h"
namespace ProcessLib
{
@@ -42,7 +42,8 @@ public:
void assemble(std::size_t const id,
NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
double const t, const GlobalVector& /*x*/,
- GlobalMatrix& /*K*/, GlobalVector& b) override
+ GlobalMatrix& /*K*/, GlobalVector& b,
+ GlobalMatrix* /*Jac*/) override
{
_local_rhs.setZero();
@@ -74,4 +75,4 @@ public:
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
-} // namespace ProcessLib
+} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h
index 693993742036af5cce389794849d4931638f251c..49ae3aacf7ee8206f61547129c18eece09e785ab 100644
--- a/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h
@@ -60,11 +60,9 @@ public:
}
void getEssentialBCValues(
- const double /*t*/,
+ const double /*t*/, GlobalVector const& /*x*/,
NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override
{
- SpatialPosition pos;
-
bc_values.ids.clear();
bc_values.values.clear();
diff --git a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h
index d995082e418215755cbd3a3d3e1e46631eb19339..3666afdcd39006e608152b6ecf6d9d5bcfbfadc7 100644
--- a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h
@@ -57,7 +57,8 @@ public:
void assemble(std::size_t const id,
NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
double const /*t*/, const GlobalVector& /*x*/,
- GlobalMatrix& /*K*/, GlobalVector& b) override
+ GlobalMatrix& /*K*/, GlobalVector& b,
+ GlobalMatrix* /*Jac*/) override
{
_local_rhs.setZero();
diff --git a/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h
index c0c26e40d9987eaeac43adf617884de960498eed..ab5fad4556ba652be1e860f72df549005efaac1f 100644
--- a/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h
@@ -65,11 +65,12 @@ void NormalTractionBoundaryCondition<
LocalAssemblerImplementation>::applyNaturalBC(const double t,
const GlobalVector& x,
GlobalMatrix& K,
- GlobalVector& b)
+ GlobalVector& b,
+ GlobalMatrix* Jac)
{
GlobalExecutor::executeMemberOnDereferenced(
&NormalTractionBoundaryConditionLocalAssemblerInterface::assemble,
- _local_assemblers, *_dof_table_boundary, t, x, K, b);
+ _local_assemblers, *_dof_table_boundary, t, x, K, b, Jac);
}
std::unique_ptr<NormalTractionBoundaryCondition<
@@ -98,4 +99,4 @@ createNormalTractionBoundaryCondition(
}
} // namespace NormalTractionBoundaryCondition
-} // ProcessLib
+} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition.h b/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition.h
index 9f849aadc890b10fc0c6a523f5a6008d778eb761..14cb196809e7b29dd583e604364f9b1655f1213d 100644
--- a/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition.h
@@ -44,10 +44,8 @@ public:
/// Calls local assemblers which calculate their contributions to the global
/// matrix and the right-hand-side.
- void applyNaturalBC(const double t,
- GlobalVector const& x,
- GlobalMatrix& K,
- GlobalVector& b) override;
+ void applyNaturalBC(const double t, GlobalVector const& x, GlobalMatrix& K,
+ GlobalVector& b, GlobalMatrix* Jac) override;
private:
MeshLib::Mesh const& _bc_mesh;
diff --git a/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h
index e296fa1071d80b2852adfbbc4f6a750fd32ea4a4..cfe97bbc449dd8988a1fdf16e79ab8efae58ae1f 100644
--- a/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h
@@ -10,8 +10,8 @@
#pragma once
-#include "MeshLib/Elements/FaceRule.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "MeshLib/Elements/FaceRule.h"
#include "NumLib/DOF/DOFTableUtil.h"
#include "ProcessLib/Parameter/Parameter.h"
@@ -45,7 +45,8 @@ public:
virtual void assemble(
std::size_t const id,
NumLib::LocalToGlobalIndexMap const& dof_table_boundary, double const t,
- const GlobalVector& /*x*/, GlobalMatrix& /*K*/, GlobalVector& b) = 0;
+ const GlobalVector& /*x*/, GlobalMatrix& /*K*/, GlobalVector& b,
+ GlobalMatrix* /*Jac*/) = 0;
virtual ~NormalTractionBoundaryConditionLocalAssemblerInterface() = default;
};
@@ -126,7 +127,8 @@ public:
void assemble(std::size_t const id,
NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
double const t, const GlobalVector& /*x*/,
- GlobalMatrix& /*K*/, GlobalVector& local_rhs)
+ GlobalMatrix& /*K*/, GlobalVector& local_rhs,
+ GlobalMatrix* /*Jac*/) override
{
_local_rhs.setZero();
diff --git a/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.cpp
index bc72d5b47334e51860f4ee502c9b6b6797c7312e..c9c70288dba6c3fc4262a6b48e6012eca0e70910 100644
--- a/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.cpp
@@ -17,7 +17,7 @@
namespace ProcessLib
{
void PhaseFieldIrreversibleDamageOracleBoundaryCondition::getEssentialBCValues(
- const double /*t*/,
+ const double /*t*/, GlobalVector const& /*x*/,
NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
{
SpatialPosition pos;
diff --git a/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.h b/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.h
index 9afa9f2426c3ceeb3016736c3ea7f5538cf87788..c14b17452717727fa872ee37070568b68ebbd8b6 100644
--- a/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/PhaseFieldIrreversibleDamageOracleBoundaryCondition.h
@@ -43,7 +43,7 @@ public:
}
void getEssentialBCValues(
- const double t,
+ const double t, const GlobalVector& x,
NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
void preTimestep(const double t, const GlobalVector& x) override;
diff --git a/ProcessLib/BoundaryCondition/Python/CMakeLists.txt b/ProcessLib/BoundaryCondition/Python/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..30cfd3828bd182894f28f97c60a02133533b06c3
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/Python/CMakeLists.txt
@@ -0,0 +1,22 @@
+add_library(ProcessLibBoundaryConditionPython
+ PythonBoundaryCondition.cpp
+ PythonBoundaryCondition.h
+ PythonBoundaryConditionLocalAssembler.h
+ PythonBoundaryConditionPythonSideInterface.h)
+
+target_compile_definitions(ProcessLibBoundaryConditionPython
+ PUBLIC OGS_USE_PYTHON)
+
+target_link_libraries(ProcessLibBoundaryConditionPython
+ PUBLIC BaseLib MathLib MeshLib NumLib logog
+ PRIVATE pybind11::pybind11)
+
+# For the embedded Python module
+add_library(ProcessLibBoundaryConditionPythonModule
+ PythonBoundaryConditionModule.cpp
+ PythonBoundaryConditionModule.h)
+
+target_link_libraries(ProcessLibBoundaryConditionPythonModule
+ PUBLIC
+ ProcessLibBoundaryConditionPython
+ pybind11::pybind11)
diff --git a/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..90989956d84d7191530db60a334e0eeaa48a4158
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.cpp
@@ -0,0 +1,239 @@
+/**
+ * \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 "PythonBoundaryCondition.h"
+
+#include <pybind11/pybind11.h>
+#include <iostream>
+
+#include "MeshLib/MeshSearch/NodeSearch.h"
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+#include "PythonBoundaryConditionLocalAssembler.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
+{
+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
+{
+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));
+
+ 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");
+
+ 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 occured 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,
+ const GlobalVector& x,
+ GlobalMatrix& K, GlobalVector& b,
+ GlobalMatrix* Jac)
+{
+ FlushStdoutGuard guard(_flush_stdout);
+
+ try
+ {
+ GlobalExecutor::executeMemberOnDereferenced(
+ &GenericNaturalBoundaryConditionLocalAssemblerInterface::assemble,
+ _local_assemblers, *_dof_table_boundary, t, x, 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)
+{
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
+ config.checkConfigParameter("type", "Python");
+
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__Python__bc_object}
+ auto const bc_object = config.getConfigParameter<std::string>("bc_object");
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__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(bc_object))
+ OGS_FATAL(
+ "Function `%s' is not defined in the python script file, or there "
+ "was no python script file specified.",
+ bc_object.c_str());
+
+ 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));
+ }
+
+ 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
diff --git a/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.h b/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.h
new file mode 100644
index 0000000000000000000000000000000000000000..f3ef475ddeeea2050621c40125266b03738d3942
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.h
@@ -0,0 +1,85 @@
+/**
+ * \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/BoundaryCondition/BoundaryCondition.h"
+#include "ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h"
+
+#include "PythonBoundaryConditionPythonSideInterface.h"
+
+namespace ProcessLib
+{
+//! Groups data used by essential and natural BCs, in particular by the
+//! local assemblers of the latter.
+struct PythonBoundaryConditionData
+{
+ //! Python object computing BC values.
+ PythonBoundaryConditionPythonSideInterface* bc_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 BC is
+ //! applied.
+ int const global_component_id;
+
+ //! The boundary mesh, i.e., the domain of this BC.
+ const MeshLib::Mesh& boundary_mesh;
+};
+
+//! A boundary condition whose values are computed by a Python script.
+class PythonBoundaryCondition final : public BoundaryCondition
+{
+public:
+ PythonBoundaryCondition(PythonBoundaryConditionData&& bc_data,
+ unsigned const integration_order,
+ unsigned const shapefunction_order,
+ unsigned const global_dim,
+ bool const flush_stdout);
+
+ void getEssentialBCValues(
+ const double t, const GlobalVector& x,
+ NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
+
+ void applyNaturalBC(const double t, const GlobalVector& x, GlobalMatrix& K,
+ GlobalVector& b, GlobalMatrix* Jac) override;
+
+private:
+ //! Auxiliary data.
+ PythonBoundaryConditionData _bc_data;
+
+ //! Local dof table for the boundary mesh.
+ std::unique_ptr<NumLib::LocalToGlobalIndexMap> _dof_table_boundary;
+
+ //! Local assemblers for all elements of the boundary mesh.
+ std::vector<
+ std::unique_ptr<GenericNaturalBoundaryConditionLocalAssemblerInterface>>
+ _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;
+};
+
+//! Creates a new PythonBoundaryCondition object.
+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);
+
+} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h
new file mode 100644
index 0000000000000000000000000000000000000000..56ba59325353e9384ccd994645e8a872d7e5a141
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h
@@ -0,0 +1,198 @@
+/**
+ * \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 "PythonBoundaryCondition.h"
+
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h"
+
+#include "PythonBoundaryConditionPythonSideInterface.h"
+
+namespace ProcessLib
+{
+//! 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 ShapeFunction, typename IntegrationMethod,
+ unsigned 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),
+ _element(e)
+ {
+ }
+
+ void assemble(std::size_t const boundary_element_id,
+ NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
+ double const t, const GlobalVector& x, GlobalMatrix& /*K*/,
+ 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 =
+ Base::_integration_method.getNumberOfPoints();
+ auto const num_var = _data.dof_table_bulk.getNumberOfVariables();
+ auto const num_nodes = _element.getNumberOfNodes();
+ auto const num_comp_total =
+ _data.dof_table_bulk.getNumberOfComponents();
+
+ auto const& bulk_node_ids_map =
+ *_data.boundary_mesh.getProperties().getPropertyVector<std::size_t>(
+ "bulk_node_ids");
+
+ // 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 = _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];
+ }
+ }
+ }
+
+ 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& sm = Base::_shape_matrices[ip];
+ auto const coords = fe.interpolateCoordinates(sm.N);
+ prim_vars =
+ sm.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);
+
+ auto const& wp = Base::_integration_method.getWeightedPoint(ip);
+ auto const w = sm.detJ * wp.getWeight() * sm.integralMeasure;
+ local_rhs.noalias() += sm.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 assignement -= takes into account the sign convention
+ // of 1st-order in time ODE systems in OpenGeoSys.
+ local_Jac.block(top, left, width, height).noalias() -=
+ sm.N.transpose() * (dFlux[comp] * w) * sm.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;
+ MeshLib::Element const& _element;
+};
+
+} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionModule.cpp b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionModule.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..45fdfed80c31c225950446f50158fa9aac82a9ff
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionModule.cpp
@@ -0,0 +1,63 @@
+/**
+ * \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 "PythonBoundaryConditionModule.h"
+
+#include <pybind11/stl.h>
+
+#include "PythonBoundaryConditionPythonSideInterface.h"
+
+namespace ProcessLib
+{
+//! Trampoline class allowing methods of class
+//! PythonBoundaryConditionPythonSideInterface to be overridden on the Python
+//! side. Cf. https://pybind11.readthedocs.io/en/stable/advanced/classes.html
+class PythonBoundaryConditionPythonSideInterfaceTrampoline
+ : public PythonBoundaryConditionPythonSideInterface
+{
+public:
+ using PythonBoundaryConditionPythonSideInterface::
+ PythonBoundaryConditionPythonSideInterface;
+
+ std::pair<bool, double> getDirichletBCValue(
+ double t, std::array<double, 3> x, std::size_t node_id,
+ std::vector<double> const& primary_variables) const override
+ {
+ using Ret = std::pair<bool, double>;
+ PYBIND11_OVERLOAD(Ret, PythonBoundaryConditionPythonSideInterface,
+ getDirichletBCValue, t, x, node_id,
+ primary_variables);
+ }
+
+ std::tuple<bool, double, std::vector<double>> getFlux(
+ double t, std::array<double, 3> x,
+ std::vector<double> const& primary_variables) const override
+ {
+ using Ret = std::tuple<bool, double, std::vector<double>>;
+ PYBIND11_OVERLOAD(Ret, PythonBoundaryConditionPythonSideInterface,
+ getFlux, t, x, primary_variables);
+ }
+};
+
+void pythonBindBoundaryCondition(pybind11::module& m)
+{
+ namespace py = pybind11;
+
+ py::class_<PythonBoundaryConditionPythonSideInterface,
+ PythonBoundaryConditionPythonSideInterfaceTrampoline>
+ pybc(m, "BoundaryCondition");
+
+ pybc.def(py::init());
+
+ pybc.def("getDirichletBCValue",
+ &PythonBoundaryConditionPythonSideInterface::getDirichletBCValue);
+ pybc.def("getFlux", &PythonBoundaryConditionPythonSideInterface::getFlux);
+}
+
+} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionModule.h b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionModule.h
new file mode 100644
index 0000000000000000000000000000000000000000..82c141a4bc3aa73734b85f315aeb0576874b6cea
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionModule.h
@@ -0,0 +1,18 @@
+/**
+ * \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>
+
+namespace ProcessLib
+{
+//! Creates Python bindings for the Python BC class.
+void pythonBindBoundaryCondition(pybind11::module& m);
+} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionPythonSideInterface.h b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionPythonSideInterface.h
new file mode 100644
index 0000000000000000000000000000000000000000..43ea47f400eebdf12c8f120c1697826e47cb937d
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionPythonSideInterface.h
@@ -0,0 +1,75 @@
+/**
+ * \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
+{
+//! Base class for boundary conditions.
+//! This class will get Python bindings and is intended to be to be derived in
+//! Python.
+class PythonBoundaryConditionPythonSideInterface
+{
+public:
+ /*!
+ * Computes Dirichlet boundary condition values for the provided arguments
+ * (time, position of the node, node id, primary variables at the node).
+ *
+ * \return a pair (is_dirichlet, value) indicating if a Dirichlet BC shall
+ * be set at that node and (if so) the value of the Dirichlet BC at that
+ * node.
+ */
+ virtual std::pair<bool, double> getDirichletBCValue(
+ double /*t*/, std::array<double, 3> /*x*/, std::size_t /*node_id*/,
+ std::vector<double> const& /*primary_variables*/) const
+ {
+ _overridden_essential = false;
+ return {false, std::numeric_limits<double>::quiet_NaN()};
+ }
+
+ /*!
+ * Computes the flux for the provided arguments (time, position, primary
+ * variables at that position).
+ *
+ * \return a pair (is_natural, flux, flux_jacobian) indicating if a natural
+ * BC shall be set at that position and (if so) the flux at that node and
+ * the derivative of the flux w.r.t. all primary variables.
+ */
+ virtual std::tuple<bool, double, std::vector<double>> getFlux(
+ double /*t*/,
+ std::array<double, 3> /*x*/,
+ std::vector<double> const& /*primary_variables*/) const
+ {
+ _overridden_natural = false;
+ return std::tuple<bool, double, std::vector<double>>{
+ false, std::numeric_limits<double>::quiet_NaN(), {}};
+ }
+
+ //! Tells if getDirichletBCValue() has been overridden in the derived class
+ //! in Python.
+ //!
+ //! \pre getDirichletBCValue() must already have been called
+ //! once.
+ bool isOverriddenEssential() const { return _overridden_essential; }
+
+ //! Tells if getFlux() has been overridden in the derived class in Python.
+ //!
+ //! \pre getFlux() must already have been called once.
+ bool isOverriddenNatural() const { return _overridden_natural; }
+
+ virtual ~PythonBoundaryConditionPythonSideInterface() = default;
+
+private:
+ //! Tells if getDirichletBCValue() has been overridden in the derived class
+ //! in Python.
+ mutable bool _overridden_essential = true;
+ //! Tells if getFlux() has been overridden in the derived class in Python.
+ mutable bool _overridden_natural = true;
+};
+} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
index f3bdc6181ee504e003ce2cc2167b42aab914b438..1b780d40b04476d9d28ab1be117f40181e11a6bd 100644
--- a/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
@@ -46,7 +46,7 @@ public:
void assemble(std::size_t const id,
NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
double const t, const GlobalVector& /*x*/, GlobalMatrix& K,
- GlobalVector& b) override
+ GlobalVector& b, GlobalMatrix* /*Jac*/) override
{
_local_K.setZero();
_local_rhs.setZero();
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index 71a28b6b06e69e8b8130ee4be412be452d0eba9b..7d432a84e61352c33031ac05c33735cf8a495afe 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/ProcessLib/CMakeLists.txt
@@ -25,6 +25,12 @@ target_link_libraries(ProcessLib
PUBLIC BaseLib MaterialLib MathLib MeshLib NumLib logog
)
+if(OGS_USE_PYTHON)
+ add_subdirectory(BoundaryCondition/Python)
+ target_link_libraries(ProcessLib
+ PUBLIC ProcessLibBoundaryConditionPython)
+endif()
+
if(OGS_INSITU)
target_link_libraries(ProcessLib InSituLib)
endif()
diff --git a/ProcessLib/ComponentTransport/Tests.cmake b/ProcessLib/ComponentTransport/Tests.cmake
index b1e59587d21004d65973a056503be1075c72233e..89070f85d7119bed5fb146d0ec8646c7934390bd 100644
--- a/ProcessLib/ComponentTransport/Tests.cmake
+++ b/ProcessLib/ComponentTransport/Tests.cmake
@@ -295,6 +295,69 @@ AddTest(
elder_pcs_0_ts_80_t_21038400.000000_reference.vtu elder_pcs_0_ts_80_t_21038400.000000.vtu conc conc 1e-1 1e-5
)
+AddTest(
+ NAME LARGE_2D_ComponentTransport_ElderPython
+ PATH Parabolic/ComponentTransport/elder
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS elder-python.prj
+ WRAPPER time
+ TESTER vtkdiff
+ REQUIREMENTS OGS_USE_PYTHON AND NOT OGS_USE_MPI
+ DIFF_DATA
+ elder_pcs_0_ts_0_t_0.000000_reference.vtu elder_python_pcs_0_ts_0_t_0.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_100_t_26298000.000000_reference.vtu elder_python_pcs_0_ts_100_t_26298000.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_120_t_31557600.000000_reference.vtu elder_python_pcs_0_ts_120_t_31557600.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_140_t_36817200.000000_reference.vtu elder_python_pcs_0_ts_140_t_36817200.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_160_t_42076800.000000_reference.vtu elder_python_pcs_0_ts_160_t_42076800.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_180_t_47336400.000000_reference.vtu elder_python_pcs_0_ts_180_t_47336400.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_200_t_52596000.000000_reference.vtu elder_python_pcs_0_ts_200_t_52596000.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_20_t_5259600.000000_reference.vtu elder_python_pcs_0_ts_20_t_5259600.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_220_t_57855600.000000_reference.vtu elder_python_pcs_0_ts_220_t_57855600.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_240_t_63115200.000000_reference.vtu elder_python_pcs_0_ts_240_t_63115200.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_260_t_68374800.000000_reference.vtu elder_python_pcs_0_ts_260_t_68374800.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_280_t_73634400.000000_reference.vtu elder_python_pcs_0_ts_280_t_73634400.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_300_t_78894000.000000_reference.vtu elder_python_pcs_0_ts_300_t_78894000.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_320_t_84153600.000000_reference.vtu elder_python_pcs_0_ts_320_t_84153600.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_340_t_89413200.000000_reference.vtu elder_python_pcs_0_ts_340_t_89413200.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_360_t_94672800.000000_reference.vtu elder_python_pcs_0_ts_360_t_94672800.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_380_t_99932400.000000_reference.vtu elder_python_pcs_0_ts_380_t_99932400.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_400_t_105192000.000000_reference.vtu elder_python_pcs_0_ts_400_t_105192000.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_40_t_10519200.000000_reference.vtu elder_python_pcs_0_ts_40_t_10519200.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_420_t_110451600.000000_reference.vtu elder_python_pcs_0_ts_420_t_110451600.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_440_t_115711200.000000_reference.vtu elder_python_pcs_0_ts_440_t_115711200.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_460_t_120970800.000000_reference.vtu elder_python_pcs_0_ts_460_t_120970800.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_480_t_126230400.000000_reference.vtu elder_python_pcs_0_ts_480_t_126230400.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_500_t_131490000.000000_reference.vtu elder_python_pcs_0_ts_500_t_131490000.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_60_t_15778800.000000_reference.vtu elder_python_pcs_0_ts_60_t_15778800.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_80_t_21038400.000000_reference.vtu elder_python_pcs_0_ts_80_t_21038400.000000.vtu pressure pressure 1e-1 1e-5
+ elder_pcs_0_ts_0_t_0.000000_reference.vtu elder_python_pcs_0_ts_0_t_0.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_100_t_26298000.000000_reference.vtu elder_python_pcs_0_ts_100_t_26298000.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_120_t_31557600.000000_reference.vtu elder_python_pcs_0_ts_120_t_31557600.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_140_t_36817200.000000_reference.vtu elder_python_pcs_0_ts_140_t_36817200.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_160_t_42076800.000000_reference.vtu elder_python_pcs_0_ts_160_t_42076800.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_180_t_47336400.000000_reference.vtu elder_python_pcs_0_ts_180_t_47336400.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_200_t_52596000.000000_reference.vtu elder_python_pcs_0_ts_200_t_52596000.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_20_t_5259600.000000_reference.vtu elder_python_pcs_0_ts_20_t_5259600.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_220_t_57855600.000000_reference.vtu elder_python_pcs_0_ts_220_t_57855600.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_240_t_63115200.000000_reference.vtu elder_python_pcs_0_ts_240_t_63115200.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_260_t_68374800.000000_reference.vtu elder_python_pcs_0_ts_260_t_68374800.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_280_t_73634400.000000_reference.vtu elder_python_pcs_0_ts_280_t_73634400.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_300_t_78894000.000000_reference.vtu elder_python_pcs_0_ts_300_t_78894000.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_320_t_84153600.000000_reference.vtu elder_python_pcs_0_ts_320_t_84153600.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_340_t_89413200.000000_reference.vtu elder_python_pcs_0_ts_340_t_89413200.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_360_t_94672800.000000_reference.vtu elder_python_pcs_0_ts_360_t_94672800.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_380_t_99932400.000000_reference.vtu elder_python_pcs_0_ts_380_t_99932400.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_400_t_105192000.000000_reference.vtu elder_python_pcs_0_ts_400_t_105192000.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_40_t_10519200.000000_reference.vtu elder_python_pcs_0_ts_40_t_10519200.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_420_t_110451600.000000_reference.vtu elder_python_pcs_0_ts_420_t_110451600.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_440_t_115711200.000000_reference.vtu elder_python_pcs_0_ts_440_t_115711200.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_460_t_120970800.000000_reference.vtu elder_python_pcs_0_ts_460_t_120970800.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_480_t_126230400.000000_reference.vtu elder_python_pcs_0_ts_480_t_126230400.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_500_t_131490000.000000_reference.vtu elder_python_pcs_0_ts_500_t_131490000.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_60_t_15778800.000000_reference.vtu elder_python_pcs_0_ts_60_t_15778800.000000.vtu conc conc 1e-1 1e-5
+ elder_pcs_0_ts_80_t_21038400.000000_reference.vtu elder_python_pcs_0_ts_80_t_21038400.000000.vtu conc conc 1e-1 1e-5
+)
+
AddTest(
NAME 2D_ComponentTransport_HeterogeneousPermeability
PATH Elliptic/square_100x100_ComponentTransport
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 7debb42fcfbdc5fea92bdcaad9f7da432aa399b3..8e4ad70fe293efab5e160756b20911ab4b861ecd 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -187,7 +187,7 @@ void Process::assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
const auto pcs_id =
(_coupled_solutions) != nullptr ? _coupled_solutions->process_id : 0;
- _boundary_conditions[pcs_id].applyNaturalBC(t, x, K, b);
+ _boundary_conditions[pcs_id].applyNaturalBC(t, x, K, b, nullptr);
_source_term_collections[pcs_id].integrateNodalSourceTerms(t, b);
}
@@ -207,7 +207,7 @@ void Process::assembleWithJacobian(const double t, GlobalVector const& x,
// TODO: apply BCs to Jacobian.
const auto pcs_id =
(_coupled_solutions) != nullptr ? _coupled_solutions->process_id : 0;
- _boundary_conditions[pcs_id].applyNaturalBC(t, x, K, b);
+ _boundary_conditions[pcs_id].applyNaturalBC(t, x, K, b, &Jac);
}
void Process::constructDofTable()
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index f4f9e86b641d76f3e7b9ad520eeb3b6b90a67e0c..b887d36dc2360751d5c0a0e9b9bc297f1be1e2d8 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -16,13 +16,13 @@
#include "NumLib/ODESolver/ODESystem.h"
#include "NumLib/ODESolver/TimeDiscretization.h"
#include "ProcessLib/BoundaryCondition/BoundaryConditionCollection.h"
-#include "ProcessLib/SourceTerms/SourceTermCollection.h"
#include "ProcessLib/Output/CachedSecondaryVariable.h"
#include "ProcessLib/Output/ExtrapolatorData.h"
#include "ProcessLib/Output/IntegrationPointWriter.h"
#include "ProcessLib/Output/SecondaryVariable.h"
#include "ProcessLib/Parameter/Parameter.h"
#include "ProcessLib/SourceTerms/NodalSourceTerm.h"
+#include "ProcessLib/SourceTerms/SourceTermCollection.h"
#include "AbstractJacobianAssembler.h"
#include "ProcessVariable.h"
@@ -103,11 +103,11 @@ public:
GlobalMatrix& Jac) final;
std::vector<NumLib::IndexValueVector<GlobalIndexType>> const*
- getKnownSolutions(double const t) const final
+ getKnownSolutions(double const t, GlobalVector const& x) const final
{
const auto pcs_id =
(_coupled_solutions) ? _coupled_solutions->process_id : 0;
- return _boundary_conditions[pcs_id].getKnownSolutions(t);
+ return _boundary_conditions[pcs_id].getKnownSolutions(t, x);
}
virtual NumLib::LocalToGlobalIndexMap const& getDOFTable(
diff --git a/ProcessLib/SmallDeformation/Tests.cmake b/ProcessLib/SmallDeformation/Tests.cmake
index 75970c97317dc22b627eb76b4aa71988b328ad07..704d03ba4abcd8d6a29ac9e87d4bf97d304ec2de 100644
--- a/ProcessLib/SmallDeformation/Tests.cmake
+++ b/ProcessLib/SmallDeformation/Tests.cmake
@@ -299,3 +299,20 @@ AddTest(
expected_two_material_gravity.vtu two_material_gravity_pcs_0_ts_1_t_1.000000.vtu sigma sigma 5e-14 0
expected_two_material_gravity.vtu two_material_gravity_pcs_0_ts_1_t_1.000000.vtu epsilon epsilon 5e-14 0
)
+
+AddTest(
+ NAME PythonBCSmallDeformationPiston
+ PATH Mechanics/Linear/PythonPiston
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS piston.prj
+ WRAPPER time
+ TESTER vtkdiff
+ REQUIREMENTS OGS_USE_PYTHON AND NOT OGS_USE_MPI
+ DIFF_DATA
+ ref_piston_pcs_0_ts_5_t_5.000000.vtu piston_pcs_0_ts_5_t_5.000000.vtu displacement displacement 1e-16 0
+ ref_piston_pcs_0_ts_5_t_5.000000.vtu piston_pcs_0_ts_5_t_5.000000.vtu epsilon epsilon 1e-14 0
+ ref_piston_pcs_0_ts_5_t_5.000000.vtu piston_pcs_0_ts_5_t_5.000000.vtu sigma sigma 1e-8 0
+ ref_piston_pcs_0_ts_10_t_10.000000.vtu piston_pcs_0_ts_10_t_10.000000.vtu displacement displacement 1e-16 0
+ ref_piston_pcs_0_ts_10_t_10.000000.vtu piston_pcs_0_ts_10_t_10.000000.vtu epsilon epsilon 1e-13 0
+ ref_piston_pcs_0_ts_10_t_10.000000.vtu piston_pcs_0_ts_10_t_10.000000.vtu sigma sigma 1e-7 0
+)
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/bcs_laplace_eq.py b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/bcs_laplace_eq.py
new file mode 100644
index 0000000000000000000000000000000000000000..e6d5c8dd7d848239b73593327590af1705e8bc45
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/bcs_laplace_eq.py
@@ -0,0 +1,52 @@
+
+import OpenGeoSys
+from math import pi, sin, cos, sinh, cosh
+
+a = 2.0*pi/3.0
+
+# analytical solution used to set the Dirichlet BCs
+def solution(x, y):
+ return sin(a*x) * sinh(a*y)
+
+# gradient of the analytical solution used to set the Neumann BCs
+def grad_solution(x, y):
+ return a * cos(a*x) * sinh(a*y), \
+ a * sin(a*x) * cosh(a*y)
+
+# Dirichlet BCs
+class BCTop(OpenGeoSys.BoundaryCondition):
+ def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+ x, y, z = coords
+ assert y == 1.0 and z == 0.0
+ value = solution(x, y)
+ return (True, value)
+
+class BCLeft(OpenGeoSys.BoundaryCondition):
+ def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+ x, y, z = coords
+ assert x == 0.0 and z == 0.0
+ value = solution(x, y)
+ return (True, value)
+
+class BCBottom(OpenGeoSys.BoundaryCondition):
+ def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+ x, y, z = coords
+ assert y == 0.0 and z == 0.0
+ value = solution(x, y)
+ return (True, value)
+
+# Neumann BC
+class BCRight(OpenGeoSys.BoundaryCondition):
+ def getFlux(self, t, coords, primary_vars):
+ x, y, z = coords
+ assert x == 1.0 and z == 0.0
+ value = grad_solution(x, y)[0]
+ Jac = [ 0.0 ] # value does not depend on primary variable
+ return (True, value, Jac)
+
+
+# instantiate BC objects referenced in OpenGeoSys' prj file
+bc_top = BCTop()
+bc_right = BCRight()
+bc_bottom = BCBottom()
+bc_left = BCLeft()
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/python_laplace_eq_ref.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/python_laplace_eq_ref.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..ca30c5172424cb4841837fb1c5b3d6c4275dab5b
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/python_laplace_eq_ref.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5afedc3c451f400d58bdc8cefe7d5989bc1a248d80470131f750b2eb4793e49f
+size 23684
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1e3_laplace_eq.prj b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1e3_laplace_eq.prj
new file mode 100644
index 0000000000000000000000000000000000000000..5e196baa339be5aa1f184de851a357748fc1dae7
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1e3_laplace_eq.prj
@@ -0,0 +1,122 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh>square_1x1_quad_1e3.vtu</mesh>
+ <geometry>square_1x1.gml</geometry>
+ <python_script>bcs_laplace_eq.py</python_script>
+ <processes>
+ <process>
+ <name>GW23</name>
+ <type>GROUNDWATER_FLOW</type>
+ <integration_order>2</integration_order>
+ <hydraulic_conductivity>K</hydraulic_conductivity>
+ <process_variables>
+ <process_variable>pressure</process_variable>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable type="static" internal_name="darcy_velocity" output_name="v"/>
+ </secondary_variables>
+
+ <jacobian_assembler>
+ <type>CentralDifferences</type>
+ </jacobian_assembler>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="GW23">
+ <nonlinear_solver>basic_newton</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstol>1.e-6</abstol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <output>
+ <variables>
+ <variable> pressure </variable>
+ <variable> v </variable>
+ </variables>
+ </output>
+ <time_stepping>
+ <type>SingleStep</type>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>square_1e3_neumann</prefix>
+ </output>
+ </time_loop>
+ <parameters>
+ <parameter>
+ <name>K</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>zero</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ </parameters>
+ <process_variables>
+ <process_variable>
+ <name>pressure</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>zero</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>left</geometry>
+ <type>Python</type>
+ <bc_object>bc_left</bc_object>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>right</geometry>
+ <type>Python</type>
+ <bc_object>bc_right</bc_object>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>top</geometry>
+ <type>Python</type>
+ <bc_object>bc_top</bc_object>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>bottom</geometry>
+ <type>Python</type>
+ <bc_object>bc_bottom</bc_object>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_newton</name>
+ <type>Newton</type>
+ <max_iter>10</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+ <eigen>
+ <solver_type>CG</solver_type>
+ <precon_type>DIAGONAL</precon_type>
+ <max_iteration_step>10000</max_iteration_step>
+ <error_tolerance>1e-16</error_tolerance>
+ </eigen>
+ <petsc>
+ <prefix>gw</prefix>
+ <parameters>-gw_ksp_type cg -gw_pc_type bjacobi -gw_ksp_rtol 1e-16 -gw_ksp_max_it 10000</parameters>
+ </petsc>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1.gml b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1.gml
new file mode 100644
index 0000000000000000000000000000000000000000..3e5bb2a80a8f6bbdc4b7776b696d7a6af0d43b51
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1.gml
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a3f52cd5a8058d967f8ead873e1afa323af39266e538bcd2c35af08683a81dad
+size 1083
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..b56982101f6be9ed1cb17945723023e95d9b9f98
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:3653db864dd917e8d5be84d68ee2de13ba3485b1ef23de4e9894fb6749a1aa11
+size 25692
diff --git a/Tests/Data/Mechanics/Linear/PythonPiston/chamber.py b/Tests/Data/Mechanics/Linear/PythonPiston/chamber.py
new file mode 100644
index 0000000000000000000000000000000000000000..2fd26893826c383eaddfaa8f44cacb13e8fd4f22
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonPiston/chamber.py
@@ -0,0 +1,19 @@
+# model of the ideal gas inside the chamber
+
+from math import pi
+
+R_piston = 0.05
+
+h_chamber_0 = 0.1
+V_chamber_0 = pi * R_piston**2 * h_chamber_0
+
+p_chamber_0 = 1e5
+
+
+def p_chamber(displacement_y):
+ V_chamber = pi * R_piston**2 * (h_chamber_0 + displacement_y)
+ return p_chamber_0 * V_chamber_0 / V_chamber
+
+def dp_chamber_du_y(displacement_y):
+ return -p_chamber_0 * V_chamber_0 / pi / R_piston**2 / \
+ (h_chamber_0 + displacement_y)**2
diff --git a/Tests/Data/Mechanics/Linear/PythonPiston/piston.gml b/Tests/Data/Mechanics/Linear/PythonPiston/piston.gml
new file mode 100644
index 0000000000000000000000000000000000000000..087abc71cf6e7f1679749c59c2d1cc29eb1fbd35
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonPiston/piston.gml
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:41db3c3cdde67b2ec37281d99840e3fa63ddf6a0cf1059c3be68126ab78fe199
+size 740
diff --git a/Tests/Data/Mechanics/Linear/PythonPiston/piston.prj b/Tests/Data/Mechanics/Linear/PythonPiston/piston.prj
new file mode 100644
index 0000000000000000000000000000000000000000..474c0517076216aab1880c8a57abbe07653dd7a7
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonPiston/piston.prj
@@ -0,0 +1,152 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh axially_symmetric="true">piston.vtu</mesh>
+ <geometry>piston.gml</geometry>
+ <python_script>piston_bc.py</python_script>
+ <processes>
+ <process>
+ <name>SD</name>
+ <type>SMALL_DEFORMATION</type>
+ <integration_order>2</integration_order>
+ <constitutive_relation>
+ <type>LinearElasticIsotropic</type>
+ <youngs_modulus>E</youngs_modulus>
+ <poissons_ratio>nu</poissons_ratio>
+ </constitutive_relation>
+ <solid_density>rho_sr</solid_density>
+ <specific_body_force>0 0</specific_body_force>
+ <process_variables>
+ <process_variable>displacement</process_variable>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+ <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
+ </secondary_variables>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="SD">
+ <nonlinear_solver>basic_newton</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>INFINITY_N</norm_type>
+ <abstol>1e-8</abstol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <output>
+ <variables>
+ <variable>displacement</variable>
+ <variable>sigma</variable>
+ <variable>epsilon</variable>
+ </variables>
+ </output>
+ <time_stepping>
+ <type>FixedTimeStepping</type>
+ <t_initial>0</t_initial>
+ <t_end>10</t_end>
+ <timesteps>
+ <pair>
+ <repeat>1</repeat>
+ <delta_t>1</delta_t>
+ </pair>
+ </timesteps>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>piston</prefix>
+ <timesteps>
+ <pair>
+ <repeat>1</repeat>
+ <each_steps>1</each_steps>
+ </pair>
+ </timesteps>
+ <output_iteration_results>false</output_iteration_results>
+ </output>
+ </time_loop>
+ <parameters>
+ <parameter>
+ <name>E</name>
+ <type>Constant</type>
+ <value>1e6</value>
+ </parameter>
+ <parameter>
+ <name>nu</name>
+ <type>Constant</type>
+ <value>.3</value>
+ </parameter>
+ <parameter>
+ <name>rho_sr</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>displacement0</name>
+ <type>Constant</type>
+ <values>0 0</values>
+ </parameter>
+ <parameter>
+ <name>zero</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ </parameters>
+ <process_variables>
+ <process_variable>
+ <name>displacement</name>
+ <components>2</components>
+ <order>1</order>
+ <initial_condition>displacement0</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>geom</geometrical_set>
+ <!-- The Python BC is set on the whole domain boundary.
+ It is determined in the Python script where the BC is set specifically. -->
+ <geometry>whole_domain_boundary</geometry>
+ <type>Python</type>
+ <component>0</component>
+ <bc_object>bc_x</bc_object>
+ <flush_stdout>false</flush_stdout> <!-- for debugging only -->
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>geom</geometrical_set>
+ <!-- The Python BC is set on the whole domain boundary.
+ It is determined in the Python script where the BC is set specifically. -->
+ <geometry>whole_domain_boundary</geometry>
+ <type>Python</type>
+ <component>1</component>
+ <bc_object>bc_y</bc_object>
+ <flush_stdout>false</flush_stdout> <!-- for debugging only -->
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_newton</name>
+ <type>Newton</type>
+ <max_iter>50</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+ <eigen>
+ <solver_type>CG</solver_type>
+ <precon_type>DIAGONAL</precon_type>
+ <max_iteration_step>10000</max_iteration_step>
+ <error_tolerance>1e-16</error_tolerance>
+ </eigen>
+ <petsc>
+ <prefix>sd</prefix>
+ <parameters>-sd_ksp_type cg -sd_pc_type bjacobi -sd_ksp_rtol 1e-16 -sd_ksp_max_it 10000</parameters>
+ </petsc>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Mechanics/Linear/PythonPiston/piston.vtu b/Tests/Data/Mechanics/Linear/PythonPiston/piston.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..9613afb36e463881d9c2e448ffe7c85e7e13b8fa
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonPiston/piston.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:9f057b8c3eb5a6f8ed954a670e57cdb6c2c9b8e8acf48d7b545cea68b763b4d0
+size 45734
diff --git a/Tests/Data/Mechanics/Linear/PythonPiston/piston_bc.py b/Tests/Data/Mechanics/Linear/PythonPiston/piston_bc.py
new file mode 100644
index 0000000000000000000000000000000000000000..e06fbf5023832ab57a0b5cb24264cb00383c7eb0
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonPiston/piston_bc.py
@@ -0,0 +1,43 @@
+import OpenGeoSys
+from chamber import *
+
+
+class BC_X(OpenGeoSys.BoundaryCondition):
+ def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+ x, y, z = coords
+
+ if x == 0.0 or x == R_piston:
+ # no x displacement at outer boundary (x==R_piston) and center (x==0.0)
+ return (True, 0.0)
+
+ return (False, 0.0)
+
+
+class BC_Y(OpenGeoSys.BoundaryCondition):
+ def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+ x, y, z = coords
+
+ if y == 0.04:
+ # upper boundary -- prescribe displacement
+ return (True, -t * h_chamber_0 / 20.0)
+
+ return (False, 0.0)
+
+
+ def getFlux(self, t, coords, primary_vars):
+ x, y, z = coords
+
+ if y == 0.0:
+ # lower boundary -- the chamber
+ ux, uy = primary_vars
+
+ p = p_chamber(uy)
+ Jac = [ 0.0, dp_chamber_du_y(uy) ]
+ return (True, p, Jac)
+
+ return (False, 0.0, [ 0.0, 0.0 ])
+
+
+# instantiate the BC objects used by OpenGeoSys
+bc_x = BC_X()
+bc_y = BC_Y()
diff --git a/Tests/Data/Mechanics/Linear/PythonPiston/post.py b/Tests/Data/Mechanics/Linear/PythonPiston/post.py
new file mode 100755
index 0000000000000000000000000000000000000000..09695a52a3b1d042f740adc8c9567c47c6780073
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonPiston/post.py
@@ -0,0 +1,99 @@
+#!/usr/bin/vtkpython
+
+from __future__ import print_function
+import vtk
+from vtk.util.numpy_support import vtk_to_numpy, numpy_to_vtk
+import numpy as np
+from scipy.interpolate import interp1d
+
+import matplotlib.pyplot as plt
+import chamber as ch
+
+pvd_file = "out/piston_pcs_0.pvd"
+
+
+### helpers ##############################################
+
+import os
+try:
+ import xml.etree.cElementTree as ET
+except:
+ import xml.etree.ElementTree as ET
+
+
+def relpathfrom(origin, relpath):
+ if os.path.isabs(relpath):
+ return relpath
+ return os.path.join(origin, relpath)
+
+def read_pvd_file(fn):
+ try:
+ path = fn.name
+ except AttributeError:
+ path = fn
+ pathroot = os.path.dirname(path)
+ pvdtree = ET.parse(fn)
+ node = pvdtree.getroot()
+ if node.tag != "VTKFile": return None, None
+ children = list(node)
+ if len(children) != 1: return None, None
+ node = children[0]
+ if node.tag != "Collection": return None, None
+
+ ts = []
+ fs = []
+
+ for child in node:
+ if child.tag != "DataSet": return None, None
+ ts.append(float(child.get("timestep")))
+ fs.append(relpathfrom(pathroot, child.get("file")))
+
+ return ts, fs
+
+### helpers end ##########################################
+
+
+ts, fns = read_pvd_file(pvd_file)
+
+reader = vtk.vtkXMLUnstructuredGridReader()
+
+
+loc_points = vtk.vtkPoints()
+loc_points.InsertNextPoint([0.0, 0.0, 0.0])
+loc = vtk.vtkPolyData()
+loc.SetPoints(loc_points)
+
+probe = vtk.vtkProbeFilter()
+probe.SetSourceConnection(reader.GetOutputPort())
+probe.SetInputData(loc)
+
+uys = np.zeros(len(ts))
+ps = np.zeros(len(ts))
+
+for i, (t, fn) in enumerate(zip(ts, fns)):
+ print("###### time", t)
+ reader.SetFileName(fn)
+ probe.Update()
+
+ grid = probe.GetOutput()
+ uy = vtk_to_numpy(grid.GetPointData().GetArray("displacement"))[0,1]
+ p = vtk_to_numpy(grid.GetPointData().GetArray("sigma"))[0,1]
+
+ uys[i] = uy
+ ps[i] = -p
+
+uys_ana = np.linspace(min(uys), max(uys), 100)
+ps_ana = ch.p_chamber(uys_ana)
+
+
+fig, ax = plt.subplots()
+ax.scatter(uys[1:], ps[1:], label="ogs")
+ax.plot(uys_ana, ps_ana, label="ref")
+ax.legend()
+
+ax.set_ylabel("pressure inside the chamber $p$ / Pa")
+ax.set_xlabel("displacement of the piston $u_y$ / m")
+
+fig.subplots_adjust(left=0.15, right=0.98)
+fig.savefig("pressure-displacement.png")
+plt.close(fig)
diff --git a/Tests/Data/Mechanics/Linear/PythonPiston/ref_piston_pcs_0_ts_10_t_10.000000.vtu b/Tests/Data/Mechanics/Linear/PythonPiston/ref_piston_pcs_0_ts_10_t_10.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..662010e3f2ff2b5f5ad5a7c7508367020f3e3264
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonPiston/ref_piston_pcs_0_ts_10_t_10.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:03fcf36250d379fee936780f926f19d0b26a783a3ec30826fc7e80f2c48d957a
+size 88956
diff --git a/Tests/Data/Mechanics/Linear/PythonPiston/ref_piston_pcs_0_ts_5_t_5.000000.vtu b/Tests/Data/Mechanics/Linear/PythonPiston/ref_piston_pcs_0_ts_5_t_5.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..fe837edcca068f6307fdbf60467f4f144361a223
--- /dev/null
+++ b/Tests/Data/Mechanics/Linear/PythonPiston/ref_piston_pcs_0_ts_5_t_5.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d7319534828115a05183c89d15464714dd1ff0e3a29ae3dd4d861f763bae0a5a
+size 88869
diff --git a/Tests/Data/Parabolic/ComponentTransport/elder/elder-python-bcs.py b/Tests/Data/Parabolic/ComponentTransport/elder/elder-python-bcs.py
new file mode 100644
index 0000000000000000000000000000000000000000..3d815d21b7bc5486b97ffaf7c6e9bb7c147cf67d
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/elder/elder-python-bcs.py
@@ -0,0 +1,34 @@
+import OpenGeoSys
+
+
+class BCPressure(OpenGeoSys.BoundaryCondition):
+ def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+ x, y, z = coords
+
+ if x == -150 and z == 75:
+ # prescribe pressure of 0
+ return (True, 0.0)
+
+ # no Dirichlet BC
+ return (False, 0.0)
+
+
+class BCConcentration(OpenGeoSys.BoundaryCondition):
+ def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+ x, y, z = coords
+
+ if z == -75:
+ # prescribe concentration of 0
+ return (True, 0.0)
+
+ if z == 75 and x >= 0:
+ # prescribe concentration of 1
+ return (True, 1.0)
+
+ # no Dirichlet BC
+ return (False, 0.0)
+
+
+# instantiate the BC objects used by OpenGeoSys
+bc_p = BCPressure()
+bc_c = BCConcentration()
diff --git a/Tests/Data/Parabolic/ComponentTransport/elder/elder-python.gml b/Tests/Data/Parabolic/ComponentTransport/elder/elder-python.gml
new file mode 100644
index 0000000000000000000000000000000000000000..d22068fc2b35f4f29064dcf89bd42163dcc610b2
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/elder/elder-python.gml
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:b4fb7c7bc4ddbc37a38f1fdbaabc23eb415fe842ae2dbd4a5c9791fa7f93dc83
+size 1282
diff --git a/Tests/Data/Parabolic/ComponentTransport/elder/elder-python.prj b/Tests/Data/Parabolic/ComponentTransport/elder/elder-python.prj
new file mode 100644
index 0000000000000000000000000000000000000000..2132e3beecfa76ea1d095e962cca295ae1b376a2
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/elder/elder-python.prj
@@ -0,0 +1,201 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<!-- vim: ft=xml sw=2
+-->
+<OpenGeoSysProject>
+ <mesh>elder.vtu</mesh>
+ <geometry>elder-python.gml</geometry>
+ <python_script>elder-python-bcs.py</python_script>
+ <processes>
+ <process>
+ <name>HC</name>
+ <type>ComponentTransport</type>
+ <integration_order>2</integration_order>
+ <process_variables>
+ <concentration>conc</concentration>
+ <pressure>pressure</pressure>
+ </process_variables>
+ <fluid>
+ <density>
+ <type>ConcentrationDependent</type>
+ <reference_density>1000</reference_density>
+ <reference_concentration>0</reference_concentration>
+ <fluid_density_difference_ratio>0.2</fluid_density_difference_ratio>
+ </density>
+ <viscosity>
+ <type>Constant</type>
+ <value>1.0e-3</value>
+ </viscosity>
+ </fluid>
+ <porous_medium>
+ <porous_medium id="0">
+ <permeability>
+ <permeability_tensor_entries>kappa1</permeability_tensor_entries>
+ <type>Constant</type>
+ </permeability>
+ <porosity>
+ <type>Constant</type>
+ <porosity_parameter>constant_porosity_parameter</porosity_parameter>
+ </porosity>
+ <storage>
+ <type>Constant</type>
+ <value>0.0</value>
+ </storage>
+ </porous_medium>
+ </porous_medium>
+ <decay_rate>decay_rate</decay_rate>
+ <fluid_reference_density>rho_fluid</fluid_reference_density>
+ <retardation_factor>retardation_factor</retardation_factor>
+ <solute_dispersivity_longitudinal>alpha_l</solute_dispersivity_longitudinal>
+ <solute_dispersivity_transverse>alpha_t</solute_dispersivity_transverse>
+ <molecular_diffusion_coefficient>Dm</molecular_diffusion_coefficient>
+ <specific_body_force>0 0 -9.81</specific_body_force>
+ <secondary_variables>
+ <secondary_variable type="static" internal_name="darcy_velocity" output_name="darcy_velocity"/>
+ </secondary_variables>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="HC">
+ <nonlinear_solver>basic_picard</nonlinear_solver>
+ <convergence_criterion>
+ <type>PerComponentDeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <reltols>1e-3 1e-3</reltols>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <time_stepping>
+ <type>FixedTimeStepping</type>
+ <t_initial> 0.0 </t_initial>
+ <t_end>1.3149e8</t_end>
+ <timesteps>
+ <pair>
+ <repeat>10</repeat>
+ <delta_t>262980</delta_t>
+ </pair>
+ </timesteps>
+ </time_stepping>
+ <output>
+ <variables>
+ <variable>conc</variable>
+ <variable>pressure</variable>
+ <variable>darcy_velocity</variable>
+ </variables>
+ </output>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>elder_python</prefix>
+ <timesteps>
+ <pair>
+ <repeat>1</repeat>
+ <each_steps>20</each_steps>
+ </pair>
+ </timesteps>
+ </output>
+ </time_loop>
+ <parameters>
+ <parameter>
+ <name>rho_fluid</name>
+ <type>Constant</type>
+ <value>1000.0</value>
+ </parameter>
+ <parameter>
+ <name>alpha_l</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>alpha_t</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>Dm</name>
+ <type>Constant</type>
+ <value>3.57e-6</value>
+ </parameter>
+ <parameter>
+ <name>decay_rate</name>
+ <type>Constant</type>
+ <value>0.0</value>
+ </parameter>
+ <parameter>
+ <name>retardation_factor</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>p_ini</name>
+ <type>MeshNode</type>
+ <field_name>p_ini</field_name>
+ </parameter>
+ <parameter>
+ <name>c_ini</name>
+ <type>MeshNode</type>
+ <field_name>c_ini</field_name>
+ </parameter>
+ <parameter>
+ <name>constant_porosity_parameter</name>
+ <type>Constant</type>
+ <value>0.1</value>
+ </parameter>
+ <parameter>
+ <name>kappa1</name>
+ <type>Constant</type>
+ <values>4.84404e-13 0 0 0 4.84404e-13 0 0 0 4.84404e-13</values>
+ </parameter>
+ </parameters>
+ <process_variables>
+ <process_variable>
+ <name>conc</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>c_ini</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>elder</geometrical_set>
+ <geometry>whole_domain_boundary</geometry>
+ <type>Python</type>
+ <bc_object>bc_c</bc_object>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ <process_variable>
+ <name>pressure</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>p_ini</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>elder</geometrical_set>
+ <geometry>whole_domain_boundary</geometry>
+ <type>Python</type>
+ <bc_object>bc_p</bc_object>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_picard</name>
+ <type>Picard</type>
+ <max_iter>300</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <eigen>
+ <precon_type>ILUT</precon_type>
+ <solver_type>BiCGSTAB</solver_type>
+ <max_iteration_step>5000</max_iteration_step>
+ <error_tolerance>1e-12</error_tolerance>
+ </eigen>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/ThirdParty/CMakeLists.txt b/ThirdParty/CMakeLists.txt
index ea686f5ee2fa134513f1f41aa989f790cf94806b..2cdae58ab87fce49132f383bc37f7bda4aa3fea1 100644
--- a/ThirdParty/CMakeLists.txt
+++ b/ThirdParty/CMakeLists.txt
@@ -45,6 +45,37 @@ if(OGS_BUILD_SWMM)
set_target_properties(SWMM SwmmInterface PROPERTIES COMPILE_FLAGS /W0)
endif()
+if(OGS_USE_PYTHON)
+ add_subdirectory(pybind11)
+
+ function(check_python_version_compatibility)
+ if(NOT EXISTS ${VTK_DIR}/VTKConfig.cmake)
+ return()
+ endif()
+
+ include(${VTK_DIR}/VTKConfig.cmake)
+
+ if(NOT ${VTK_WRAP_PYTHON})
+ return()
+ endif()
+ if(NOT EXISTS ${VTK_MODULES_DIR}/vtkPython.cmake)
+ return()
+ endif()
+
+ include(${VTK_MODULES_DIR}/vtkPython.cmake)
+
+ if (NOT "${vtkPython_LIBRARIES}" STREQUAL "${PYTHON_LIBRARIES}")
+ message(SEND_ERROR "Mismatch between VTK's and OpenGeoSys' Python "
+ "libraries: ${vtkPython_LIBRARIES} vs. ${PYTHON_LIBRARIES}. "
+ "This will lead to compilation or linking errors. "
+ "You can fix this error by using the same Python version for "
+ "OpenGeoSys as VTK is built with.")
+ endif()
+ endfunction()
+
+ check_python_version_compatibility()
+endif()
+
if(OGS_BUILD_UTILS)
include(${PROJECT_SOURCE_DIR}/scripts/cmake/MetisSetup.cmake)
endif()
diff --git a/ThirdParty/pybind11 b/ThirdParty/pybind11
new file mode 160000
index 0000000000000000000000000000000000000000..8edc147d67ca85a93ed1f53628004528dc36a04d
--- /dev/null
+++ b/ThirdParty/pybind11
@@ -0,0 +1 @@
+Subproject commit 8edc147d67ca85a93ed1f53628004528dc36a04d
diff --git a/netlify.toml b/netlify.toml
index 87d97a136355f0b8fa7a8127a56a65762870985f..286b6ae44489035f8e6e2fbef1410a64ae108333 100644
--- a/netlify.toml
+++ b/netlify.toml
@@ -9,7 +9,7 @@
hugo"""
[context.production.environment]
- HUGO_VERSION = "0.42.1"
+ HUGO_VERSION = "0.47"
HUGO_ENV = "production"
HUGO_ENABLEGITINFO = "true"
@@ -21,7 +21,7 @@
hugo -b $DEPLOY_PRIME_URL"""
[context.deploy-preview.environment]
- HUGO_VERSION = "0.42.1"
+ HUGO_VERSION = "0.47"
[context.branch-deploy]
command = """
@@ -32,4 +32,4 @@
node_modules/.bin/hugo-algolia --toml -s"""
[context.branch-deploy.environment]
- HUGO_VERSION = "0.42.1"
+ HUGO_VERSION = "0.47"
diff --git a/scripts/cmake/SubmoduleSetup.cmake b/scripts/cmake/SubmoduleSetup.cmake
index 09d8d084c1f9ae92197c939cbc2f69bab8629af3..a07dbb57ab1c8ee4c946f65b8a03fb3cc0c2e5a5 100644
--- a/scripts/cmake/SubmoduleSetup.cmake
+++ b/scripts/cmake/SubmoduleSetup.cmake
@@ -19,6 +19,9 @@ endif()
if(OGS_BUILD_SWMM)
list(APPEND REQUIRED_SUBMODULES ThirdParty/SwmmInterface)
endif()
+if(OGS_USE_PYTHON)
+ list(APPEND REQUIRED_SUBMODULES ThirdParty/pybind11)
+endif()
# Sync submodules, which is required when a submodule changed its URL
if(OGS_SYNC_SUBMODULES)
diff --git a/web/config.toml b/web/config.toml
index 29f2e488c2674f2ed39c7f89404a7ca2353c2d8a..32f601f5c16f816740afcade65298e8130c9a630 100644
--- a/web/config.toml
+++ b/web/config.toml
@@ -93,6 +93,10 @@ staticDir = ["dist", "static"]
name = "Two-phase Flow"
identifier = "two-phase-flow"
weight = 8
+[[menu.benchmarks]]
+ name = "Python Boundary Conditions"
+ identifier = "python-bc"
+ weight = 9
# Quickstart sidebar top-level categories
[[menu.quickstart]]
diff --git a/web/content/docs/benchmarks/python-bc/elder-benchmark/elder.md b/web/content/docs/benchmarks/python-bc/elder-benchmark/elder.md
new file mode 100644
index 0000000000000000000000000000000000000000..04b5ba1dde5c0fe7fb7f3219e34caba7d9303309
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/elder-benchmark/elder.md
@@ -0,0 +1,27 @@
++++
+project = "Parabolic/ComponentTransport/elder/elder-python.prj"
+author = "Christoph Lehmann"
+date = "2018-08-16T09:18:00+02:00"
+title = "Saturated Variable-Density Flow and Mass Transport (Elder) with Python BC"
+weight = 3
+
+[menu]
+ [menu.benchmarks]
+ parent = "python-bc"
+
++++
+
+{{< data-link >}}
+
+## Motivation of this test case
+
+The aim of this test is:
+
+* to show that it is possible to prescribe BCs only on parts of a given geometry
+* to assert that Python BCs work with processes involving more than one physical
+ field.
+
+## Details
+
+This test is a copy of [this test case]({{< ref "../../hydro-component/elder.md" >}}).
+Please check the original test case for any details.
diff --git a/web/content/docs/benchmarks/python-bc/laplace-equation/python-laplace-eq.md b/web/content/docs/benchmarks/python-bc/laplace-equation/python-laplace-eq.md
new file mode 100644
index 0000000000000000000000000000000000000000..0e68c1d6c6a1e218ced82b290e056b4c2fe0f5ba
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/laplace-equation/python-laplace-eq.md
@@ -0,0 +1,95 @@
++++
+project = "https://github.com/ufz/ogs-data/blob/master/Elliptic/square_1x1_GroundWaterFlow_Python/square_1e3_laplace_eq.prj"
+author = "Christoph Lehmann"
+date = "2018-06-01T14:16:55+02:00"
+title = "Manufactured Solution for Laplace's Equation with Python"
+weight = 1
+
+[menu]
+ [menu.benchmarks]
+ parent = "python-bc"
+
++++
+
+{{< data-link >}}
+
+## Motivation of this test case
+
+The aim of this test is:
+
+* to provide a simple introductory example of how Python BCs can be used
+* to show that Python BCs can be used to easily prescribe BCs based on
+ analytical formulas
+* to check whether both essential and natural BCs
+ are implemented correctly in OpenGeoSys' Python BC.
+
+## Problem description
+
+We solve Laplace's Equation in 2D on a $1 \times 1$ square domain.
+ The boundary conditions will be specified below.
+
+## Weak form
+
+Laplace's equation is
+$$
+\begin{equation}
+- \mathop{\mathrm{div}} (a \mathop{\mathrm{grad}} u) = 0
+\end{equation}
+$$
+The weak form is derived as usual by multiplying with a test function $v$ and
+integrating over the domain $\Omega$:
+$$
+\begin{equation}
+- \int\_{\Omega} v \mathop{\mathrm{div}} (a \mathop{\mathrm{grad}} u) \, \mathrm{d}\Omega = 0
+\,,
+\end{equation}
+$$
+which can be transformed further to
+$$
+\begin{equation}
+\int\_{\Omega} a \mathop{\mathrm{grad}} v \cdot \mathop{\mathrm{grad}} u \, \mathrm{d}\Omega = \int\_{\Omega} \mathop{\mathrm{div}} (v a \mathop{\mathrm{grad}} u) \, \mathrm{d}\Omega = \int\_{\Gamma\_{\mathrm{N}}} v a \mathop{\mathrm{grad}} u \cdot n \, \mathrm{d}\Gamma \,,
+\end{equation}
+$$
+where in the second equality Gauss's theorem has been applied.
+As usual, the domain boundary $\partial\Omega = \Gamma\_{\mathrm{D}} \cup \Gamma\_{\mathrm{N}}$ is subdivided
+into the dirichlet and the Neumann boundary and $v$ vanishes on
+$\Gamma\_{\mathrm{D}}$.
+The r.h.s. of the above equation is the total flux associated with $u$ flowing
+**into** the domain $\Omega$ through $\Gamma\_{\mathrm{N}}$:
+$-a \mathop{\mathrm{grad}} u$ is the flux density and $-n$ is the inwards directed surface
+normal.
+
+The weak form just derived is implemented (after FEM discretization) in the
+groundwater flow process in OpenGeoSys.
+Note that for the application of Neumann boundary conditions, it is necessary to
+know whether the flux has to be applied with a positive or a negative sign!
+
+## Analytial solution
+
+The coefficient $a$ of Laplace's equation is taken to be unity.
+By differentiation it can be easily checked that
+$$
+\begin{equation}
+u(x, y) = \sin(bx) \sinh(by)
+\end{equation}
+$$
+solves Laplace's equation inside $\Omega$ for any $b$.
+In this example we set $b = \tfrac 23 \pi$.
+
+As boundary conditions we apply Dirichlet BCs at the top, left and bottom of the
+domain with values from $u(x,y)|\_{\Gamma\_{\mathrm{D}}}$.
+On the right boundary of the domain a Neumann BC is applied.
+There $n = (1, 0)$, which implies that $a \mathop{\mathrm{grad}} u \cdot n
+= a \, \partial u / \partial x$.
+
+
+## Results
+
+The numerical result obtained from OpenGeoSys is:
+
+{{< img src="../python_laplace_eq_solution.png" >}}
+
+The absolute difference between the analytical and numerical solutions is
+smaller than $4 \cdot 10^{-4}$:
+
+{{< img src="../python_laplace_eq_diff.png" >}}
diff --git a/web/content/docs/benchmarks/python-bc/laplace-equation/python_laplace_eq_diff.png b/web/content/docs/benchmarks/python-bc/laplace-equation/python_laplace_eq_diff.png
new file mode 100644
index 0000000000000000000000000000000000000000..5c4d28f5451b43ec2d55f2ad98dc69388d216e75
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/laplace-equation/python_laplace_eq_diff.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:37b86d2f172fa864d6bfecabf99bd203ddaded8210808956325b4a97b3020f25
+size 34437
diff --git a/web/content/docs/benchmarks/python-bc/laplace-equation/python_laplace_eq_solution.png b/web/content/docs/benchmarks/python-bc/laplace-equation/python_laplace_eq_solution.png
new file mode 100644
index 0000000000000000000000000000000000000000..ad1bab5c5da68cb72014ef2e44fa6c08dea1a18e
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/laplace-equation/python_laplace_eq_solution.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d6aa873d3c8cb0f6b4e7cbd30cf32c07d1ae261e313d503a15d6e505902626fc
+size 39593
diff --git a/web/content/docs/benchmarks/python-bc/piston/load-steps.png b/web/content/docs/benchmarks/python-bc/piston/load-steps.png
new file mode 100644
index 0000000000000000000000000000000000000000..c86f2e8257628ea21ca2785128da08f174346bd4
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/piston/load-steps.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:b287e1f87187e948495d595c4d1a1d0b0aaf2f24975d965c251883180b85cd86
+size 24953
diff --git a/web/content/docs/benchmarks/python-bc/piston/piston.md b/web/content/docs/benchmarks/python-bc/piston/piston.md
new file mode 100644
index 0000000000000000000000000000000000000000..c0894d69265c9355dadd3b6422a2ca4d72ebe0f8
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/piston/piston.md
@@ -0,0 +1,48 @@
++++
+project = "Mechanics/Linear/PythonPiston/piston.prj"
+author = "Christoph Lehmann"
+date = "2018-08-16T12:18:00+02:00"
+title = "Ideal Gas Compressed by an Elastic Piston"
+weight = 2
+
+[menu]
+ [menu.benchmarks]
+ parent = "python-bc"
+
++++
+
+{{< data-link >}}
+
+## Motivation of this test case
+
+The aim of this test is:
+
+* to show that it is possible to prescribe BCs only on parts of a given geometry
+* to assert that the computation of fluxes that depend nonlinearly on the
+ primary variables is implemented correctly.
+
+## Details
+
+A chamber filled with ideal gas is sealed tightly with a movable, elastic
+piston. The position of the piston is varied between different load steps.
+Friction between the piston and the chamber wall is neglected.
+For simplicitly, also initially the elastic piston is in an unstressed state.
+
+{{<img src="../sketch-piston.png" >}}
+
+
+## Results
+
+{{<img src="../load-steps.png" >}}
+
+The figure above shows that the piston is being compressed
+($y$ displacement has larger negative values at the top)
+by the forces acting on it.
+The initial position of the top part of the piston is indicated as a wireframe.
+
+{{<img src="../pressure-displacement.png" >}}
+
+The plot shows that the relation between the stress in the piston and its
+displacement coincides with the pressure-volume relation of the chamber.
+That indicates that the Jacobian computation inside the Python BC classes was
+correctly implemented.
diff --git a/web/content/docs/benchmarks/python-bc/piston/pressure-displacement.png b/web/content/docs/benchmarks/python-bc/piston/pressure-displacement.png
new file mode 100644
index 0000000000000000000000000000000000000000..a666eb785b20faa51632f4b43ae48177b1607b81
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/piston/pressure-displacement.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5f531b45437626649d49c146c830868dad04d6ce51385fb1493f0736de1374bf
+size 25290
diff --git a/web/content/docs/benchmarks/python-bc/piston/sketch-piston.png b/web/content/docs/benchmarks/python-bc/piston/sketch-piston.png
new file mode 100644
index 0000000000000000000000000000000000000000..1d047c42cab5f708e22d38b81dc257c91044d8a9
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/piston/sketch-piston.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:bd7a9e1658f925dff4f8da98abc106541b61524d9f97d914bfb965ec7cce80de
+size 22296
diff --git a/web/content/docs/benchmarks/python-bc/piston/sketch-piston.svg b/web/content/docs/benchmarks/python-bc/piston/sketch-piston.svg
new file mode 100644
index 0000000000000000000000000000000000000000..8dbfdaa991928a121eb5ec5d2ce98ed8152c46e2
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/piston/sketch-piston.svg
@@ -0,0 +1,589 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<!-- Created with Inkscape (http://www.inkscape.org/) -->
+
+<svg
+ xmlns:dc="http://purl.org/dc/elements/1.1/"
+ xmlns:cc="http://creativecommons.org/ns#"
+ xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
+ xmlns:svg="http://www.w3.org/2000/svg"
+ xmlns="http://www.w3.org/2000/svg"
+ xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
+ xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
+ width="359.36136"
+ height="539.42145"
+ viewBox="0 0 95.081026 142.72192"
+ version="1.1"
+ id="svg8"
+ inkscape:version="0.92.2 2405546, 2018-03-11"
+ sodipodi:docname="sketch-piston.svg">
+ <defs
+ id="defs2">
+ <marker
+ inkscape:isstock="true"
+ style="overflow:visible"
+ id="marker5183"
+ refX="0"
+ refY="0"
+ orient="auto"
+ inkscape:stockid="Arrow1Mstart">
+ <path
+ transform="matrix(0.4,0,0,0.4,4,0)"
+ style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.00000003pt;stroke-opacity:1"
+ d="M 0,0 5,-5 -12.5,0 5,5 Z"
+ id="path5181"
+ inkscape:connector-curvature="0" />
+ </marker>
+ <marker
+ inkscape:stockid="Arrow1Mend"
+ orient="auto"
+ refY="0"
+ refX="0"
+ id="Arrow1Mend"
+ style="overflow:visible"
+ inkscape:isstock="true">
+ <path
+ id="path4533"
+ d="M 0,0 5,-5 -12.5,0 5,5 Z"
+ style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.00000003pt;stroke-opacity:1"
+ transform="matrix(-0.4,0,0,-0.4,-4,0)"
+ inkscape:connector-curvature="0" />
+ </marker>
+ </defs>
+ <sodipodi:namedview
+ id="base"
+ pagecolor="#ffffff"
+ bordercolor="#666666"
+ borderopacity="1.0"
+ inkscape:pageopacity="0.0"
+ inkscape:pageshadow="2"
+ inkscape:zoom="0.7"
+ inkscape:cx="305.61598"
+ inkscape:cy="278.35678"
+ inkscape:document-units="mm"
+ inkscape:current-layer="layer1"
+ showgrid="false"
+ units="px"
+ inkscape:window-width="1600"
+ inkscape:window-height="875"
+ inkscape:window-x="0"
+ inkscape:window-y="25"
+ inkscape:window-maximized="1"
+ fit-margin-top="0"
+ fit-margin-left="0"
+ fit-margin-right="0"
+ fit-margin-bottom="0" />
+ <metadata
+ id="metadata5">
+ <rdf:RDF>
+ <cc:Work
+ rdf:about="">
+ <dc:format>image/svg+xml</dc:format>
+ <dc:type
+ rdf:resource="http://purl.org/dc/dcmitype/StillImage" />
+ <dc:title></dc:title>
+ </cc:Work>
+ </rdf:RDF>
+ </metadata>
+ <g
+ inkscape:label="Ebene 1"
+ inkscape:groupmode="layer"
+ id="layer1"
+ transform="translate(-25.589826,-142.51819)">
+ <rect
+ style="opacity:1;fill:#e6e6e6;fill-opacity:1;fill-rule:evenodd;stroke:none;stroke-width:0.26458335;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ id="rect6468"
+ width="85.963516"
+ height="60.540501"
+ x="32.922787"
+ y="218.64417" />
+ <path
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:0.26458332;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ d="M 120.53856,153.209 V 280.69656 H 31.270535 V 153.209"
+ id="rect4518"
+ inkscape:connector-curvature="0"
+ sodipodi:nodetypes="cccc" />
+ <rect
+ y="193.83394"
+ x="32.874153"
+ height="24.054256"
+ width="42.898102"
+ id="rect5168"
+ style="opacity:1;fill:#cccccc;fill-opacity:1;fill-rule:evenodd;stroke:none;stroke-width:0.26458332;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1" />
+ <rect
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:0.26458332;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ id="rect4520"
+ width="86.060783"
+ height="24.054256"
+ x="32.874153"
+ y="193.83394" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:0.26458332;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:0.26458332, 0.52916663;stroke-dashoffset:0;stroke-opacity:1"
+ d="M 75.904545,142.51819 V 285.24012"
+ id="path5170"
+ inkscape:connector-curvature="0" />
+ <path
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:0.26458332;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1;marker-start:url(#marker5183)"
+ id="path5173"
+ sodipodi:type="arc"
+ sodipodi:cx="75.903915"
+ sodipodi:cy="150.36656"
+ sodipodi:rx="10.824416"
+ sodipodi:ry="4.0090427"
+ sodipodi:start="5.1836395"
+ sodipodi:end="4.342102"
+ d="m 80.81821,146.7945 a 10.824416,4.0090427 0 0 1 5.543899,4.60606 10.824416,4.0090427 0 0 1 -11.004397,2.96994 10.824416,4.0090427 0 0 1 -10.140123,-3.36561 10.824416,4.0090427 0 0 1 6.769154,-4.37565"
+ sodipodi:open="true" />
+ <g
+ id="g6322"
+ transform="matrix(0,-0.24041841,0.23377061,0,45.747789,206.26434)"
+ style="stroke-width:4.21814156">
+ <path
+ transform="matrix(0.73936612,0,0,0.82399862,-6.1239083,26.230728)"
+ inkscape:transform-center-y="-3.1008631"
+ d="m -11.910359,148.04422 -26.072139,0 13.03607,-22.57914 z"
+ inkscape:randomized="0"
+ inkscape:rounded="0"
+ inkscape:flatsided="true"
+ sodipodi:arg2="1.5707963"
+ sodipodi:arg1="0.52359878"
+ sodipodi:r2="7.5263782"
+ sodipodi:r1="15.052756"
+ sodipodi:cy="140.51784"
+ sodipodi:cx="-24.946428"
+ sodipodi:sides="3"
+ id="path6290"
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.42985046;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ sodipodi:type="star" />
+ <circle
+ r="3.7797618"
+ cy="152.27068"
+ cx="-29.453543"
+ id="path6292"
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1" />
+ <circle
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ id="circle6294"
+ cx="-19.683363"
+ cy="152.27068"
+ r="3.7797618" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6296"
+ d="m -34.423952,156.31501 h 19.710999"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6298"
+ d="m -16.68344,156.45039 -1.992943,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -19.734346,156.45039 -1.992942,2.82356"
+ id="path6300"
+ inkscape:connector-curvature="0" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6302"
+ d="m -22.785251,156.45039 -1.992942,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -25.836157,156.45039 -1.992942,2.82356"
+ id="path6304"
+ inkscape:connector-curvature="0" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6306"
+ d="m -28.887062,156.45039 -1.992942,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -31.937966,156.45039 -1.992942,2.82356"
+ id="path6308"
+ inkscape:connector-curvature="0" />
+ </g>
+ <g
+ style="stroke-width:4.21814156"
+ transform="matrix(0,-0.24041841,-0.23377061,0,62.898619,194.68627)"
+ id="g6344">
+ <path
+ sodipodi:type="star"
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.42985046;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ id="path6324"
+ sodipodi:sides="3"
+ sodipodi:cx="-24.946428"
+ sodipodi:cy="140.51784"
+ sodipodi:r1="15.052756"
+ sodipodi:r2="7.5263782"
+ sodipodi:arg1="0.52359878"
+ sodipodi:arg2="1.5707963"
+ inkscape:flatsided="true"
+ inkscape:rounded="0"
+ inkscape:randomized="0"
+ d="m -11.910359,148.04422 -26.072139,0 13.03607,-22.57914 z"
+ inkscape:transform-center-y="-3.1008631"
+ transform="matrix(0.73936612,0,0,0.82399862,-6.1239083,26.230728)" />
+ <circle
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ id="circle6326"
+ cx="-29.453543"
+ cy="152.27068"
+ r="3.7797618" />
+ <circle
+ r="3.7797618"
+ cy="152.27068"
+ cx="-19.683363"
+ id="circle6328"
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -34.423952,156.31501 h 19.710999"
+ id="path6330"
+ inkscape:connector-curvature="0" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -16.68344,156.45039 -1.992943,2.82356"
+ id="path6332"
+ inkscape:connector-curvature="0" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6334"
+ d="m -19.734346,156.45039 -1.992942,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -22.785251,156.45039 -1.992942,2.82356"
+ id="path6336"
+ inkscape:connector-curvature="0" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6338"
+ d="m -25.836157,156.45039 -1.992942,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -28.887062,156.45039 -1.992942,2.82356"
+ id="path6340"
+ inkscape:connector-curvature="0" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6342"
+ d="m -31.937966,156.45039 -1.992942,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ </g>
+ <g
+ style="stroke-width:4.21814156"
+ transform="matrix(0,-0.24041841,0.23377061,0,45.747789,194.68627)"
+ id="g6366">
+ <path
+ sodipodi:type="star"
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.42985046;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ id="path6346"
+ sodipodi:sides="3"
+ sodipodi:cx="-24.946428"
+ sodipodi:cy="140.51784"
+ sodipodi:r1="15.052756"
+ sodipodi:r2="7.5263782"
+ sodipodi:arg1="0.52359878"
+ sodipodi:arg2="1.5707963"
+ inkscape:flatsided="true"
+ inkscape:rounded="0"
+ inkscape:randomized="0"
+ d="m -11.910359,148.04422 -26.072139,0 13.03607,-22.57914 z"
+ inkscape:transform-center-y="-3.1008631"
+ transform="matrix(0.73936612,0,0,0.82399862,-6.1239083,26.230728)" />
+ <circle
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ id="circle6348"
+ cx="-29.453543"
+ cy="152.27068"
+ r="3.7797618" />
+ <circle
+ r="3.7797618"
+ cy="152.27068"
+ cx="-19.683363"
+ id="circle6350"
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -34.423952,156.31501 h 19.710999"
+ id="path6352"
+ inkscape:connector-curvature="0" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -16.68344,156.45039 -1.992943,2.82356"
+ id="path6354"
+ inkscape:connector-curvature="0" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6356"
+ d="m -19.734346,156.45039 -1.992942,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -22.785251,156.45039 -1.992942,2.82356"
+ id="path6358"
+ inkscape:connector-curvature="0" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6360"
+ d="m -25.836157,156.45039 -1.992942,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -28.887062,156.45039 -1.992942,2.82356"
+ id="path6362"
+ inkscape:connector-curvature="0" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6364"
+ d="m -31.937966,156.45039 -1.992942,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ </g>
+ <g
+ id="g6388"
+ transform="matrix(0,-0.24041841,-0.23377061,0,62.898619,206.26434)"
+ style="stroke-width:4.21814156">
+ <path
+ transform="matrix(0.73936612,0,0,0.82399862,-6.1239083,26.230728)"
+ inkscape:transform-center-y="-3.1008631"
+ d="m -11.910359,148.04422 -26.072139,0 13.03607,-22.57914 z"
+ inkscape:randomized="0"
+ inkscape:rounded="0"
+ inkscape:flatsided="true"
+ sodipodi:arg2="1.5707963"
+ sodipodi:arg1="0.52359878"
+ sodipodi:r2="7.5263782"
+ sodipodi:r1="15.052756"
+ sodipodi:cy="140.51784"
+ sodipodi:cx="-24.946428"
+ sodipodi:sides="3"
+ id="path6368"
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.42985046;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ sodipodi:type="star" />
+ <circle
+ r="3.7797618"
+ cy="152.27068"
+ cx="-29.453543"
+ id="circle6370"
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1" />
+ <circle
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ id="circle6372"
+ cx="-19.683363"
+ cy="152.27068"
+ r="3.7797618" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6374"
+ d="m -34.423952,156.31501 h 19.710999"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6376"
+ d="m -16.68344,156.45039 -1.992943,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -19.734346,156.45039 -1.992942,2.82356"
+ id="path6378"
+ inkscape:connector-curvature="0" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6380"
+ d="m -22.785251,156.45039 -1.992942,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -25.836157,156.45039 -1.992942,2.82356"
+ id="path6382"
+ inkscape:connector-curvature="0" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6384"
+ d="m -28.887062,156.45039 -1.992942,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -31.937966,156.45039 -1.992942,2.82356"
+ id="path6386"
+ inkscape:connector-curvature="0" />
+ </g>
+ <g
+ id="g6410"
+ transform="matrix(0.24041841,0,0,-0.23377061,49.374542,223.8584)"
+ style="stroke-width:4.21814156">
+ <path
+ transform="matrix(0.73936612,0,0,0.82399862,-6.1239083,26.230728)"
+ inkscape:transform-center-y="-3.1008631"
+ d="m -11.910359,148.04422 -26.072139,0 13.03607,-22.57914 z"
+ inkscape:randomized="0"
+ inkscape:rounded="0"
+ inkscape:flatsided="true"
+ sodipodi:arg2="1.5707963"
+ sodipodi:arg1="0.52359878"
+ sodipodi:r2="7.5263782"
+ sodipodi:r1="15.052756"
+ sodipodi:cy="140.51784"
+ sodipodi:cx="-24.946428"
+ sodipodi:sides="3"
+ id="path6390"
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.42985046;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ sodipodi:type="star" />
+ <circle
+ r="3.7797618"
+ cy="152.27068"
+ cx="-29.453543"
+ id="circle6392"
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1" />
+ <circle
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ id="circle6394"
+ cx="-19.683363"
+ cy="152.27068"
+ r="3.7797618" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6396"
+ d="m -34.423952,156.31501 h 19.710999"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6398"
+ d="m -16.68344,156.45039 -1.992943,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -19.734346,156.45039 -1.992942,2.82356"
+ id="path6400"
+ inkscape:connector-curvature="0" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6402"
+ d="m -22.785251,156.45039 -1.992942,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -25.836157,156.45039 -1.992942,2.82356"
+ id="path6404"
+ inkscape:connector-curvature="0" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6406"
+ d="m -28.887062,156.45039 -1.992942,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -31.937966,156.45039 -1.992942,2.82356"
+ id="path6408"
+ inkscape:connector-curvature="0" />
+ </g>
+ <g
+ style="stroke-width:4.21814156"
+ transform="matrix(0.24041841,0,0,-0.23377061,69.629943,223.8584)"
+ id="g6432">
+ <path
+ sodipodi:type="star"
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.42985046;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ id="path6412"
+ sodipodi:sides="3"
+ sodipodi:cx="-24.946428"
+ sodipodi:cy="140.51784"
+ sodipodi:r1="15.052756"
+ sodipodi:r2="7.5263782"
+ sodipodi:arg1="0.52359878"
+ sodipodi:arg2="1.5707963"
+ inkscape:flatsided="true"
+ inkscape:rounded="0"
+ inkscape:randomized="0"
+ d="m -11.910359,148.04422 -26.072139,0 13.03607,-22.57914 z"
+ inkscape:transform-center-y="-3.1008631"
+ transform="matrix(0.73936612,0,0,0.82399862,-6.1239083,26.230728)" />
+ <circle
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ id="circle6414"
+ cx="-29.453543"
+ cy="152.27068"
+ r="3.7797618" />
+ <circle
+ r="3.7797618"
+ cy="152.27068"
+ cx="-19.683363"
+ id="circle6416"
+ style="opacity:1;fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -34.423952,156.31501 h 19.710999"
+ id="path6418"
+ inkscape:connector-curvature="0" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -16.68344,156.45039 -1.992943,2.82356"
+ id="path6420"
+ inkscape:connector-curvature="0" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6422"
+ d="m -19.734346,156.45039 -1.992942,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -22.785251,156.45039 -1.992942,2.82356"
+ id="path6424"
+ inkscape:connector-curvature="0" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6426"
+ d="m -25.836157,156.45039 -1.992942,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ <path
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
+ d="m -28.887062,156.45039 -1.992942,2.82356"
+ id="path6428"
+ inkscape:connector-curvature="0" />
+ <path
+ inkscape:connector-curvature="0"
+ id="path6430"
+ d="m -31.937966,156.45039 -1.992942,2.82356"
+ style="fill:none;fill-rule:evenodd;stroke:#000000;stroke-width:1.11605px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1" />
+ </g>
+ <text
+ xml:space="preserve"
+ style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-size:2.82222223px;line-height:6.61458302px;font-family:'Latin Modern Sans';-inkscape-font-specification:'Latin Modern Sans';font-variant-ligatures:normal;font-variant-position:normal;font-variant-caps:normal;font-variant-numeric:normal;font-variant-alternates:normal;font-feature-settings:normal;text-indent:0;text-align:start;text-decoration:none;text-decoration-line:none;text-decoration-style:solid;text-decoration-color:#000000;letter-spacing:0px;word-spacing:0px;text-transform:none;writing-mode:lr-tb;direction:ltr;text-orientation:mixed;dominant-baseline:auto;baseline-shift:baseline;text-anchor:start;white-space:normal;shape-padding:0;opacity:1;vector-effect:none;fill:#000000;fill-opacity:1;stroke:none;stroke-width:0.26458332px;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ x="40.921589"
+ y="207.0658"
+ id="text6436"><tspan
+ sodipodi:role="line"
+ id="tspan6434"
+ x="40.921589"
+ y="207.0658"
+ style="font-size:3.52777791px;stroke-width:0.26458332px">simulation domain</tspan></text>
+ <text
+ xml:space="preserve"
+ style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-size:2.82222223px;line-height:4px;font-family:'Latin Modern Sans';-inkscape-font-specification:'Latin Modern Sans';font-variant-ligatures:normal;font-variant-position:normal;font-variant-caps:normal;font-variant-numeric:normal;font-variant-alternates:normal;font-feature-settings:normal;text-indent:0;text-align:start;text-decoration:none;text-decoration-line:none;text-decoration-style:solid;text-decoration-color:#000000;letter-spacing:0px;word-spacing:0px;text-transform:none;writing-mode:lr-tb;direction:ltr;text-orientation:mixed;dominant-baseline:auto;baseline-shift:baseline;text-anchor:start;white-space:normal;shape-padding:0;opacity:1;vector-effect:none;fill:#000000;fill-opacity:1;stroke:none;stroke-width:0.26458332px;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ x="93.009796"
+ y="204.63506"
+ id="text6440"><tspan
+ sodipodi:role="line"
+ id="tspan6438"
+ x="93.009796"
+ y="204.63506"
+ style="font-size:3.52777791px;line-height:4px;stroke-width:0.26458332px">elastic</tspan><tspan
+ sodipodi:role="line"
+ x="93.009796"
+ y="208.85095"
+ style="font-size:3.52777791px;line-height:4px;stroke-width:0.26458332px"
+ id="tspan6442">piston</tspan></text>
+ <text
+ xml:space="preserve"
+ style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-size:2.82222223px;line-height:6.61458302px;font-family:'Latin Modern Sans';-inkscape-font-specification:'Latin Modern Sans';font-variant-ligatures:normal;font-variant-position:normal;font-variant-caps:normal;font-variant-numeric:normal;font-variant-alternates:normal;font-feature-settings:normal;text-indent:0;text-align:start;text-decoration:none;text-decoration-line:none;text-decoration-style:solid;text-decoration-color:#000000;letter-spacing:0px;word-spacing:0px;text-transform:none;writing-mode:lr-tb;direction:ltr;text-orientation:mixed;dominant-baseline:auto;baseline-shift:baseline;text-anchor:start;white-space:normal;shape-padding:0;opacity:1;vector-effect:none;fill:#000000;fill-opacity:1;stroke:none;stroke-width:0.26458332px;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0;stroke-opacity:1"
+ x="54.790409"
+ y="249.7946"
+ id="text6472"><tspan
+ sodipodi:role="line"
+ id="tspan6470"
+ x="54.790409"
+ y="249.7946"
+ style="font-size:3.52777767px;stroke-width:0.26458332px">chamber filled with ideal gas</tspan></text>
+ </g>
+</svg>