Detailed Description

The ConstraintDirichletBoundaryCondition class describes a Dirichlet-type boundary condition that is constant in space and time where the domain can shrink and grow within the simulation. The expected parameter in the passed configuration is "value" which, when not present defaults to zero.

Definition at line 25 of file ConstraintDirichletBoundaryCondition.h.

#include <ConstraintDirichletBoundaryCondition.h>

Inheritance diagram for ProcessLib::ConstraintDirichletBoundaryCondition:

[legend]

Collaboration diagram for ProcessLib::ConstraintDirichletBoundaryCondition:

[legend]

Public Member Functions
	ConstraintDirichletBoundaryCondition (ParameterLib::Parameter< double > const &parameter, NumLib::LocalToGlobalIndexMap const &dof_table_bulk, int const variable_id, int const component_id, MeshLib::Mesh const &bc_mesh, unsigned const integration_order, MeshLib::Mesh const &bulk_mesh, double const constraint_threshold, bool const lower, std::function< Eigen::Vector3d(std::size_t const, MathLib::Point3d const &, double const, std::vector< GlobalVector * > const &)> getFlux)

void	preTimestep (double const t, std::vector< GlobalVector * > const &x, int const process_id) override

void	getEssentialBCValues (const double t, const GlobalVector &x, NumLib::IndexValueVector< GlobalIndexType > &bc_values) const override
	Writes the values of essential BCs to `bc_values`.

Public Member Functions inherited from ProcessLib::BoundaryCondition
virtual void	applyNaturalBC (const double, std::vector< GlobalVector * > const &, int const, GlobalMatrix , GlobalVector &, GlobalMatrix )

virtual void	postTimestep (const double, std::vector< GlobalVector * > const &, int const)

virtual	~BoundaryCondition ()=default

Private Attributes
ParameterLib::Parameter< double > const &	_parameter

std::unique_ptr< NumLib::LocalToGlobalIndexMap >	_dof_table_boundary

int const	_variable_id

int const	_component_id

MeshLib::Mesh const &	_bc_mesh

unsigned const	_integration_order
	Integration order for integration over the lower-dimensional elements.

std::vector< std::pair< std::size_t, unsigned > >	_bulk_ids

std::vector< double >	_flux_values
	Stores the results of the flux computations per boundary element.

std::vector< std::unique_ptr< ConstraintDirichletBoundaryConditionLocalAssemblerInterface > >	_local_assemblers
	Local assemblers for each boundary element.

double const	_constraint_threshold

bool const	_lower

MeshLib::Mesh const &	_bulk_mesh

std::function< Eigen::Vector3d(std::size_t const, MathLib::Point3d const &, double const, std::vector< GlobalVector * > const &)>	_getFlux
	The function _getFlux calculates the flux through the boundary element.

Constructor & Destructor Documentation

◆ ConstraintDirichletBoundaryCondition()

ProcessLib::ConstraintDirichletBoundaryCondition::ConstraintDirichletBoundaryCondition	(	ParameterLib::Parameter< double > const &	parameter,
		NumLib::LocalToGlobalIndexMap const &	dof_table_bulk,
		int const	variable_id,
		int const	component_id,
		MeshLib::Mesh const &	bc_mesh,
		unsigned const	integration_order,
		MeshLib::Mesh const &	bulk_mesh,
		double const	constraint_threshold,
		bool const	lower,
		std::function< Eigen::Vector3d(std::size_t const, MathLib::Point3d const &, double const, std::vector< GlobalVector * > const &)>	getFlux )

Parameters

parameter	Used for setting the values for the boundary condition.
dof_table_bulk	The bulk local to global index map is used to derive the local to global index map for the boundary.
variable_id	The variable id is needed to determine the global index.
component_id	The component id is needed to determine the global index.
bc_mesh	Lower dimensional mesh the boundary condition is defined on. The bc_mesh must have the two PropertyVector objects 'bulk_element_ids' and 'bulk_node_ids' containing the corresponding information.
integration_order	Order the order of integration used to compute the constraint.
bulk_mesh	The FE mesh for the simulation.
constraint_threshold	The threshold value used for the switch off/on decision.
lower	Boolean value used for the calculation of the constraint criterion, i.e., if lower is set to true the criterion 'calculated_value < constraint_threshold' is evaluated to switch on/off the boundary condition, else 'calculated_value > constraint_threshold' is evaluated.
getFlux	The function used for the flux calculation.

