diff --git a/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp
index e3518dd4c0b5b8edc069192d506f336c9f93b088..c4f2b4cb3dfbe01b9e6b065e96528fea49090e4e 100644
--- a/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp
@@ -15,8 +15,6 @@
 #include "DirichletBoundaryCondition.h"
 #include "DirichletBoundaryConditionWithinTimeInterval.h"
 #include "NeumannBoundaryCondition.h"
-#include "NonuniformDirichletBoundaryCondition.h"
-#include "NonuniformNeumannBoundaryCondition.h"
 #include "NonuniformVariableDependentNeumannBoundaryCondition.h"
 #include "NormalTractionBoundaryCondition.h"
 #include "PhaseFieldIrreversibleDamageOracleBoundaryCondition.h"
@@ -78,20 +76,6 @@ std::unique_ptr<BoundaryCondition> createBoundaryCondition(
             *config.component_id, integration_order, shapefunction_order,
             bulk_mesh.getDimension(), parameters);
     }
-    if (type == "NonuniformDirichlet")
-    {
-        return ProcessLib::createNonuniformDirichletBoundaryCondition(
-            config.config, config.boundary_mesh, dof_table, variable_id,
-            *config.component_id);
-    }
-    if (type == "NonuniformNeumann")
-    {
-        return ProcessLib::createNonuniformNeumannBoundaryCondition(
-            config.config, config.boundary_mesh, dof_table, variable_id,
-            *config.component_id, integration_order, shapefunction_order,
-            bulk_mesh);
-    }
-
     if (type == "NonuniformVariableDependentNeumann")
     {
         return ProcessLib::
diff --git a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp
index fc757a7b5c0db70d132d9622b4dca63c233d8539..c714a531dec83f6ec4c75e076bfd18b546aa4d0d 100644
--- a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp
@@ -67,7 +67,16 @@ std::unique_ptr<DirichletBoundaryCondition> createDirichletBoundaryCondition(
     auto const param_name = config.getConfigParameter<std::string>("parameter");
     DBUG("Using parameter %s", param_name.c_str());
 
-    auto& param = findParameter<double>(param_name, parameters, 1);
+    auto& parameter = findParameter<double>(param_name, parameters, 1);
+
+    if (parameter.mesh() && *parameter.mesh() != bc_mesh)
+    {
+        OGS_FATAL(
+            "The boundary condition is defined on a different domain than the "
+            "parameter: boundary condition is defined on mesh '%s', parameter "
+            "is defined on mesh '%s'.",
+            bc_mesh.getName().c_str(), parameter.mesh()->getName().c_str());
+    }
 
 // In case of partitioned mesh the boundary could be empty, i.e. there is no
 // boundary condition.
@@ -84,7 +93,7 @@ std::unique_ptr<DirichletBoundaryCondition> createDirichletBoundaryCondition(
 #endif  // USE_PETSC
 
     return std::make_unique<DirichletBoundaryCondition>(
-        param, bc_mesh, dof_table_bulk, variable_id, component_id);
+        parameter, bc_mesh, dof_table_bulk, variable_id, component_id);
 }
 
 }  // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
index 529313734041734872b3440f6a7e5bf55632ec94..20907cf23c499b3d303cf3729ec74d9f604f6db2 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
@@ -26,6 +26,7 @@ class NeumannBoundaryConditionLocalAssembler final
     using Base = GenericNaturalBoundaryConditionLocalAssembler<
         ShapeFunction, IntegrationMethod, GlobalDim>;
     using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+    using NodalVectorType = typename Base::NodalVectorType;
 
 public:
     /// The neumann_bc_value factor is directly integrated into the local
@@ -61,6 +62,11 @@ public:
         FemType fe(
              static_cast<const typename ShapeFunction::MeshElement&>(Base::_element));
 
+        // Get element nodes for the interpolation from nodes to
+        // integration point.
+        NodalVectorType parameter_node_values =
+            _neumann_bc_parameter.getNodalValuesOnElement(Base::_element, t);
+
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
             pos.setIntegrationPoint(ip);
@@ -70,7 +76,7 @@ public:
             MathLib::TemplatePoint<double, 3> coords_ip(fe.interpolateCoordinates(sm.N));
             pos.setCoordinates(coords_ip);
 
-            _local_rhs.noalias() += sm.N * _neumann_bc_parameter(t, pos)[0] *
+            _local_rhs.noalias() += sm.N * parameter_node_values.dot(sm.N) *
                                     sm.detJ * wp.getWeight() *
                                     sm.integralMeasure;
         }
@@ -81,7 +87,7 @@ public:
 
 private:
     Parameter<double> const& _neumann_bc_parameter;
