Adapt Robin / Neumann boundary conditions to user needs.

Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/Neumann/t_area_parameter.md

@@ -0,0 +1,8 @@
+Typically boundary conditions are set on domains with dimension exactly one
+lower than the bulk mesh dimension. In some cases the boundary condition should
+be set on entities that have lower dimension than one lower. For instance
+boundary conditions on points in 2d or 3d domains or boundary conditions on
+lines in 3d domains.
+
+In order not to integrate over a null set the area tag can be used to specify
+the area the boundary condition should be effective.

Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/Robin/t_area_parameter.md

@@ -0,0 +1,8 @@
+Typically boundary conditions are set on domains with dimension exactly one
+lower than the bulk mesh dimension. In some cases the boundary condition should
+be set on entities that have lower dimension than one lower. For instance
+boundary conditions on points in 2d or 3d domains or boundary conditions on
+lines in 3d domains.
+
+In order not to integrate over a null set the area tag can be used to specify
+the area the boundary condition should be effective.

ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h

@@ -43,14 +43,6 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
             dof_table_bulk.getNumberOfVariableComponents(variable_id));
     }
 
-    if (_bc_mesh.getDimension() + 1 != global_dim)
-    {
-        OGS_FATAL(
-            "The dimension (%d) of the given boundary mesh '%s' is not by one "
-            "lower than the bulk dimension (%d).",
-            _bc_mesh.getDimension(), _bc_mesh.getName().c_str(), global_dim);
-    }
-
     if (!_bc_mesh.getProperties().template existsPropertyVector<std::size_t>(
             "bulk_node_ids"))
     {

ProcessLib/BoundaryCondition/HCNonAdvectiveFreeComponentFlowBoundaryCondition.cpp

@@ -30,6 +30,14 @@ createHCNonAdvectiveFreeComponentFlowBoundaryCondition(
     config.checkConfigParameter("type",
                                 "HCNonAdvectiveFreeComponentFlowBoundary");
 
+    if (bc_mesh.getDimension() + 1 != global_dim)
+    {
+        OGS_FATAL(
+            "The dimension (%d) of the given boundary mesh '%s' is not by one "
+            "lower than the bulk dimension (%d).",
+            bc_mesh.getDimension(), bc_mesh.getName().c_str(), global_dim);
+    }
+
     auto const boundary_permeability_name =
         //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__HCNonAdvectiveFreeComponentFlowBoundary__parameter}
         config.getConfigParameter<std::string>("parameter");

ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp

@@ -31,8 +31,8 @@ std::unique_ptr<NeumannBoundaryCondition> createNeumannBoundaryCondition(
     auto const& param = ParameterLib::findParameter<double>(
         param_name, parameters, 1, &bc_mesh);
 
-    // In case of partitioned mesh the boundary could be empty, i.e. there is no
-    // boundary condition.
+    // 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.
@@ -45,9 +45,29 @@ std::unique_ptr<NeumannBoundaryCondition> createNeumannBoundaryCondition(
     }
 #endif  // USE_PETSC
 
+    ParameterLib::Parameter<double> const* integral_measure(nullptr);
+    if (global_dim - bc_mesh.getDimension() != 1)
+    {
+        //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__Neumann__area_parameter}
+        auto const area_parameter_name =
+            config.getConfigParameter<std::string>("area_parameter");
+        DBUG("area parameter name '%s'", area_parameter_name.c_str());
+        integral_measure = &ParameterLib::findParameter<double>(
+            area_parameter_name, parameters, 1, &bc_mesh);
+    }
+
+    if (bc_mesh.getDimension() >= global_dim)
+    {
+        OGS_FATAL(
+            "The dimension (%d) of the given boundary mesh '%s' is not lower "
+            "than the bulk dimension (%d).",
+            bc_mesh.getDimension(), bc_mesh.getName().c_str(), global_dim);
+    }
+
     return std::make_unique<NeumannBoundaryCondition>(
         integration_order, shapefunction_order, dof_table, variable_id,
-        component_id, global_dim, bc_mesh, param);
+        component_id, global_dim, bc_mesh,
+        NeumannBoundaryConditionData{param, integral_measure});
 }
 
 }  // namespace ProcessLib

ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h