Note: The function has to be stored by value, else the process value is not captured properly.

Definition at line 24 of file ConstraintDirichletBoundaryCondition.cpp.

    : _parameter(parameter),
      _variable_id(variable_id),
      _component_id(component_id),
      _bc_mesh(bc_mesh),
      _integration_order(integration_order),
      _constraint_threshold(constraint_threshold),
      _lower(lower),
      _bulk_mesh(bulk_mesh),
      _getFlux(getFlux)
{
    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));
    }
 
    std::vector<MeshLib::Node*> const& bc_nodes = _bc_mesh.getNodes();
    DBUG(
        "Found {:d} nodes for constraint Dirichlet BCs for the variable {:d} "
        "and component {:d}",
        bc_nodes.size(), variable_id, component_id);
 
    MeshLib::MeshSubset bc_mesh_subset{_bc_mesh, bc_nodes};
 
    // Create local DOF table from intersected mesh subsets for the given
    // variable and component ids.
    _dof_table_boundary = dof_table_bulk.deriveBoundaryConstrainedMap(
        variable_id, {component_id}, std::move(bc_mesh_subset));
 
    auto const& bc_elements(_bc_mesh.getElements());
    _local_assemblers.resize(bc_elements.size());
    _flux_values.resize(bc_elements.size());
    // create _bulk_ids vector
    auto const* bulk_element_ids = MeshLib::bulkElementIDs(_bc_mesh);
    auto const* bulk_node_ids = MeshLib::bulkNodeIDs(_bc_mesh);
    auto const& bulk_nodes = bulk_mesh.getNodes();
 
    auto get_bulk_element_face_id =
        [&](auto const bulk_element_id, MeshLib::Element const* bc_elem)
    {
        auto const* bulk_elem = _bulk_mesh.getElement(bulk_element_id);
        std::array<MeshLib::Node const*, 3> nodes{
            {bulk_nodes[(*bulk_node_ids)[bc_elem->getNode(0)->getID()]],
             bulk_nodes[(*bulk_node_ids)[bc_elem->getNode(1)->getID()]],
             bulk_nodes[(*bulk_node_ids)[bc_elem->getNode(2)->getID()]]}};
        return bulk_elem->identifyFace(nodes.data());
    };
 
    _bulk_ids.reserve(bc_elements.size());
    std::transform(
        begin(bc_elements), end(bc_elements), std::back_inserter(_bulk_ids),
        [&](auto const* bc_element)
        {
            auto const bulk_element_id =
                (*bulk_element_ids)[bc_element->getID()];
            return std::make_pair(
                bulk_element_id,
                get_bulk_element_face_id(bulk_element_id, bc_element));
        });
 
    const int shape_function_order = 1;
 
    BoundaryConditionAndSourceTerm::createLocalAssemblers<
        ConstraintDirichletBoundaryConditionLocalAssembler>(
        _bulk_mesh.getDimension(), _bc_mesh.getElements(), *_dof_table_boundary,
        shape_function_order, _local_assemblers,
        NumLib::IntegrationOrder{_integration_order},
        _bc_mesh.isAxiallySymmetric(), bulk_mesh, _bulk_ids);
}

References _bc_mesh, _bulk_ids, _bulk_mesh, _dof_table_boundary, _flux_values, _integration_order, _local_assemblers, MeshLib::bulkElementIDs(), MeshLib::bulkNodeIDs(), ProcessLib::BoundaryConditionAndSourceTerm::createLocalAssemblers(), DBUG(), NumLib::LocalToGlobalIndexMap::deriveBoundaryConstrainedMap(), MeshLib::Mesh::getDimension(), MeshLib::Mesh::getElement(), MeshLib::Mesh::getElements(), MeshLib::Mesh::getNodes(), NumLib::LocalToGlobalIndexMap::getNumberOfVariableComponents(), NumLib::LocalToGlobalIndexMap::getNumberOfVariables(), MeshLib::Element::identifyFace(), MeshLib::Mesh::isAxiallySymmetric(), and OGS_FATAL.

Member Function Documentation

◆ getEssentialBCValues()

void ProcessLib::ConstraintDirichletBoundaryCondition::getEssentialBCValues	(	const double	,
		const GlobalVector &	,
		NumLib::IndexValueVector< GlobalIndexType > &	) const

overridevirtual