-    typename Base::NodalVectorType _local_rhs;
+    NodalVectorType _local_rhs;
 
 public:
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
diff --git a/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.cpp
deleted file mode 100644
index 5786e0947c6f2c3840e12746fe58fdfcdb1ec01b..0000000000000000000000000000000000000000
--- a/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.cpp
+++ /dev/null
@@ -1,63 +0,0 @@
-/**
- * \file
- *
- * \copyright
- * Copyright (c) 2012-2019, 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 "NonuniformDirichletBoundaryCondition.h"
-
-#include "MeshLib/IO/readMeshFromFile.h"
-#include "ProcessLib/Utils/ProcessUtils.h"
-
-namespace ProcessLib
-{
-std::unique_ptr<NonuniformDirichletBoundaryCondition>
-createNonuniformDirichletBoundaryCondition(
-    BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh,
-    NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
-    int const component_id)
-{
-    DBUG("Constructing NonuniformDirichlet BC from config.");
-    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
-    config.checkConfigParameter("type", "NonuniformDirichlet");
-
-    // TODO finally use ProcessLib::Parameter here
-    auto const field_name =
-        //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__NonuniformDirichlet__field_name}
-        config.getConfigParameter<std::string>("field_name");
-
-    auto const* const property =
-        boundary_mesh.getProperties().getPropertyVector<double>(
-            field_name, MeshLib::MeshItemType::Node, 1);
-
-    if (!property)
-    {
-        OGS_FATAL("A property with name `%s' does not exist in `%s'.",
-                  field_name.c_str(), boundary_mesh.getName().c_str());
-    }
-
-    // In case of partitioned mesh the boundary could be empty, i.e. there is no
-    // boundary condition.
-#ifdef USE_PETSC
-    // This can be extracted to createBoundaryCondition() but then the config
-    // parameters are not read and will cause an error.
-    // TODO (naumov): Add a function to ConfigTree for skipping the tags of the
-    // subtree and move the code up in createBoundaryCondition().
-    if (boundary_mesh.getDimension() == 0 &&
-        boundary_mesh.getNumberOfNodes() == 0 &&
-        boundary_mesh.getNumberOfElements() == 0)
-    {
-        return nullptr;
-    }
-#endif  // USE_PETSC
-
-    return std::make_unique<NonuniformDirichletBoundaryCondition>(
-        boundary_mesh, *property, dof_table, variable_id, component_id);
-}
-
-}  // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h
deleted file mode 100644
index ba67cb886f3236a3089ce0f421a5140486f857f8..0000000000000000000000000000000000000000
--- a/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h
+++ /dev/null
@@ -1,129 +0,0 @@
-/**
- * \file
- *
- * \copyright
- * Copyright (c) 2012-2019, 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/PropertyVector.h"
-#include "NumLib/DOF/DOFTableUtil.h"
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "NumLib/IndexValueVector.h"
-#include "ProcessLib/Parameter/Parameter.h"
-#include "ProcessLib/Utils/ProcessUtils.h"
-
-namespace ProcessLib
-{
-class NonuniformDirichletBoundaryCondition final : public BoundaryCondition
-{
-public:
-    NonuniformDirichletBoundaryCondition(
-        // int const bulk_mesh_dimension,
-        MeshLib::Mesh const& boundary_mesh,
-        MeshLib::PropertyVector<double> const& values,
-        NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
-        int const variable_id,
-        int const component_id)
-        : _values(values),
-          _boundary_mesh(boundary_mesh),
-          _variable_id(variable_id),
-          _component_id(component_id)
-    {
-        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));
-        }
-
-        if (!_boundary_mesh.getProperties().existsPropertyVector<std::size_t>(
-                "bulk_node_ids"))
-        {
-            OGS_FATAL(
-                "The required bulk node ids map does not exist in the boundary "
-                "mesh '%s'.",
-                _boundary_mesh.getName().c_str());
-        }
-
-        std::vector<MeshLib::Node*> const& bc_nodes = _boundary_mesh.getNodes();
-        DBUG(
-            "Found %d nodes for Natural BCs for the variable %d and component "
-            "%d",
-            bc_nodes.size(), variable_id, component_id);
-
-        MeshLib::MeshSubset boundary_mesh_subset(_boundary_mesh, bc_nodes);
-
-        // Create local DOF table from the BC mesh subset for the given variable
-        // and component id.
-        _dof_table_boundary.reset(dof_table_bulk.deriveBoundaryConstrainedMap(
-            variable_id, {component_id}, std::move(boundary_mesh_subset)));
-    }
-
-    void getEssentialBCValues(
-        const double /*t*/, GlobalVector const& /*x*/,
-        NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override
-    {
-        bc_values.ids.clear();
-        bc_values.values.clear();
-
-        // Convert mesh node ids to global index for the given component.
-        bc_values.ids.reserve(_values.size());
-        bc_values.values.reserve(_values.size());
-
-        // Map boundary dof indices to bulk dof indices and the corresponding
-        // values.
-        for (auto const* const node : _boundary_mesh.getNodes())
-        {
-            auto const node_id = node->getID();
-            auto const global_index = _dof_table_boundary->getGlobalIndex(
-                {_boundary_mesh.getID(), MeshLib::MeshItemType::Node, node_id},
-                _variable_id, _component_id);
-            if (global_index == NumLib::MeshComponentMap::nop)
-            {
-                continue;
-            }
-            // For the DDC approach (e.g. with PETSc option), the negative index
-            // of global_index 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 (global_index >= 0)
-            {
-                bc_values.ids.emplace_back(global_index);
-                bc_values.values.push_back(_values[node_id]);
-            }
-        }
-    }
-
-private:
-    MeshLib::PropertyVector<double> const& _values;
-    MeshLib::Mesh const& _boundary_mesh;
-    std::unique_ptr<NumLib::LocalToGlobalIndexMap const> _dof_table_boundary;
-    int const _variable_id;
-    int const _component_id;
-};
-
-std::unique_ptr<NonuniformDirichletBoundaryCondition>
-createNonuniformDirichletBoundaryCondition(
-    BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh,
-    NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
-    int const component_id);
-
-}  // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp
deleted file mode 100644
index bbe963d3cae71e39349c59936f8ab09d2a3d15c0..0000000000000000000000000000000000000000
--- a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp
+++ /dev/null
@@ -1,92 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2019, 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, MeshLib::Mesh const& boundary_mesh,
-    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 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(), boundary_mesh.getName().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 const mapping_to_bulk_nodes_property = "bulk_node_ids";
-    auto const* const mapping_to_bulk_nodes =
-        boundary_mesh.getProperties().getPropertyVector<std::size_t>(
-            mapping_to_bulk_nodes_property);
-
-    if (mapping_to_bulk_nodes == nullptr ||
-        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());
-    }
-
-    // In case of partitioned mesh the boundary could be empty, i.e. there is no
-    // boundary condition.
-#ifdef USE_PETSC
-    // This can be extracted to createBoundaryCondition() but then the config
-    // parameters are not read and will cause an error.
-    // TODO (naumov): Add a function to ConfigTree for skipping the tags of the
-    // subtree and move the code up in createBoundaryCondition().
-    if (boundary_mesh.getDimension() == 0 &&
-        boundary_mesh.getNumberOfNodes() == 0 &&
-        boundary_mesh.getNumberOfElements() == 0)
-    {
-        return nullptr;
-    }
-#endif  // USE_PETSC
-
-    return std::make_unique<NonuniformNeumannBoundaryCondition>(
-        integration_order, shapefunction_order, dof_table, variable_id,
-        component_id, bulk_mesh.getDimension(), boundary_mesh,
-        NonuniformNeumannBoundaryConditionData{
-            *property, bulk_mesh.getID(), *mapping_to_bulk_nodes, dof_table,
-            variable_id, component_id});
-}
-
-}  // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h
deleted file mode 100644
index 50cb24c44bbb2fecef3e63215f83fc149e864a40..0000000000000000000000000000000000000000
--- a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h
+++ /dev/null
@@ -1,30 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2019, 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, MeshLib::Mesh const& boundary_mesh,
-    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);
-
-}  // namespace ProcessLib
diff --git a/ProcessLib/Parameter/MeshNodeParameter.cpp b/ProcessLib/Parameter/MeshNodeParameter.cpp
index 48fc8f5a281de41edd95e91a6145056f32da9d3e..11b9ecbabe2186d6c82605503bf7ad8560c262ea 100644
--- a/ProcessLib/Parameter/MeshNodeParameter.cpp
+++ b/ProcessLib/Parameter/MeshNodeParameter.cpp
@@ -32,7 +32,7 @@ std::unique_ptr<ParameterBase> createMeshNodeParameter(
                   field_name.c_str());
     }
 
-    return std::make_unique<MeshNodeParameter<double>>(name, *property);
+    return std::make_unique<MeshNodeParameter<double>>(name, mesh, *property);
 }
 
 }  // namespace ProcessLib
diff --git a/ProcessLib/Parameter/MeshNodeParameter.h b/ProcessLib/Parameter/MeshNodeParameter.h
index 5c013c22dcf87e60395b2489022091d78c3d2aea..5ea4da35a38b86bfbb3924d83d1c3459d2521a0f 100644
--- a/ProcessLib/Parameter/MeshNodeParameter.h
+++ b/ProcessLib/Parameter/MeshNodeParameter.h
@@ -27,9 +27,9 @@ namespace ProcessLib
 template <typename T>
 struct MeshNodeParameter final : public Parameter<T> {
     MeshNodeParameter(std::string const& name_,
+                      MeshLib::Mesh const& mesh,
                       MeshLib::PropertyVector<T> const& property)
-        : Parameter<T>(name_),
-          _property(property)
+        : Parameter<T>(name_, &mesh), _property(property)
     {
     }