DOLFINx 0.9.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
DirichletBC.h
1// Copyright (C) 2007-2021 Michal Habera, Anders Logg, Garth N. Wells
2// and Jørgen S.Dokken
3//
4// This file is part of DOLFINx (https://www.fenicsproject.org)
5//
6// SPDX-License-Identifier: LGPL-3.0-or-later
7
8#pragma once
9
10#include "Constant.h"
11#include "DofMap.h"
12#include "Function.h"
13#include "FunctionSpace.h"
14#include <algorithm>
15#include <array>
16#include <concepts>
17#include <dolfinx/common/types.h>
18#include <functional>
19#include <memory>
20#include <optional>
21#include <span>
22#include <type_traits>
23#include <utility>
24#include <variant>
25#include <vector>
26
27namespace dolfinx::fem
28{
29
55std::vector<std::int32_t>
56locate_dofs_topological(const mesh::Topology& topology, const DofMap& dofmap,
57 int dim, std::span<const std::int32_t> entities,
58 bool remote = true);
59
86std::array<std::vector<std::int32_t>, 2> locate_dofs_topological(
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);
90
102template <std::floating_point T, typename U>
103std::vector<std::int32_t> locate_dofs_geometrical(const FunctionSpace<T>& V,
104 U marker_fn)
105{
106 // FIXME: Calling V.tabulate_dof_coordinates() is very expensive,
107 // especially when we usually want the boundary dofs only. Add
108 // interface that computes dofs coordinates only for specified cell.
109
110 assert(V.element());
111 if (V.element()->is_mixed())
112 {
113 throw std::runtime_error(
114 "Cannot locate dofs geometrically for mixed space. Use subspaces.");
115 }
116
117 // Compute dof coordinates
118 const std::vector<T> dof_coordinates = V.tabulate_dof_coordinates(true);
119
120 using cmdspan3x_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
121 const T,
122 MDSPAN_IMPL_STANDARD_NAMESPACE::extents<
123 std::size_t, 3, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent>>;
124
125 // Compute marker for each dof coordinate
126 cmdspan3x_t x(dof_coordinates.data(), 3, dof_coordinates.size() / 3);
127 const std::vector<std::int8_t> marked_dofs = marker_fn(x);
128
129 std::vector<std::int32_t> dofs;
130 dofs.reserve(std::count(marked_dofs.begin(), marked_dofs.end(), true));
131 for (std::size_t i = 0; i < marked_dofs.size(); ++i)
132 {
133 if (marked_dofs[i])
134 dofs.push_back(i);
135 }
136
137 return dofs;
138}
139
153template <std::floating_point T, typename U>
154std::array<std::vector<std::int32_t>, 2> locate_dofs_geometrical(
155 const std::array<std::reference_wrapper<const FunctionSpace<T>>, 2>& V,
156 U marker_fn)
157{
158 // FIXME: Calling V.tabulate_dof_coordinates() is very expensive,
159 // especially when we usually want the boundary dofs only. Add
160 // interface that computes dofs coordinates only for specified cell.
161
162 // Get function spaces
163 const FunctionSpace<T>& V0 = V.at(0).get();
164 const FunctionSpace<T>& V1 = V.at(1).get();
165
166 // Get mesh
167 auto mesh = V0.mesh();
168 assert(mesh);
169 assert(V1.mesh());
170 if (mesh != V1.mesh())
171 throw std::runtime_error("Meshes are not the same.");
172 const int tdim = mesh->topology()->dim();
173
174 assert(V0.element());
175 assert(V1.element());
176 if (*V0.element() != *V1.element())
177 throw std::runtime_error("Function spaces must have the same element.");
178
179 // Compute dof coordinates
180 const std::vector<T> dof_coordinates = V1.tabulate_dof_coordinates(true);
181
182 using cmdspan3x_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
183 const T,
184 MDSPAN_IMPL_STANDARD_NAMESPACE::extents<
185 std::size_t, 3, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent>>;
186
187 // Evaluate marker for each dof coordinate
188 cmdspan3x_t x(dof_coordinates.data(), 3, dof_coordinates.size() / 3);
189 const std::vector<std::int8_t> marked_dofs = marker_fn(x);
190
191 // Get dofmaps
192 std::shared_ptr<const DofMap> dofmap0 = V0.dofmap();
193 assert(dofmap0);
194 const int bs0 = dofmap0->bs();
195 std::shared_ptr<const DofMap> dofmap1 = V1.dofmap();
196 assert(dofmap1);
197 const int bs1 = dofmap1->bs();
198
199 const int element_bs = dofmap0->element_dof_layout().block_size();
200 assert(element_bs == dofmap1->element_dof_layout().block_size());
201
202 // Iterate over cells
203 auto topology = mesh->topology();
204 assert(topology);
205 std::vector<std::array<std::int32_t, 2>> bc_dofs;
206 for (int c = 0; c < topology->connectivity(tdim, 0)->num_nodes(); ++c)
207 {
208 // Get cell dofmaps
209 auto cell_dofs0 = dofmap0->cell_dofs(c);
210 auto cell_dofs1 = dofmap1->cell_dofs(c);
211
212 // Loop over cell dofs and add to bc_dofs if marked.
213 for (std::size_t i = 0; i < cell_dofs1.size(); ++i)
214 {
215 if (marked_dofs[cell_dofs1[i]])
216 {
217 // Unroll over blocks
218 for (int k = 0; k < element_bs; ++k)
219 {
220 const int local_pos = element_bs * i + k;
221 const std::div_t pos0 = std::div(local_pos, bs0);
222 const std::div_t pos1 = std::div(local_pos, bs1);
223 const std::int32_t dof_index0
224 = bs0 * cell_dofs0[pos0.quot] + pos0.rem;
225 const std::int32_t dof_index1
226 = bs1 * cell_dofs1[pos1.quot] + pos1.rem;
227 bc_dofs.push_back({dof_index0, dof_index1});
228 }
229 }
230 }
231 }
232
233 // Remove duplicates
234 std::ranges::sort(bc_dofs);
235 auto [unique_end, range_end] = std::ranges::unique(bc_dofs);
236 bc_dofs.erase(unique_end, range_end);
237
238 // Copy to separate array
239 std::array dofs = {std::vector<std::int32_t>(bc_dofs.size()),
240 std::vector<std::int32_t>(bc_dofs.size())};
241 std::ranges::transform(bc_dofs, dofs[0].begin(),
242 [](auto dof) { return dof[0]; });
243 std::ranges::transform(bc_dofs, dofs[1].begin(),
244 [](auto dof) { return dof[1]; });
245
246 return dofs;
247}
248
257template <dolfinx::scalar T,
258 std::floating_point U = dolfinx::scalar_value_type_t<T>>
259class DirichletBC
260{
261private:
264 std::size_t num_owned(const DofMap& dofmap,
265 std::span<const std::int32_t> dofs)
266 {
267 int bs = dofmap.index_map_bs();
268 std::int32_t map_size = dofmap.index_map->size_local();
269 std::int32_t owned_size = bs * map_size;
270 auto it = std::ranges::lower_bound(dofs, owned_size);
271 return std::distance(dofs.begin(), it);
272 }
273
275 static std::vector<std::int32_t>
276 unroll_dofs(std::span<const std::int32_t> dofs, int bs)
277 {
278 std::vector<std::int32_t> dofs_unrolled(bs * dofs.size());
279 for (std::size_t i = 0; i < dofs.size(); ++i)
280 for (int k = 0; k < bs; ++k)
281 dofs_unrolled[bs * i + k] = bs * dofs[i] + k;
282 return dofs_unrolled;
283 }
284
285public:
303 template <typename S, typename X,
304 typename
305 = std::enable_if_t<std::is_convertible_v<S, T>
306 or std::is_convertible_v<S, std::span<const T>>>>
307 requires std::is_convertible_v<std::remove_cvref_t<X>,
308 std::vector<std::int32_t>>
309 DirichletBC(const S& g, X&& dofs, std::shared_ptr<const FunctionSpace<U>> V)
310 : DirichletBC(std::make_shared<Constant<T>>(g), dofs, V)
311 {
312 }
313
329 template <typename X>
330 requires std::is_convertible_v<std::remove_cvref_t<X>,
331 std::vector<std::int32_t>>
332 DirichletBC(std::shared_ptr<const Constant<T>> g, X&& dofs,
333 std::shared_ptr<const FunctionSpace<U>> V)
334 : _function_space(V), _g(g), _dofs0(std::forward<X>(dofs)),
335 _owned_indices0(num_owned(*V->dofmap(), _dofs0))
336 {
337 assert(g);
338 assert(V);
339 if (g->shape.size() != V->value_shape().size())
340 {
341 throw std::runtime_error(
342 "Rank mis-match between Constant and function space in DirichletBC");
343 }
344
345 if (g->value.size() != _function_space->dofmap()->bs())
346 {
347 throw std::runtime_error(
348 "Creating a DirichletBC using a Constant is not supported when the "
349 "Constant size is not equal to the block size of the constrained "
350 "(sub-)space. Use a fem::Function to create the fem::DirichletBC.");
351 }
352
353 if (!V->element()->interpolation_ident())
354 {
355 throw std::runtime_error(
356 "Constant can be used only with point-evaluation elements");
357 }
358
359 // Unroll _dofs0 if dofmap block size > 1
360 if (const int bs = V->dofmap()->bs(); bs > 1)
361 {
362 _owned_indices0 *= bs;
363 _dofs0 = unroll_dofs(_dofs0, bs);
364 }
365 }
366
379 template <typename X>
380 requires std::is_convertible_v<std::remove_cvref_t<X>,
381 std::vector<std::int32_t>>
382 DirichletBC(std::shared_ptr<const Function<T, U>> g, X&& dofs)
383 : _function_space(g->function_space()), _g(g),
384 _dofs0(std::forward<X>(dofs)),
385 _owned_indices0(num_owned(*_function_space->dofmap(), _dofs0))
386 {
387 assert(_function_space);
388
389 // Unroll _dofs0 if dofmap block size > 1
390 if (const int bs = _function_space->dofmap()->bs(); bs > 1)
391 {
392 _owned_indices0 *= bs;
393 _dofs0 = unroll_dofs(_dofs0, bs);
394 }
395 }
396
418 template <typename X>
419 DirichletBC(std::shared_ptr<const Function<T, U>> g, X&& V_g_dofs,
420 std::shared_ptr<const FunctionSpace<U>> V)
421 : _function_space(V), _g(g),
422 _dofs0(std::forward<typename X::value_type>(V_g_dofs[0])),
423 _dofs1_g(std::forward<typename X::value_type>(V_g_dofs[1])),
424 _owned_indices0(num_owned(*_function_space->dofmap(), _dofs0))
425 {
426 }
427
430 DirichletBC(const DirichletBC& bc) = default;
431
434 DirichletBC(DirichletBC&& bc) = default;
435
437 ~DirichletBC() = default;
438
441 DirichletBC& operator=(const DirichletBC& bc) = default;
442
445
448 std::shared_ptr<const FunctionSpace<U>> function_space() const
449 {
450 return _function_space;
451 }
452
455 std::variant<std::shared_ptr<const Function<T, U>>,
456 std::shared_ptr<const Constant<T>>>
457 value() const
458 {
459 return _g;
460 }
461
468 std::pair<std::span<const std::int32_t>, std::int32_t> dof_indices() const
469 {
470 return {_dofs0, _owned_indices0};
471 }
472
495 void set(std::span<T> x, std::optional<std::span<const T>> x0,
496 T alpha = 1) const
497 {
498 std::int32_t x_size = x.size();
499 if (alpha == T(0)) // Optimisation for when alpha == 0
500 {
501 for (std::int32_t idx : _dofs0)
502 {
503 if (idx < x_size)
504 x[idx] = 0;
505 }
506 }
507 else
508 {
509 if (std::holds_alternative<std::shared_ptr<const Function<T, U>>>(_g))
510 {
511 auto g = std::get<std::shared_ptr<const Function<T, U>>>(_g);
512 assert(g);
513 auto dofs1_g
514 = _dofs1_g.empty() ? std::span(_dofs0) : std::span(_dofs1_g);
515 std::span<const T> values = g->x()->array();
516 if (x0.has_value())
517 {
518 std::span<const T> _x0 = x0.value();
519 assert(x.size() <= _x0.size());
520 for (std::size_t i = 0; i < _dofs0.size(); ++i)
521 {
522 if (_dofs0[i] < x_size)
523 {
524 assert(dofs1_g[i] < (std::int32_t)values.size());
525 x[_dofs0[i]] = alpha * (values[dofs1_g[i]] - _x0[_dofs0[i]]);
526 }
527 }
528 }
529 else
530 {
531 for (std::size_t i = 0; i < _dofs0.size(); ++i)
532 {
533 if (_dofs0[i] < x_size)
534 {
535 assert(dofs1_g[i] < (std::int32_t)values.size());
536 x[_dofs0[i]] = alpha * values[dofs1_g[i]];
537 }
538 }
539 }
540 }
541 else if (std::holds_alternative<std::shared_ptr<const Constant<T>>>(_g))
542 {
543 auto g = std::get<std::shared_ptr<const Constant<T>>>(_g);
544 const std::vector<T>& value = g->value;
545 std::int32_t bs = _function_space->dofmap()->bs();
546 if (x0.has_value())
547 {
548 assert(x.size() <= x0.value().size());
549 std::ranges::for_each(
550 _dofs0,
551 [x_size, &x, x0 = x0.value(), &value, alpha, bs](auto dof)
552 {
553 if (dof < x_size)
554 x[dof] = alpha * (value[dof % bs] - x0[dof]);
555 });
556 }
557 else
558 {
559 std::ranges::for_each(_dofs0,
560 [x_size, bs, alpha, &value, &x](auto dof)
561 {
562 if (dof < x_size)
563 x[dof] = alpha * value[dof % bs];
564 });
565 }
566 }
567 }
568 }
569
579 void mark_dofs(std::span<std::int8_t> markers) const
580 {
581 for (std::int32_t idx : _dofs0)
582 {
583 assert(idx < (std::int32_t)markers.size());
584 markers[idx] = true;
585 }
586 }
587
588private:
589 // The function space (possibly a sub function space)
590 std::shared_ptr<const FunctionSpace<U>> _function_space;
591
592 // The function
593 std::variant<std::shared_ptr<const Function<T, U>>,
594 std::shared_ptr<const Constant<T>>>
595 _g;
596
597 // Dof indices (_dofs0) in _function_space and (_dofs1_g) in the space
598 // of _g. _dofs1_g may be empty if _dofs0 can be re-used
599 std::vector<std::int32_t> _dofs0, _dofs1_g;
600
601 // The first _owned_indices in _dofs are owned by this process
602 std::int32_t _owned_indices0 = -1;
603};
604} // namespace dolfinx::fem
Degree-of-freedom map representations and tools.
Constant value which can be attached to a Form.
Definition Form.h:29
Definition petsc.h:32
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:332
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:495
std::variant< std::shared_ptr< const Function< T, U > >, std::shared_ptr< const Constant< T > > > value() const
Definition DirichletBC.h:457
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:468
std::shared_ptr< const FunctionSpace< U > > function_space() const
Definition DirichletBC.h:448
~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:382
DirichletBC & operator=(DirichletBC &&bc)=default
Move assignment operator.
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:309
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:579
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:419
This class represents a finite element function space defined by a mesh, a finite element,...
Definition vtk_utils.h:32
std::shared_ptr< const DofMap > dofmap() const
The dofmap.
Definition FunctionSpace.h:335
std::shared_ptr< const mesh::Mesh< geometry_type > > mesh() const
The mesh.
Definition FunctionSpace.h:323
std::shared_ptr< const FiniteElement< geometry_type > > element() const
The finite element.
Definition FunctionSpace.h:329
std::vector< geometry_type > tabulate_dof_coordinates(bool transpose) const
Tabulate the physical coordinates of all dofs on this process.
Definition FunctionSpace.h:190
Definition XDMFFile.h:29
Definition types.h:20
Finite element method functionality.
Definition assemble_matrix_impl.h:26
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