Writes the values of essential BCs to bc_values.

Reimplemented from ProcessLib::BoundaryCondition.

Definition at line 127 of file ConstraintDirichletBoundaryCondition.cpp.

{
    ParameterLib::SpatialPosition pos;
 
    bc_values.ids.clear();
    bc_values.values.clear();
 
    std::vector<std::pair<GlobalIndexType, double>> tmp_bc_values;
 
    auto isFlux = [&](const std::size_t element_id)
    {
        return _lower ? _flux_values[element_id] < _constraint_threshold
                      : _flux_values[element_id] > _constraint_threshold;
    };
 
    for (auto const* boundary_element : _bc_mesh.getElements())
    {
        // check if the boundary element is active
        if (isFlux(boundary_element->getID()))
        {
            continue;
        }
 
        // loop over the boundary element nodes and set the dirichlet value
        unsigned const number_nodes = boundary_element->getNumberOfNodes();
        for (unsigned i = 0; i < number_nodes; ++i)
        {
            auto const id = boundary_element->getNode(i)->getID();
            pos.setAll(id, boundary_element->getID(), {}, {});
 
            MeshLib::Location l(_bc_mesh.getID(), MeshLib::MeshItemType::Node,
                                id);
            // TODO: that might be slow, but only done once
            const auto g_idx = _dof_table_boundary->getGlobalIndex(
                l, _variable_id, _component_id);
            if (g_idx == NumLib::MeshComponentMap::nop)
            {
                continue;
            }
            // For the DDC approach (e.g. with PETSc option), the negative
            // index of g_idx 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 (g_idx >= 0)
            {
                tmp_bc_values.emplace_back(g_idx, _parameter(t, pos).front());
            }
        }
    }
 
    if (tmp_bc_values.empty())
    {
        DBUG("The domain on which the boundary is defined is empty");
        return;
    }
 
    // store unique pairs of node id and value with respect to the node id in
    // the bc_values vector. The values related to the same node id are
    // averaged.
    // first: sort the (node id, value) pairs according to the node id
    std::sort(tmp_bc_values.begin(), tmp_bc_values.end(),
              [](std::pair<GlobalIndexType, double> const& a,
                 std::pair<GlobalIndexType, double> const& b)
              { return a.first < b.first; });
    // second: average the values over equal node id ranges
    unsigned cnt = 1;
    GlobalIndexType current_id = tmp_bc_values.begin()->first;
    double sum = tmp_bc_values.begin()->second;
    for (auto const& tmp_bc_value : tmp_bc_values)
    {
        if (tmp_bc_value.first == current_id)
        {
            cnt++;
            sum += tmp_bc_value.second;
        }
        else
        {
            bc_values.ids.emplace_back(current_id);
            bc_values.values.emplace_back(sum / cnt);
            cnt = 1;
            sum = tmp_bc_value.second;
            current_id = tmp_bc_value.first;
        }
    }
    bc_values.ids.emplace_back(current_id);
    bc_values.values.emplace_back(sum / cnt);
 
    DBUG("Found {:d} constraint dirichlet boundary condition values.",
         bc_values.ids.size());
}

References _bc_mesh, _component_id, _constraint_threshold, _dof_table_boundary, _flux_values, _lower, _parameter, _variable_id, DBUG(), MeshLib::Mesh::getElements(), MeshLib::Mesh::getID(), NumLib::IndexValueVector< typename >::ids, MeshLib::Node, NumLib::MeshComponentMap::nop, ParameterLib::SpatialPosition::setAll(), and NumLib::IndexValueVector< typename >::values.

◆ preTimestep()

void ProcessLib::ConstraintDirichletBoundaryCondition::preTimestep	(	double const	t,
		std::vector< GlobalVector * > const &	x,
		int const	process_id )

overridevirtual

Reimplemented from ProcessLib::BoundaryCondition.

Definition at line 110 of file ConstraintDirichletBoundaryCondition.cpp.

{
    DBUG(
        "ConstraintDirichletBoundaryCondition::preTimestep: computing flux "
        "constraints");
    for (auto const id : _bc_mesh.getElements() | MeshLib::views::ids)
    {
        _flux_values[id] = _local_assemblers[id]->integrate(
            x, t,
            [this](std::size_t const element_id, MathLib::Point3d const& pnt,
                   double const t, std::vector<GlobalVector*> const& x)
            { return _getFlux(element_id, pnt, t, x); });
    }
}

