DOLFINx
DOLFINx C++ interface
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 "Function.h"
12#include "FunctionSpace.h"
13#include <array>
14#include <functional>
15#include <memory>
16#include <span>
17#include <type_traits>
18#include <utility>
19#include <variant>
20#include <vector>
21#include <xtensor/xtensor.hpp>
22
23namespace dolfinx::fem
24{
25
47std::vector<std::int32_t>
48locate_dofs_topological(const FunctionSpace& V, int dim,
49 const std::span<const std::int32_t>& entities,
50 bool remote = true);
51
75std::array<std::vector<std::int32_t>, 2> locate_dofs_topological(
76 const std::array<std::reference_wrapper<const FunctionSpace>, 2>& V,
77 int dim, const std::span<const std::int32_t>& entities, bool remote = true);
78
90std::vector<std::int32_t> locate_dofs_geometrical(
91 const FunctionSpace& V,
92 const std::function<xt::xtensor<bool, 1>(const xt::xtensor<double, 2>&)>&
93 marker_fn);
94
108std::array<std::vector<std::int32_t>, 2> locate_dofs_geometrical(
109 const std::array<std::reference_wrapper<const FunctionSpace>, 2>& V,
110 const std::function<xt::xtensor<bool, 1>(const xt::xtensor<double, 2>&)>&
111 marker_fn);
112
123template <typename T>
125{
126private:
129 std::size_t num_owned(const FunctionSpace& V,
130 std::span<const std::int32_t> dofs)
131 {
132 int bs = V.dofmap()->index_map_bs();
133 std::int32_t map_size = V.dofmap()->index_map->size_local();
134 std::int32_t owned_size = bs * map_size;
135 auto it = std::lower_bound(dofs.begin(), dofs.end(), owned_size);
136 return std::distance(dofs.begin(), it);
137 }
138
140 std::vector<std::int32_t> unroll_dofs(std::span<const std::int32_t> dofs,
141 int bs)
142 {
143 std::vector<std::int32_t> dofs_unrolled(bs * dofs.size());
144 for (std::size_t i = 0; i < dofs.size(); ++i)
145 for (int k = 0; k < bs; ++k)
146 dofs_unrolled[bs * i + k] = bs * dofs[i] + k;
147 return dofs_unrolled;
148 }
149
150public:
169 template <typename S, typename U,
170 typename = std::enable_if_t<
171 std::is_convertible_v<
172 S, T> or std::is_convertible_v<S, std::span<const T>>>>
173 DirichletBC(const S& g, U&& dofs,
174 const std::shared_ptr<const FunctionSpace>& V)
175 : DirichletBC(std::make_shared<Constant<T>>(g), dofs, V)
176 {
177 }
178
195 template <typename U>
196 DirichletBC(const std::shared_ptr<const Constant<T>>& g, U&& dofs,
197 const std::shared_ptr<const FunctionSpace>& V)
198 : _function_space(V), _g(g), _dofs0(std::forward<U>(dofs)),
199 _owned_indices0(num_owned(*V, _dofs0))
200 {
201 assert(g);
202 assert(V);
203 if (g->shape.size() != V->element()->value_shape().size())
204 {
205 throw std::runtime_error(
206 "Rank mis-match between Constant and function space in DirichletBC");
207 }
208
209 if (g->value.size() != _function_space->dofmap()->bs())
210 {
211 throw std::runtime_error(
212 "Creating a DirichletBC using a Constant is not supported when the "
213 "Constant size is not equal to the block size of the constrained "
214 "(sub-)space. Use a fem::Function to create the fem::DirichletBC.");
215 }
216
217 if (!V->element()->interpolation_ident())
218 {
219 throw std::runtime_error(
220 "Constant can be used only with point-evaluation elements");
221 }
222
223 // Unroll _dofs0 if dofmap block size > 1
224 if (const int bs = V->dofmap()->bs(); bs > 1)
225 {
226 _owned_indices0 *= bs;
227 _dofs0 = unroll_dofs(_dofs0, bs);
228 }
229 }
230
244 template <typename U>
245 DirichletBC(const std::shared_ptr<const Function<T>>& g, U&& dofs)
246 : _function_space(g->function_space()), _g(g),
247 _dofs0(std::forward<U>(dofs)),
248 _owned_indices0(num_owned(*_function_space, _dofs0))
249 {
250 assert(_function_space);
251
252 // Unroll _dofs0 if dofmap block size > 1
253 if (const int bs = _function_space->dofmap()->bs(); bs > 1)
254 {
255 _owned_indices0 *= bs;
256 _dofs0 = unroll_dofs(_dofs0, bs);
257 }
258 }
259
281 template <typename U>
282 DirichletBC(const std::shared_ptr<const Function<T>>& g, U&& V_g_dofs,
283 const std::shared_ptr<const FunctionSpace>& V)
284 : _function_space(V), _g(g),
285 _dofs0(std::forward<typename U::value_type>(V_g_dofs[0])),
286 _dofs1_g(std::forward<typename U::value_type>(V_g_dofs[1])),
287 _owned_indices0(num_owned(*_function_space, _dofs0))
288 {
289 }
290
293 DirichletBC(const DirichletBC& bc) = default;
294
297 DirichletBC(DirichletBC&& bc) = default;
298
300 ~DirichletBC() = default;
301
304 DirichletBC& operator=(const DirichletBC& bc) = default;
305
308
311 std::shared_ptr<const FunctionSpace> function_space() const
312 {
313 return _function_space;
314 }
315
318 std::variant<std::shared_ptr<const Function<T>>,
319 std::shared_ptr<const Constant<T>>>
320 value() const
321 {
322 return _g;
323 }
324
331 std::pair<std::span<const std::int32_t>, std::int32_t> dof_indices() const
332 {
333 return {_dofs0, _owned_indices0};
334 }
335
347 void set(std::span<T> x, double scale = 1.0) const
348 {
349 if (std::holds_alternative<std::shared_ptr<const Function<T>>>(_g))
350 {
351 auto g = std::get<std::shared_ptr<const Function<T>>>(_g);
352 assert(g);
353 std::span<const T> values = g->x()->array();
354 auto dofs1_g = _dofs1_g.empty() ? std::span(_dofs0) : std::span(_dofs1_g);
355 std::int32_t x_size = x.size();
356 for (std::size_t i = 0; i < _dofs0.size(); ++i)
357 {
358 if (_dofs0[i] < x_size)
359 {
360 assert(dofs1_g[i] < (std::int32_t)values.size());
361 x[_dofs0[i]] = scale * values[dofs1_g[i]];
362 }
363 }
364 }
365 else if (std::holds_alternative<std::shared_ptr<const Constant<T>>>(_g))
366 {
367 auto g = std::get<std::shared_ptr<const Constant<T>>>(_g);
368 std::vector<T> value = g->value;
369 int bs = _function_space->dofmap()->bs();
370 std::int32_t x_size = x.size();
371 std::for_each(_dofs0.cbegin(), _dofs0.cend(),
372 [x_size, bs, scale, &value, &x](auto dof)
373 {
374 if (dof < x_size)
375 x[dof] = scale * value[dof % bs];
376 });
377 }
378 }
379
384 void set(std::span<T> x, const std::span<const T>& x0,
385 double scale = 1.0) const
386 {
387 if (std::holds_alternative<std::shared_ptr<const Function<T>>>(_g))
388 {
389 auto g = std::get<std::shared_ptr<const Function<T>>>(_g);
390 assert(g);
391 std::span<const T> values = g->x()->array();
392 assert(x.size() <= x0.size());
393 auto dofs1_g = _dofs1_g.empty() ? std::span(_dofs0) : std::span(_dofs1_g);
394 std::int32_t x_size = x.size();
395 for (std::size_t i = 0; i < _dofs0.size(); ++i)
396 {
397 if (_dofs0[i] < x_size)
398 {
399 assert(dofs1_g[i] < (std::int32_t)values.size());
400 x[_dofs0[i]] = scale * (values[dofs1_g[i]] - x0[_dofs0[i]]);
401 }
402 }
403 }
404 else if (std::holds_alternative<std::shared_ptr<const Constant<T>>>(_g))
405 {
406 auto g = std::get<std::shared_ptr<const Constant<T>>>(_g);
407 const std::vector<T>& value = g->value;
408 std::int32_t bs = _function_space->dofmap()->bs();
409 std::for_each(_dofs0.cbegin(), _dofs0.cend(),
410 [&x, &x0, &value, scale, bs](auto dof)
411 {
412 if (dof < (std::int32_t)x.size())
413 x[dof] = scale * (value[dof % bs] - x0[dof]);
414 });
415 }
416 }
417
426 void dof_values(std::span<T> values) const
427 {
428 if (std::holds_alternative<std::shared_ptr<const Function<T>>>(_g))
429 {
430 auto g = std::get<std::shared_ptr<const Function<T>>>(_g);
431 assert(g);
432 std::span<const T> g_values = g->x()->array();
433 auto dofs1_g = _dofs1_g.empty() ? std::span(_dofs0) : std::span(_dofs1_g);
434 for (std::size_t i = 0; i < dofs1_g.size(); ++i)
435 values[_dofs0[i]] = g_values[dofs1_g[i]];
436 }
437 else if (std::holds_alternative<std::shared_ptr<const Constant<T>>>(_g))
438 {
439 auto g = std::get<std::shared_ptr<const Constant<T>>>(_g);
440 assert(g);
441 const std::vector<T>& g_value = g->value;
442 const std::int32_t bs = _function_space->dofmap()->bs();
443 for (std::size_t i = 0; i < _dofs0.size(); ++i)
444 values[_dofs0[i]] = g_value[_dofs0[i] % bs];
445 }
446 }
447
454 void mark_dofs(const std::span<std::int8_t>& markers) const
455 {
456 for (std::size_t i = 0; i < _dofs0.size(); ++i)
457 {
458 assert(_dofs0[i] < (std::int32_t)markers.size());
459 markers[_dofs0[i]] = true;
460 }
461 }
462
463private:
464 // The function space (possibly a sub function space)
465 std::shared_ptr<const FunctionSpace> _function_space;
466
467 // The function
468 std::variant<std::shared_ptr<const Function<T>>,
469 std::shared_ptr<const Constant<T>>>
470 _g;
471
472 // Dof indices (_dofs0) in _function_space and (_dofs1_g) in the
473 // space of _g. _dofs1_g may be empty if _dofs0 can be re-used
474 std::vector<std::int32_t> _dofs0, _dofs1_g;
475
476 // The first _owned_indices in _dofs are owned by this process
477 std::int32_t _owned_indices0 = -1;
478};
479} // namespace dolfinx::fem
Constant value which can be attached to a Form. Constants may be scalar (rank 0), vector (rank 1),...
Definition: Constant.h:20
Object for setting (strong) Dirichlet boundary conditions.
Definition: DirichletBC.h:125
DirichletBC(const std::shared_ptr< const Function< T > > &g, U &&dofs)
Create a representation of a Dirichlet boundary condition where the space being constrained is the sa...
Definition: DirichletBC.h:245
void dof_values(std::span< T > values) const
Definition: DirichletBC.h:426
void set(std::span< T > x, double scale=1.0) const
Set bc entries in x to scale * x_bc
Definition: DirichletBC.h:347
void set(std::span< T > x, const std::span< const T > &x0, double scale=1.0) const
Set bc entries in x to scale * (x0 - x_bc)
Definition: DirichletBC.h:384
DirichletBC & operator=(const DirichletBC &bc)=default
Assignment operator.
DirichletBC(const DirichletBC &bc)=default
Copy constructor.
DirichletBC(DirichletBC &&bc)=default
Move constructor.
std::variant< std::shared_ptr< const Function< T > >, std::shared_ptr< const Constant< T > > > value() const
Return boundary value function g.
Definition: DirichletBC.h:320
std::pair< std::span< const std::int32_t >, std::int32_t > dof_indices() const
Access dof indices (local indices, unrolled), including ghosts, to which a Dirichlet condition is app...
Definition: DirichletBC.h:331
DirichletBC(const std::shared_ptr< const Constant< T > > &g, U &&dofs, const std::shared_ptr< const FunctionSpace > &V)
Create a representation of a Dirichlet boundary condition constrained by a fem::Constant.
Definition: DirichletBC.h:196
std::shared_ptr< const FunctionSpace > function_space() const
The function space to which boundary conditions are applied.
Definition: DirichletBC.h:311
DirichletBC(const std::shared_ptr< const Function< T > > &g, U &&V_g_dofs, const std::shared_ptr< const FunctionSpace > &V)
Create a representation of a Dirichlet boundary condition where the space being constrained and the f...
Definition: DirichletBC.h:282
~DirichletBC()=default
Destructor.
DirichletBC & operator=(DirichletBC &&bc)=default
Move assignment operator.
DirichletBC(const S &g, U &&dofs, const std::shared_ptr< const FunctionSpace > &V)
Create a representation of a Dirichlet boundary condition constrained by a scalar- or vector-valued c...
Definition: DirichletBC.h:173
void mark_dofs(const std::span< std::int8_t > &markers) const
Set markers[i] = true if dof i has a boundary condition applied. Value of markers[i] is not changed o...
Definition: DirichletBC.h:454
This class represents a finite element function space defined by a mesh, a finite element,...
Definition: FunctionSpace.h:31
std::shared_ptr< const DofMap > dofmap() const
The dofmap.
Definition: FunctionSpace.cpp:253
This class represents a function in a finite element function space , given by.
Definition: Function.h:45
Finite element method functionality.
Definition: assemble_matrix_impl.h:25
std::vector< std::int32_t > locate_dofs_topological(const FunctionSpace &V, int dim, const std::span< const std::int32_t > &entities, bool remote=true)
Find degrees-of-freedom which belong to the provided mesh entities (topological). Note that degrees-o...
Definition: DirichletBC.cpp:182
std::vector< std::int32_t > locate_dofs_geometrical(const FunctionSpace &V, const std::function< xt::xtensor< bool, 1 >(const xt::xtensor< double, 2 > &)> &marker_fn)
Finds degrees of freedom whose geometric coordinate is true for the provided marking function.
Definition: DirichletBC.cpp:452