From 25433246e123ddb2b223788ed12eb4b8723e9819 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <github@naumov.de>
Date: Thu, 21 Feb 2019 00:31:55 +0100
Subject: [PATCH] [PL] Replace non-uniform BCs.
---
.../CreateBoundaryCondition.cpp | 16 ---
.../DirichletBoundaryCondition.cpp | 13 +-
.../NeumannBoundaryConditionLocalAssembler.h | 10 +-
.../NonuniformDirichletBoundaryCondition.cpp | 63 ---------
.../NonuniformDirichletBoundaryCondition.h | 129 ------------------
.../NonuniformNeumannBoundaryCondition.cpp | 92 -------------
.../NonuniformNeumannBoundaryCondition.h | 30 ----
ProcessLib/Parameter/MeshNodeParameter.cpp | 2 +-
ProcessLib/Parameter/MeshNodeParameter.h | 4 +-
9 files changed, 22 insertions(+), 337 deletions(-)
delete mode 100644 ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.cpp
delete mode 100644 ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h
delete mode 100644 ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp
delete mode 100644 ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h
diff --git a/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp
index e3518dd4c0b..c4f2b4cb3df 100644
--- a/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/CreateBoundaryCondition.cpp
@@ -15,8 +15,6 @@
#include "DirichletBoundaryCondition.h"
#include "DirichletBoundaryConditionWithinTimeInterval.h"
#include "NeumannBoundaryCondition.h"
-#include "NonuniformDirichletBoundaryCondition.h"
-#include "NonuniformNeumannBoundaryCondition.h"
#include "NonuniformVariableDependentNeumannBoundaryCondition.h"
#include "NormalTractionBoundaryCondition.h"
#include "PhaseFieldIrreversibleDamageOracleBoundaryCondition.h"
@@ -78,20 +76,6 @@ std::unique_ptr<BoundaryCondition> createBoundaryCondition(
*config.component_id, integration_order, shapefunction_order,
bulk_mesh.getDimension(), parameters);
}
- if (type == "NonuniformDirichlet")
- {
- return ProcessLib::createNonuniformDirichletBoundaryCondition(
- config.config, config.boundary_mesh, dof_table, variable_id,
- *config.component_id);
- }
- if (type == "NonuniformNeumann")
- {
- return ProcessLib::createNonuniformNeumannBoundaryCondition(
- config.config, config.boundary_mesh, dof_table, variable_id,
- *config.component_id, integration_order, shapefunction_order,
- bulk_mesh);
- }
-
if (type == "NonuniformVariableDependentNeumann")
{
return ProcessLib::
diff --git a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp
index fc757a7b5c0..c714a531dec 100644
--- a/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/DirichletBoundaryCondition.cpp
@@ -67,7 +67,16 @@ std::unique_ptr<DirichletBoundaryCondition> createDirichletBoundaryCondition(
auto const param_name = config.getConfigParameter<std::string>("parameter");
DBUG("Using parameter %s", param_name.c_str());
- auto& param = findParameter<double>(param_name, parameters, 1);
+ auto& parameter = findParameter<double>(param_name, parameters, 1);
+
+ if (parameter.mesh() && *parameter.mesh() != bc_mesh)
+ {
+ OGS_FATAL(
+ "The boundary condition is defined on a different domain than the "
+ "parameter: boundary condition is defined on mesh '%s', parameter "
+ "is defined on mesh '%s'.",
+ bc_mesh.getName().c_str(), parameter.mesh()->getName().c_str());
+ }
// In case of partitioned mesh the boundary could be empty, i.e. there is no
// boundary condition.
@@ -84,7 +93,7 @@ std::unique_ptr<DirichletBoundaryCondition> createDirichletBoundaryCondition(
#endif // USE_PETSC
return std::make_unique<DirichletBoundaryCondition>(
- param, bc_mesh, dof_table_bulk, variable_id, component_id);
+ parameter, bc_mesh, dof_table_bulk, variable_id, component_id);
}
} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
index 52931373404..20907cf23c4 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
@@ -26,6 +26,7 @@ class NeumannBoundaryConditionLocalAssembler final
using Base = GenericNaturalBoundaryConditionLocalAssembler<
ShapeFunction, IntegrationMethod, GlobalDim>;
using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+ using NodalVectorType = typename Base::NodalVectorType;
public:
/// The neumann_bc_value factor is directly integrated into the local
@@ -61,6 +62,11 @@ public:
FemType fe(
static_cast<const typename ShapeFunction::MeshElement&>(Base::_element));
+ // Get element nodes for the interpolation from nodes to
+ // integration point.
+ NodalVectorType parameter_node_values =
+ _neumann_bc_parameter.getNodalValuesOnElement(Base::_element, t);
+
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
pos.setIntegrationPoint(ip);
@@ -70,7 +76,7 @@ public:
MathLib::TemplatePoint<double, 3> coords_ip(fe.interpolateCoordinates(sm.N));
pos.setCoordinates(coords_ip);
- _local_rhs.noalias() += sm.N * _neumann_bc_parameter(t, pos)[0] *
+ _local_rhs.noalias() += sm.N * parameter_node_values.dot(sm.N) *
sm.detJ * wp.getWeight() *
sm.integralMeasure;
}
@@ -81,7 +87,7 @@ public:
private:
Parameter<double> const& _neumann_bc_parameter;
- typename Base::NodalVectorType _local_rhs;
+ NodalVectorType _local_rhs;
public:
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
diff --git a/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.cpp
deleted file mode 100644
index 5786e0947c6..00000000000
--- a/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.cpp
+++ /dev/null
@@ -1,63 +0,0 @@
-/**
- * \file
- *
- * \copyright
- * Copyright (c) 2012-2019, 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 "NonuniformDirichletBoundaryCondition.h"
-
-#include "MeshLib/IO/readMeshFromFile.h"
-#include "ProcessLib/Utils/ProcessUtils.h"
-
-namespace ProcessLib
-{
-std::unique_ptr<NonuniformDirichletBoundaryCondition>
-createNonuniformDirichletBoundaryCondition(
- BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh,
- NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
- int const component_id)
-{
- DBUG("Constructing NonuniformDirichlet BC from config.");
- //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
- config.checkConfigParameter("type", "NonuniformDirichlet");
-
- // TODO finally use ProcessLib::Parameter here
- auto const field_name =
- //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__NonuniformDirichlet__field_name}
- config.getConfigParameter<std::string>("field_name");
-
- auto const* const property =
- boundary_mesh.getProperties().getPropertyVector<double>(
- field_name, MeshLib::MeshItemType::Node, 1);
-
- if (!property)
- {
- OGS_FATAL("A property with name `%s' does not exist in `%s'.",
- field_name.c_str(), boundary_mesh.getName().c_str());
- }
-
- // 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.
- // TODO (naumov): Add a function to ConfigTree for skipping the tags of the
- // subtree and move the code up in createBoundaryCondition().
- if (boundary_mesh.getDimension() == 0 &&
- boundary_mesh.getNumberOfNodes() == 0 &&
- boundary_mesh.getNumberOfElements() == 0)
- {
- return nullptr;
- }
-#endif // USE_PETSC
-
- return std::make_unique<NonuniformDirichletBoundaryCondition>(
- boundary_mesh, *property, dof_table, variable_id, component_id);
-}
-
-} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h
deleted file mode 100644
index ba67cb886f3..00000000000
--- a/ProcessLib/BoundaryCondition/NonuniformDirichletBoundaryCondition.h
+++ /dev/null
@@ -1,129 +0,0 @@
-/**
- * \file
- *
- * \copyright
- * Copyright (c) 2012-2019, 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/PropertyVector.h"
-#include "NumLib/DOF/DOFTableUtil.h"
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "NumLib/IndexValueVector.h"
-#include "ProcessLib/Parameter/Parameter.h"
-#include "ProcessLib/Utils/ProcessUtils.h"
-
-namespace ProcessLib
-{
-class NonuniformDirichletBoundaryCondition final : public BoundaryCondition
-{
-public:
- NonuniformDirichletBoundaryCondition(
- // int const bulk_mesh_dimension,
- MeshLib::Mesh const& boundary_mesh,
- MeshLib::PropertyVector<double> const& values,
- NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
- int const variable_id,
- int const component_id)
- : _values(values),
- _boundary_mesh(boundary_mesh),
- _variable_id(variable_id),
- _component_id(component_id)
- {
- if (_variable_id >=
- static_cast<int>(dof_table_bulk.getNumberOfVariables()) ||
- _component_id >=
- dof_table_bulk.getNumberOfVariableComponents(_variable_id))
- {
- OGS_FATAL(
- "Variable id or component id too high. Actual values: (%d, "
- "%d), "
- "maximum values: (%d, %d).",
- _variable_id, _component_id,
- dof_table_bulk.getNumberOfVariables(),
- dof_table_bulk.getNumberOfVariableComponents(_variable_id));
- }
-
- if (!_boundary_mesh.getProperties().existsPropertyVector<std::size_t>(
- "bulk_node_ids"))
- {
- OGS_FATAL(
- "The required bulk node ids map does not exist in the boundary "
- "mesh '%s'.",
- _boundary_mesh.getName().c_str());
- }
-
- std::vector<MeshLib::Node*> const& bc_nodes = _boundary_mesh.getNodes();
- DBUG(
- "Found %d nodes for Natural BCs for the variable %d and component "
- "%d",
- bc_nodes.size(), variable_id, component_id);
-
- MeshLib::MeshSubset boundary_mesh_subset(_boundary_mesh, bc_nodes);
-
- // Create local DOF table from the BC mesh subset for the given variable
- // and component id.
- _dof_table_boundary.reset(dof_table_bulk.deriveBoundaryConstrainedMap(
- variable_id, {component_id}, std::move(boundary_mesh_subset)));
- }
-
- void getEssentialBCValues(
- const double /*t*/, GlobalVector const& /*x*/,
- NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override
- {
- bc_values.ids.clear();
- bc_values.values.clear();
-
- // Convert mesh node ids to global index for the given component.
- bc_values.ids.reserve(_values.size());
- bc_values.values.reserve(_values.size());
-
- // Map boundary dof indices to bulk dof indices and the corresponding
- // values.
- for (auto const* const node : _boundary_mesh.getNodes())
- {
- auto const node_id = node->getID();
- auto const global_index = _dof_table_boundary->getGlobalIndex(
- {_boundary_mesh.getID(), MeshLib::MeshItemType::Node, node_id},
- _variable_id, _component_id);
- if (global_index == NumLib::MeshComponentMap::nop)
- {
- continue;
- }
- // For the DDC approach (e.g. with PETSc option), the negative index
- // of global_index means that the entry by that index is a ghost
- // one, which should be dropped. Especially for PETSc routines
- // MatZeroRows and MatZeroRowsColumns, which are called to apply the
- // Dirichlet BC, the negative index is not accepted like other
- // matrix or vector PETSc routines. Therefore, the following
- // if-condition is applied.
- if (global_index >= 0)
- {
- bc_values.ids.emplace_back(global_index);
- bc_values.values.push_back(_values[node_id]);
- }
- }
- }
-
-private:
- MeshLib::PropertyVector<double> const& _values;
- MeshLib::Mesh const& _boundary_mesh;
- std::unique_ptr<NumLib::LocalToGlobalIndexMap const> _dof_table_boundary;
- int const _variable_id;
- int const _component_id;
-};
-
-std::unique_ptr<NonuniformDirichletBoundaryCondition>
-createNonuniformDirichletBoundaryCondition(
- BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh,
- NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
- int const component_id);
-
-} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp
deleted file mode 100644
index bbe963d3cae..00000000000
--- a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.cpp
+++ /dev/null
@@ -1,92 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2019, 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 "MeshLib/IO/readMeshFromFile.h"
-#include "ProcessLib/Utils/ProcessUtils.h"
-
-namespace ProcessLib
-{
-std::unique_ptr<NonuniformNeumannBoundaryCondition>
-createNonuniformNeumannBoundaryCondition(
- BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh,
- NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
- 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");
-
- // TODO finally use ProcessLib::Parameter here
- auto const field_name =
- //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__NonuniformNeumann__field_name}
- config.getConfigParameter<std::string>("field_name");
-
- auto const* const property =
- boundary_mesh.getProperties().getPropertyVector<double>(field_name);
-
- if (!property)
- {
- OGS_FATAL("A property with name `%s' does not exist in `%s'.",
- field_name.c_str(), boundary_mesh.getName().c_str());
- }
-
- if (property->getMeshItemType() != MeshLib::MeshItemType::Node)
- {
- OGS_FATAL(
- "Only nodal fields are supported for 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());
- }
-
- std::string const mapping_to_bulk_nodes_property = "bulk_node_ids";
- auto const* const mapping_to_bulk_nodes =
- boundary_mesh.getProperties().getPropertyVector<std::size_t>(
- mapping_to_bulk_nodes_property);
-
- if (mapping_to_bulk_nodes == nullptr ||
- mapping_to_bulk_nodes->getMeshItemType() !=
- MeshLib::MeshItemType::Node ||
- mapping_to_bulk_nodes->getNumberOfComponents() != 1)
- {
- OGS_FATAL("Field `%s' is not set up properly.",
- mapping_to_bulk_nodes_property.c_str());
- }
-
- // 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.
- // TODO (naumov): Add a function to ConfigTree for skipping the tags of the
- // subtree and move the code up in createBoundaryCondition().
- if (boundary_mesh.getDimension() == 0 &&
- boundary_mesh.getNumberOfNodes() == 0 &&
- boundary_mesh.getNumberOfElements() == 0)
- {
- return nullptr;
- }
-#endif // USE_PETSC
-
- return std::make_unique<NonuniformNeumannBoundaryCondition>(
- integration_order, shapefunction_order, dof_table, variable_id,
- component_id, bulk_mesh.getDimension(), boundary_mesh,
- NonuniformNeumannBoundaryConditionData{
- *property, bulk_mesh.getID(), *mapping_to_bulk_nodes, dof_table,
- variable_id, component_id});
-}
-
-} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h b/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h
deleted file mode 100644
index 50cb24c44bb..00000000000
--- a/ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryCondition.h
+++ /dev/null
@@ -1,30 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2019, 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 "MeshLib/PropertyVector.h"
-#include "NonuniformNeumannBoundaryConditionLocalAssembler.h"
-
-namespace ProcessLib
-{
-using NonuniformNeumannBoundaryCondition =
- GenericNonuniformNaturalBoundaryCondition<
- NonuniformNeumannBoundaryConditionData,
- NonuniformNeumannBoundaryConditionLocalAssembler>;
-
-std::unique_ptr<NonuniformNeumannBoundaryCondition>
-createNonuniformNeumannBoundaryCondition(
- BaseLib::ConfigTree const& config, MeshLib::Mesh const& boundary_mesh,
- NumLib::LocalToGlobalIndexMap const& dof_table, int const variable_id,
- int const component_id, unsigned const integration_order,
- unsigned const shapefunction_order, const MeshLib::Mesh& bulk_mesh);
-
-} // namespace ProcessLib
diff --git a/ProcessLib/Parameter/MeshNodeParameter.cpp b/ProcessLib/Parameter/MeshNodeParameter.cpp
index 48fc8f5a281..11b9ecbabe2 100644
--- a/ProcessLib/Parameter/MeshNodeParameter.cpp
+++ b/ProcessLib/Parameter/MeshNodeParameter.cpp
@@ -32,7 +32,7 @@ std::unique_ptr<ParameterBase> createMeshNodeParameter(
field_name.c_str());
}
- return std::make_unique<MeshNodeParameter<double>>(name, *property);
+ return std::make_unique<MeshNodeParameter<double>>(name, mesh, *property);
}
} // namespace ProcessLib
diff --git a/ProcessLib/Parameter/MeshNodeParameter.h b/ProcessLib/Parameter/MeshNodeParameter.h
index 5c013c22dcf..5ea4da35a38 100644
--- a/ProcessLib/Parameter/MeshNodeParameter.h
+++ b/ProcessLib/Parameter/MeshNodeParameter.h
@@ -27,9 +27,9 @@ namespace ProcessLib
template <typename T>
struct MeshNodeParameter final : public Parameter<T> {
MeshNodeParameter(std::string const& name_,
+ MeshLib::Mesh const& mesh,
MeshLib::PropertyVector<T> const& property)
- : Parameter<T>(name_),
- _property(property)
+ : Parameter<T>(name_, &mesh), _property(property)
{
}
--
GitLab