@@ -17,7 +17,7 @@
 namespace ProcessLib
 {
 using NeumannBoundaryCondition =
-    GenericNaturalBoundaryCondition<ParameterLib::Parameter<double> const&,
+    GenericNaturalBoundaryCondition<NeumannBoundaryConditionData,
                                     NeumannBoundaryConditionLocalAssembler>;
 
 std::unique_ptr<NeumannBoundaryCondition> createNeumannBoundaryCondition(

ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h

@@ -18,6 +18,12 @@
 
 namespace ProcessLib
 {
+struct NeumannBoundaryConditionData final
+{
+    ParameterLib::Parameter<double> const& neumann_bc_parameter;
+    ParameterLib::Parameter<double> const* const integral_measure;
+};
+
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 class NeumannBoundaryConditionLocalAssembler final
@@ -37,9 +43,9 @@ class NeumannBoundaryConditionLocalAssembler final
         std::size_t const local_matrix_size,
         bool const is_axially_symmetric,
         unsigned const integration_order,
-        ParameterLib::Parameter<double> const& neumann_bc_parameter)
+        NeumannBoundaryConditionData const& data)
         : Base(e, is_axially_symmetric, integration_order),
-          _neumann_bc_parameter(neumann_bc_parameter),
+          _data(data),
           _local_rhs(local_matrix_size)
     {
     }
@@ -58,24 +64,36 @@ class NeumannBoundaryConditionLocalAssembler final
         // Get element nodes for the interpolation from nodes to integration
         // point.
         NodalVectorType parameter_node_values =
-            _neumann_bc_parameter.getNodalValuesOnElement(Base::_element, t)
+            _data.neumann_bc_parameter
+                .getNodalValuesOnElement(Base::_element, t)
                 .template topRows<ShapeFunction::MeshElement::n_all_nodes>();
 
+        ParameterLib::SpatialPosition position;
+        position.setElementID(Base::_element.getID());
+
+        double integral_measure = 1.0;
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
+            position.setIntegrationPoint(ip);
             auto const& ip_data = Base::_ns_and_weights[ip];
+            auto const& N = ip_data.N;
+            auto const& w = ip_data.weight;
 
-            _local_rhs.noalias() += ip_data.N *
-                                    parameter_node_values.dot(ip_data.N) *
-                                    ip_data.weight;
+            if (_data.integral_measure)
+            {
+                integral_measure = (*_data.integral_measure)(t, position)[0];
+            }
+            _local_rhs.noalias() +=
+                N * parameter_node_values.dot(N) * w * integral_measure;
         }
 
         auto const indices = NumLib::getIndices(id, dof_table_boundary);
         b.add(indices, _local_rhs);
     }
 
 private:
-    ParameterLib::Parameter<double> const& _neumann_bc_parameter;
+    NeumannBoundaryConditionData const& _data;
+
     NodalVectorType _local_rhs;
 
 public:

ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp

@@ -25,17 +25,35 @@ std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
     //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
     config.checkConfigParameter("type", "Robin");
 
+    if (bc_mesh.getDimension() >= global_dim)
+    {
+        OGS_FATAL(
+            "The dimension (%d) of the given boundary mesh '%s' is not lower "
+            "than the bulk dimension (%d).",
+            bc_mesh.getDimension(), bc_mesh.getName().c_str(), global_dim);
+    }
+
     //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__Robin__alpha}
     auto const alpha_name = config.getConfigParameter<std::string>("alpha");
     //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__Robin__u_0}
     auto const u_0_name = config.getConfigParameter<std::string>("u_0");
 
     auto const& alpha = ParameterLib::findParameter<double>(
         alpha_name, parameters, 1, &bc_mesh);
-
     auto const& u_0 =
         ParameterLib::findParameter<double>(u_0_name, parameters, 1, &bc_mesh);
 
+    ParameterLib::Parameter<double> const* integral_measure(nullptr);
+    if (global_dim - bc_mesh.getDimension() != 1)
+    {
+        //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__Robin__area_parameter}
+        auto const area_parameter_name =
+            config.getConfigParameter<std::string>("area_parameter");
+        DBUG("area parameter name '%s'", area_parameter_name.c_str());
+        integral_measure = &ParameterLib::findParameter<double>(
+            area_parameter_name, parameters, 1, &bc_mesh);
+    }
+
     // In case of partitioned mesh the boundary could be empty, i.e. there is no
     // boundary condition.
 #ifdef USE_PETSC
@@ -53,7 +71,7 @@ std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
     return std::make_unique<RobinBoundaryCondition>(
         integration_order, shapefunction_order, dof_table, variable_id,
         component_id, global_dim, bc_mesh,
-        RobinBoundaryConditionData{alpha, u_0});
+        RobinBoundaryConditionData{alpha, u_0, integral_measure});
 }
 
 }  // namespace ProcessLib

ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h

@@ -20,6 +20,7 @@ struct RobinBoundaryConditionData final
 {
     ParameterLib::Parameter<double> const& alpha;
     ParameterLib::Parameter<double> const& u_0;
+    ParameterLib::Parameter<double> const* const integral_measure;
 };
 
 template <typename ShapeFunction, typename IntegrationMethod,
@@ -64,17 +65,29 @@ class RobinBoundaryConditionLocalAssembler final
             _data.u_0.getNodalValuesOnElement(Base::_element, t)
                 .template topRows<ShapeFunction::MeshElement::n_all_nodes>();
 
+        ParameterLib::SpatialPosition position;
+        position.setElementID(Base::_element.getID());
+
         for (unsigned ip = 0; ip < n_integration_points; ++ip)
         {
+            position.setIntegrationPoint(ip);
             auto const& ip_data = Base::_ns_and_weights[ip];
             auto const& N = ip_data.N;
             auto const& w = ip_data.weight;
 
+            double integral_measure = 1.0;
+            if (_data.integral_measure)
+            {
+                integral_measure = (*_data.integral_measure)(t, position)[0];
+            }
+
             // 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() += N * alpha.dot(N) * w;
-            _local_rhs.noalias() += N * alpha.dot(N) * u_0.dot(N) * w;
+            _local_K.diagonal().noalias() +=
+                N * alpha.dot(N) * w * integral_measure;
+            _local_rhs.noalias() +=
+                N * alpha.dot(N) * u_0.dot(N) * w * integral_measure;
         }
 
         auto const indices = NumLib::getIndices(id, dof_table_boundary);

ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryCondition.cpp

