13#include "FunctionSpace.h"
17#include <dolfinx/common/types.h>
55std::vector<std::int32_t>
57 int dim, std::span<const std::int32_t> entities,
87 const mesh::Topology& topology,
88 std::array<std::reference_wrapper<const DofMap>, 2> dofmaps,
int dim,
89 std::span<const std::int32_t> entities,
bool remote =
true);
102template <std::
floating_po
int T,
typename U>
111 if (V.element()->is_mixed())
113 throw std::runtime_error(
114 "Cannot locate dofs geometrically for mixed space. Use subspaces.");
118 const std::vector<T> dof_coordinates = V.tabulate_dof_coordinates(
true);
121 = md::mdspan<const T, md::extents<std::size_t, 3, md::dynamic_extent>>;
124 cmdspan3x_t x(dof_coordinates.data(), 3, dof_coordinates.size() / 3);
125 const std::vector<std::int8_t> marked_dofs = marker_fn(x);
127 std::vector<std::int32_t> dofs;
128 dofs.reserve(std::count(marked_dofs.begin(), marked_dofs.end(),
true));
129 for (std::size_t i = 0; i < marked_dofs.size(); ++i)
151template <std::
floating_po
int T,
typename U>
165 auto mesh = V0.mesh();
169 throw std::runtime_error(
"Meshes are not the same.");
170 const int tdim =
mesh->topology()->dim();
172 assert(V0.element());
174 if (*V0.element() != *V1.
element())
175 throw std::runtime_error(
"Function spaces must have the same element.");
181 = md::mdspan<const T, md::extents<std::size_t, 3, md::dynamic_extent>>;
184 cmdspan3x_t x(dof_coordinates.data(), 3, dof_coordinates.size() / 3);
185 const std::vector<std::int8_t> marked_dofs = marker_fn(x);
188 std::shared_ptr<const DofMap> dofmap0 = V0.dofmap();
190 const int bs0 = dofmap0->bs();
191 std::shared_ptr<const DofMap> dofmap1 = V1.
dofmap();
193 const int bs1 = dofmap1->bs();
195 const int element_bs = dofmap0->element_dof_layout().block_size();
196 assert(element_bs == dofmap1->element_dof_layout().block_size());
199 auto topology =
mesh->topology();
201 std::vector<std::array<std::int32_t, 2>> bc_dofs;
202 for (
int c = 0; c < topology->connectivity(tdim, 0)->num_nodes(); ++c)
205 auto cell_dofs0 = dofmap0->cell_dofs(c);
206 auto cell_dofs1 = dofmap1->cell_dofs(c);
209 for (std::size_t i = 0; i < cell_dofs1.size(); ++i)
211 if (marked_dofs[cell_dofs1[i]])
214 for (
int k = 0; k < element_bs; ++k)
216 const int local_pos = element_bs * i + k;
217 const std::div_t pos0 = std::div(local_pos, bs0);
218 const std::div_t pos1 = std::div(local_pos, bs1);
219 const std::int32_t dof_index0
220 = bs0 * cell_dofs0[pos0.quot] + pos0.rem;
221 const std::int32_t dof_index1
222 = bs1 * cell_dofs1[pos1.quot] + pos1.rem;
223 bc_dofs.push_back({dof_index0, dof_index1});
230 std::ranges::sort(bc_dofs);
231 auto [unique_end, range_end] = std::ranges::unique(bc_dofs);
232 bc_dofs.erase(unique_end, range_end);
235 std::array dofs = {std::vector<std::int32_t>(bc_dofs.size()),
236 std::vector<std::int32_t>(bc_dofs.size())};
237 std::ranges::transform(bc_dofs, dofs[0].begin(),
238 [](
auto dof) {
return dof[0]; });
239 std::ranges::transform(bc_dofs, dofs[1].begin(),
240 [](
auto dof) {
return dof[1]; });
253template <dolfinx::scalar T, std::
floating_po
int U = dolfinx::scalar_value_t<T>>
259 std::size_t num_owned(
const DofMap& dofmap,
260 std::span<const std::int32_t> dofs)
263 std::int32_t map_size = dofmap.
index_map->size_local();
264 std::int32_t owned_size = bs * map_size;
265 auto it = std::ranges::lower_bound(dofs, owned_size);
266 return std::distance(dofs.begin(), it);
270 static std::vector<std::int32_t>
271 unroll_dofs(std::span<const std::int32_t> dofs,
int bs)
273 std::vector<std::int32_t> dofs_unrolled(bs * dofs.size());
274 for (std::size_t i = 0; i < dofs.size(); ++i)
275 for (
int k = 0; k < bs; ++k)
276 dofs_unrolled[bs * i + k] = bs * dofs[i] + k;
277 return dofs_unrolled;
298 template <
typename S,
typename X>
299 requires(std::is_convertible_v<S, T>
300 || std::is_convertible_v<S, std::span<const T>>)
301 && std::is_convertible_v<std::remove_cvref_t<X>,
302 std::vector<std::int32_t>>
323 template <
typename X>
324 requires std::is_convertible_v<std::remove_cvref_t<X>,
325 std::vector<std::int32_t>>
328 : _function_space(V), _g(g), _dofs0(std::forward<X>(dofs)),
329 _owned_indices0(num_owned(*V->dofmap(), _dofs0))
333 if (g->shape.size() != V->element()->value_shape().size())
335 throw std::runtime_error(
336 "Rank mismatch between Constant and function space in DirichletBC");
339 if (g->value.size() != (std::size_t)_function_space->dofmap()->bs())
341 throw std::runtime_error(
342 "Creating a DirichletBC using a Constant is not supported when the "
343 "Constant size is not equal to the block size of the constrained "
344 "(sub-)space. Use a fem::Function to create the fem::DirichletBC.");
347 if (!V->element()->interpolation_ident())
349 throw std::runtime_error(
350 "Constant can be used only with point-evaluation elements");
354 if (
const int bs = V->dofmap()->bs(); bs > 1)
356 _owned_indices0 *= bs;
357 _dofs0 = unroll_dofs(_dofs0, bs);
373 template <
typename X>
374 requires std::is_convertible_v<std::remove_cvref_t<X>,
375 std::vector<std::int32_t>>
378 _dofs0(std::forward<X>(dofs)),
379 _owned_indices0(num_owned(*_function_space->dofmap(), _dofs0))
381 assert(_function_space);
384 if (
const int bs = _function_space->dofmap()->bs(); bs > 1)
386 _owned_indices0 *= bs;
387 _dofs0 = unroll_dofs(_dofs0, bs);
412 template <
typename X>
415 : _function_space(V), _g(g),
416 _dofs0(std::forward<typename std::remove_reference_t<X>::value_type>(
418 _dofs1_g(std::forward<typename std::remove_reference_t<X>::value_type>(
420 _owned_indices0(num_owned(*_function_space->dofmap(), _dofs0))
446 return _function_space;
451 std::variant<std::shared_ptr<const Function<T, U>>,
452 std::shared_ptr<const Constant<T>>>
464 std::pair<std::span<const std::int32_t>, std::int32_t>
dof_indices()
const
466 return {_dofs0, _owned_indices0};
491 void set(std::span<T> x, std::optional<std::span<const T>> x0,
496 auto apply = [&](std::invocable<std::int32_t>
auto set_fn)
499 std::is_same_v<std::invoke_result_t<
decltype(set_fn), std::int32_t>,
502 std::int32_t x_size = x.size();
503 for (std::size_t i = 0; i < _dofs0.size(); ++i)
505 if (_dofs0[i] < x_size)
506 x[_dofs0[i]] = set_fn(i);
512 apply([](std::int32_t) -> T {
return 0; });
516 if (std::holds_alternative<std::shared_ptr<
const Function<T, U>>>(_g))
518 auto g = std::get<std::shared_ptr<const Function<T, U>>>(_g);
520 std::span<const T> values = g->x()->array();
525 auto dofs_g = _dofs1_g.empty() ? std::span(_dofs0) : std::span(_dofs1_g);
529 assert(x.size() <= x0->size());
531 [dofs_g, x0 = *x0, alpha, values,
532 &dofs0 = this->_dofs0](std::int32_t i) -> T
534 assert(dofs_g[i] <
static_cast<std::int32_t
>(values.size()));
535 return alpha * (values[dofs_g[i]] - x0[dofs0[i]]);
541 [dofs_g, values, alpha](std::int32_t i) -> T
543 assert(dofs_g[i] <
static_cast<std::int32_t
>(values.size()));
544 return alpha * values[dofs_g[i]];
548 else if (std::holds_alternative<std::shared_ptr<
const Constant<T>>>(_g))
550 auto g = std::get<std::shared_ptr<const Constant<T>>>(_g);
551 const std::vector<T>&
value = g->value;
552 std::int32_t bs = _function_space->dofmap()->bs();
555 assert(x.size() <= x0->size());
557 [x0 = *x0, alpha, bs, &
value, &dofs0 = _dofs0](std::int32_t i) -> T
560 return alpha * (
value[dof % bs] - x0[dof]);
565 apply([alpha, bs, &
value, &dofs0 = _dofs0](std::int32_t i) -> T
566 {
return alpha *
value[dofs0[i] % bs]; });
587 for (std::int32_t idx : _dofs0)
589 assert(idx < (std::int32_t)markers.size());
596 std::shared_ptr<const FunctionSpace<U>> _function_space;
599 std::variant<std::shared_ptr<const Function<T, U>>,
600 std::shared_ptr<const Constant<T>>>
605 std::vector<std::int32_t> _dofs0, _dofs1_g;
608 std::int32_t _owned_indices0 = -1;
Degree-of-freedom map representations and tools.
Constant value which can be attached to a Form.
Definition Constant.h:23
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
DirichletBC & operator=(const DirichletBC &bc)=default
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
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:491
std::variant< std::shared_ptr< const Function< T, U > >, std::shared_ptr< const Constant< T > > > value() const
Definition DirichletBC.h:453
DirichletBC(const DirichletBC &bc)=default
DirichletBC(DirichletBC &&bc)=default
std::pair< std::span< const std::int32_t >, std::int32_t > dof_indices() const
Definition DirichletBC.h:464
std::shared_ptr< const FunctionSpace< U > > function_space() const
Definition DirichletBC.h:444
~DirichletBC()=default
Destructor.
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:376
DirichletBC & operator=(DirichletBC &&bc)=default
Move assignment operator.
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:585
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:413
Degree-of-freedom map.
Definition DofMap.h:73
std::shared_ptr< const common::IndexMap > index_map
Index map that describes the parallel distribution of the dofmap.
Definition DofMap.h:164
int index_map_bs() const
Block size associated with the index_map.
Definition DofMap.cpp:276
This class represents a finite element function space defined by a mesh, a finite element,...
Definition FunctionSpace.h:34
std::shared_ptr< const DofMap > dofmap() const
The dofmap.
Definition FunctionSpace.h:361
std::shared_ptr< const mesh::Mesh< geometry_type > > mesh() const
The mesh.
Definition FunctionSpace.h:336
std::shared_ptr< const FiniteElement< geometry_type > > element() const
The finite element.
Definition FunctionSpace.h:342
std::vector< geometry_type > tabulate_dof_coordinates(bool transpose) const
Tabulate the physical coordinates of all dofs on this process.
Definition FunctionSpace.h:208
Finite element method functionality.
Definition assemble_expression_impl.h:23
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
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:187
Mesh data structures and algorithms on meshes.
Definition DofMap.h:32