ufz · endJunction · Aug 12, 2016 · Aug 3, 2016 · Aug 3, 2016 · Aug 3, 2016
diff --git a/Applications/CLI/Tests.cmake b/Applications/CLI/Tests.cmake
@@ -134,6 +134,30 @@ if(NOT OGS_USE_MPI)
             DIFF_DATA
             line_1_line_${mesh_size}.vtu line_${mesh_size}_neumann_pcs_0_ts_1_t_1.000000.vtu D1_left_N1_right pressure
         )
+
+        AddTest(
+            NAME GroundWaterFlowProcess_line_1_Robin_Right_Picard_${mesh_size}
+            PATH Elliptic/line_1_GroundWaterFlow
+            EXECUTABLE ogs
+            EXECUTABLE_ARGS line_${mesh_size}_robin_right_picard.prj
+            WRAPPER time
+            TESTER vtkdiff
+            ABSTOL 1e-14 RELTOL 1e-14
+            DIFF_DATA
+            line_1_line_${mesh_size}.vtu line_${mesh_size}_robin_right_picard_pcs_0_ts_1_t_1.000000.vtu D1_left_N1_right pressure
+        )
+
+        AddTest(
+            NAME GroundWaterFlowProcess_line_1_Robin_Left_Picard_${mesh_size}
+            PATH Elliptic/line_1_GroundWaterFlow
+            EXECUTABLE ogs
+            EXECUTABLE_ARGS line_${mesh_size}_robin_left_picard.prj
+            WRAPPER time
+            TESTER vtkdiff
+            ABSTOL 1e-14 RELTOL 1e-14
+            DIFF_DATA
+            line_1_line_${mesh_size}.vtu line_${mesh_size}_robin_left_picard_pcs_0_ts_1_t_1.000000.vtu D1_left_N1_right pressure
+        )
     endforeach()
 
     # Some Neumann BC tests

diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
@@ -13,6 +13,23 @@
 #include "BoundaryConditionConfig.h"
 #include "UniformDirichletBoundaryCondition.h"
 #include "UniformNeumannBoundaryCondition.h"
