diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
index 55edb6fda627ceb6b2ce6162d397d9c3d9c30840..5ad6afe6310d0516f691b6796f32cf6ba9cb3c2b 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
@@ -15,6 +15,7 @@
 #include "MeshGeoToolsLib/MeshNodeSearcher.h"
 #include "MeshGeoToolsLib/SearchLength.h"
 #include "NeumannBoundaryCondition.h"
+#include "NonuniformDirichletBoundaryCondition.h"
 #include "NonuniformNeumannBoundaryCondition.h"
 #include "PressureBoundaryCondition.h"
 #include "RobinBoundaryCondition.h"
@@ -50,6 +51,11 @@ BoundaryConditionBuilder::createBoundaryCondition(
                     config, dof_table, mesh, variable_id,
                     integration_order, shapefunction_order, parameters);
     }
+    if (type == "NonuniformDirichlet")
+    {
+        return createNonuniformDirichletBoundaryCondition(config, dof_table,
+                                                          mesh, variable_id);
+    }
     if (type == "NonuniformNeumann")
     {
         return createNonuniformNeumannBoundaryCondition(
@@ -175,6 +181,16 @@ BoundaryConditionBuilder::createRobinBoundaryCondition(
         parameters);
 }
 
+std::unique_ptr<BoundaryCondition>
+BoundaryConditionBuilder::createNonuniformDirichletBoundaryCondition(
+    const BoundaryConditionConfig& config,
+    const NumLib::LocalToGlobalIndexMap& dof_table, const MeshLib::Mesh& mesh,
+    const int variable_id)
+{
+    return ProcessLib::createNonuniformDirichletBoundaryCondition(
+        config.config, dof_table, variable_id, *config.component_id, mesh);
+}
+
 std::unique_ptr<BoundaryCondition>
 BoundaryConditionBuilder::createNonuniformNeumannBoundaryCondition(
     const BoundaryConditionConfig& config,
diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.h b/ProcessLib/BoundaryCondition/BoundaryCondition.h
index 8c6a28a1b0b183185cf3abddf4938c200c557360..dc9088f9cf21a64cd5aba5bf051df80b1c88bf10 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.h
@@ -107,6 +107,12 @@ protected:
         const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>&
             parameters);
 
+    virtual std::unique_ptr<BoundaryCondition>
+    createNonuniformDirichletBoundaryCondition(
+        const BoundaryConditionConfig& config,
+        const NumLib::LocalToGlobalIndexMap& dof_table,
+        const MeshLib::Mesh& mesh, const int variable_id);
+
     virtual std::unique_ptr<BoundaryCondition>
     createNonuniformNeumannBoundaryCondition(
         const BoundaryConditionConfig& config,
diff --git a/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c7802df5f86d8c150f515b9ee6c831bc5a390d88
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.cpp
@@ -0,0 +1,93 @@
+/**
+ * \file
+ *
+ * \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 "NonuniformDirichletBoundaryCondition.h"
+
+#include "MeshLib/IO/readMeshFromFile.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+namespace ProcessLib
+{
+std::unique_ptr<NonuniformDirichletBoundaryCondition>
+createNonuniformDirichletBoundaryCondition(
+    BaseLib::ConfigTree const& config,
+    NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
+    int const component_id, MeshLib::Mesh const& bulk_mesh)
+{
+    DBUG("Constructing NonuniformDirichlet BC from config.");
+    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
+    config.checkConfigParameter("type", "NonuniformDirichlet");
+
+    // TODO handle paths correctly
+    //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__NonuniformDirichlet__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());
+    }
+
+    // The axial symmetry is not used in the Dirichlet BC but kept here for
+    // consistency.
+    //
+    // 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__NonuniformDirichlet__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<NonuniformDirichletBoundaryCondition>(
+        // bulk_mesh.getDimension(),
+        std::move(boundary_mesh), *property, bulk_mesh.getID(),
+        *mapping_to_bulk_nodes, dof_table, variable_id, component_id);
+}
+
+}  // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h
new file mode 100644
index 0000000000000000000000000000000000000000..6d23ea8abb0399c58ac4b8336da662ca1bb25fa6
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h
@@ -0,0 +1,119 @@
+/**
+ * \file
+ *
+ * \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/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,
+        std::unique_ptr<MeshLib::Mesh>
+            boundary_mesh,
+        MeshLib::PropertyVector<double> const& values,
+        std::size_t const 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)
+        : _bulk_mesh_id(bulk_mesh_id),
+          _values(values),
+          _boundary_mesh(std::move(boundary_mesh)),
+          _mapping_to_bulk_nodes(mapping_to_bulk_nodes),
+          _dof_table_bulk(dof_table_bulk),
+          _variable_id_bulk(variable_id_bulk),
+          _component_id_bulk(component_id_bulk)
+    {
+        if (_variable_id_bulk >=
+                static_cast<int>(_dof_table_bulk.getNumberOfVariables()) ||
+            _component_id_bulk >= _dof_table_bulk.getNumberOfVariableComponents(
+                                      _variable_id_bulk))
+        {
+            OGS_FATAL(
+                "Variable id or component id too high. Actual values: (%d, "
+                "%d), "
+                "maximum values: (%d, %d).",
+                _variable_id_bulk, _component_id_bulk,
+                _dof_table_bulk.getNumberOfVariables(),
+                _dof_table_bulk.getNumberOfVariableComponents(
+                    _variable_id_bulk));
+        }
+    }
+
+    void preTimestep(const double /*t*/) override {}
+
+    void getEssentialBCValues(
+        const double /*t*/,
+        NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override
+    {
+        // TODO: Reenable when fixed ;)
+        // if (_already_computed)
+        // return;
+
+        _already_computed = true;
+
+        SpatialPosition pos;
+
+        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 (std::size_t i = 0; i < _boundary_mesh->getNumberOfNodes(); ++i)
+        {
+            auto const bulk_node_id = _mapping_to_bulk_nodes.getComponent(i, 0);
+
+            MeshLib::Location const l{
+                _bulk_mesh_id, MeshLib::MeshItemType::Node, bulk_node_id};
+
+            auto const global_index = _dof_table_bulk.getGlobalIndex(
+                l, _variable_id_bulk, _component_id_bulk);
+            assert(global_index != NumLib::MeshComponentMap::nop);
+
+            bc_values.ids.push_back(global_index);
+            bc_values.values.push_back(_values.getComponent(i, 0));
+        }
+    }
+
+private:
+    std::size_t _bulk_mesh_id;
+    MeshLib::PropertyVector<double> const& _values;
+    std::unique_ptr<MeshLib::Mesh> _boundary_mesh;
+    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;
+
+    mutable bool _already_computed = false;
+};
+
+std::unique_ptr<NonuniformDirichletBoundaryCondition>
+createNonuniformDirichletBoundaryCondition(
+    BaseLib::ConfigTree const& config,
+    NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
+    int const component_id, const MeshLib::Mesh& bulk_mesh);
+
+}  // ProcessLib