diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/Neumann/t_area_parameter.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/Neumann/t_area_parameter.md
new file mode 100644
index 0000000000000000000000000000000000000000..302b1130109f71035b28c2c62f4269669d5209be
--- /dev/null
+++ b/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.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/Robin/t_area_parameter.md b/Documentation/ProjectFile/prj/process_variables/process_variable/boundary_conditions/boundary_condition/Robin/t_area_parameter.md
new file mode 100644
index 0000000000000000000000000000000000000000..302b1130109f71035b28c2c62f4269669d5209be
--- /dev/null
+++ b/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.
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
index bc716873285ea35bca49501ba9018b63587c491a..a3dd7eee7254c70cfaa88f531356ebc1d2bf368b 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
+++ b/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"))
{
diff --git a/ProcessLib/BoundaryCondition/HCNonAdvectiveFreeComponentFlowBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/HCNonAdvectiveFreeComponentFlowBoundaryCondition.cpp
index df566bfa3cc6ec4286216a0951ac1f0247528d3a..fee37c8c8cd8b8238caa5f0f388794e1f76bb3c9 100644
--- a/ProcessLib/BoundaryCondition/HCNonAdvectiveFreeComponentFlowBoundaryCondition.cpp
+++ b/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");
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp
index 3bec7dca1cbfc9be01e2e2aaca9c2bb838b654cf..0e51abfbe2f0f1b3e3f9a155d45bbaed1f17f24e 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.cpp
+++ b/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
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h b/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h
index 71ab4e441e4704755bd40679863563fbd6b097a0..1b42a8f7d22345abc78e831e0ea73b1cf1151ee1 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryCondition.h
+++ b/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(
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
index 11939130fa9d3a520edce452c4de0c75adaacdd4..1afc028a6742b44d221c04d07cc9160232b44bfe 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
+++ b/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 @@ public:
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,16 +64,27 @@ public:
// 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);
@@ -75,7 +92,8 @@ public:
}
private:
- ParameterLib::Parameter<double> const& _neumann_bc_parameter;
+ NeumannBoundaryConditionData const& _data;
+
NodalVectorType _local_rhs;
public:
diff --git a/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp
index 2d90c71b0444b10183d7f51694563120e43e0b7e..ebb59a87d69b151df4882ab82eb5cedfeaee0888 100644
--- a/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/RobinBoundaryCondition.cpp
@@ -25,6 +25,14 @@ 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}
@@ -32,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
@@ -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
diff --git a/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
index 64d9190bc2c520bae87f7ad6c137685d8e360d5d..60ba7d4a3a3dd3971b807101c8b0bd8d551770ca 100644
--- a/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
+++ b/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);
diff --git a/ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryCondition.cpp
index 583821c38a0490575bee3f37c6d7852be659cdfb..2f3bdfe81a66ed18a979da13a5976c4f37ec4544 100644
--- a/ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryCondition.cpp
+++ b/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");
diff --git a/ProcessLib/ComponentTransport/Tests.cmake b/ProcessLib/ComponentTransport/Tests.cmake
index 9ec09fc1d8f789f6424e15e9604de4301cfce7a9..6bf1362f711d7062d5efe2a04442a4a6b5862b58 100644
--- a/ProcessLib/ComponentTransport/Tests.cmake
+++ b/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(
diff --git a/ProcessLib/GroundwaterFlow/Tests.cmake b/ProcessLib/GroundwaterFlow/Tests.cmake
index 19916a48a90ba652bdd269e28a6015cd50214a0b..54408f819e49fef7f3250e46a4c267d82964a21d 100644
--- a/ProcessLib/GroundwaterFlow/Tests.cmake
+++ b/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
diff --git a/Tests/Data/Elliptic/quarter_circle/quarter_circle_neumann.prj b/Tests/Data/Elliptic/quarter_circle/quarter_circle_neumann.prj
new file mode 100644
index 0000000000000000000000000000000000000000..a777723bac1cde8b2c095640091b22be0c86030b
--- /dev/null
+++ b/Tests/Data/Elliptic/quarter_circle/quarter_circle_neumann.prj
@@ -0,0 +1,120 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <meshes>
+ <mesh>quarter_circle_r_1.vtu</mesh>
+ <mesh>quarter_circle_r_1_circle_line.vtu</mesh>
+ <mesh>quarter_circle_r_1_center_point.vtu</mesh>
+ </meshes>
+ <processes>
+ <process>
+ <name>GW23</name>
+ <type>GROUNDWATER_FLOW</type>
+ <integration_order>2</integration_order>
+ <hydraulic_conductivity>K</hydraulic_conductivity>
+ <process_variables>
+ <process_variable>pressure</process_variable>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable internal_name="darcy_velocity" output_name="v"/>
+ </secondary_variables>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="GW23">
+ <nonlinear_solver>basic_picard</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstol>1.e-6</abstol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <time_stepping>
+ <type>SingleStep</type>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>neumann_in_center</prefix>
+ <variables>
+ <variable> pressure </variable>
+ <variable> v </variable>
+ </variables>
+ </output>
+ </time_loop>
+ <parameters>
+ <parameter>
+ <name>K</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>p0</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>inflow</name>
+ <type>Constant</type>
+ <value>0.125</value>
+ </parameter>
+ <parameter>
+ <name>p_Dirichlet_barrel</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>one</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ </parameters>
+ <process_variables>
+ <process_variable>
+ <name>pressure</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>p0</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <mesh>quarter_circle_r_1_circle_line</mesh>
+ <type>Dirichlet</type>
+ <parameter>p_Dirichlet_barrel</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <mesh>quarter_circle_r_1_center_point</mesh>
+ <type>Neumann</type>
+ <parameter>inflow</parameter>
+ <area_parameter>one</area_parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_picard</name>
+ <type>Picard</type>
+ <max_iter>10</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+ <eigen>
+ <solver_type>CG</solver_type>
+ <precon_type>DIAGONAL</precon_type>
+ <max_iteration_step>10000</max_iteration_step>
+ <error_tolerance>1e-16</error_tolerance>
+ </eigen>
+ <petsc>
+ <prefix>gw</prefix>
+ <parameters>-gw_ksp_type cg -gw_pc_type bjacobi -gw_ksp_rtol 1e-16 -gw_ksp_max_it 10000</parameters>
+ </petsc>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Elliptic/quarter_circle/quarter_circle_nodal_source_term.prj b/Tests/Data/Elliptic/quarter_circle/quarter_circle_nodal_source_term.prj
new file mode 100644
index 0000000000000000000000000000000000000000..9f14e245fd26d8ad2285946c2f649a79780c9c14
--- /dev/null
+++ b/Tests/Data/Elliptic/quarter_circle/quarter_circle_nodal_source_term.prj
@@ -0,0 +1,116 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <meshes>
+ <mesh>quarter_circle_r_1.vtu</mesh>
+ <mesh>quarter_circle_r_1_circle_line.vtu</mesh>
+ <mesh>quarter_circle_r_1_center_point.vtu</mesh>
+ </meshes>
+ <processes>
+ <process>
+ <name>GW23</name>
+ <type>GROUNDWATER_FLOW</type>
+ <integration_order>2</integration_order>
+ <hydraulic_conductivity>K</hydraulic_conductivity>
+ <process_variables>
+ <process_variable>pressure</process_variable>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable internal_name="darcy_velocity" output_name="v"/>
+ </secondary_variables>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="GW23">
+ <nonlinear_solver>basic_picard</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstol>1.e-6</abstol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <time_stepping>
+ <type>SingleStep</type>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>nodal_source_term_in_center</prefix>
+ <variables>
+ <variable> pressure </variable>
+ <variable> v </variable>
+ </variables>
+ </output>
+ </time_loop>
+ <parameters>
+ <parameter>
+ <name>K</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>p0</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>inflow</name>
+ <type>Constant</type>
+ <value>0.125</value>
+ </parameter>
+ <parameter>
+ <name>p_Dirichlet_barrel</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ </parameters>
+ <process_variables>
+ <process_variable>
+ <name>pressure</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>p0</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <mesh>quarter_circle_r_1_circle_line</mesh>
+ <type>Dirichlet</type>
+ <parameter>p_Dirichlet_barrel</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ <source_terms>
+ <source_term>
+ <mesh>quarter_circle_r_1_center_point</mesh>
+ <type>Nodal</type>
+ <parameter>inflow</parameter>
+ </source_term>
+ </source_terms>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_picard</name>
+ <type>Picard</type>
+ <max_iter>10</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+ <eigen>
+ <solver_type>CG</solver_type>
+ <precon_type>DIAGONAL</precon_type>
+ <max_iteration_step>10000</max_iteration_step>
+ <error_tolerance>1e-16</error_tolerance>
+ </eigen>
+ <petsc>
+ <prefix>gw</prefix>
+ <parameters>-gw_ksp_type cg -gw_pc_type bjacobi -gw_ksp_rtol 1e-16 -gw_ksp_max_it 10000</parameters>
+ </petsc>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Elliptic/quarter_circle/quarter_circle_r_1.vtu b/Tests/Data/Elliptic/quarter_circle/quarter_circle_r_1.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..01989dd392117b05855ae028b86abd1068fcb5f3
--- /dev/null
+++ b/Tests/Data/Elliptic/quarter_circle/quarter_circle_r_1.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:303f677330f145b500d48c9f87716881a4ae33952b1966b043dd8184ef56a4cb
+size 6412
diff --git a/Tests/Data/Elliptic/quarter_circle/quarter_circle_r_1_center_point.vtu b/Tests/Data/Elliptic/quarter_circle/quarter_circle_r_1_center_point.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..9b78fb0812d16195a6292413d7ae03ab0116965b
--- /dev/null
+++ b/Tests/Data/Elliptic/quarter_circle/quarter_circle_r_1_center_point.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:0467be9935d5e71213e863db4546de7a3c512eacb05b699e7da9db266c737723
+size 1786
diff --git a/Tests/Data/Elliptic/quarter_circle/quarter_circle_r_1_circle_line.vtu b/Tests/Data/Elliptic/quarter_circle/quarter_circle_r_1_circle_line.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..9d9d4dc79e25a50e188bb0f5783d533db3ff2c7b
--- /dev/null
+++ b/Tests/Data/Elliptic/quarter_circle/quarter_circle_r_1_circle_line.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d2353f9fe025c1c24cf34d6314f70390ae6482adf79cb3f68eaa65c02430dcb1
+size 3564
diff --git a/Tests/Data/Elliptic/quarter_disc/quarter_disc_neumann.prj b/Tests/Data/Elliptic/quarter_disc/quarter_disc_neumann.prj
new file mode 100644
index 0000000000000000000000000000000000000000..79fcd1c3af86cb920d7f67869c08bb072eaf1a0d
--- /dev/null
+++ b/Tests/Data/Elliptic/quarter_disc/quarter_disc_neumann.prj
@@ -0,0 +1,120 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <meshes>
+ <mesh>quarter_disc_r_1.vtu</mesh>
+ <mesh>quarter_disc_r_1_line_in_center.vtu</mesh>
+ <mesh>quarter_disc_r_1_outer_surface.vtu</mesh>
+ </meshes>
+ <processes>
+ <process>
+ <name>GW23</name>
+ <type>GROUNDWATER_FLOW</type>
+ <integration_order>2</integration_order>
+ <hydraulic_conductivity>K</hydraulic_conductivity>
+ <process_variables>
+ <process_variable>pressure</process_variable>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable internal_name="darcy_velocity" output_name="v"/>
+ </secondary_variables>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="GW23">
+ <nonlinear_solver>basic_picard</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstol>1.e-6</abstol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <time_stepping>
+ <type>SingleStep</type>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>neumann_along_line</prefix>
+ <variables>
+ <variable> pressure </variable>
+ <variable> v </variable>
+ </variables>
+ </output>
+ </time_loop>
+ <parameters>
+ <parameter>
+ <name>K</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>p0</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>inflow</name>
+ <type>Constant</type>
+ <value>0.125</value>
+ </parameter>
+ <parameter>
+ <name>p_Dirichlet_barrel</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>one</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ </parameters>
+ <process_variables>
+ <process_variable>
+ <name>pressure</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>p0</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <mesh>quarter_disc_r_1_outer_surface</mesh>
+ <type>Dirichlet</type>
+ <parameter>p_Dirichlet_barrel</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <mesh>quarter_disc_r_1_line_in_center</mesh>
+ <type>Neumann</type>
+ <parameter>inflow</parameter>
+ <area_parameter>one</area_parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_picard</name>
+ <type>Picard</type>
+ <max_iter>10</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+ <eigen>
+ <solver_type>CG</solver_type>
+ <precon_type>DIAGONAL</precon_type>
+ <max_iteration_step>10000</max_iteration_step>
+ <error_tolerance>1e-16</error_tolerance>
+ </eigen>
+ <petsc>
+ <prefix>gw</prefix>
+ <parameters>-gw_ksp_type cg -gw_pc_type bjacobi -gw_ksp_rtol 1e-16 -gw_ksp_max_it 10000</parameters>
+ </petsc>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Elliptic/quarter_disc/quarter_disc_r_1.vtu b/Tests/Data/Elliptic/quarter_disc/quarter_disc_r_1.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..c3d8ddb70ef07539b59ea4d656748b2fa8af54a9
--- /dev/null
+++ b/Tests/Data/Elliptic/quarter_disc/quarter_disc_r_1.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:69176df7b59a265e99a206730ba9c63d671a244489fc595721d88d1c99afd4b1
+size 2301717
diff --git a/Tests/Data/Elliptic/quarter_disc/quarter_disc_r_1_line_in_center.vtu b/Tests/Data/Elliptic/quarter_disc/quarter_disc_r_1_line_in_center.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..d73c1bf13d7b07a4f57323a870228d7ce661a04e
--- /dev/null
+++ b/Tests/Data/Elliptic/quarter_disc/quarter_disc_r_1_line_in_center.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a45d38ddbf7f1822c0de5d1dbcba93d9da372b2043330a62345d706506866392
+size 1683
diff --git a/Tests/Data/Elliptic/quarter_disc/quarter_disc_r_1_outer_surface.vtu b/Tests/Data/Elliptic/quarter_disc/quarter_disc_r_1_outer_surface.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..7e2ac534fc01d7d577a465bee518fc55665769b5
--- /dev/null
+++ b/Tests/Data/Elliptic/quarter_disc/quarter_disc_r_1_outer_surface.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:448f143d9e5aec970f421291a191fe3bf300382c2caa3a47d9b04fbb301cb542
+size 26910