diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
index 8df5b5e7bf7354726d678fb06b8969f5e76ebe87..3b0875555d438a3131721567f659fab023dcd51f 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
@@ -8,13 +8,14 @@
*/
#include "BoundaryCondition.h"
+#include "BoundaryConditionConfig.h"
+#include "DirichletBoundaryCondition.h"
#include "MeshGeoToolsLib/BoundaryElementsSearcher.h"
+#include "MeshGeoToolsLib/CreateSearchLength.h"
#include "MeshGeoToolsLib/MeshNodeSearcher.h"
#include "MeshGeoToolsLib/SearchLength.h"
-#include "MeshGeoToolsLib/CreateSearchLength.h"
-#include "BoundaryConditionConfig.h"
-#include "DirichletBoundaryCondition.h"
#include "NeumannBoundaryCondition.h"
+#include "NonuniformNeumannBoundaryCondition.h"
#include "RobinBoundaryCondition.h"
namespace ProcessLib
@@ -48,6 +49,12 @@ BoundaryConditionBuilder::createBoundaryCondition(
config, dof_table, mesh, variable_id,
integration_order, shapefunction_order, parameters);
}
+ if (type == "NonuniformNeumann")
+ {
+ return createNonuniformNeumannBoundaryCondition(
+ config, dof_table, mesh, variable_id, integration_order,
+ shapefunction_order, parameters);
+ }
OGS_FATAL("Unknown boundary condition type: `%s'.", type.c_str());
}
@@ -158,6 +165,32 @@ BoundaryConditionBuilder::createRobinBoundaryCondition(
parameters);
}
+std::unique_ptr<BoundaryCondition>
+BoundaryConditionBuilder::createNonuniformNeumannBoundaryCondition(
+ const BoundaryConditionConfig& config,
+ const NumLib::LocalToGlobalIndexMap& dof_table, const MeshLib::Mesh& mesh,
+ const int variable_id, const unsigned integration_order,
+ const unsigned shapefunction_order,
+ const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>& parameters)
+{
+ std::unique_ptr<MeshGeoToolsLib::SearchLength> search_length_algorithm =
+ MeshGeoToolsLib::createSearchLengthAlgorithm(config.config, mesh);
+
+ MeshGeoToolsLib::MeshNodeSearcher const& mesh_node_searcher =
+ MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher(
+ mesh, std::move(search_length_algorithm));
+
+ MeshGeoToolsLib::BoundaryElementsSearcher boundary_element_searcher(
+ mesh, mesh_node_searcher);
+
+ return ProcessLib::createNonuniformNeumannBoundaryCondition(
+ config.config,
+ getClonedElements(boundary_element_searcher, config.geometry),
+ dof_table, variable_id, config.component_id, mesh.isAxiallySymmetric(),
+ integration_order, shapefunction_order, mesh.getDimension(),
+ parameters);
+}
+
std::vector<MeshLib::Element*> BoundaryConditionBuilder::getClonedElements(
MeshGeoToolsLib::BoundaryElementsSearcher& boundary_element_searcher,
GeoLib::GeoObject const& geometry)
diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.h b/ProcessLib/BoundaryCondition/BoundaryCondition.h
index e72923933438ca38cd5bdb42769f3e3051c260cc..fce30ffe566cda1aba4ca838ce90491d137d9da9 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.h
@@ -107,6 +107,15 @@ protected:
const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>&
parameters);
+ virtual std::unique_ptr<BoundaryCondition>
+ createNonuniformNeumannBoundaryCondition(
+ const BoundaryConditionConfig& config,
+ const NumLib::LocalToGlobalIndexMap& dof_table,
+ const MeshLib::Mesh& mesh, const int variable_id,
+ const unsigned integration_order, const unsigned shapefunction_order,
+ const std::vector<std::unique_ptr<ProcessLib::ParameterBase>>&
+ parameters);
+
static std::vector<MeshLib::Element*> getClonedElements(
MeshGeoToolsLib::BoundaryElementsSearcher& boundary_element_searcher,
GeoLib::GeoObject const& geometry);
diff --git a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..3a08a06be3b0ec07d4fa61bd7a8d0f07b8acdaec
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h
@@ -0,0 +1,90 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "GenericNonuniformNaturalBoundaryCondition.h"
+#include "GenericNonuniformNaturalBoundaryConditionLocalAssembler.h"
+#include "MeshLib/MeshSearch/NodeSearch.h"
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+
+namespace ProcessLib
+{
+template <typename BoundaryConditionData,
+ template <typename, typename, unsigned>
+ class LocalAssemblerImplementation>
+template <typename Data>
+GenericNonuniformNaturalBoundaryCondition<BoundaryConditionData,
+ LocalAssemblerImplementation>::
+ GenericNonuniformNaturalBoundaryCondition(
+ typename std::enable_if<
+ std::is_same<typename std::decay<BoundaryConditionData>::type,
+ typename std::decay<Data>::type>::value,
+ bool>::type is_axially_symmetric,
+ unsigned const integration_order, unsigned const shapefunction_order,
+ NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
+ int const variable_id, int const component_id,
+ unsigned const global_dim, std::vector<MeshLib::Element*>&& elements,
+ Data&& data)
+ : _data(std::forward<Data>(data)),
+ _elements(std::move(elements)),
+ _integration_order(integration_order)
+{
+ assert(component_id <
+ static_cast<int>(dof_table_bulk.getNumberOfComponents()));
+
+ std::vector<MeshLib::Node*> nodes = MeshLib::getUniqueNodes(_elements);
+ DBUG("Found %d nodes for Natural BCs for the variable %d and component %d",
+ nodes.size(), variable_id, component_id);
+
+ auto const& mesh_subsets =
+ dof_table_bulk.getMeshSubsets(variable_id, component_id);
+
+ // TODO extend the node intersection to all parts of mesh_subsets, i.e.
+ // to each of the MeshSubset in the mesh_subsets.
+ _mesh_subset_all_nodes.reset(
+ mesh_subsets.getMeshSubset(0).getIntersectionByNodes(nodes));
+ MeshLib::MeshSubsets all_mesh_subsets{_mesh_subset_all_nodes.get()};
+
+ // Create local DOF table from intersected mesh subsets for the given
+ // variable and component ids.
+ _dof_table_boundary.reset(dof_table_bulk.deriveBoundaryConstrainedMap(
+ variable_id, {component_id}, std::move(all_mesh_subsets), _elements));
+
+ createLocalAssemblers<LocalAssemblerImplementation>(
+ global_dim, _elements, *_dof_table_boundary, shapefunction_order,
+ _local_assemblers, is_axially_symmetric, _integration_order, _data);
+}
+
+template <typename BoundaryConditionData,
+ template <typename, typename, unsigned>
+ class LocalAssemblerImplementation>
+GenericNonuniformNaturalBoundaryCondition<
+ BoundaryConditionData,
+ LocalAssemblerImplementation>::~GenericNonuniformNaturalBoundaryCondition()
+{
+ for (auto e : _elements)
+ delete e;
+}
+
+template <typename BoundaryConditionData,
+ template <typename, typename, unsigned>
+ class LocalAssemblerImplementation>
+void GenericNonuniformNaturalBoundaryCondition<
+ BoundaryConditionData,
+ LocalAssemblerImplementation>::applyNaturalBC(const double t,
+ const GlobalVector& x,
+ GlobalMatrix& K,
+ GlobalVector& b)
+{
+ GlobalExecutor::executeMemberOnDereferenced(
+ &GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface::
+ assemble,
+ _local_assemblers, *_dof_table_boundary, t, x, K, b);
+}
+
+} // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h
new file mode 100644
index 0000000000000000000000000000000000000000..495b361adc9ba69e7a5610a97a4d9af8864fb239
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h
@@ -0,0 +1,74 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "BoundaryCondition.h"
+#include "MeshLib/MeshSubset.h"
+
+namespace ProcessLib
+{
+class GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface;
+
+template <typename BoundaryConditionData,
+ template <typename, typename, unsigned>
+ class LocalAssemblerImplementation>
+class GenericNonuniformNaturalBoundaryCondition final : public BoundaryCondition
+{
+public:
+ /// Create a boundary condition process from given config,
+ /// DOF-table, and a mesh subset for a given variable and its component.
+ /// A local DOF-table, a subset of the given one, is constructed.
+ template <typename Data>
+ GenericNonuniformNaturalBoundaryCondition(
+ typename std::enable_if<
+ std::is_same<typename std::decay<BoundaryConditionData>::type,
+ typename std::decay<Data>::type>::value,
+ bool>::type is_axially_symmetric,
+ unsigned const integration_order, unsigned const shapefunction_order,
+ NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
+ int const variable_id, int const component_id,
+ unsigned const global_dim, std::vector<MeshLib::Element*>&& elements,
+ Data&& data);
+
+ ~GenericNonuniformNaturalBoundaryCondition() override;
+
+ /// Calls local assemblers which calculate their contributions to the global
+ /// matrix and the right-hand-side.
+ void applyNaturalBC(const double t,
+ GlobalVector const& x,
+ GlobalMatrix& K,
+ GlobalVector& b) override;
+
+private:
+ /// Data used in the assembly of the specific boundary condition.
+ BoundaryConditionData _data;
+
+ /// Vector of lower-dimensional elements on which the boundary condition is
+ /// defined.
+ std::vector<MeshLib::Element*> _elements;
+
+ std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_all_nodes;
+
+ /// Local dof table, a subset of the global one restricted to the
+ /// participating #_elements of the boundary condition.
+ std::unique_ptr<NumLib::LocalToGlobalIndexMap> _dof_table_boundary;
+
+ /// Integration order for integration over the lower-dimensional elements
+ unsigned const _integration_order;
+
+ /// Local assemblers for each element of #_elements.
+ std::vector<std::unique_ptr<
+ GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface>>
+ _local_assemblers;
+};
+
+} // ProcessLib
+
+#include "GenericNonuniformNaturalBoundaryCondition-impl.h"
diff --git a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryConditionLocalAssembler.h
new file mode 100644
index 0000000000000000000000000000000000000000..73eb6f9c14c8a49e5a1648bc0247d69ae8839499
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryConditionLocalAssembler.h
@@ -0,0 +1,58 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+namespace ProcessLib
+{
+class GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface
+{
+public:
+ virtual ~GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface() =
+ default;
+
+ virtual void assemble(
+ std::size_t const id,
+ NumLib::LocalToGlobalIndexMap const& dof_table_boundary, double const t,
+ const GlobalVector& x, GlobalMatrix& K, GlobalVector& b) = 0;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+class GenericNonuniformNaturalBoundaryConditionLocalAssembler
+ : public GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface
+{
+protected:
+ using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+ using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
+ using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
+
+public:
+ GenericNonuniformNaturalBoundaryConditionLocalAssembler(
+ MeshLib::Element const& e, bool is_axially_symmetric,
+ unsigned const integration_order)
+ : _integration_method(integration_order),
+ _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
+ IntegrationMethod, GlobalDim>(
+ e, is_axially_symmetric, _integration_method))
+ {
+ }
+
+protected:
+ IntegrationMethod const _integration_method;
+ std::vector<typename ShapeMatricesType::ShapeMatrices,
+ Eigen::aligned_allocator<
+ typename ShapeMatricesType::ShapeMatrices>> const
+ _shape_matrices;
+};
+
+} // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7f50ca64070b9dca929910f7f1a82b2c18cbbebf
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp
@@ -0,0 +1,40 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "NonuniformNeumannBoundaryCondition.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+namespace ProcessLib
+{
+std::unique_ptr<NonuniformNeumannBoundaryCondition>
+createNonuniformNeumannBoundaryCondition(
+ BaseLib::ConfigTree const& config,
+ std::vector<MeshLib::Element*>&& elements,
+ NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
+ int const component_id, bool is_axially_symmetric,
+ unsigned const integration_order, unsigned const shapefunction_order,
+ unsigned const global_dim,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters)
+{
+ DBUG("Constructing NonuniformNeumann BC from config.");
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
+ config.checkConfigParameter("type", "NonuniformNeumann");
+
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__NonuniformNeumann__parameter}
+ auto const param_name = config.getConfigParameter<std::string>("parameter");
+ DBUG("Using parameter %s", param_name.c_str());
+
+ auto const& param = findParameter<double>(param_name, parameters, 1);
+
+ return std::make_unique<NonuniformNeumannBoundaryCondition>(
+ is_axially_symmetric, integration_order, shapefunction_order, dof_table,
+ variable_id, component_id, global_dim, std::move(elements), param);
+}
+
+} // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h
new file mode 100644
index 0000000000000000000000000000000000000000..c391a8dfc859270ce25d1e74e41e3ab8c7f20786
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h
@@ -0,0 +1,33 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "GenericNonuniformNaturalBoundaryCondition.h"
+#include "NonuniformNeumannBoundaryConditionLocalAssembler.h"
+#include "ProcessLib/Parameter/Parameter.h"
+
+namespace ProcessLib
+{
+using NonuniformNeumannBoundaryCondition =
+ GenericNonuniformNaturalBoundaryCondition<
+ Parameter<double> const&,
+ NonuniformNeumannBoundaryConditionLocalAssembler>;
+
+std::unique_ptr<NonuniformNeumannBoundaryCondition>
+createNonuniformNeumannBoundaryCondition(
+ BaseLib::ConfigTree const& config,
+ std::vector<MeshLib::Element*>&& elements,
+ NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
+ int const component_id, bool is_axially_symmetric,
+ unsigned const integration_order, unsigned const shapefunction_order,
+ unsigned const global_dim,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters);
+
+} // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h
new file mode 100644
index 0000000000000000000000000000000000000000..4d736f8a60011ef746ba1af7985ca972b4baff87
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h
@@ -0,0 +1,77 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "GenericNonuniformNaturalBoundaryConditionLocalAssembler.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "ProcessLib/Parameter/Parameter.h"
+
+namespace ProcessLib
+{
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+class NonuniformNeumannBoundaryConditionLocalAssembler final
+ : public GenericNonuniformNaturalBoundaryConditionLocalAssembler<
+ ShapeFunction, IntegrationMethod, GlobalDim>
+{
+ using Base = GenericNonuniformNaturalBoundaryConditionLocalAssembler<
+ ShapeFunction, IntegrationMethod, GlobalDim>;
+
+public:
+ /// The neumann_bc_value factor is directly integrated into the local
+ /// element matrix.
+ NonuniformNeumannBoundaryConditionLocalAssembler(
+ MeshLib::Element const& e,
+ std::size_t const local_matrix_size,
+ bool is_axially_symmetric,
+ unsigned const integration_order,
+ Parameter<double> const& neumann_bc_parameter)
+ : Base(e, is_axially_symmetric, integration_order),
+ _neumann_bc_parameter(neumann_bc_parameter),
+ _local_rhs(local_matrix_size)
+ {
+ }
+
+ void assemble(std::size_t const id,
+ NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
+ double const t, const GlobalVector& /*x*/,
+ GlobalMatrix& /*K*/, GlobalVector& b) override
+ {
+ _local_rhs.setZero();
+
+ unsigned const n_integration_points =
+ Base::_integration_method.getNumberOfPoints();
+
+ SpatialPosition pos;
+ pos.setElementID(id);
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ pos.setIntegrationPoint(ip);
+ auto const& sm = Base::_shape_matrices[ip];
+ auto const& wp = Base::_integration_method.getWeightedPoint(ip);
+ _local_rhs.noalias() += sm.N * _neumann_bc_parameter(t, pos)[0] *
+ sm.detJ * wp.getWeight() *
+ sm.integralMeasure;
+ }
+
+ auto const indices = NumLib::getIndices(id, dof_table_boundary);
+ b.add(indices, _local_rhs);
+ }
+
+private:
+ Parameter<double> const& _neumann_bc_parameter;
+ typename Base::NodalVectorType _local_rhs;
+
+public:
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+} // namespace ProcessLib