diff --git a/Applications/CLI/Tests.cmake b/Applications/CLI/Tests.cmake
index 4c5d737e915c60882cb119e48fabbdf9b0f0002b..daf0ae6c3e41c640757bfeb6c7d59856b9a250de 100644
--- a/Applications/CLI/Tests.cmake
+++ b/Applications/CLI/Tests.cmake
@@ -14,6 +14,15 @@ foreach(mesh_size 1e0 1e1 1e2 1e3)
DIFF_DATA cube_${mesh_size}_result.dat
DATA cube_${mesh_size}.prj cube_1x1x1_hex_${mesh_size}.vtu cube_1x1x1.gml
)
+
+ AddTest(
+ NAME GroundWaterFlowProcess_cube_1x1x1_Neumann_${mesh_size}
+ PATH Elliptic/cube_1x1x1_GroundWaterFlow
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS cube_${mesh_size}_neumann.prj
+ WRAPPER time
+ DATA cube_${mesh_size}_neumann.prj cube_1x1x1_hex_${mesh_size}.vtu cube_1x1x1.gml
+ )
endforeach()
foreach(mesh_size 1e4 2e4 3e4 4e4 5e4 1e5 1e6)
@@ -27,6 +36,15 @@ foreach(mesh_size 1e4 2e4 3e4 4e4 5e4 1e5 1e6)
DIFF_DATA cube_${mesh_size}_result.dat
DATA cube_${mesh_size}.prj cube_1x1x1_hex_${mesh_size}.vtu cube_1x1x1.gml
)
+
+ AddTest(
+ NAME LARGE_GroundWaterFlowProcess_cube_1x1x1_Neumann_${mesh_size}
+ PATH Elliptic/cube_1x1x1_GroundWaterFlow
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS cube_${mesh_size}_neumann.prj
+ WRAPPER time
+ DATA cube_${mesh_size}_neumann.prj cube_1x1x1_hex_${mesh_size}.vtu cube_1x1x1.gml
+ )
endforeach()
# SQUARE 1x1 GROUNDWATER FLOW TESTS
@@ -41,6 +59,15 @@ foreach(mesh_size 1e0 1e1 1e2 1e3 1e4)
DIFF_DATA square_${mesh_size}_result.dat
DATA square_${mesh_size}.prj square_1x1_quad_${mesh_size}.vtu square_1x1.gml
)
+
+ AddTest(
+ NAME GroundWaterFlowProcess_square_1x1_Neumann_${mesh_size}
+ PATH Elliptic/square_1x1_GroundWaterFlow
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS square_${mesh_size}_neumann.prj
+ WRAPPER time
+ DATA square_${mesh_size}_neumann.prj square_1x1_quad_${mesh_size}.vtu square_1x1.gml
+ )
endforeach()
foreach(mesh_size 1e5 1e6)
@@ -54,4 +81,13 @@ foreach(mesh_size 1e5 1e6)
DIFF_DATA square_${mesh_size}_result.dat
DATA square_${mesh_size}.prj square_1x1_quad_${mesh_size}.vtu square_1x1.gml
)
+
+ AddTest(
+ NAME LARGE_GroundWaterFlowProcess_square_1x1_Neumann_${mesh_size}
+ PATH Elliptic/square_1x1_GroundWaterFlow
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS square_${mesh_size}_neumann.prj
+ WRAPPER time
+ DATA square_${mesh_size}_neumann.prj square_1x1_quad_${mesh_size}.vtu square_1x1.gml
+ )
endforeach()
diff --git a/MathLib/ConstantFunction.h b/MathLib/ConstantFunction.h
new file mode 100644
index 0000000000000000000000000000000000000000..707ecbfa3e52f2fb7b155efaa670627bac807b45
--- /dev/null
+++ b/MathLib/ConstantFunction.h
@@ -0,0 +1,46 @@
+/**
+ * \author Norihiro Watanabe
+ * \date 2013-08-13
+ *
+ * \copyright
+ * Copyright (c) 2012-2015, 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 MATHLIB_CONSTANTFUNCTION_H_
+#define MATHLIB_CONSTANTFUNCTION_H_
+
+#include <array>
+#include <numeric>
+
+namespace MathLib
+{
+
+/**
+ * A constant function.
+ * \f[
+ * f(x_1,...,x_k)=a_0
+ * \f]
+ */
+template <typename T_TYPE>
+class ConstantFunction
+{
+public:
+ explicit ConstantFunction(T_TYPE const& value)
+ : _value(value)
+ {}
+
+ T_TYPE operator()() const
+ {
+ return _value;
+ }
+private:
+ T_TYPE const _value;
+};
+
+}
+
+#endif // MATHLIB_CONSTANTFUNCTION_H_
diff --git a/ProcessLib/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlowProcess.h
index bb089f92c19a12dd3838e87af0774c6f37d771fd..59ed1e35b12b52bb2583eeb5b30d3482abbebbd6 100644
--- a/ProcessLib/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlowProcess.h
@@ -28,10 +28,13 @@
#include "MeshLib/MeshSubset.h"
#include "MeshLib/MeshSubsets.h"
#include "MeshLib/NodeAdjacencyTable.h"
+#include "MeshLib/MeshSearcher.h"
#include "MeshGeoToolsLib/MeshNodeSearcher.h"
-#include "BoundaryCondition.h"
+#include "UniformDirichletBoundaryCondition.h"
#include "GroundwaterFlowFEM.h"
+#include "NeumannBcAssembler.h"
+#include "NeumannBc.h"
#include "ProcessVariable.h"
namespace ProcessLib
@@ -55,7 +58,7 @@ public:
// Find the corresponding process variable.
std::string const name = config.get<std::string>("process_variable");
- auto const& variable = std::find_if(variables.cbegin(), variables.cend(),
+ auto variable = std::find_if(variables.cbegin(), variables.cend(),
[&name](ProcessVariable const& v) {
return v.getName() == name;
});
@@ -66,7 +69,7 @@ public:
DBUG("Associate hydraulic_head with process variable \'%s\'.",
name.c_str());
- _hydraulic_head = &*variable;
+ _hydraulic_head = const_cast<ProcessVariable*>(&*variable);
}
@@ -130,15 +133,32 @@ public:
MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher(
_hydraulic_head->getMesh());
- using BCCI = ProcessVariable::BoundaryConditionCI;
- for (BCCI bc = _hydraulic_head->beginBoundaryConditions();
- bc != _hydraulic_head->endBoundaryConditions(); ++bc)
- {
- (*bc)->initialize(
+ _hydraulic_head->initializeDirichletBCs(
hydraulic_head_mesh_node_searcher,
_dirichlet_bc.global_ids, _dirichlet_bc.values);
+
+ //
+ // Neumann boundary conditions.
+ //
+ {
+ // Find mesh nodes.
+ MeshGeoToolsLib::BoundaryElementsSearcher hydraulic_head_mesh_element_searcher(
+ _hydraulic_head->getMesh(), hydraulic_head_mesh_node_searcher);
+
+ // Create a neumann BC for the hydraulic head storing them in the
+ // _neumann_bcs vector.
+ _hydraulic_head->createNeumannBcs(
+ std::back_inserter(_neumann_bcs),
+ hydraulic_head_mesh_element_searcher,
+ _global_setup,
+ _integration_order,
+ *_local_to_global_index_map,
+ *_mesh_subset_all_nodes);
}
+ for (auto bc : _neumann_bcs)
+ bc->initialize(_global_setup, *_A, *_rhs);
+
}
void solve()
@@ -152,6 +172,10 @@ public:
// Call global assembler for each local assembly item.
_global_setup.execute(*_global_assembler, _local_assemblers);
+ // Call global assembler for each Neumann boundary local assembler.
+ for (auto bc : _neumann_bcs)
+ bc->integrate(_global_setup);
+
// Apply known values from the Dirichlet boundary conditions.
MathLib::applyKnownSolution(*_A, *_rhs, _dirichlet_bc.global_ids, _dirichlet_bc.values);
@@ -179,9 +203,9 @@ public:
}
private:
- ProcessVariable const* _hydraulic_head = nullptr;
+ ProcessVariable* _hydraulic_head = nullptr;
- double const _hydraulic_conductivity = 1e-6;
+ double const _hydraulic_conductivity = 1;
MeshLib::MeshSubset const* _mesh_subset_all_nodes = nullptr;
std::vector<MeshLib::MeshSubsets*> _all_mesh_subsets;
@@ -214,6 +238,8 @@ private:
std::vector<double> values;
} _dirichlet_bc;
+ std::vector<NeumannBc<GlobalSetup>*> _neumann_bcs;
+
MeshLib::NodeAdjacencyTable _node_adjacency_table;
};
diff --git a/ProcessLib/NeumannBc.h b/ProcessLib/NeumannBc.h
new file mode 100644
index 0000000000000000000000000000000000000000..995ee16ab2e4107acd8554c9aa7eb275c2669e39
--- /dev/null
+++ b/ProcessLib/NeumannBc.h
@@ -0,0 +1,179 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, 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 PROCESS_LIB_NEUMANN_BC_H_
+#define PROCESS_LIB_NEUMANN_BC_H_
+
+#include <memory>
+#include <vector>
+
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+
+#include "AssemblerLib/VectorMatrixAssembler.h"
+#include "AssemblerLib/LocalToGlobalIndexMap.h"
+#include "AssemblerLib/LocalDataInitializer.h"
+#include "AssemblerLib/LocalAssemblerBuilder.h"
+#include "MeshLib/MeshSubset.h"
+
+#include "NeumannBcConfig.h"
+#include "NeumannBcAssembler.h"
+#include "MeshLib/MeshSearcher.h"
+
+namespace ProcessLib
+{
+
+/// The NeumannBc class is a process which is integrating a single, Neumann type
+/// boundary condition in to the global matrix and the right-hand-side.
+///
+/// The process operates on a set of elements and a subset of the DOF-table (the
+/// local to global index map). For each element a local assembler is created
+/// (NeumannBcAssembler).
+///
+/// In the construction phase the NeumannBcConfig together with global DOF-table
+/// and mesh subset are used.
+/// The creation of the local assemblers and binding to the global matrix and
+/// right-hand-sides happen in the initialize() function.
+/// The integration() function provides calls then the actual integration of the
+/// Neumann boundary condition.
+template <typename GlobalSetup_>
+class NeumannBc
+{
+public:
+ /// Create a Neumann boundary condition process from given config,
+ /// DOF-table, and a mesh subset.
+ /// A local DOF-table, a subset of the given one, is constructed.
+ NeumannBc(
+ NeumannBcConfig const& bc,
+ unsigned const integration_order,
+ AssemblerLib::LocalToGlobalIndexMap const& local_to_global_index_map,
+ MeshLib::MeshSubset const& mesh_subset_all_nodes
+ )
+ :
+ _function(*bc.getFunction()),
+ _integration_order(integration_order)
+ {
+ // deep copy because the neumann bc config destroys the elements.
+ std::transform(bc.elementsBegin(), bc.elementsEnd(),
+ std::back_inserter(_elements),
+ std::mem_fn(&MeshLib::Element::clone));
+
+ std::vector<MeshLib::Node*> nodes = MeshLib::selectNodes(_elements);
+
+ _mesh_subset_all_nodes =
+ mesh_subset_all_nodes.getIntersectionByNodes(nodes);
+ _all_mesh_subsets.push_back(new MeshLib::MeshSubsets(_mesh_subset_all_nodes));
+
+ _local_to_global_index_map.reset(
+ local_to_global_index_map.deriveBoundaryConstrainedMap(
+ _all_mesh_subsets, _elements));
+ }
+
+ ~NeumannBc()
+ {
+ for (auto e : _elements)
+ delete e;
+
+ for (auto p : _local_assemblers)
+ delete p;
+ }
+
+ /// Allocates the local assemblers for each element and stores references to
+ /// global matrix and the right-hand-side.
+ template <typename GlobalSetup>
+ void
+ initialize(GlobalSetup const& global_setup,
+ typename GlobalSetup::MatrixType& A,
+ typename GlobalSetup::VectorType& rhs)
+ {
+ // Shape matrices initializer
+ using LocalDataInitializer = AssemblerLib::LocalDataInitializer<
+ LocalNeumannBcAsmDataInterface,
+ LocalNeumannBcAsmData,
+ typename GlobalSetup::MatrixType,
+ typename GlobalSetup::VectorType>;
+
+ LocalDataInitializer initializer;
+
+ using LocalAssemblerBuilder =
+ AssemblerLib::LocalAssemblerBuilder<
+ MeshLib::Element,
+ LocalDataInitializer>;
+
+ // Populate the vector of local assemblers.
+ _local_assemblers.resize(_elements.size());
+ LocalAssemblerBuilder local_asm_builder(
+ initializer, *_local_to_global_index_map);
+
+ auto elementValueLookup = [this](MeshLib::Element const& e)
+ {
+ return _function();
+ };
+
+ DBUG("Calling local Neumann assembler builder for Neumann boundary elements.");
+ global_setup.execute(
+ local_asm_builder,
+ _elements,
+ _local_assemblers,
+ elementValueLookup,
+ _integration_order);
+
+ DBUG("Create global assembler.");
+ _global_assembler.reset(
+ new GlobalAssembler(A, rhs, *_local_to_global_index_map));
+ }
+
+ /// Calls local assemblers which calculate their contributions to the global
+ /// matrix and the right-hand-side.
+ template <typename GlobalSetup>
+ void
+ integrate(GlobalSetup const& global_setup)
+ {
+ global_setup.execute(*_global_assembler, _local_assemblers);
+ }
+
+
+private:
+ /// The right-hand-side function of the Neumann boundary condition given as
+ /// \f$ \alpha(x) \, \partial u(x) / \partial n = \text{_function}(x)\f$.
+ MathLib::ConstantFunction<double> const _function;
+
+ /// Vector of lower-dimensional elements on which the boundary condition is
+ /// defined.
+ std::vector<MeshLib::Element const*> _elements;
+
+ MeshLib::MeshSubset const* _mesh_subset_all_nodes = nullptr;
+ std::vector<MeshLib::MeshSubsets*> _all_mesh_subsets;
+
+ /// Local dof table, a subset of the global one restricted to the
+ /// participating #_elements of the boundary condition.
+ std::unique_ptr<AssemblerLib::LocalToGlobalIndexMap> _local_to_global_index_map;
+
+ /// Integration order for integration over the lower-dimensional elements of
+ /// the #_function.
+ unsigned const _integration_order;
+
+ using GlobalAssembler =
+ AssemblerLib::VectorMatrixAssembler<
+ typename GlobalSetup_::MatrixType,
+ typename GlobalSetup_::VectorType>;
+
+ std::unique_ptr<GlobalAssembler> _global_assembler;
+
+ using LocalAssembler = LocalNeumannBcAsmDataInterface<
+ typename GlobalSetup_::MatrixType, typename GlobalSetup_::VectorType>;
+
+ /// Local assemblers for each element of #_elements.
+ std::vector<LocalAssembler*> _local_assemblers;
+
+};
+
+} // namespace ProcessLib
+
+#endif // PROCESS_LIB_NEUMANN_BC_H_
diff --git a/ProcessLib/NeumannBcAssembler.h b/ProcessLib/NeumannBcAssembler.h
new file mode 100644
index 0000000000000000000000000000000000000000..4694ce3fab1570c2c2cca2015ed84ca8a77a3cef
--- /dev/null
+++ b/ProcessLib/NeumannBcAssembler.h
@@ -0,0 +1,122 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, 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 PROCESS_LIB_NEUMANN_BC_ASSEMBLER_H_
+#define PROCESS_LIB_NEUMANN_BC_ASSEMBLER_H_
+
+#include <memory>
+#include <vector>
+
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+
+namespace ProcessLib
+{
+
+template <typename GlobalMatrix, typename GlobalVector>
+class LocalNeumannBcAsmDataInterface
+{
+public:
+ virtual ~LocalNeumannBcAsmDataInterface() = default;
+
+ virtual void init(MeshLib::Element const& e,
+ std::size_t const local_matrix_size,
+ std::function<double (MeshLib::Element const&)> const& value_lookup,
+ unsigned const integration_order) = 0;
+
+ virtual void assemble() = 0;
+
+ virtual void addToGlobal(GlobalMatrix& A, GlobalVector& rhs,
+ AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const&) const = 0;
+};
+
+template <typename ShapeFunction_,
+ typename IntegrationMethod_,
+ typename GlobalMatrix,
+ typename GlobalVector>
+class LocalNeumannBcAsmData : public LocalNeumannBcAsmDataInterface<GlobalMatrix, GlobalVector>
+{
+public:
+ using ShapeFunction = ShapeFunction_;
+ using NodalMatrixType = typename ShapeMatrixPolicyType<ShapeFunction>::NodalMatrixType;
+ using NodalVectorType = typename ShapeMatrixPolicyType<ShapeFunction>::NodalVectorType;
+
+ using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction>;
+ using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+
+ /// The neumann_bc_value factor is directly integrated into the local
+ /// element matrix.
+ void
+ init(MeshLib::Element const& e,
+ std::size_t const local_matrix_size,
+ std::function<double (MeshLib::Element const&)> const& value_lookup,
+ unsigned const integration_order)
+ {
+ using FemType = NumLib::TemplateIsoparametric<
+ ShapeFunction, ShapeMatricesType>;
+
+ FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(&e));
+
+
+ _integration_order = integration_order;
+ IntegrationMethod_ integration_method(_integration_order);
+ std::size_t const n_integration_points = integration_method.getNPoints();
+
+ _shape_matrices.resize(n_integration_points);
+ for (std::size_t ip(0); ip < n_integration_points; ip++) {
+ _shape_matrices[ip].resize(ShapeFunction::DIM, ShapeFunction::NPOINTS);
+ fe.computeShapeFunctions(
+ integration_method.getWeightedPoint(ip).getCoords(),
+ _shape_matrices[ip]);
+ }
+
+ _neumann_bc_value = value_lookup(e);
+
+ _localA.reset(new NodalMatrixType(local_matrix_size, local_matrix_size));
+ _localRhs.reset(new NodalVectorType(local_matrix_size));
+ }
+
+ void
+ assemble()
+ {
+ _localA->setZero();
+ _localRhs->setZero();
+
+ IntegrationMethod_ integration_method(_integration_order);
+ std::size_t const n_integration_points = integration_method.getNPoints();
+
+ for (std::size_t ip(0); ip < n_integration_points; ip++) {
+ auto const& sm = _shape_matrices[ip];
+ auto const& wp = integration_method.getWeightedPoint(ip);
+ _localRhs->noalias() += sm.N * _neumann_bc_value
+ * sm.detJ * wp.getWeight();
+ }
+ }
+
+ void
+ addToGlobal(GlobalMatrix& A, GlobalVector& rhs,
+ AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const& indices) const
+ {
+ A.add(indices, *_localA);
+ rhs.add(indices.rows, *_localRhs);
+ }
+
+private:
+ std::vector<ShapeMatrices> _shape_matrices;
+ double _neumann_bc_value;
+
+ std::unique_ptr<NodalMatrixType> _localA;
+ std::unique_ptr<NodalVectorType> _localRhs;
+
+ unsigned _integration_order = 2;
+};
+
+} // namespace ProcessLib
+
+#endif // PROCESS_LIB_NEUMANN_BC_ASSEMBLER_H_
diff --git a/ProcessLib/NeumannBcConfig.h b/ProcessLib/NeumannBcConfig.h
new file mode 100644
index 0000000000000000000000000000000000000000..ac4ae385d4907b0d506800966a1d82a7c1518246
--- /dev/null
+++ b/ProcessLib/NeumannBcConfig.h
@@ -0,0 +1,112 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, 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 PROCESS_LIB_NEUMANN_BC_CONFIG_H_
+#define PROCESS_LIB_NEUMANN_BC_CONFIG_H_
+
+#include <boost/property_tree/ptree.hpp>
+#include "logog/include/logog.hpp"
+
+#include "MathLib/ConstantFunction.h"
+#include "MeshGeoToolsLib/BoundaryElementsSearcher.h"
+#include "MeshLib/Elements/Element.h"
+
+namespace GeoLib
+{
+ class GeoObject;
+}
+
+namespace ProcessLib
+{
+
+class BoundaryConditionConfig
+{
+public:
+ BoundaryConditionConfig(GeoLib::GeoObject const* const geometry)
+ : _geometry(geometry)
+ { }
+
+ virtual ~BoundaryConditionConfig() = default;
+
+protected:
+ GeoLib::GeoObject const* const _geometry;
+};
+
+
+/// Configuration of a Neumann type boundary condition read from input file.
+class NeumannBcConfig : public BoundaryConditionConfig
+{
+ using ConfigTree = boost::property_tree::ptree;
+public:
+ NeumannBcConfig(GeoLib::GeoObject const* const geometry,
+ ConfigTree const& config)
+ : BoundaryConditionConfig(geometry)
+ {
+ DBUG("Constructing NeumannBcConfig from config.");
+
+ double const value = config.get<double>("value", 0);
+ DBUG("Using value %g", value);
+
+ _function = new MathLib::ConstantFunction<double>(value);
+ }
+
+ ~NeumannBcConfig()
+ {
+ for (auto e : _elements)
+ delete e;
+
+ delete _function;
+ }
+
+ /// Initialize Neumann type boundary conditions.
+ /// Fills in elements of the particular geometry of the boundary condition.
+ /// The elements are appended to the \c elements vector.
+ void initialize(MeshGeoToolsLib::BoundaryElementsSearcher& searcher)
+ {
+ std::vector<MeshLib::Element*> elems = searcher.getBoundaryElements(*_geometry);
+
+ // deep copy because the searcher destroys the elements.
+ std::transform(elems.cbegin(), elems.cend(),
+ std::back_inserter(_elements),
+ std::mem_fn(&MeshLib::Element::clone));
+ }
+
+ /// Iterator over elements of the boundary condition.
+ std::vector<MeshLib::Element*>::const_iterator
+ elementsBegin() const
+ {
+ return _elements.cbegin();
+ }
+
+ /// Past the end iterator over elements of the boundary condition.
+ /// \sa elementsBegin().
+ std::vector<MeshLib::Element*>::const_iterator
+ elementsEnd() const
+ {
+ return _elements.cend();
+ }
+
+ /// Returns the right-hand-side function of the boundary condition.
+ MathLib::ConstantFunction<double>*
+ getFunction() const
+ {
+ return _function;
+ }
+
+private:
+ /// Domain of the boundary condition.
+ std::vector<MeshLib::Element*> _elements;
+
+ /// A function given on the domain of the boundary condition.
+ MathLib::ConstantFunction<double>* _function = nullptr;
+};
+
+} // namespace ProcessLib
+
+#endif // PROCESS_LIB_NEUMANN_BC_CONFIG_H_
diff --git a/ProcessLib/ProcessVariable.cpp b/ProcessLib/ProcessVariable.cpp
index 78d33361bcc870de7c1737944d4a35899627310a..c60f01ead86c9b2f79b30e2ba21cf54433783357 100644
--- a/ProcessLib/ProcessVariable.cpp
+++ b/ProcessLib/ProcessVariable.cpp
@@ -13,7 +13,7 @@
#include "GeoLib/GEOObjects.h"
#include "MeshLib/Mesh.h"
-#include "BoundaryCondition.h"
+#include "UniformDirichletBoundaryCondition.h"
#include "InitialCondition.h"
#include "ProcessVariable.h"
@@ -76,10 +76,15 @@ ProcessVariable::ProcessVariable(
if (type == "UniformDirichlet")
{
- _boundary_conditions.emplace_back(
+ _dirichlet_bcs.emplace_back(
new UniformDirichletBoundaryCondition(
geometry, bc_config));
}
+ else if (type == "UniformNeumann")
+ {
+ _neumann_bc_configs.emplace_back(
+ new NeumannBcConfig(geometry, bc_config));
+ }
else
{
ERR("Unknown type \'%s\' of the boundary condition.",
@@ -94,7 +99,10 @@ ProcessVariable::~ProcessVariable()
{
delete _initial_condition;
- for(auto p : _boundary_conditions)
+ for(auto p : _dirichlet_bcs)
+ delete p;
+
+ for(auto p : _neumann_bc_configs)
delete p;
}
@@ -108,4 +116,12 @@ MeshLib::Mesh const& ProcessVariable::getMesh() const
return _mesh;
}
+void ProcessVariable::initializeDirichletBCs(
+ MeshGeoToolsLib::MeshNodeSearcher& searcher,
+ std::vector<std::size_t>& global_ids, std::vector<double>& values)
+{
+ for (UniformDirichletBoundaryCondition* bc : _dirichlet_bcs)
+ bc->initialize(searcher, global_ids, values);
+}
+
} // namespace ProcessLib
diff --git a/ProcessLib/ProcessVariable.h b/ProcessLib/ProcessVariable.h
index 7ac4c81f6d96f9d1ed8a8dca49820897250cf8a3..e17a572dba3bb12adfa4538b730e5dd95075e72b 100644
--- a/ProcessLib/ProcessVariable.h
+++ b/ProcessLib/ProcessVariable.h
@@ -15,10 +15,16 @@
#include "GeoLib/GEOObjects.h"
#include "MeshLib/Mesh.h"
+#include "MeshGeoToolsLib/MeshNodeSearcher.h"
+#include "MeshGeoToolsLib/BoundaryElementsSearcher.h"
+
+#include "NeumannBcConfig.h"
+#include "NeumannBc.h"
+
namespace ProcessLib
{
- class BoundaryCondition;
class InitialCondition;
+ class UniformDirichletBoundaryCondition;
}
namespace ProcessLib
@@ -40,28 +46,28 @@ public:
/// Returns a mesh on which the process variable is defined.
MeshLib::Mesh const& getMesh() const;
- /// Const iterator over boundary conditions of the process variable.
- using BoundaryConditionCI = std::vector<BoundaryCondition*>::const_iterator;
-
- /// Returns a BoundaryConditionCI iterator to the beginning.
- BoundaryConditionCI
- beginBoundaryConditions() const
- {
- return _boundary_conditions.cbegin();
- }
+ void initializeDirichletBCs(MeshGeoToolsLib::MeshNodeSearcher& searcher,
+ std::vector<std::size_t>& global_ids, std::vector<double>& values);
- /// Returns a past-the-end BoundaryConditionCI iterator.
- BoundaryConditionCI
- endBoundaryConditions() const
+ template <typename OutputIterator, typename GlobalSetup, typename ...Args>
+ void createNeumannBcs(OutputIterator bcs,
+ MeshGeoToolsLib::BoundaryElementsSearcher& searcher,
+ GlobalSetup const& global_setup,
+ Args&&... args)
{
- return _boundary_conditions.cend();
+ for (NeumannBcConfig* config : _neumann_bc_configs)
+ {
+ config->initialize(searcher);
+ bcs = new NeumannBc<GlobalSetup>(*config, std::forward<Args>(args)...);
+ }
}
private:
std::string const _name;
MeshLib::Mesh const& _mesh;
InitialCondition* _initial_condition;
- std::vector<BoundaryCondition*> _boundary_conditions;
+ std::vector<UniformDirichletBoundaryCondition*> _dirichlet_bcs;
+ std::vector<NeumannBcConfig*> _neumann_bc_configs;
};
} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition.h b/ProcessLib/UniformDirichletBoundaryCondition.h
similarity index 76%
rename from ProcessLib/BoundaryCondition.h
rename to ProcessLib/UniformDirichletBoundaryCondition.h
index b460359e089f36252eff1463c7ab7e58d7b2ec43..2d0154fd733bee90f1659aabf2f2ab2869a51bd6 100644
--- a/ProcessLib/BoundaryCondition.h
+++ b/ProcessLib/UniformDirichletBoundaryCondition.h
@@ -26,38 +26,17 @@ namespace GeoLib
namespace ProcessLib
{
-class BoundaryCondition
-{
-public:
- BoundaryCondition(GeoLib::GeoObject const* const geometry)
- : _geometry(geometry)
- { }
-
- virtual ~BoundaryCondition() = default;
-
- /// Initialization interface of Dirichlet type boundary conditions.
- /// Fills in global_ids of the particular geometry of the boundary condition
- /// and the corresponding values.
- /// The ids are appended to the global_ids and also the values.
- virtual void initialize(
- MeshGeoToolsLib::MeshNodeSearcher& searcher,
- std::vector<std::size_t>& global_ids, std::vector<double>& values) = 0;
-
-protected:
- GeoLib::GeoObject const* const _geometry;
-};
-
/// The UniformDirichletBoundaryCondition class describes a constant in space
/// and time Dirichlet boundary condition.
/// The expected parameter in the passed configuration is "value" which, when
/// not present defaults to zero.
-class UniformDirichletBoundaryCondition : public BoundaryCondition
+class UniformDirichletBoundaryCondition
{
using ConfigTree = boost::property_tree::ptree;
public:
UniformDirichletBoundaryCondition(GeoLib::GeoObject const* const geometry,
ConfigTree const& config)
- : BoundaryCondition(geometry)
+ : _geometry(geometry)
{
DBUG("Constructing UniformDirichletBoundaryCondition from config.");
@@ -88,6 +67,7 @@ public:
private:
double _value;
+ GeoLib::GeoObject const* const _geometry;
};
diff --git a/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e0_neumann.prj.md5 b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e0_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..484a6f50d65cf95a2bb0363ed8a6c31422d561ac
--- /dev/null
+++ b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e0_neumann.prj.md5
@@ -0,0 +1 @@
+7b04a1ecb48658cd44cacf44e7c12a5f
diff --git a/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e1_neumann.prj.md5 b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e1_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..926108db4cbc6bbf1d9e46b42f9eb8ab42be95fd
--- /dev/null
+++ b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e1_neumann.prj.md5
@@ -0,0 +1 @@
+1d136ab1ffa5ea419e26978ca7d06e80
diff --git a/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e2_neumann.prj.md5 b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e2_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..eaf59bc98167cbfa68b87766243c112ca6e6352d
--- /dev/null
+++ b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e2_neumann.prj.md5
@@ -0,0 +1 @@
+8cae23bf96bc9189a479ab067c17beb0
diff --git a/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e3_neumann.prj.md5 b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e3_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..63e0d3493c4142a401774118c8252d6886d8df6d
--- /dev/null
+++ b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e3_neumann.prj.md5
@@ -0,0 +1 @@
+0ee43123514a4920ee069da96473121f
diff --git a/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e4_neumann.prj.md5 b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e4_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..32ff4708667cca087283176da95a5883b276fa33
--- /dev/null
+++ b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e4_neumann.prj.md5
@@ -0,0 +1 @@
+e6c2e2317e153268396e302939fd3c37
diff --git a/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e5_neumann.prj.md5 b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e5_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..daf6445ad1905a1954b3d3cbcc4eb630c1c4a4c5
--- /dev/null
+++ b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e5_neumann.prj.md5
@@ -0,0 +1 @@
+f3ff20ea4d49a9dd17ebbdf12e9ff544
diff --git a/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e6_neumann.prj.md5 b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e6_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..5bffdf58bae3d30b54f22c123d3bedaed33f2702
--- /dev/null
+++ b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_1e6_neumann.prj.md5
@@ -0,0 +1 @@
+e9a7983b85ab94a9cd906072f9b08c1f
diff --git a/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_2e4_neumann.prj.md5 b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_2e4_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..112267ed840dc0ce6fcd915c9076e90ecb691420
--- /dev/null
+++ b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_2e4_neumann.prj.md5
@@ -0,0 +1 @@
+e7dcf3f96b40cb2c3a6c9bb9e3eb67a9
diff --git a/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_3e4_neumann.prj.md5 b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_3e4_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..f5282e3934eab928b6be99b5f4011c38f4fef9e9
--- /dev/null
+++ b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_3e4_neumann.prj.md5
@@ -0,0 +1 @@
+54c43279079a2b333c9fe89e7668cd86
diff --git a/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_4e4_neumann.prj.md5 b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_4e4_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..1edf255c550ab8bd27062ba0c8fb0cf0055b220a
--- /dev/null
+++ b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_4e4_neumann.prj.md5
@@ -0,0 +1 @@
+847f7a7edb6669b21fc9f71580d4f654
diff --git a/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_5e4_neumann.prj.md5 b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_5e4_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..4461a8cfafca756dc2f0ea29219f0be1479bdf14
--- /dev/null
+++ b/Tests/Data/Elliptic/cube_1x1x1_GroundWaterFlow/cube_5e4_neumann.prj.md5
@@ -0,0 +1 @@
+180b76953459ae33e63ec051e26b1c01
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e0_neumann.prj.md5 b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e0_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..bef579a5eedff6228bbf7a151a19a2f38b94e2c0
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e0_neumann.prj.md5
@@ -0,0 +1 @@
+dbb942f5c176305be5b40a5b047e46de
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e1_neumann.prj.md5 b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e1_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..1b5f03d762e3e8c3a1f26245e7c3cb63f4287451
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e1_neumann.prj.md5
@@ -0,0 +1 @@
+ba884f3d14dd0d1fe47486eabf587112
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e2_neumann.prj.md5 b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e2_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..fb0870befcd33ead321a58edf9df451d3a35e28f
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e2_neumann.prj.md5
@@ -0,0 +1 @@
+dc9a86cc8b98420f0700bbd0844a1183
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e3_neumann.prj.md5 b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e3_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..4542619249cbd3f080647fb221e1af964805c2ad
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e3_neumann.prj.md5
@@ -0,0 +1 @@
+4324cd1de69be7c0311f2befc3797ddb
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e4_neumann.prj.md5 b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e4_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..e75b25c7953c93a876f5e8014bc9718a963864a0
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e4_neumann.prj.md5
@@ -0,0 +1 @@
+a2232351400427dab309c89fa4e038cf
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e5_neumann.prj.md5 b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e5_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..e2302e869b154a6a2814b051616238ca7b9873eb
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e5_neumann.prj.md5
@@ -0,0 +1 @@
+d6ae99cd82ac619c79f19d6d30fe078a
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e6_neumann.prj.md5 b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e6_neumann.prj.md5
new file mode 100644
index 0000000000000000000000000000000000000000..875c7f1c5766be2ff9b0f6c4a684eca6c08a457e
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e6_neumann.prj.md5
@@ -0,0 +1 @@
+7ce4190aaad81a77b1d2b0cbd7c1e216