Merge pull request #1952 from endJunction/NonuniformDirichlet

Nonuniform Dirichlet boundary conditions

Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformDirichlet/c_NonuniformDirichlet.md

+Declares a Dirichlet boundary condition that is non-uniform in space.
+The values are read from the given boundary mesh.

Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformDirichlet/t_field_name.md

+This specifies a nodal field defined on the surface mesh (see the \c mesh
+parameter of the nonuniform Dirichlet BC) and the values used for the boundary
+condition.

Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/NonuniformDirichlet/t_mesh.md

+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.

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,

ProcessLib/BoundaryCondition/BoundaryCondition.h

@@ -62,8 +62,6 @@ public:
         // Therefore there is nothing to do here.
     }
 
-    virtual void preTimestep(const double /*t*/) {}
-
     virtual ~BoundaryCondition() = default;
 };
 
@@ -107,6 +105,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,

ProcessLib/BoundaryCondition/BoundaryConditionCollection.cpp

@@ -43,10 +43,4 @@ void BoundaryConditionCollection::addBCsForProcessVariables(
     // object if needed.
     _dirichlet_bcs.resize(_boundary_conditions.size());
 }
-
-void BoundaryConditionCollection::preTimestep(const double t)
-{
-    for (auto& bc : _boundary_conditions)
-        bc->preTimestep(t);
-}
 }

ProcessLib/BoundaryCondition/BoundaryConditionCollection.h

@@ -44,8 +44,6 @@ public:
         NumLib::LocalToGlobalIndexMap const& dof_table,
         unsigned const integration_order);
 
-    void preTimestep(const double t);
-
 private:
     mutable std::vector<NumLib::IndexValueVector<GlobalIndexType>> _dirichlet_bcs;
     std::vector<std::unique_ptr<BoundaryCondition>> _boundary_conditions;

ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp

@@ -16,21 +16,9 @@
 
 namespace ProcessLib
 {
-void DirichletBoundaryCondition::preTimestep(const double /*t*/)
-{
-    if (_parameter.isTimeDependent())
-        _already_computed = false;
-}
-
 void DirichletBoundaryCondition::getEssentialBCValues(
     const double t, NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
 {
-    // TODO: Reenable when fixed ;)
-    //if (_already_computed)
-        //return;
-
-    _already_computed = true;
-
     SpatialPosition pos;
 
     bc_values.ids.clear();

ProcessLib/BoundaryCondition/DirichletBoundaryCondition.h

@@ -49,8 +49,6 @@ public:
         }
     }
 
-    void preTimestep(const double t) override;
-
     void getEssentialBCValues(
         const double t,
         NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
@@ -63,7 +61,6 @@ private:
     std::size_t const _mesh_id;
     int const _variable_id;
     int const _component_id;
-    mutable bool _already_computed = false;
 };
 
 std::unique_ptr<DirichletBoundaryCondition> createDirichletBoundaryCondition(

ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.cpp

+/**
+ * \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 const 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

ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h

+/**
+ * \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 getEssentialBCValues(
+        const double /*t*/,
+        NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override
+    {
+        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;
+};
+
+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

ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp

@@ -67,7 +67,7 @@ createNonuniformNeumannBoundaryCondition(
         OGS_FATAL("`%s' is not a one-component field.", field_name.c_str());
     }
 
-    std::string mapping_to_bulk_nodes_property = "OriginalSubsurfaceNodeIDs";
+    std::string const 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);

ProcessLib/GroundwaterFlow/Tests.cmake

@@ -585,3 +585,15 @@ AddTest(
     DIFF_DATA
     inhomogeneous_permeability.vtu neumann_nonuniform_pcs_0_ts_1_t_1.000000.vtu mass_flux_ref mass_flux
 )
+
+AddTest(
+    NAME GroundWaterFlowProcess_Dirichlet_nonuniform
+    PATH Elliptic/nonuniform_bc_Groundwaterflow
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS dirichlet_nonuniform.prj
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    ABSTOL 1e-14 RELTOL 1e-16
+    DIFF_DATA
+    expected_dirichlet_nonuniform_pcs_0_ts_1_t_1.000000.vtu dirichlet_nonuniform_pcs_0_ts_1_t_1.000000.vtu pressure pressure
+)

ProcessLib/Process.cpp

@@ -245,7 +245,6 @@ void Process::preTimestep(GlobalVector const& x, const double t,
 
     MathLib::LinAlg::setLocalAccessibleVector(x);
     preTimestepConcreteProcess(x, t, delta_t);
-    _boundary_conditions.preTimestep(t);
 }
 
 void Process::postTimestep(GlobalVector const& x)

Data @ accf02ea

-Subproject commit 094bd3a2516375c774efd69c99f36404120bfc73
+Subproject commit accf02eaaf1c4c0581f1a4aeb91fc86e714a2966