References _bc_mesh, _flux_values, _getFlux, _local_assemblers, DBUG(), MeshLib::Mesh::getElements(), and MeshLib::views::ids.

Member Data Documentation

◆ _bc_mesh

MeshLib::Mesh const& ProcessLib::ConstraintDirichletBoundaryCondition::_bc_mesh

private

Vector of (lower-dimensional) boundary elements on which the boundary condition is defined.

Definition at line 82 of file ConstraintDirichletBoundaryCondition.h.

Referenced by ConstraintDirichletBoundaryCondition(), getEssentialBCValues(), and preTimestep().

◆ _bulk_ids

std::vector<std::pair<std::size_t, unsigned> > ProcessLib::ConstraintDirichletBoundaryCondition::_bulk_ids

private

The first item of the pair is the element id in the bulk mesh, the second item is the face id of the bulk element that is part of the boundary

Definition at line 90 of file ConstraintDirichletBoundaryCondition.h.

Referenced by ConstraintDirichletBoundaryCondition().

◆ _bulk_mesh

MeshLib::Mesh const& ProcessLib::ConstraintDirichletBoundaryCondition::_bulk_mesh

private

The mesh _bulk_mesh is the discretized domain the process(es) are defined on. It is needed to get values for the constraint calculation.

Definition at line 112 of file ConstraintDirichletBoundaryCondition.h.

Referenced by ConstraintDirichletBoundaryCondition().

◆ _component_id

int const ProcessLib::ConstraintDirichletBoundaryCondition::_component_id

private

Definition at line 78 of file ConstraintDirichletBoundaryCondition.h.

Referenced by getEssentialBCValues().

◆ _constraint_threshold

double const ProcessLib::ConstraintDirichletBoundaryCondition::_constraint_threshold

private

The threshold value used to the switch off/on the Dirichlet-type boundary condition.

Definition at line 102 of file ConstraintDirichletBoundaryCondition.h.

Referenced by getEssentialBCValues().

◆ _dof_table_boundary

std::unique_ptr<NumLib::LocalToGlobalIndexMap> ProcessLib::ConstraintDirichletBoundaryCondition::_dof_table_boundary

private

Local dof table, a subset of the global one restricted to the participating number of elements of the boundary condition.

Definition at line 75 of file ConstraintDirichletBoundaryCondition.h.

Referenced by ConstraintDirichletBoundaryCondition(), and getEssentialBCValues().

◆ _flux_values

std::vector<double> ProcessLib::ConstraintDirichletBoundaryCondition::_flux_values

private

Stores the results of the flux computations per boundary element.

Definition at line 93 of file ConstraintDirichletBoundaryCondition.h.

Referenced by ConstraintDirichletBoundaryCondition(), getEssentialBCValues(), and preTimestep().

◆ _getFlux

std::function<Eigen::Vector3d(std::size_t const, MathLib::Point3d const&, double const, std::vector<GlobalVector*> const&)> ProcessLib::ConstraintDirichletBoundaryCondition::_getFlux

private

The function _getFlux calculates the flux through the boundary element.

Definition at line 118 of file ConstraintDirichletBoundaryCondition.h.

Referenced by preTimestep().

◆ _integration_order

unsigned const ProcessLib::ConstraintDirichletBoundaryCondition::_integration_order

private

Integration order for integration over the lower-dimensional elements.

Definition at line 85 of file ConstraintDirichletBoundaryCondition.h.

Referenced by ConstraintDirichletBoundaryCondition().

◆ _local_assemblers

std::vector<std::unique_ptr< ConstraintDirichletBoundaryConditionLocalAssemblerInterface> > ProcessLib::ConstraintDirichletBoundaryCondition::_local_assemblers

private

Local assemblers for each boundary element.

Definition at line 98 of file ConstraintDirichletBoundaryCondition.h.

Referenced by ConstraintDirichletBoundaryCondition(), and preTimestep().

◆ _lower

bool const ProcessLib::ConstraintDirichletBoundaryCondition::_lower

private

The boolean value lower is used for the calculation of the constraint criterion, i.e., if lower is set to true the criterion 'calculated_value < constraint_threshold' is evaluated to switch on/off the boundary condition, else 'calculated_value > constraint_threshold' is evaluated.

Definition at line 108 of file ConstraintDirichletBoundaryCondition.h.

Referenced by getEssentialBCValues().