Loading...
Searching...
No Matches
DirichletBC.h
// Copyright (C) 2007-2024 Michal Habera, Anders Logg, Garth N. Wells, Jørgen
// S.Dokken and Paul T. Kühner
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier: LGPL-3.0-or-later
#pragma once
#include "Constant.h"
#include "DofMap.h"
#include "Function.h"
#include "FunctionSpace.h"
#include <algorithm>
#include <array>
#include <concepts>
#include <dolfinx/common/types.h>
#include <functional>
#include <memory>
#include <optional>
#include <span>
#include <type_traits>
#include <utility>
#include <variant>
#include <vector>
namespace dolfinx::fem
{
std::vector<std::int32_t>
locate_dofs_topological(const mesh::Topology& topology, const DofMap& dofmap,
int dim, std::span<const std::int32_t> entities,
bool remote = true);
std::array<std::vector<std::int32_t>, 2> locate_dofs_topological(
const mesh::Topology& topology,
std::array<std::reference_wrapper<const DofMap>, 2> dofmaps, int dim,
std::span<const std::int32_t> entities, bool remote = true);
template <std::floating_point T, typename U>
std::vector<std::int32_t> locate_dofs_geometrical(const FunctionSpace<T>& V,
U marker_fn)
{
// FIXME: Calling V.tabulate_dof_coordinates() is very expensive,
// especially when we usually want the boundary dofs only. Add
// interface that computes dofs coordinates only for specified cell.
assert(V.element());
if (V.element()->is_mixed())
{
throw std::runtime_error(
"Cannot locate dofs geometrically for mixed space. Use subspaces.");
}
// Compute dof coordinates
const std::vector<T> dof_coordinates = V.tabulate_dof_coordinates(true);
using cmdspan3x_t
= md::mdspan<const T, md::extents<std::size_t, 3, md::dynamic_extent>>;
// Compute marker for each dof coordinate
cmdspan3x_t x(dof_coordinates.data(), 3, dof_coordinates.size() / 3);
const std::vector<std::int8_t> marked_dofs = marker_fn(x);
std::vector<std::int32_t> dofs;
dofs.reserve(std::count(marked_dofs.begin(), marked_dofs.end(), true));
for (std::size_t i = 0; i < marked_dofs.size(); ++i)
{
if (marked_dofs[i])
dofs.push_back(i);
}
return dofs;
}
template <std::floating_point T, typename U>
std::array<std::vector<std::int32_t>, 2> locate_dofs_geometrical(
std::array<std::reference_wrapper<const FunctionSpace<T>>, 2> V,
U marker_fn)
{
// FIXME: Calling V.tabulate_dof_coordinates() is very expensive,
// especially when we usually want the boundary dofs only. Add
// interface that computes dofs coordinates only for specified cell.
// Get function spaces
const FunctionSpace<T>& V0 = V.at(0).get();
const FunctionSpace<T>& V1 = V.at(1).get();
// Get mesh
auto mesh = V0.mesh();
assert(mesh);
assert(V1.mesh());
if (mesh != V1.mesh())
throw std::runtime_error("Meshes are not the same.");
const int tdim = mesh->topology()->dim();
assert(V0.element());
assert(V1.element());
if (*V0.element() != *V1.element())
throw std::runtime_error("Function spaces must have the same element.");
// Compute dof coordinates
const std::vector<T> dof_coordinates = V1.tabulate_dof_coordinates(true);
using cmdspan3x_t
= md::mdspan<const T, md::extents<std::size_t, 3, md::dynamic_extent>>;
// Evaluate marker for each dof coordinate
cmdspan3x_t x(dof_coordinates.data(), 3, dof_coordinates.size() / 3);
const std::vector<std::int8_t> marked_dofs = marker_fn(x);
// Get dofmaps
std::shared_ptr<const DofMap> dofmap0 = V0.dofmap();
assert(dofmap0);
const int bs0 = dofmap0->bs();
std::shared_ptr<const DofMap> dofmap1 = V1.dofmap();
assert(dofmap1);
const int bs1 = dofmap1->bs();
const int element_bs = dofmap0->element_dof_layout().block_size();
assert(element_bs == dofmap1->element_dof_layout().block_size());
// Iterate over cells
auto topology = mesh->topology();
assert(topology);
std::vector<std::array<std::int32_t, 2>> bc_dofs;
for (int c = 0; c < topology->connectivity(tdim, 0)->num_nodes(); ++c)
{
// Get cell dofmaps
auto cell_dofs0 = dofmap0->cell_dofs(c);
auto cell_dofs1 = dofmap1->cell_dofs(c);
// Loop over cell dofs and add to bc_dofs if marked.
for (std::size_t i = 0; i < cell_dofs1.size(); ++i)
{
if (marked_dofs[cell_dofs1[i]])
{
// Unroll over blocks
for (int k = 0; k < element_bs; ++k)
{
const int local_pos = element_bs * i + k;
const std::div_t pos0 = std::div(local_pos, bs0);
const std::div_t pos1 = std::div(local_pos, bs1);
const std::int32_t dof_index0
= bs0 * cell_dofs0[pos0.quot] + pos0.rem;
const std::int32_t dof_index1
= bs1 * cell_dofs1[pos1.quot] + pos1.rem;
bc_dofs.push_back({dof_index0, dof_index1});
}
}
}
}
// Remove duplicates
std::ranges::sort(bc_dofs);
auto [unique_end, range_end] = std::ranges::unique(bc_dofs);
bc_dofs.erase(unique_end, range_end);
// Copy to separate array
std::array dofs = {std::vector<std::int32_t>(bc_dofs.size()),
std::vector<std::int32_t>(bc_dofs.size())};
std::ranges::transform(bc_dofs, dofs[0].begin(),
[](auto dof) { return dof[0]; });
std::ranges::transform(bc_dofs, dofs[1].begin(),
[](auto dof) { return dof[1]; });
return dofs;
}
template <dolfinx::scalar T, std::floating_point U = dolfinx::scalar_value_t<T>>
class DirichletBC
{
private:
std::size_t num_owned(const DofMap& dofmap,
std::span<const std::int32_t> dofs)
{
int bs = dofmap.index_map_bs();
std::int32_t map_size = dofmap.index_map->size_local();
std::int32_t owned_size = bs * map_size;
auto it = std::ranges::lower_bound(dofs, owned_size);
return std::distance(dofs.begin(), it);
}
static std::vector<std::int32_t>
unroll_dofs(std::span<const std::int32_t> dofs, int bs)
{
std::vector<std::int32_t> dofs_unrolled(bs * dofs.size());
for (std::size_t i = 0; i < dofs.size(); ++i)
for (int k = 0; k < bs; ++k)
dofs_unrolled[bs * i + k] = bs * dofs[i] + k;
return dofs_unrolled;
}
public:
template <typename S, typename X>
requires(std::is_convertible_v<S, T>
|| std::is_convertible_v<S, std::span<const T>>)
&& std::is_convertible_v<std::remove_cvref_t<X>,
std::vector<std::int32_t>>
DirichletBC(const S& g, X&& dofs, std::shared_ptr<const FunctionSpace<U>> V)
: DirichletBC(std::make_shared<Constant<T>>(g), dofs, V)
{
}
template <typename X>
requires std::is_convertible_v<std::remove_cvref_t<X>,
std::vector<std::int32_t>>
DirichletBC(std::shared_ptr<const Constant<T>> g, X&& dofs,
std::shared_ptr<const FunctionSpace<U>> V)
: _function_space(V), _g(g), _dofs0(std::forward<X>(dofs))
{
assert(g);
assert(V);
if (g->shape.size() != V->element()->value_shape().size())
{
throw std::runtime_error(
"Rank mismatch between Constant and function space in DirichletBC");
}
if (g->value.size() != (std::size_t)_function_space->dofmap()->bs())
{
throw std::runtime_error(
"Creating a DirichletBC using a Constant is not supported when the "
"Constant size is not equal to the block size of the constrained "
"(sub-)space. Use a fem::Function to create the fem::DirichletBC.");
}
if (!V->element()->interpolation_ident())
{
throw std::runtime_error(
"Constant can be used only with point-evaluation elements");
}
// Unroll _dofs0 if dofmap block size > 1
if (const int bs = V->dofmap()->bs(); bs > 1)
_dofs0 = unroll_dofs(_dofs0, bs);
_owned_indices0 = num_owned(*_function_space->dofmap(), _dofs0);
}
template <typename X>
requires std::is_convertible_v<std::remove_cvref_t<X>,
std::vector<std::int32_t>>
DirichletBC(std::shared_ptr<const Function<T, U>> g, X&& dofs)
: _function_space(g->function_space()), _g(g),
_dofs0(std::forward<X>(dofs))
{
assert(_function_space);
// Unroll _dofs0 if dofmap block size > 1
if (const int bs = _function_space->dofmap()->bs(); bs > 1)
_dofs0 = unroll_dofs(_dofs0, bs);
_owned_indices0 = num_owned(*_function_space->dofmap(), _dofs0);
}
template <typename X>
DirichletBC(std::shared_ptr<const Function<T, U>> g, X&& V_g_dofs,
std::shared_ptr<const FunctionSpace<U>> V)
: _function_space(V), _g(g),
_dofs0(std::forward<typename std::remove_reference_t<X>::value_type>(
V_g_dofs[0])),
_dofs1_g(std::forward<typename std::remove_reference_t<X>::value_type>(
V_g_dofs[1])),
_owned_indices0(num_owned(*_function_space->dofmap(), _dofs0))
{
}
DirichletBC(const DirichletBC& bc) = default;
DirichletBC(DirichletBC&& bc) = default;
~DirichletBC() = default;
DirichletBC& operator=(const DirichletBC& bc) = default;
DirichletBC& operator=(DirichletBC&& bc) = default;
std::shared_ptr<const FunctionSpace<U>> function_space() const
{
return _function_space;
}
std::variant<std::shared_ptr<const Function<T, U>>,
std::shared_ptr<const Constant<T>>>
value() const
{
return _g;
}
std::pair<std::span<const std::int32_t>, std::int32_t> dof_indices() const
{
return {_dofs0, _owned_indices0};
}
void set(std::span<T> x, std::optional<std::span<const T>> x0,
T alpha = 1) const
{
// set_fn is a lambda which gets evaluated for every index in [0,
// _dofs0.size()) and its result is assigned to x[_dofs0[i]].
auto apply = [&](std::invocable<std::int32_t> auto set_fn)
{
static_assert(
std::is_same_v<std::invoke_result_t<decltype(set_fn), std::int32_t>,
T>);
std::int32_t x_size = x.size();
for (std::size_t i = 0; i < _dofs0.size(); ++i)
{
if (_dofs0[i] < x_size)
x[_dofs0[i]] = set_fn(i);
}
};
if (alpha == T(0)) // Optimisation for when alpha == 0
{
apply([](std::int32_t) -> T { return 0; });
return;
}
if (std::holds_alternative<std::shared_ptr<const Function<T, U>>>(_g))
{
auto g = std::get<std::shared_ptr<const Function<T, U>>>(_g);
assert(g);
std::span<const T> values = g->x()->array();
// Extract degrees of freedom associated with g. If g is in a collapsed
// sub-space, get the dofs in this space, otherwise the degrees of g is
// the same as for x
auto dofs_g = _dofs1_g.empty() ? std::span(_dofs0) : std::span(_dofs1_g);
if (x0)
{
assert(x.size() <= x0->size());
apply(
[dofs_g, x0 = *x0, alpha, values,
&dofs0 = this->_dofs0](std::int32_t i) -> T
{
assert(dofs_g[i] < static_cast<std::int32_t>(values.size()));
return alpha * (values[dofs_g[i]] - x0[dofs0[i]]);
});
}
else
{
apply(
[dofs_g, values, alpha](std::int32_t i) -> T
{
assert(dofs_g[i] < static_cast<std::int32_t>(values.size()));
return alpha * values[dofs_g[i]];
});
}
}
else if (std::holds_alternative<std::shared_ptr<const Constant<T>>>(_g))
{
auto g = std::get<std::shared_ptr<const Constant<T>>>(_g);
const std::vector<T>& value = g->value;
std::int32_t bs = _function_space->dofmap()->bs();
if (x0)
{
assert(x.size() <= x0->size());
apply(
[x0 = *x0, alpha, bs, &value, &dofs0 = _dofs0](std::int32_t i) -> T
{
auto dof = dofs0[i];
return alpha * (value[dof % bs] - x0[dof]);
});
}
else
{
apply([alpha, bs, &value, &dofs0 = _dofs0](std::int32_t i) -> T
{ return alpha * value[dofs0[i] % bs]; });
}
}
else
{
// replace with std::unreachable once C++23 is supported
assert(false);
}
}
void mark_dofs(std::span<std::int8_t> markers) const
{
for (std::int32_t idx : _dofs0)
{
assert(idx < (std::int32_t)markers.size());
markers[idx] = true;
}
}
private:
// The function space (possibly a sub function space)
std::shared_ptr<const FunctionSpace<U>> _function_space;
// The function
std::variant<std::shared_ptr<const Function<T, U>>,
std::shared_ptr<const Constant<T>>>
_g;
// Dof indices (_dofs0) in _function_space and (_dofs1_g) in the space
// of _g. _dofs1_g may be empty if _dofs0 can be re-used
std::vector<std::int32_t> _dofs0, _dofs1_g;
// The first _owned_indices in _dofs are owned by this process
std::int32_t _owned_indices0 = -1;
};
} // namespace dolfinx::fem
DofMap.h
Degree-of-freedom map representations and tools.
dolfinx::fem::Constant
Constant (in space) value which can be attached to a Form.
Definition Constant.h:22
dolfinx::fem::DirichletBC::DirichletBC
DirichletBC(const S &g, X &&dofs, std::shared_ptr< const FunctionSpace< U > > V)
Create a representation of a Dirichlet boundary condition constrained by a scalar- or vector-valued c...
Definition DirichletBC.h:303
dolfinx::fem::DirichletBC::operator=
DirichletBC & operator=(const DirichletBC &bc)=default
dolfinx::fem::DirichletBC::DirichletBC
DirichletBC(std::shared_ptr< const Constant< T > > g, X &&dofs, std::shared_ptr< const FunctionSpace< U > > V)
Create a representation of a Dirichlet boundary condition constrained by a fem::Constant.
Definition DirichletBC.h:326
dolfinx::fem::DirichletBC::set
void set(std::span< T > x, std::optional< std::span< const T > > x0, T alpha=1) const
Set entries in an array that are constrained by Dirichlet boundary conditions.
Definition DirichletBC.h:487
dolfinx::fem::DirichletBC::value
std::variant< std::shared_ptr< const Function< T, U > >, std::shared_ptr< const Constant< T > > > value() const
Definition DirichletBC.h:449
dolfinx::fem::DirichletBC::DirichletBC
DirichletBC(const DirichletBC &bc)=default
dolfinx::fem::DirichletBC::DirichletBC
DirichletBC(DirichletBC &&bc)=default
dolfinx::fem::DirichletBC::dof_indices
std::pair< std::span< const std::int32_t >, std::int32_t > dof_indices() const
Definition DirichletBC.h:460
dolfinx::fem::DirichletBC::function_space
std::shared_ptr< const FunctionSpace< U > > function_space() const
Definition DirichletBC.h:440
dolfinx::fem::DirichletBC::~DirichletBC
~DirichletBC()=default
Destructor.
dolfinx::fem::DirichletBC::DirichletBC
DirichletBC(std::shared_ptr< const Function< T, U > > g, X &&dofs)
Create a representation of a Dirichlet boundary condition where the space being constrained is the sa...
Definition DirichletBC.h:374
dolfinx::fem::DirichletBC::operator=
DirichletBC & operator=(DirichletBC &&bc)=default
Move assignment operator.
dolfinx::fem::DirichletBC::mark_dofs
void mark_dofs(std::span< std::int8_t > markers) const
Set markers[i] = true if dof i has a boundary condition applied.
Definition DirichletBC.h:581
dolfinx::fem::DirichletBC::DirichletBC
DirichletBC(std::shared_ptr< const Function< T, U > > g, X &&V_g_dofs, std::shared_ptr< const FunctionSpace< U > > V)
Create a representation of a Dirichlet boundary condition where the space being constrained and the f...
Definition DirichletBC.h:409
dolfinx::fem::DofMap
Degree-of-freedom map.
Definition DofMap.h:73
dolfinx::fem::DofMap::index_map
std::shared_ptr< const common::IndexMap > index_map
Index map that describes the parallel distribution of the dofmap.
Definition DofMap.h:164
dolfinx::fem::DofMap::index_map_bs
int index_map_bs() const
Block size associated with the index_map.
Definition DofMap.cpp:276
dolfinx::fem::FunctionSpace
This class represents a finite element function space defined by a mesh, a finite element,...
Definition FunctionSpace.h:34
dolfinx::fem::FunctionSpace::dofmap
std::shared_ptr< const DofMap > dofmap() const
The dofmap.
Definition FunctionSpace.h:361
dolfinx::fem::FunctionSpace::mesh
std::shared_ptr< const mesh::Mesh< geometry_type > > mesh() const
The mesh.
Definition FunctionSpace.h:336
dolfinx::fem::FunctionSpace::element
std::shared_ptr< const FiniteElement< geometry_type > > element() const
The finite element.
Definition FunctionSpace.h:342
dolfinx::fem::FunctionSpace::tabulate_dof_coordinates
std::vector< geometry_type > tabulate_dof_coordinates(bool transpose) const
Tabulate the physical coordinates of all dofs on this process.
Definition FunctionSpace.h:208
dolfinx::fem::Function
Definition Function.h:47
dolfinx::fem
Finite element method functionality.
Definition assemble_expression_impl.h:23
dolfinx::fem::locate_dofs_geometrical
std::vector< std::int32_t > locate_dofs_geometrical(const FunctionSpace< T > &V, U marker_fn)
Find degrees of freedom whose geometric coordinate is true for the provided marking function.
Definition DirichletBC.h:103
dolfinx::fem::locate_dofs_topological
std::vector< std::int32_t > locate_dofs_topological(const mesh::Topology &topology, const DofMap &dofmap, int dim, std::span< const std::int32_t > entities, bool remote=true)
Find degrees-of-freedom which belong to the provided mesh entities (topological).
Definition DirichletBC.cpp:188
dolfinx::mesh
Mesh data structures and algorithms on meshes.
Definition DofMap.h:32