diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformNeumann/c_NonuniformNeumann.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformNeumann/c_NonuniformNeumann.md
new file mode 100644
index 0000000000000000000000000000000000000000..cf46bf6546ca6b7816fb337f8920904f8173952c
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformNeumann/c_NonuniformNeumann.md
@@ -0,0 +1 @@
+Declares a Neumann boundary condition that is non-uniform in space.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformNeumann/t_field_name.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformNeumann/t_field_name.md
new file mode 100644
index 0000000000000000000000000000000000000000..7f36683addd9a72514cca3010606b6cd24bab004
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformNeumann/t_field_name.md
@@ -0,0 +1,5 @@
+This specifies a nodal field defined on the surface mesh (see the \c mesh
+parameter of the nonuniform Neumann BC).
+
+From that nodal field the flux values, which are applied via this Neumann BC,
+are read.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformNeumann/t_mesh.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformNeumann/t_mesh.md
new file mode 100644
index 0000000000000000000000000000000000000000..b448b3494b3b8e6e94b50a58ec0e10733171e06e
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformNeumann/t_mesh.md
@@ -0,0 +1,8 @@
+Path to the surface mesh vtu file.
+
+The surface mesh must contain a nodal integer-valued field (unsigned 64bit)
+named \c OriginalSubsurfaceNodeIDs. That field establishes the mapping between
+the nodes of the surface mesh to some notes in the bulk mesh.
+\warning It is not checked if the surface mesh and the bulk mesh correspond to each
+other; in particular it is not checked if surface and bulk nodes coincide and if
+surface elements coincide with the faces of bulk elements.
diff --git a/MeshLib/IO/VtkIO/VtuInterface.cpp b/MeshLib/IO/VtkIO/VtuInterface.cpp
index 7fcd2e2c487d0080f159c20fca7147dbe39a2916..75d2290de7fea6e83964349c346d25d6c6d9ecdc 100644
--- a/MeshLib/IO/VtkIO/VtuInterface.cpp
+++ b/MeshLib/IO/VtkIO/VtuInterface.cpp
@@ -61,7 +61,10 @@ MeshLib::Mesh* VtuInterface::readVTUFile(std::string const &file_name)
vtkUnstructuredGrid* vtkGrid = reader->GetOutput();
if (vtkGrid->GetNumberOfPoints() == 0)
+ {
+ ERR("Mesh \"%s\" contains zero points.", file_name.c_str());
return nullptr;
+ }
std::string const mesh_name (BaseLib::extractBaseNameWithoutExtension(file_name));
return MeshLib::VtkMeshConverter::convertUnstructuredGrid(vtkGrid, mesh_name);
diff --git a/MeshLib/Mesh.h b/MeshLib/Mesh.h
index bb307ed8acb8c66ffc2d07a1b7abaa73e411ccef..9937e9bb863ee8eaceb23a97ea36af7afafb1f27 100644
--- a/MeshLib/Mesh.h
+++ b/MeshLib/Mesh.h
@@ -133,11 +133,6 @@ public:
bool isAxiallySymmetric() const { return _is_axially_symmetric; }
void setAxiallySymmetric(bool is_axial_symmetric) {
_is_axially_symmetric = is_axial_symmetric;
- if (_is_axially_symmetric && getDimension() != 2) {
- OGS_FATAL(
- "Axial symmetry is implemented only for two-dimensional "
- "meshes.");
- }
}
protected:
diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
index 8df5b5e7bf7354726d678fb06b8969f5e76ebe87..4a164b8deca7e682cfa4f3c6de88cf0b00d96b4e 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
@@ -8,13 +8,14 @@
*/
#include "BoundaryCondition.h"
+#include "BoundaryConditionConfig.h"
+#include "DirichletBoundaryCondition.h"
#include "MeshGeoToolsLib/BoundaryElementsSearcher.h"
+#include "MeshGeoToolsLib/CreateSearchLength.h"
#include "MeshGeoToolsLib/MeshNodeSearcher.h"
#include "MeshGeoToolsLib/SearchLength.h"
-#include "MeshGeoToolsLib/CreateSearchLength.h"
-#include "BoundaryConditionConfig.h"
-#include "DirichletBoundaryCondition.h"
#include "NeumannBoundaryCondition.h"
+#include "NonuniformNeumannBoundaryCondition.h"
#include "RobinBoundaryCondition.h"
namespace ProcessLib
@@ -48,6 +49,12 @@ BoundaryConditionBuilder::createBoundaryCondition(
config, dof_table, mesh, variable_id,
integration_order, shapefunction_order, parameters);
}
+ if (type == "NonuniformNeumann")
+ {
+ return createNonuniformNeumannBoundaryCondition(
+ config, dof_table, mesh, variable_id, integration_order,
+ shapefunction_order);
+ }
OGS_FATAL("Unknown boundary condition type: `%s'.", type.c_str());
}
@@ -158,6 +165,18 @@ BoundaryConditionBuilder::createRobinBoundaryCondition(
parameters);
}
+std::unique_ptr<BoundaryCondition>
+BoundaryConditionBuilder::createNonuniformNeumannBoundaryCondition(
+ const BoundaryConditionConfig& config,
+ const NumLib::LocalToGlobalIndexMap& dof_table, const MeshLib::Mesh& mesh,
+ const int variable_id, const unsigned integration_order,
+ const unsigned shapefunction_order)
+{
+ return ProcessLib::createNonuniformNeumannBoundaryCondition(
+ config.config, dof_table, variable_id, config.component_id,
+ integration_order, shapefunction_order, mesh);
+}
+
std::vector<MeshLib::Element*> BoundaryConditionBuilder::getClonedElements(
MeshGeoToolsLib::BoundaryElementsSearcher& boundary_element_searcher,
GeoLib::GeoObject const& geometry)
diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.h b/ProcessLib/BoundaryCondition/BoundaryCondition.h
index e72923933438ca38cd5bdb42769f3e3051c260cc..c7dda2d155414cf94b9872092e3dd74f9a8a42c7 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.h
@@ -107,6 +107,13 @@ protected:
const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>&
parameters);
+ virtual std::unique_ptr<BoundaryCondition>
+ createNonuniformNeumannBoundaryCondition(
+ const BoundaryConditionConfig& config,
+ const NumLib::LocalToGlobalIndexMap& dof_table,
+ const MeshLib::Mesh& mesh, const int variable_id,
+ const unsigned integration_order, const unsigned shapefunction_order);
+
static std::vector<MeshLib::Element*> getClonedElements(
MeshGeoToolsLib::BoundaryElementsSearcher& boundary_element_searcher,
GeoLib::GeoObject const& geometry);
diff --git a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.h
index ce485b3b7f0676138d1d540d064eeae260951a3b..0654e5605ea6884a3639f369aacdd3c937a3b344 100644
--- a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.h
@@ -36,6 +36,17 @@ public:
_variable_id(variable_id),
_component_id(component_id)
{
+ 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));
+ }
}
void preTimestep(const double t) override;
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
index d2f0ea3ea4c5086546f44ce24c59bc1b5701324b..b15f8a80c590d0e7ea6b7741c1a2bc3ad1369758 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
@@ -35,8 +35,18 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
_elements(std::move(elements)),
_integration_order(integration_order)
{
- assert(component_id <
- static_cast<int>(dof_table_bulk.getNumberOfComponents()));
+ // check basic data consistency
+ if (variable_id >=
+ static_cast<int>(dof_table_bulk.getNumberOfVariables()) ||
+ component_id >=
+ dof_table_bulk.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_bulk.getNumberOfVariables(),
+ dof_table_bulk.getNumberOfVariableComponents(variable_id));
+ }
std::vector<MeshLib::Node*> nodes = MeshLib::getUniqueNodes(_elements);
DBUG("Found %d nodes for Natural BCs for the variable %d and component %d",
diff --git a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..6450b93aa5dcebefed8973310d33c78af099542b
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h
@@ -0,0 +1,101 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "GenericNonuniformNaturalBoundaryCondition.h"
+#include "GenericNonuniformNaturalBoundaryConditionLocalAssembler.h"
+#include "MeshLib/MeshSearch/NodeSearch.h"
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+
+namespace ProcessLib
+{
+template <typename BoundaryConditionData,
+ template <typename, typename, unsigned>
+ class LocalAssemblerImplementation>
+template <typename Data>
+GenericNonuniformNaturalBoundaryCondition<BoundaryConditionData,
+ LocalAssemblerImplementation>::
+ GenericNonuniformNaturalBoundaryCondition(
+ unsigned const integration_order, unsigned const shapefunction_order,
+ unsigned const global_dim,
+ std::unique_ptr<MeshLib::Mesh>&& boundary_mesh, Data&& data)
+ : _data(std::forward<Data>(data)), _boundary_mesh(std::move(boundary_mesh))
+{
+ static_assert(std::is_same<typename std::decay<BoundaryConditionData>::type,
+ typename std::decay<Data>::type>::value,
+ "Type mismatch between declared and passed BC data.");
+
+ // check basic data consistency
+ if (_data.variable_id_bulk >=
+ static_cast<int>(_data.dof_table_bulk.getNumberOfVariables()) ||
+ _data.component_id_bulk >=
+ _data.dof_table_bulk.getNumberOfVariableComponents(
+ _data.variable_id_bulk))
+ {
+ OGS_FATAL(
+ "Variable id or component id too high. Actual values: (%d, %d), "
+ "maximum values: (%d, %d).",
+ _data.variable_id_bulk, _data.component_id_bulk,
+ _data.dof_table_bulk.getNumberOfVariables(),
+ _data.dof_table_bulk.getNumberOfVariableComponents(
+ _data.variable_id_bulk));
+ }
+
+ if (_boundary_mesh->getDimension() + 1 != global_dim)
+ {
+ OGS_FATAL(
+ "The dimension of the given boundary mesh (%d) is not by one lower "
+ "than the bulk dimension (%d).",
+ _boundary_mesh->getDimension(), global_dim);
+ }
+
+ constructDofTable();
+
+ createLocalAssemblers<LocalAssemblerImplementation>(
+ global_dim, _boundary_mesh->getElements(), *_dof_table_boundary,
+ shapefunction_order, _local_assemblers,
+ _boundary_mesh->isAxiallySymmetric(), integration_order, _data);
+}
+
+template <typename BoundaryConditionData,
+ template <typename, typename, unsigned>
+ class LocalAssemblerImplementation>
+void GenericNonuniformNaturalBoundaryCondition<
+ BoundaryConditionData, LocalAssemblerImplementation>::constructDofTable()
+{
+ // construct one-component DOF-table for the surface mesh
+ _mesh_subset_all_nodes.reset(
+ new MeshLib::MeshSubset(*_boundary_mesh, &_boundary_mesh->getNodes()));
+
+ std::vector<MeshLib::MeshSubsets> all_mesh_subsets{
+ _mesh_subset_all_nodes.get()};
+
+ std::vector<unsigned> vec_var_n_components{1};
+
+ _dof_table_boundary = std::make_unique<NumLib::LocalToGlobalIndexMap>(
+ std::move(all_mesh_subsets), vec_var_n_components,
+ NumLib::ComponentOrder::BY_LOCATION);
+}
+
+template <typename BoundaryConditionData,
+ template <typename, typename, unsigned>
+ class LocalAssemblerImplementation>
+void GenericNonuniformNaturalBoundaryCondition<
+ BoundaryConditionData,
+ LocalAssemblerImplementation>::applyNaturalBC(const double t,
+ const GlobalVector& x,
+ GlobalMatrix& K,
+ GlobalVector& b)
+{
+ GlobalExecutor::executeMemberOnDereferenced(
+ &GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface::
+ assemble,
+ _local_assemblers, *_dof_table_boundary, t, x, K, b);
+}
+
+} // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h
new file mode 100644
index 0000000000000000000000000000000000000000..4a4f75ce52a94936ca4701a4e91bf573d230f0b9
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h
@@ -0,0 +1,62 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "BoundaryCondition.h"
+#include "MeshLib/MeshSubset.h"
+
+namespace ProcessLib
+{
+class GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface;
+
+template <typename BoundaryConditionData,
+ template <typename, typename, unsigned>
+ class LocalAssemblerImplementation>
+class GenericNonuniformNaturalBoundaryCondition final : public BoundaryCondition
+{
+public:
+ /// Create a boundary condition process from given config,
+ /// DOF-table, and a mesh subset for a given variable and its component.
+ /// A local DOF-table, a subset of the given one, is constructed.
+ template <typename Data>
+ GenericNonuniformNaturalBoundaryCondition(
+ unsigned const integration_order, unsigned const shapefunction_order,
+ unsigned const global_dim,
+ std::unique_ptr<MeshLib::Mesh>&& boundary_mesh, Data&& data);
+
+ /// 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;
+
+private:
+ void constructDofTable();
+
+ /// Data used in the assembly of the specific boundary condition.
+ BoundaryConditionData _data;
+
+ std::unique_ptr<MeshLib::Mesh> _boundary_mesh;
+
+ std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_all_nodes;
+
+ /// DOF-table (one single-component variable) of the boundary mesh.
+ std::unique_ptr<NumLib::LocalToGlobalIndexMap> _dof_table_boundary;
+
+ /// Local assemblers for each element of the boundary mesh.
+ std::vector<std::unique_ptr<
+ GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface>>
+ _local_assemblers;
+};
+
+} // ProcessLib
+
+#include "GenericNonuniformNaturalBoundaryCondition-impl.h"
diff --git a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryConditionLocalAssembler.h
new file mode 100644
index 0000000000000000000000000000000000000000..873d05f838c11650e2913182977aed254099cac8
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryConditionLocalAssembler.h
@@ -0,0 +1,90 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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/Fem/ShapeMatrixPolicy.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+namespace ProcessLib
+{
+class GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface
+{
+public:
+ virtual ~GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface() =
+ default;
+
+ 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;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+class GenericNonuniformNaturalBoundaryConditionLocalAssembler
+ : public GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface
+{
+protected:
+ using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+ using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
+ using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
+
+ struct NAndWeight
+ {
+ NAndWeight(typename ShapeMatricesType::ShapeMatrices::ShapeType&& N_,
+ double const weight_)
+ : N(std::move(N_)), weight(weight_)
+ {
+ }
+ typename ShapeMatricesType::ShapeMatrices::ShapeType const N;
+ double const weight;
+ };
+
+private:
+ static std::vector<NAndWeight, Eigen::aligned_allocator<NAndWeight>>
+ initNsAndWeights(MeshLib::Element const& e, bool is_axially_symmetric,
+ unsigned const integration_order)
+ {
+ IntegrationMethod const integration_method(integration_order);
+ std::vector<NAndWeight, Eigen::aligned_allocator<NAndWeight>>
+ ns_and_weights;
+ ns_and_weights.reserve(integration_method.getNumberOfPoints());
+
+ auto sms = initShapeMatrices<ShapeFunction, ShapeMatricesType,
+ IntegrationMethod, GlobalDim>(
+ e, is_axially_symmetric, integration_method);
+ for (unsigned ip = 0; ip < sms.size(); ++ip)
+ {
+ auto& sm = sms[ip];
+ double const w =
+ sm.detJ * sm.integralMeasure *
+ integration_method.getWeightedPoint(ip).getWeight();
+
+ ns_and_weights.emplace_back(std::move(sm.N), w);
+ }
+
+ return ns_and_weights;
+ }
+
+public:
+ GenericNonuniformNaturalBoundaryConditionLocalAssembler(
+ MeshLib::Element const& e, bool is_axially_symmetric,
+ unsigned const integration_order)
+ : _ns_and_weights(
+ initNsAndWeights(e, is_axially_symmetric, integration_order))
+ {
+ }
+
+protected:
+ std::vector<NAndWeight, Eigen::aligned_allocator<NAndWeight>> const
+ _ns_and_weights;
+};
+
+} // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dc8c456eed4ad4ce0ef356722eec77bbd859e03e
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp
@@ -0,0 +1,92 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "NonuniformNeumannBoundaryCondition.h"
+
+#include "MeshLib/IO/readMeshFromFile.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+namespace ProcessLib
+{
+std::unique_ptr<NonuniformNeumannBoundaryCondition>
+createNonuniformNeumannBoundaryCondition(
+ BaseLib::ConfigTree const& config,
+ NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
+ int const component_id, unsigned const integration_order,
+ unsigned const shapefunction_order, MeshLib::Mesh const& bulk_mesh)
+{
+ DBUG("Constructing NonuniformNeumann BC from config.");
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
+ config.checkConfigParameter("type", "NonuniformNeumann");
+
+ // TODO handle paths correctly
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__NonuniformNeumann__mesh}
+ auto const mesh_file = config.getConfigParameter<std::string>("mesh");
+
+ std::unique_ptr<MeshLib::Mesh> boundary_mesh(
+ MeshLib::IO::readMeshFromFile(mesh_file));
+
+ if (!boundary_mesh)
+ {
+ OGS_FATAL("Error reading mesh `%s'", mesh_file.c_str());
+ }
+
+ // Surface mesh and bulk mesh must have equal axial symmetry flags!
+ boundary_mesh->setAxiallySymmetric(bulk_mesh.isAxiallySymmetric());
+
+ // TODO finally use ProcessLib::Parameter here
+ auto const field_name =
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__NonuniformNeumann__field_name}
+ config.getConfigParameter<std::string>("field_name");
+
+ auto const* const property =
+ boundary_mesh->getProperties().getPropertyVector<double>(field_name);
+
+ if (!property)
+ {
+ OGS_FATAL("A property with name `%s' does not exist in `%s'.",
+ field_name.c_str(), mesh_file.c_str());
+ }
+
+ if (property->getMeshItemType() != MeshLib::MeshItemType::Node)
+ {
+ OGS_FATAL(
+ "Only nodal fields are supported for non-uniform fields. Field "
+ "`%s' is not nodal.",
+ field_name.c_str());
+ }
+
+ if (property->getNumberOfComponents() != 1)
+ {
+ OGS_FATAL("`%s' is not a one-component field.", field_name.c_str());
+ }
+
+ std::string mapping_to_bulk_nodes_property = "OriginalSubsurfaceNodeIDs";
+ auto const* const mapping_to_bulk_nodes =
+ boundary_mesh->getProperties().getPropertyVector<std::size_t>(
+ mapping_to_bulk_nodes_property);
+
+ if (!(mapping_to_bulk_nodes &&
+ mapping_to_bulk_nodes->getMeshItemType() ==
+ MeshLib::MeshItemType::Node) &&
+ mapping_to_bulk_nodes->getNumberOfComponents() == 1)
+ {
+ OGS_FATAL("Field `%s' is not set up properly.",
+ mapping_to_bulk_nodes_property.c_str());
+ }
+
+ return std::make_unique<NonuniformNeumannBoundaryCondition>(
+ integration_order, shapefunction_order, bulk_mesh.getDimension(),
+ std::move(boundary_mesh),
+ NonuniformNeumannBoundaryConditionData{
+ *property, bulk_mesh.getID(), *mapping_to_bulk_nodes, dof_table,
+ variable_id, component_id});
+}
+
+} // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h
new file mode 100644
index 0000000000000000000000000000000000000000..2073e7eb9f305c46c6712b246eb76972ea2a0f29
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h
@@ -0,0 +1,30 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "GenericNonuniformNaturalBoundaryCondition.h"
+#include "MeshLib/PropertyVector.h"
+#include "NonuniformNeumannBoundaryConditionLocalAssembler.h"
+
+namespace ProcessLib
+{
+using NonuniformNeumannBoundaryCondition =
+ GenericNonuniformNaturalBoundaryCondition<
+ NonuniformNeumannBoundaryConditionData,
+ NonuniformNeumannBoundaryConditionLocalAssembler>;
+
+std::unique_ptr<NonuniformNeumannBoundaryCondition>
+createNonuniformNeumannBoundaryCondition(
+ BaseLib::ConfigTree const& config,
+ NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
+ int const component_id, unsigned const integration_order,
+ unsigned const shapefunction_order, const MeshLib::Mesh& bulk_mesh);
+
+} // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h
new file mode 100644
index 0000000000000000000000000000000000000000..12a2b88a89dd425823bec9ec3f3e36755757826a
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h
@@ -0,0 +1,110 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "MeshLib/PropertyVector.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/Function/Interpolation.h"
+
+#include "GenericNonuniformNaturalBoundaryConditionLocalAssembler.h"
+
+namespace ProcessLib
+{
+struct NonuniformNeumannBoundaryConditionData
+{
+ /* TODO use Parameter */
+ MeshLib::PropertyVector<double> const& values;
+
+ // Used for mapping boundary nodes to bulk nodes.
+ std::size_t bulk_mesh_id;
+ MeshLib::PropertyVector<std::size_t> const& mapping_to_bulk_nodes;
+ NumLib::LocalToGlobalIndexMap const& dof_table_bulk;
+ int const variable_id_bulk;
+ int const component_id_bulk;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+class NonuniformNeumannBoundaryConditionLocalAssembler final
+ : public GenericNonuniformNaturalBoundaryConditionLocalAssembler<
+ ShapeFunction, IntegrationMethod, GlobalDim>
+{
+ using Base = GenericNonuniformNaturalBoundaryConditionLocalAssembler<
+ ShapeFunction, IntegrationMethod, GlobalDim>;
+
+public:
+ /// The neumann_bc_value factor is directly integrated into the local
+ /// element matrix.
+ NonuniformNeumannBoundaryConditionLocalAssembler(
+ MeshLib::Element const& e,
+ std::size_t const local_matrix_size,
+ bool is_axially_symmetric,
+ unsigned const integration_order,
+ NonuniformNeumannBoundaryConditionData const& data)
+ : Base(e, is_axially_symmetric, integration_order),
+ _data(data),
+ _local_rhs(local_matrix_size)
+ {
+ }
+
+ void assemble(std::size_t const id,
+ NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
+ double const /*t*/, const GlobalVector& /*x*/,
+ GlobalMatrix& /*K*/, GlobalVector& b) override
+ {
+ _local_rhs.setZero();
+
+ auto indices = NumLib::getIndices(id, dof_table_boundary);
+
+ // TODO lots of heap allocations.
+ std::vector<double> neumann_param_nodal_values_local;
+ neumann_param_nodal_values_local.reserve(indices.size());
+ for (auto i : indices)
+ {
+ neumann_param_nodal_values_local.push_back(
+ _data.values.getComponent(i, 0));
+ }
+
+ auto const n_integration_points = Base::_ns_and_weights.size();
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto const& n_and_weight = Base::_ns_and_weights[ip];
+ double flux;
+ NumLib::shapeFunctionInterpolate(neumann_param_nodal_values_local,
+ n_and_weight.N, flux);
+ _local_rhs.noalias() +=
+ n_and_weight.N * (n_and_weight.weight * flux);
+ }
+
+ // map boundary dof indices to bulk dof indices
+ for (auto& i : indices)
+ {
+ auto const bulk_node_id =
+ _data.mapping_to_bulk_nodes.getComponent(i, 0);
+
+ MeshLib::Location const l{
+ _data.bulk_mesh_id, MeshLib::MeshItemType::Node, bulk_node_id};
+
+ i = _data.dof_table_bulk.getGlobalIndex(l, _data.variable_id_bulk,
+ _data.component_id_bulk);
+ assert(i != NumLib::MeshComponentMap::nop);
+ }
+ b.add(indices, _local_rhs);
+ }
+
+private:
+ NonuniformNeumannBoundaryConditionData const& _data;
+ typename Base::NodalVectorType _local_rhs;
+
+public:
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+} // namespace ProcessLib
diff --git a/ProcessLib/GroundwaterFlow/Tests.cmake b/ProcessLib/GroundwaterFlow/Tests.cmake
index 76d21f130281533aa25c3b71150d8ed671d4858f..a155dc97a64482a861f099572a826b859163e230 100644
--- a/ProcessLib/GroundwaterFlow/Tests.cmake
+++ b/ProcessLib/GroundwaterFlow/Tests.cmake
@@ -573,3 +573,15 @@ foreach(mesh_size 1e1)
line_1_line_${mesh_size}.vtu line_${mesh_size}_neumann_pcs_0_ts_1_t_1_000000_0.vtu D1_left_N1_right pressure
)
endforeach()
+
+AddTest(
+ NAME GroundWaterFlowProcess_Neumann_nonuniform
+ PATH Elliptic/nonuniform_bc_Groundwaterflow
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS neumann_nonuniform.prj
+ TESTER vtkdiff
+ REQUIREMENTS NOT OGS_USE_MPI
+ ABSTOL 4e-2 RELTOL 1e-16
+ DIFF_DATA
+ inhomogeneous_permeability.vtu neumann_nonuniform_pcs_0_ts_1_t_1.000000.vtu mass_flux_ref mass_flux
+)
diff --git a/Tests/Data b/Tests/Data
index 656a73bb38263daf2ece092092906c6ffcab5d85..54154c90a718bc33443156fcf08b3e00237d2316 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit 656a73bb38263daf2ece092092906c6ffcab5d85
+Subproject commit 54154c90a718bc33443156fcf08b3e00237d2316