Commit 86e72910 authored by renchao.lu's avatar renchao.lu

[PL/BC] Implement solution dependent dirichlet boundary condition.

parent 2c68e01c
......@@ -21,6 +21,7 @@
#include "PhaseFieldIrreversibleDamageOracleBoundaryCondition.h"
#include "PrimaryVariableConstraintDirichletBoundaryCondition.h"
#include "RobinBoundaryCondition.h"
#include "SolutionDependentDirichletBoundaryCondition.h"
#include "VariableDependentNeumannBoundaryCondition.h"
#include "BaseLib/TimeInterval.h"
......@@ -114,6 +115,12 @@ std::unique_ptr<BoundaryCondition> createBoundaryCondition(
config.config, config.boundary_mesh, dof_table, variable_id,
*config.component_id, parameters);
}
if (type == "SolutionDependentDirichlet")
{
return ProcessLib::createSolutionDependentDirichletBoundaryCondition(
config.config, config.boundary_mesh, dof_table, variable_id,
*config.component_id, parameters);
}
if (type == "HCNonAdvectiveFreeComponentFlowBoundary")
{
return createHCNonAdvectiveFreeComponentFlowBoundaryCondition(
......
/**
* \file
* \copyright
* Copyright (c) 2012-2020, 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 "SolutionDependentDirichletBoundaryCondition.h"
#include <memory>
#include <vector>
#include "BaseLib/ConfigTree.h"
#include "BaseLib/Logging.h"
#include "DirichletBoundaryConditionAuxiliaryFunctions.h"
#include "MeshLib/Mesh.h"
#include "NumLib/DOF/LocalToGlobalIndexMap.h"
#include "ParameterLib/Utils.h"
namespace ProcessLib
{
SolutionDependentDirichletBoundaryCondition::
SolutionDependentDirichletBoundaryCondition(
ParameterLib::Parameter<double> const& parameter,
MeshLib::Mesh const& bc_mesh,
NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
int const variable_id, int const component_id)
: _bc_mesh(bc_mesh), _variable_id(variable_id), _component_id(component_id)
{
checkParametersOfDirichletBoundaryCondition(_bc_mesh, dof_table_bulk,
_variable_id, _component_id);
std::vector<MeshLib::Node*> const& bc_nodes = bc_mesh.getNodes();
MeshLib::MeshSubset bc_mesh_subset(_bc_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(bc_mesh_subset)));
std::string const property_name = "solution_dependent_bc";
if (bc_mesh.getProperties().existsPropertyVector<double>(property_name))
{
OGS_FATAL(
"Found mesh property '{:s}' which is the built-in property of the "
"class SolutionDependentDirichletBoundaryCondition.",
property_name);
}
auto& solution_dependent_bc = *MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(bc_mesh), property_name,
MeshLib::MeshItemType::Node, 1);
solution_dependent_bc.resize(bc_mesh.getNumberOfNodes());
ParameterLib::SpatialPosition pos;
auto const& nodes = bc_mesh.getNodes();
for (std::size_t i = 0; i < solution_dependent_bc.size(); ++i)
{
auto const id = nodes[i]->getID();
pos.setNodeID(id);
solution_dependent_bc[i] = parameter(0, pos)[0];
}
_parameter = std::make_unique<ParameterLib::MeshNodeParameter<double>>(
property_name, bc_mesh, solution_dependent_bc);
}
void SolutionDependentDirichletBoundaryCondition::getEssentialBCValues(
double const t, GlobalVector const& x,
NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
{
getEssentialBCValuesLocal(*_parameter, _bc_mesh, _bc_mesh.getNodes(),
*_dof_table_boundary, _variable_id, _component_id,
t, x, bc_values);
}
void SolutionDependentDirichletBoundaryCondition::postTimestep(
double const /*t*/, std::vector<GlobalVector*> const& x,
int const process_id)
{
auto& solution_dependent_bc =
*const_cast<MeshLib::Properties&>(_bc_mesh.getProperties())
.getPropertyVector<double>("solution_dependent_bc");
auto const& nodes = _bc_mesh.getNodes();
for (std::size_t i = 0; i < solution_dependent_bc.size(); ++i)
{
auto const id = nodes[i]->getID();
auto const global_index = _dof_table_boundary->getGlobalIndex(
{_bc_mesh.getID(), MeshLib::MeshItemType::Node, id}, _variable_id,
_component_id);
assert(global_index >= 0);
solution_dependent_bc[i] = x[process_id]->get(global_index);
}
}
std::unique_ptr<SolutionDependentDirichletBoundaryCondition>
createSolutionDependentDirichletBoundaryCondition(
BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
int const component_id,
std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters)
{
DBUG(
"Constructing SolutionDependentDirichletBoundaryCondition from "
"config.");
//! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
config.checkConfigParameter("type", "SolutionDependentDirichlet");
auto& initial_value_parameter = ParameterLib::findParameter<double>(
//! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__SolutionDependentDirichlet__initial_value_parameter}
config.getConfigParameter<std::string>("initial_value_parameter"),
parameters, 1, &bc_mesh);
// 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 (bc_mesh.getDimension() == 0 && bc_mesh.getNumberOfNodes() == 0 &&
bc_mesh.getNumberOfElements() == 0)
{
return nullptr;
}
#endif // USE_PETSC
return std::make_unique<SolutionDependentDirichletBoundaryCondition>(
initial_value_parameter, bc_mesh, dof_table_bulk, variable_id,
component_id);
}
} // namespace ProcessLib
/**
* \file
* \copyright
* Copyright (c) 2012-2020, 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 "ParameterLib/MeshNodeParameter.h"
namespace ProcessLib
{
/// The SolutionDependentDirichletBoundaryCondition belongs to the
/// Dirichlet-type boundary condition.
///
/// This class is a special category of Dirichlet boundary condition,
/// applied in the situation where the value assigned for the boundary condition
/// is dependent on the process solution of the last time step. This particular
/// boundary condition is widely used in the reactive transport problems and has
/// the potential to be used in other processes.
class SolutionDependentDirichletBoundaryCondition final
: public BoundaryCondition
{
public:
SolutionDependentDirichletBoundaryCondition(
ParameterLib::Parameter<double> const& parameter,
MeshLib::Mesh const& bc_mesh,
NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
int const variable_id, int const component_id);
void getEssentialBCValues(
double const t, GlobalVector const& x,
NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
/// Renchao: The original idea to place the update of the boundary condition
/// value at the preTimestep stage. The update could be achieved within the
/// class function "Process::preTimestep". Whereas, I find it not doable to
/// implement in this way. The class function "Process::preTimestep" is
/// called in a row by the functions "preTimestepForAllProcesses" and
/// "TimeLoop::outputSolutions". These two functions are called when
/// initializing and subsequently looping over the "TimeLoop". Actually, it
/// is not intended to make the boundary condition value updated in the
/// first loop. Instead, the update is intended to start from the second
/// loop. For these two reasons, I think it more proper to do the
/// implementation at the postTimestep stage.
void postTimestep(double const /*t*/,
std::vector<GlobalVector*> const& x,
int const process_id) override;
private:
MeshLib::Mesh const& _bc_mesh;
int const _variable_id;
int const _component_id;
std::unique_ptr<NumLib::LocalToGlobalIndexMap const> _dof_table_boundary;
std::unique_ptr<ParameterLib::MeshNodeParameter<double>> _parameter;
};
std::unique_ptr<SolutionDependentDirichletBoundaryCondition>
createSolutionDependentDirichletBoundaryCondition(
BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
int const component_id,
std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
parameters);
} // namespace ProcessLib
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment