[PL/BC] Enable 'integral_measure' in Robin-BC.

ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp

@@ -25,11 +25,11 @@ 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() + 1 != global_dim)
+    if (bc_mesh.getDimension() >= global_dim)
     {
         OGS_FATAL(
-            "The dimension (%d) of the given boundary mesh '%s' is not by one "
-            "lower than the bulk dimension (%d).",
+            "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);
     }
 
@@ -40,10 +40,20 @@ std::unique_ptr<RobinBoundaryCondition> createRobinBoundaryCondition(
 
     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
@@ -61,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 @@ public:
             _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);