@@ -33,6 +33,14 @@ createVariableDependentNeumannBoundaryCondition(
     }
     assert(variable_id == 0 || variable_id == 1);
 
+    if (bc_mesh.getDimension() + 1 != global_dim)
+    {
+        OGS_FATAL(
+            "The dimension (%d) of the given boundary mesh '%s' is not by one "
+            "lower than the bulk dimension (%d).",
+            bc_mesh.getDimension(), bc_mesh.getName().c_str(), global_dim);
+    }
+
     auto const constant_name =
         //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__VariableDependentNeumann__constant_name}
         config.getConfigParameter<std::string>("constant_name");

ProcessLib/ComponentTransport/Tests.cmake

@@ -651,13 +651,13 @@ AddTest(
     REQUIREMENTS OGS_USE_MPI
     DIFF_DATA
     TracerSimulation_pcs_1_ts_100_t_100000_000000_0_expected.vtu TracerSimulation_pcs_1_ts_100_t_100000_000000_0.vtu Cs Cs 1e-10 1e-16
-    TracerSimulation_pcs_1_ts_100_t_100000_000000_0_expected.vtu TracerSimulation_pcs_1_ts_100_t_100000_000000_0.vtu pressure pressure 1e-6 1e-10
+    TracerSimulation_pcs_1_ts_100_t_100000_000000_0_expected.vtu TracerSimulation_pcs_1_ts_100_t_100000_000000_0.vtu pressure pressure 2.5e-5 1.7e-9
     TracerSimulation_pcs_1_ts_100_t_100000_000000_1_expected.vtu TracerSimulation_pcs_1_ts_100_t_100000_000000_1.vtu Cs Cs 1e-10 1e-16
-    TracerSimulation_pcs_1_ts_100_t_100000_000000_1_expected.vtu TracerSimulation_pcs_1_ts_100_t_100000_000000_1.vtu pressure pressure 1e-6 1e-10
+    TracerSimulation_pcs_1_ts_100_t_100000_000000_1_expected.vtu TracerSimulation_pcs_1_ts_100_t_100000_000000_1.vtu pressure pressure 2.5e-5 1.7e-9
     TracerSimulation_pcs_1_ts_100_t_100000_000000_2_expected.vtu TracerSimulation_pcs_1_ts_100_t_100000_000000_2.vtu Cs Cs 1e-10 1e-16
-    TracerSimulation_pcs_1_ts_100_t_100000_000000_2_expected.vtu TracerSimulation_pcs_1_ts_100_t_100000_000000_2.vtu pressure pressure 1e-6 1e-10
+    TracerSimulation_pcs_1_ts_100_t_100000_000000_2_expected.vtu TracerSimulation_pcs_1_ts_100_t_100000_000000_2.vtu pressure pressure 2.5e-5 1.7e-9
     TracerSimulation_pcs_1_ts_100_t_100000_000000_3_expected.vtu TracerSimulation_pcs_1_ts_100_t_100000_000000_3.vtu Cs Cs 1e-10 1e-16
-    TracerSimulation_pcs_1_ts_100_t_100000_000000_3_expected.vtu TracerSimulation_pcs_1_ts_100_t_100000_000000_3.vtu pressure pressure 1e-6 1e-10
+    TracerSimulation_pcs_1_ts_100_t_100000_000000_3_expected.vtu TracerSimulation_pcs_1_ts_100_t_100000_000000_3.vtu pressure pressure 2.5e-5 1.7e-9
 )
 
 AddTest(

ProcessLib/GroundwaterFlow/Tests.cmake

@@ -708,6 +708,30 @@ AddTest(
     square_1x1_quad_1e3_volumetricsourceterm_analytical_solution.vtu square_1e3_volumetricsourceterm_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure 1e-10 1e-11
 )
 
+AddTest(
+    NAME GroundWaterFlowProcess_NeumannBC_Along_Line_in_3D_domain
+    PATH Elliptic/quarter_disc
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS quarter_disc_neumann.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    DIFF_DATA
+    quarter_disc_r_1.vtu neumann_along_line_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure 6e-5 0
+)
+
+AddTest(
+    NAME GroundWaterFlowProcess_NeumannBC_In_Center_Point_2D_domain
+    PATH Elliptic/quarter_circle
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS quarter_circle_neumann.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    DIFF_DATA
+    quarter_circle_r_1.vtu neumann_in_center_pcs_0_ts_1_t_1.000000.vtu pressure_nodal_source_term pressure 1e-14 1e-14
+)
+
 AddTest(
     NAME GroundWaterFlowProcess_VolumetricSourceTerm_sin_x_sin_y_square_1e3
     PATH Elliptic/square_1x1_GroundWaterFlow