diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
index 6e90534adf2333067594174a34f1c220c7082cf3..55edb6fda627ceb6b2ce6162d397d9c3d9c30840 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
@@ -16,6 +16,7 @@
#include "MeshGeoToolsLib/SearchLength.h"
#include "NeumannBoundaryCondition.h"
#include "NonuniformNeumannBoundaryCondition.h"
+#include "PressureBoundaryCondition.h"
#include "RobinBoundaryCondition.h"
namespace ProcessLib
@@ -56,6 +57,15 @@ BoundaryConditionBuilder::createBoundaryCondition(
shapefunction_order);
}
+ //
+ // Special boundary conditions
+ //
+ if (type == "Pressure")
+ {
+ return createPressureBoundaryCondition(config, dof_table, mesh,
+ variable_id, integration_order,
+ shapefunction_order, parameters);
+ }
OGS_FATAL("Unknown boundary condition type: `%s'.", type.c_str());
}
@@ -177,6 +187,31 @@ BoundaryConditionBuilder::createNonuniformNeumannBoundaryCondition(
integration_order, shapefunction_order, mesh);
}
+std::unique_ptr<BoundaryCondition>
+BoundaryConditionBuilder::createPressureBoundaryCondition(
+ 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::PressureBoundaryCondition::createPressureBoundaryCondition(
+ config.config,
+ getClonedElements(boundary_element_searcher, config.geometry),
+ dof_table, variable_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 c7dda2d155414cf94b9872092e3dd74f9a8a42c7..8c6a28a1b0b183185cf3abddf4938c200c557360 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.h
@@ -114,6 +114,14 @@ protected:
const MeshLib::Mesh& mesh, const int variable_id,
const unsigned integration_order, const unsigned shapefunction_order);
+ virtual std::unique_ptr<BoundaryCondition> createPressureBoundaryCondition(
+ 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/PressureBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/PressureBoundaryCondition-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..d1b62717f20bf44577dd95154fc266fa5fc329e4
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/PressureBoundaryCondition-impl.h
@@ -0,0 +1,115 @@
+/**
+ * \file
+ * \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 "MeshLib/MeshSearch/NodeSearch.h"
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+#include "PressureBoundaryConditionLocalAssembler.h"
+
+namespace ProcessLib
+{
+namespace PressureBoundaryCondition
+{
+template <template <typename, typename, unsigned>
+ class LocalAssemblerImplementation>
+PressureBoundaryCondition<LocalAssemblerImplementation>::
+ PressureBoundaryCondition(
+ bool const is_axially_symmetric, unsigned const integration_order,
+ unsigned const shapefunction_order,
+ NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
+ int const variable_id, unsigned const global_dim,
+ std::vector<MeshLib::Element*>&& elements,
+ Parameter<double> const& pressure)
+ : _elements(std::move(elements)),
+ _integration_order(integration_order),
+ _pressure(pressure)
+{
+ std::vector<MeshLib::Node*> const nodes =
+ MeshLib::getUniqueNodes(_elements);
+ DBUG("Found %d nodes for Natural BCs for the variable %d", nodes.size(),
+ variable_id);
+
+ // Assume that the mesh subsets are equal for all components of the
+ // variable.
+ auto const& mesh_subsets = dof_table_bulk.getMeshSubsets(variable_id, 0);
+
+ // 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 component ids vector for the current variable.
+ auto const& number_of_components =
+ dof_table_bulk.getNumberOfVariableComponents(variable_id);
+ std::vector<int> component_ids(number_of_components);
+ std::iota(std::begin(component_ids), std::end(component_ids), 0);
+
+ // 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_ids, std::move(all_mesh_subsets), _elements));
+
+ createLocalAssemblers<LocalAssemblerImplementation>(
+ global_dim, _elements, *_dof_table_boundary, shapefunction_order,
+ _local_assemblers, is_axially_symmetric, _integration_order, _pressure);
+}
+
+template <template <typename, typename, unsigned>
+ class LocalAssemblerImplementation>
+PressureBoundaryCondition<
+ LocalAssemblerImplementation>::~PressureBoundaryCondition()
+{
+ for (auto e : _elements)
+ delete e;
+}
+
+template <template <typename, typename, unsigned>
+ class LocalAssemblerImplementation>
+void PressureBoundaryCondition<LocalAssemblerImplementation>::applyNaturalBC(
+ const double t, const GlobalVector& x, GlobalMatrix& K, GlobalVector& b)
+{
+ GlobalExecutor::executeMemberOnDereferenced(
+ &PressureBoundaryConditionLocalAssemblerInterface::assemble,
+ _local_assemblers, *_dof_table_boundary, t, x, K, b);
+}
+
+std::unique_ptr<
+ PressureBoundaryCondition<PressureBoundaryConditionLocalAssembler>>
+createPressureBoundaryCondition(
+ BaseLib::ConfigTree const& config,
+ std::vector<MeshLib::Element*>&& elements,
+ NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_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 PressureBoundaryCondition from config.");
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
+ config.checkConfigParameter("type", "Pressure");
+
+ auto const parameter_name =
+ //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__Pressure__parameter}
+ config.getConfigParameter<std::string>("parameter");
+ DBUG("Using parameter %s", parameter_name.c_str());
+
+ auto const& pressure = findParameter<double>(parameter_name, parameters, 1);
+ return std::unique_ptr<
+ PressureBoundaryCondition<PressureBoundaryConditionLocalAssembler>>{
+ new PressureBoundaryCondition<PressureBoundaryConditionLocalAssembler>{
+ is_axially_symmetric, integration_order, shapefunction_order,
+ dof_table, variable_id, global_dim, std::move(elements), pressure}};
+}
+
+} // namespace PressureBoundaryCondition
+} // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/PressureBoundaryCondition.h b/ProcessLib/BoundaryCondition/PressureBoundaryCondition.h
new file mode 100644
index 0000000000000000000000000000000000000000..2feeb227b3424c5548a9814c1dec2d906bd3ad0e
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/PressureBoundaryCondition.h
@@ -0,0 +1,83 @@
+/**
+ * \file
+ * \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 "MeshLib/MeshSubset.h"
+#include "BoundaryCondition.h"
+#include "PressureBoundaryConditionLocalAssembler.h"
+
+namespace ProcessLib
+{
+namespace PressureBoundaryCondition
+{
+class PressureBoundaryConditionLocalAssemblerInterface;
+
+template <template <typename, typename, unsigned>
+ class LocalAssemblerImplementation>
+class PressureBoundaryCondition 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.
+ PressureBoundaryCondition(
+ bool const is_axially_symmetric, unsigned const integration_order,
+ unsigned const shapefunction_order,
+ NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
+ int const variable_id, unsigned const global_dim,
+ std::vector<MeshLib::Element*>&& elements,
+ Parameter<double> const& pressure);
+
+ ~PressureBoundaryCondition() 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:
+ /// 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<PressureBoundaryConditionLocalAssemblerInterface>>
+ _local_assemblers;
+
+ Parameter<double> const& _pressure;
+};
+
+std::unique_ptr<
+ PressureBoundaryCondition<PressureBoundaryConditionLocalAssembler>>
+createPressureBoundaryCondition(
+ BaseLib::ConfigTree const& config,
+ std::vector<MeshLib::Element*>&& elements,
+ NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_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);
+
+} // namespace PressureBoundaryCondition
+} // namespace ProcessLib
+
+#include "PressureBoundaryCondition-impl.h"
diff --git a/ProcessLib/BoundaryCondition/PressureBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/PressureBoundaryConditionLocalAssembler.h
new file mode 100644
index 0000000000000000000000000000000000000000..59cdb469afb697c2645133c901c6637e4909911a
--- /dev/null
+++ b/ProcessLib/BoundaryCondition/PressureBoundaryConditionLocalAssembler.h
@@ -0,0 +1,174 @@
+/**
+ * \file
+ * \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 "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "ProcessLib/Parameter/Parameter.h"
+
+#include "GenericNaturalBoundaryConditionLocalAssembler.h"
+
+namespace ProcessLib
+{
+namespace PressureBoundaryCondition
+{
+template <typename ShapeMatricesTypeDisplacement, int GlobalDim, int NPoints>
+struct IntegrationPointData final
+{
+ IntegrationPointData(
+ typename ShapeMatricesTypeDisplacement::template VectorType<
+ NPoints * GlobalDim> const& Nu_times_n_,
+ double const integration_weight_)
+ : Nu_times_n(Nu_times_n_), integration_weight(integration_weight_)
+ {
+ }
+
+ // Shape matrix (for displacement) times element's normal vector.
+ typename ShapeMatricesTypeDisplacement::template VectorType<
+ NPoints * GlobalDim> const Nu_times_n;
+
+ double const integration_weight;
+};
+
+class PressureBoundaryConditionLocalAssemblerInterface
+{
+public:
+ 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;
+ virtual ~PressureBoundaryConditionLocalAssemblerInterface() = default;
+};
+
+template <typename ShapeFunctionDisplacement, typename IntegrationMethod,
+ unsigned GlobalDim>
+class PressureBoundaryConditionLocalAssembler final
+ : public PressureBoundaryConditionLocalAssemblerInterface
+{
+public:
+ using ShapeMatricesTypeDisplacement =
+ ShapeMatrixPolicyType<ShapeFunctionDisplacement, GlobalDim>;
+ using GlobalDimVectorType =
+ typename ShapeMatrixPolicyType<ShapeFunctionDisplacement,
+ GlobalDim>::GlobalDimVectorType;
+
+ PressureBoundaryConditionLocalAssembler(MeshLib::Element const& e,
+ std::size_t const local_matrix_size,
+ bool is_axially_symmetric,
+ unsigned const integration_order,
+ Parameter<double> const& pressure)
+ : _integration_method(integration_order), _pressure(pressure)
+ {
+ _local_rhs.setZero(local_matrix_size);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ _ip_data.reserve(n_integration_points);
+
+ auto const shape_matrices_u =
+ initShapeMatrices<ShapeFunctionDisplacement,
+ ShapeMatricesTypeDisplacement, IntegrationMethod,
+ GlobalDim>(e, is_axially_symmetric,
+ _integration_method);
+
+ GlobalDimVectorType element_normal(GlobalDim);
+
+ // TODO Extend to rotated 2d meshes and line elements.
+ if (e.getGeomType() == MeshLib::MeshElemType::LINE)
+ {
+ auto v1 = (*e.getNode(1)) - (*e.getNode(0));
+ element_normal[0] = -v1[1];
+ element_normal[1] = v1[0];
+ element_normal.normalize();
+ }
+ else
+ {
+ auto const element_normal_vector =
+ MeshLib::FaceRule::getSurfaceNormal(&e).getNormalizedVector();
+
+ std::copy_n(element_normal_vector.getCoords(), GlobalDim,
+ element_normal.data());
+ }
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ typename ShapeMatricesTypeDisplacement::template MatrixType<
+ GlobalDim, displacement_size>
+ N_u = ShapeMatricesTypeDisplacement::template MatrixType<
+ GlobalDim, displacement_size>::Zero(GlobalDim,
+ displacement_size);
+ for (int i = 0; i < static_cast<int>(GlobalDim); ++i)
+ N_u.template block<1, displacement_size / GlobalDim>(
+ i, i * displacement_size / GlobalDim)
+ .noalias() = shape_matrices_u[ip].N;
+
+ double const integration_weight =
+ _integration_method.getWeightedPoint(ip).getWeight() *
+ shape_matrices_u[ip].integralMeasure *
+ shape_matrices_u[ip].detJ;
+
+ _ip_data.emplace_back(N_u.transpose() * element_normal,
+ integration_weight);
+ }
+ }
+
+ void assemble(std::size_t const id,
+ NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
+ double const t, const GlobalVector& /*x*/,
+ GlobalMatrix& /*K*/, GlobalVector& local_rhs)
+ {
+ _local_rhs.setZero();
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ SpatialPosition pos;
+ pos.setElementID(id);
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ pos.setIntegrationPoint(ip);
+
+ auto const& w = _ip_data[ip].integration_weight;
+ auto const& Nu_times_n = _ip_data[ip].Nu_times_n;
+
+ _local_rhs.noalias() -=
+ Nu_times_n.transpose() * _pressure(t, pos)[0] * w;
+ }
+
+ auto const indices = NumLib::getIndices(id, dof_table_boundary);
+ local_rhs.add(indices, _local_rhs);
+ }
+
+private:
+ IntegrationMethod const _integration_method;
+ Parameter<double> const& _pressure;
+
+ static const int displacement_size =
+ ShapeFunctionDisplacement::NPOINTS * GlobalDim;
+ std::vector<IntegrationPointData<ShapeMatricesTypeDisplacement, GlobalDim,
+ ShapeFunctionDisplacement::NPOINTS>,
+ Eigen::aligned_allocator<IntegrationPointData<
+ ShapeMatricesTypeDisplacement, GlobalDim,
+ ShapeFunctionDisplacement::NPOINTS>>>
+ _ip_data;
+
+ typename ShapeMatricesTypeDisplacement::template VectorType<
+ displacement_size>
+ _local_rhs;
+
+public:
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+} // namespace PressureBoundaryCondition
+} // namespace ProcessLib