+#include "UniformRobinBoundaryCondition.h"
+
+static std::vector<MeshLib::Element*> getClonedElements(
+    MeshGeoToolsLib::BoundaryElementsSearcher& boundary_element_searcher,
+    GeoLib::GeoObject const& geometry)
+{
+    std::vector<MeshLib::Element*> elements =
+        boundary_element_searcher.getBoundaryElements(geometry);
+
+    // Deep copy all the elements, because the searcher might destroy the
+    // originals. Store pointers to the copies in the elements vector (i.e.,
+    // in-place modification).
+    for (auto& e : elements)
+        e = e->clone();
+
+    return elements;
+}
 
 namespace ProcessLib
 {
@@ -44,18 +61,18 @@ std::unique_ptr<BoundaryCondition> createBoundaryCondition(
     }
     else if (type == "UniformNeumann")
     {
-        std::vector<MeshLib::Element*> elements =
-            boundary_element_searcher.getBoundaryElements(config.geometry);
-
-        // Deep copy all the elements, because the searcher might destroy the
-        // originals. Store pointers to the copies in the elements vector (i.e.,
-        // in-place modification).
-        for (auto& e : elements)
-            e = e->clone();
-
         return createUniformNeumannBoundaryCondition(
-            config.config, std::move(elements), dof_table, variable_id,
-            config.component_id, integration_order, mesh.getDimension());
+            config.config,
+            getClonedElements(boundary_element_searcher, config.geometry),
+            dof_table, variable_id, config.component_id, integration_order,
+            mesh.getDimension());
+    }
+    else if (type == "UniformRobin") {
+        return createUniformRobinBoundaryCondition(
+            config.config,
+            getClonedElements(boundary_element_searcher, config.geometry),
+            dof_table, variable_id, config.component_id, integration_order,
+            mesh.getDimension());
     }
     else
     {

diff --git a/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.cpp
@@ -0,0 +1,38 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, 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 "UniformRobinBoundaryCondition.h"
+
+namespace ProcessLib
+{
+std::unique_ptr<UniformRobinBoundaryCondition>
+createUniformRobinBoundaryCondition(
+    BaseLib::ConfigTree const& config,
+    std::vector<MeshLib::Element*>&& elements,
+    NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
+    int const component_id, unsigned const integration_order,
+    unsigned const global_dim)
+{
+    DBUG("Constructing RobinBcConfig from config.");
+    //! \ogs_file_param{boundary_condition__type}
+    config.checkConfigParameter("type", "UniformRobin");
+
+    //! \ogs_file_param{boundary_condition__UniformRobin__alpha}
+    double const alpha = config.getConfigParameter<double>("alpha");
+    //! \ogs_file_param{boundary_condition__UniformRobin__u_0}
+    double const u_0 = config.getConfigParameter<double>("u_0");
+
+    return std::unique_ptr<UniformRobinBoundaryCondition>(
+        new UniformRobinBoundaryCondition(
+            integration_order, dof_table, variable_id, component_id, global_dim,
+            std::move(elements),
+            UniformRobinBoundaryConditionData{alpha, u_0}));
+}
+
+}  // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.h b/ProcessLib/BoundaryCondition/UniformRobinBoundaryCondition.h
@@ -0,0 +1
5D32
,48 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_UNIFORMROBINBOUNDARYCONDITION_H
+#define PROCESSLIB_UNIFORMROBINBOUNDARYCONDITION_H
+
+#include "GenericNaturalBoundaryCondition.h"
+#include "UniformRobinBoundaryConditionLocalAssembler.h"
+
+namespace MeshGeoToolsLib
+{
+class BoundaryElementsSearcher;
+}
+
+namespace ProcessLib
+{
+using UniformRobinBoundaryCondition = GenericNaturalBoundaryCondition<
+    UniformRobinBoundaryConditionData,
+    UniformRobinBoundaryConditionLocalAssembler>;
+
+/*! Creates a new uniform Robin boundary condition from the given data.
+ *
+ * The Robin boundary condition is given in the form
+ * \f$ \alpha \cdot [ u_0 - u(x) ] \f$,
+ * where the coefficients \f$ \alpha \f$ and \f$ u_0 \f$ are obtained from the
+ * \c config, and \f$ u \f$ is the unknown to which the boundary condition is
+ * applied.
+ *
+ * The value \f$ \alpha \cdot [ u_0 - u(x) ] \f$ is a flux. It replaces the
+ * integrand in the boundary integral for the variable \f$ u \f$.
+ */
+std::unique_ptr<UniformRobinBoundaryCondition>
+createUniformRobinBoundaryCondition(
+    BaseLib::ConfigTree const& config,
+    std::vector<MeshLib::Element*>&& elements,
+    NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
+    int const component_id, unsi
AE96
gned const integration_order,
+    unsigned const global_dim);
+
+}  // ProcessLib
+
+#endif  // PROCESSLIB_UNIFORMROBINBOUNDARYCONDITION_H
diff --git a/ProcessLib/BoundaryCondition/UniformRobinBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/UniformRobinBoundaryConditionLocalAssembler.h
@@ -0,0 +1,86 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_UNIFORMROBINBOUNDARYCONDITIONLOCALASSEMBLER_H
+#define PROCESSLIB_UNIFORMROBINBOUNDARYCONDITIONLOCALASSEMBLER_H
+
+#include "NumLib/DOF/DOFTableUtil.h"
+
+namespace ProcessLib
+{
+struct UniformRobinBoundaryConditionData final {
+    double const alpha;
+    double const u_0;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+          unsigned GlobalDim>
+class UniformRobinBoundaryConditionLocalAssembler final
+    : public GenericNaturalBoundaryConditionLocalAssembler<
+          ShapeFunction, IntegrationMethod, GlobalDim>
+{
+    using Base = GenericNaturalBoundaryConditionLocalAssembler<
+        ShapeFunction, IntegrationMethod, GlobalDim>;
+
+public:
+    UniformRobinBoundaryConditionLocalAssembler(
+        MeshLib::Element const& e, std::size_t const local_matrix_size,
+        unsigned const integration_order,
+        UniformRobinBoundaryConditionData const& data)
+        : Base(e, integration_order),
+          _data(data),
+          _local_K(local_matrix_size, local_matrix_size),
+          _local_rhs(local_matrix_size)
+    {
+    }
+
+    // TODO also implement derivative for Jacobian in Newton scheme.
+    void assemble(std::size_t const id,
+                  NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
+                  double const t, const GlobalVector& /*x*/, GlobalMatrix& K,
+                  GlobalVector& b) override
+    {
+        (void)t;  // TODO time-dependent Robin BCs
+
+        _local_K.setZero();
+        _local_rhs.setZero();
+
+        IntegrationMethod integration_method(Base::_integration_order);
+        std::size_t const n_integration_points =
+            integration_method.getNumberOfPoints();
+
+        for (std::size_t ip = 0; ip < n_integration_points; ++ip) {
+            auto const& sm = Base::_shape_matrices[ip];
+            auto const& wp = integration_method.getWeightedPoint(ip);
+
+            // flux = alpha * ( u_0 - u )
+            // adding a alpha term to the diagonal of the stiffness matrix
+            // and a alpha * u_0 term to the rhs vector
+            _local_K.diagonal().noalias() +=
+                sm.N * _data.alpha * sm.detJ * wp.getWeight();
+            _local_rhs.noalias() +=
+                sm.N * _data.alpha * _data.u_0 * sm.detJ * wp.getWeight();
+        }
+
+        auto const indices = NumLib::getIndices(id, dof_table_boundary);
+        K.add(NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices),
+              _local_K);
+        b.add(indices, _local_rhs);
+    }
+
+private:
+    UniformRobinBoundaryConditionData const& _data;
+
+    typename Base::NodalMatrixType _local_K;
+    typename Base::NodalVectorType _local_rhs;
+};
+
+}  // ProcessLib
+
+#endif  // PROCESSLIB_UNIFORMROBINBOUNDARYCONDITIONLOCALASSEMBLER_H
diff --git a/Tests/Data b/Tests/Data
+107 −0		Elliptic/line_1_GroundWaterFlow/line_1e1_robin_left_picard.prj
+107 −0		Elliptic/line_1_GroundWaterFlow/line_1e1_robin_right_picard.prj