From 2b84db3f4ba7f29af9fdba4a8ac225fb0882814e Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Fri, 28 Jul 2017 09:08:29 +0200
Subject: [PATCH] [PL] implemented nonuniform BC functionality
---
.../BoundaryCondition/BoundaryCondition.cpp | 22 +-----
.../BoundaryCondition/BoundaryCondition.h | 4 +-
...cNonuniformNaturalBoundaryCondition-impl.h | 76 ++++++++++---------
...enericNonuniformNaturalBoundaryCondition.h | 23 ++----
.../NonuniformNeumannBoundaryCondition.cpp | 65 +++++++++++++---
.../NonuniformNeumannBoundaryCondition.h | 11 +--
...rmNeumannBoundaryConditionLocalAssembler.h | 64 ++++++++++++----
7 files changed, 160 insertions(+), 105 deletions(-)
diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
index 3b0875555d4..4a164b8deca 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.cpp
@@ -53,7 +53,7 @@ BoundaryConditionBuilder::createBoundaryCondition(
{
return createNonuniformNeumannBoundaryCondition(
config, dof_table, mesh, variable_id, integration_order,
- shapefunction_order, parameters);
+ shapefunction_order);
}
OGS_FATAL("Unknown boundary condition type: `%s'.", type.c_str());
@@ -170,25 +170,11 @@ 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)
+ const unsigned shapefunction_order)
{
- 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);
+ config.config, dof_table, variable_id, config.component_id,
+ integration_order, shapefunction_order, mesh);
}
std::vector<MeshLib::Element*> BoundaryConditionBuilder::getClonedElements(
diff --git a/ProcessLib/BoundaryCondition/BoundaryCondition.h b/ProcessLib/BoundaryCondition/BoundaryCondition.h
index fce30ffe566..c7dda2d1554 100644
--- a/ProcessLib/BoundaryCondition/BoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/BoundaryCondition.h
@@ -112,9 +112,7 @@ protected:
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);
+ const unsigned integration_order, const unsigned shapefunction_order);
static std::vector<MeshLib::Element*> getClonedElements(
MeshGeoToolsLib::BoundaryElementsSearcher& boundary_element_searcher,
diff --git a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h
index 3a08a06be3b..a2673341111 100644
--- a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h
@@ -21,54 +21,58 @@ 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)
+ bool is_axially_symmetric, unsigned const integration_order,
+ unsigned const shapefunction_order, unsigned const global_dim,
+ std::unique_ptr<MeshLib::Mesh>&& boundary_mesh, Data&& data)
+ : _data(std::forward<Data>(data)), _boundary_mesh(std::move(boundary_mesh))
{
- assert(component_id <
- static_cast<int>(dof_table_bulk.getNumberOfComponents()));
+ static_assert(std::is_same<typename std::decay<BoundaryConditionData>::type,
+ typename std::decay<Data>::type>::value,
+ "Type mismatch between declared and passed BC data.");
- 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);
+#if 0
+ // TODO fix/improve check!
+ if (component_id >=
+ static_cast<int>(dof_table_bulk.getNumberOfComponents()))
+ {
+ OGS_FATAL(""); // TODO better error message.
+ }
+#endif
- auto const& mesh_subsets =
- dof_table_bulk.getMeshSubsets(variable_id, component_id);
+ if (_boundary_mesh->getDimension() + 1 != global_dim)
+ {
+ OGS_FATAL(
+ "The dimension of the given boundary mesh (%d) is not by one lower "
+ "than the bulk dimension (%d).",
+ _boundary_mesh->getDimension(), global_dim);
+ }
- // 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));
+ constructDofTable();
createLocalAssemblers<LocalAssemblerImplementation>(
- global_dim, _elements, *_dof_table_boundary, shapefunction_order,
- _local_assemblers, is_axially_symmetric, _integration_order, _data);
+ global_dim, _boundary_mesh->getElements(), *_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()
+void GenericNonuniformNaturalBoundaryCondition<
+ BoundaryConditionData, LocalAssemblerImplementation>::constructDofTable()
{
- for (auto e : _elements)
- delete e;
+ // construct one-component dof table for the surface mesh
+ _mesh_subset_all_nodes.reset(
+ new MeshLib::MeshSubset(*_boundary_mesh, &_boundary_mesh->getNodes()));
+
+ std::vector<MeshLib::MeshSubsets> all_mesh_subsets{
+ _mesh_subset_all_nodes.get()};
+
+ std::vector<unsigned> vec_var_n_components{1};
+
+ _dof_table_boundary = std::make_unique<NumLib::LocalToGlobalIndexMap>(
+ std::move(all_mesh_subsets), vec_var_n_components,
+ NumLib::ComponentOrder::BY_LOCATION);
}
template <typename BoundaryConditionData,
diff --git a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h
index 495b361adc9..be3a369e444 100644
--- a/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h
@@ -27,17 +27,9 @@ public:
/// 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;
+ bool is_axially_symmetric, unsigned const integration_order,
+ unsigned const shapefunction_order, unsigned const global_dim,
+ std::unique_ptr<MeshLib::Mesh>&& boundary_mesh, Data&& data);
/// Calls local assemblers which calculate their contributions to the global
/// matrix and the right-hand-side.
@@ -47,12 +39,12 @@ public:
GlobalVector& b) override;
private:
+ void constructDofTable();
+
/// 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::Mesh> _boundary_mesh;
std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_all_nodes;
@@ -60,9 +52,6 @@ private:
/// 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>>
diff --git a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp
index 7f50ca64070..64a584bc7f2 100644
--- a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp
@@ -8,6 +8,8 @@
*/
#include "NonuniformNeumannBoundaryCondition.h"
+
+#include "MeshLib/IO/readMeshFromFile.h"
#include "ProcessLib/Utils/ProcessUtils.h"
namespace ProcessLib
@@ -15,26 +17,67 @@ 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)
+ int const component_id, unsigned const integration_order,
+ unsigned const shapefunction_order, MeshLib::Mesh const& bulk_mesh)
{
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());
+ // TODO handle paths correctly
+ auto const mesh_file = config.getConfigParameter<std::string>("mesh");
+
+ std::unique_ptr<MeshLib::Mesh> surface_mesh(
+ MeshLib::IO::readMeshFromFile(mesh_file));
+
+ // Surface mesh and bulk mesh must have equal axial symmetry flags!
+ surface_mesh->setAxiallySymmetric(bulk_mesh.isAxiallySymmetric());
+
+ // TODO add field type?
+ auto const field_name =
+ config.getConfigParameter<std::string>("field_name");
+
+ auto const* const property =
+ surface_mesh->getProperties().getPropertyVector<double>(field_name);
+
+ if (!property)
+ {
+ OGS_FATAL("A property with name `%s' does not exist in `%s'.",
+ field_name.c_str(), mesh_file.c_str());
+ }
+
+ if (property->getMeshItemType() != MeshLib::MeshItemType::Node)
+ {
+ OGS_FATAL(
+ "Only nodal fields are supported fur non-uniform fields. Field "
+ "`%s' is not nodal.",
+ field_name.c_str());
+ }
+
+ if (property->getNumberOfComponents() != 1)
+ {
+ OGS_FATAL("`%s' is not a one-component field.", field_name.c_str());
+ }
+
+ auto const* const mapping_to_bulk_nodes =
+ surface_mesh->getProperties().getPropertyVector<long>(
+ "bulk_mesh_node_ids");
- auto const& param = findParameter<double>(param_name, parameters, 1);
+ if (!(mapping_to_bulk_nodes &&
+ mapping_to_bulk_nodes->getMeshItemType() ==
+ MeshLib::MeshItemType::Node) &&
+ mapping_to_bulk_nodes->getNumberOfComponents() == 1)
+ {
+ OGS_FATAL("Field `bulk_mesh_node_ids' is not set up properly.");
+ }
return std::make_unique<NonuniformNeumannBoundaryCondition>(
- is_axially_symmetric, integration_order, shapefunction_order, dof_table,
- variable_id, component_id, global_dim, std::move(elements), param);
+ bulk_mesh.isAxiallySymmetric(), integration_order, shapefunction_order,
+ bulk_mesh.getDimension(), std::move(surface_mesh),
+ NonuniformNeumannBoundaryConditionData{
+ *property, bulk_mesh.getID(), *mapping_to_bulk_nodes, dof_table,
+ variable_id, component_id});
}
} // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h
index c391a8dfc85..2073e7eb9f3 100644
--- a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h
@@ -10,24 +10,21 @@
#pragma once
#include "GenericNonuniformNaturalBoundaryCondition.h"
+#include "MeshLib/PropertyVector.h"
#include "NonuniformNeumannBoundaryConditionLocalAssembler.h"
-#include "ProcessLib/Parameter/Parameter.h"
namespace ProcessLib
{
using NonuniformNeumannBoundaryCondition =
GenericNonuniformNaturalBoundaryCondition<
- Parameter<double> const&,
+ NonuniformNeumannBoundaryConditionData,
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);
+ int const component_id, unsigned const integration_order,
+ unsigned const shapefunction_order, const MeshLib::Mesh& bulk_mesh);
} // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h
index 4d736f8a600..c53c3de6f03 100644
--- a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h
@@ -9,12 +9,27 @@
#pragma once
-#include "GenericNonuniformNaturalBoundaryConditionLocalAssembler.h"
+#include "MeshLib/PropertyVector.h"
#include "NumLib/DOF/DOFTableUtil.h"
-#include "ProcessLib/Parameter/Parameter.h"
+#include "NumLib/Function/Interpolation.h"
+
+#include "GenericNonuniformNaturalBoundaryConditionLocalAssembler.h"
namespace ProcessLib
{
+struct NonuniformNeumannBoundaryConditionData
+{
+ /* TODO use Parameter */
+ MeshLib::PropertyVector<double> const& values;
+
+ // Used for mapping boundary nodes to bulk nodes.
+ std::size_t bulk_mesh_id;
+ MeshLib::PropertyVector<long> const& mapping_to_bulk_nodes;
+ NumLib::LocalToGlobalIndexMap const& dof_table_bulk;
+ int const variable_id_bulk;
+ int const component_id_bulk;
+};
+
template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
class NonuniformNeumannBoundaryConditionLocalAssembler final
@@ -32,16 +47,16 @@ public:
std::size_t const local_matrix_size,
bool is_axially_symmetric,
unsigned const integration_order,
- Parameter<double> const& neumann_bc_parameter)
+ NonuniformNeumannBoundaryConditionData const& data)
: Base(e, is_axially_symmetric, integration_order),
- _neumann_bc_parameter(neumann_bc_parameter),
+ _data(data),
_local_rhs(local_matrix_size)
{
}
void assemble(std::size_t const id,
NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
- double const t, const GlobalVector& /*x*/,
+ double const /*t*/, const GlobalVector& /*x*/,
GlobalMatrix& /*K*/, GlobalVector& b) override
{
_local_rhs.setZero();
@@ -49,25 +64,48 @@ public:
unsigned const n_integration_points =
Base::_integration_method.getNumberOfPoints();
- SpatialPosition pos;
- pos.setElementID(id);
+ auto indices = NumLib::getIndices(id, dof_table_boundary);
+
+ // TODO lots of heap allocations.
+ std::vector<double> neumann_param_nodal_values_local;
+ neumann_param_nodal_values_local.reserve(indices.size());
+ for (auto i : indices)
+ {
+ neumann_param_nodal_values_local.push_back(
+ _data.values.getComponent(i, 0));
+ }
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
- pos.setIntegrationPoint(ip);
auto const& sm = Base::_shape_matrices[ip];
+ double flux;
+ NumLib::shapeFunctionInterpolate(neumann_param_nodal_values_local,
+ sm.N, flux);
+
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;
+ _local_rhs.noalias() +=
+ sm.N * flux * sm.detJ * wp.getWeight() * sm.integralMeasure;
}
- auto const indices = NumLib::getIndices(id, dof_table_boundary);
+ // map boundary dof indices to bulk dof indices
+ for (auto& i : indices)
+ {
+ auto const bulk_node_id =
+ _data.mapping_to_bulk_nodes.getComponent(i, 0);
+
+ MeshLib::Location const l{_data.bulk_mesh_id,
+ MeshLib::MeshItemType::Node,
+ static_cast<std::size_t>(bulk_node_id)};
+
+ i = _data.dof_table_bulk.getGlobalIndex(l, _data.variable_id_bulk,
+ _data.component_id_bulk);
+ assert(i != NumLib::MeshComponentMap::nop);
+ }
b.add(indices, _local_rhs);
}
private:
- Parameter<double> const& _neumann_bc_parameter;
+ NonuniformNeumannBoundaryConditionData const& _data;
typename Base::NodalVectorType _local_rhs;
public:
--
GitLab