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 <type_traits>
17#include <utility>
18#include <variant>
19#include <vector>
20#include <xtensor/xtensor.hpp>
21#include <xtl/xspan.hpp>
22
23namespace dolfinx::fem
24{
25
47std::vector<std::int32_t>
48locate_dofs_topological(const FunctionSpace& V, int dim,
49 const xtl::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 xtl::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:
127 template <typename U>
128 DirichletBC(const std::variant<std::shared_ptr<const Function<T>>,
129 std::shared_ptr<const Constant<T>>>& g,
130 U&& dofs, const std::shared_ptr<const FunctionSpace>& V, void*)
131 : _function_space(V), _g(g), _dofs0(std::forward<U>(dofs))
132 {
133 assert(_function_space);
134 if (auto c = std::get_if<std::shared_ptr<const Constant<T>>>(&_g))
135 {
136 assert(*c);
137 if ((*c)->value.size() != _function_space->dofmap()->bs())
138 {
139 throw std::runtime_error(
140 "Creating a DirichletBC using a Constant is not supported when the "
141 "Constant size is not equal to the block size of the constrained "
142 "(sub-)space. Use a Function to create the DirichletBC.");
143 }
144 }
145
146 // Compute number of owned dofs indices in the full space (will
147 // contain 'gaps' for sub-spaces)
148 const int map0_bs = _function_space->dofmap()->index_map_bs();
149 const int map0_size = _function_space->dofmap()->index_map->size_local();
150 const int owned_size0 = map0_bs * map0_size;
151
152 // Find number of owned indices in _dofs0
153 auto it0 = std::lower_bound(_dofs0.begin(), _dofs0.end(), owned_size0);
154 _owned_indices0 = std::distance(_dofs0.begin(), it0);
155
156 // Unroll _dofs0 for dofmap block size > 1
157 if (const int bs = _function_space->dofmap()->bs(); bs > 1)
158 {
159 _owned_indices0 *= bs;
160 const std::vector<std::int32_t> dof_tmp = _dofs0;
161 _dofs0.resize(bs * dof_tmp.size());
162 for (std::size_t i = 0; i < dof_tmp.size(); ++i)
163 {
164 for (int k = 0; k < bs; ++k)
165 _dofs0[bs * i + k] = bs * dof_tmp[i] + k;
166 }
167 }
168 }
169
170public:
186 template <typename S, typename U,
187 typename = std::enable_if_t<
188 std::is_convertible_v<
189 S, T> or std::is_convertible_v<S, xt::xarray<T>>>>
190 DirichletBC(const S& g, U&& dofs,
191 const std::shared_ptr<const FunctionSpace>& V)
192 : DirichletBC(std::make_shared<Constant<T>>(g), dofs, V)
193 {
194 }
195
211 template <typename U>
212 DirichletBC(const std::shared_ptr<const Constant<T>>& g, U&& dofs,
213 const std::shared_ptr<const FunctionSpace>& V)
214 : DirichletBC(g, dofs, V, nullptr)
215 {
216 assert(g);
217 assert(V);
218 if (V->element()->value_shape().size() != g->shape.size())
219 {
220 throw std::runtime_error(
221 "Rank mis-match between Constant and function space in DirichletBC");
222 }
223
224 if (!V->element()->interpolation_ident())
225 {
226 throw std::runtime_error(
227 "Constant can be used only with point-evaluation elements");
228 }
229 }
230
243 template <typename U>
244 DirichletBC(const std::shared_ptr<const Function<T>>& g, U&& dofs)
245 : DirichletBC(g, dofs, g->function_space(), nullptr)
246 {
247 }
248
268 template <typename U>
269 DirichletBC(const std::shared_ptr<const Function<T>>& g, U&& V_g_dofs,
270 const std::shared_ptr<const FunctionSpace>& V)
271 : _function_space(V), _g(g),
272 _dofs0(std::forward<typename U::value_type>(V_g_dofs[0])),
273 _dofs1_g(std::forward<typename U::value_type>(V_g_dofs[1]))
274 {
275 assert(_dofs0.size() == _dofs1_g.size());
276 assert(_function_space);
277 const int map0_bs = _function_space->dofmap()->index_map_bs();
278 const int map0_size = _function_space->dofmap()->index_map->size_local();
279 const int owned_size0 = map0_bs * map0_size;
280 auto it0 = std::lower_bound(_dofs0.begin(), _dofs0.end(), owned_size0);
281 _owned_indices0 = std::distance(_dofs0.begin(), it0);
282 }
283
286 DirichletBC(const DirichletBC& bc) = default;
287
290 DirichletBC(DirichletBC&& bc) = default;
291
293 ~DirichletBC() = default;
294
297 DirichletBC& operator=(const DirichletBC& bc) = default;
298
301
304 std::shared_ptr<const fem::FunctionSpace> function_space() const
305 {
306 return _function_space;
307 }
308
311 std::variant<std::shared_ptr<const Function<T>>,
312 std::shared_ptr<const Constant<T>>>
313 value() const
314 {
315 return _g;
316 }
317
324 std::pair<xtl::span<const std::int32_t>, std::int32_t> dof_indices() const
325 {
326 return {_dofs0, _owned_indices0};
327 }
328
340 void set(xtl::span<T> x, double scale = 1.0) const
341 {
342 if (std::holds_alternative<std::shared_ptr<const Function<T>>>(_g))
343 {
344 auto g = std::get<std::shared_ptr<const Function<T>>>(_g);
345 assert(g);
346 xtl::span<const T> values = g->x()->array();
347 auto dofs1_g = _dofs1_g.empty() ? xtl::span(_dofs0) : xtl::span(_dofs1_g);
348 std::int32_t x_size = x.size();
349 for (std::size_t i = 0; i < _dofs0.size(); ++i)
350 {
351 if (_dofs0[i] < x_size)
352 {
353 assert(dofs1_g[i] < (std::int32_t)values.size());
354 x[_dofs0[i]] = scale * values[dofs1_g[i]];
355 }
356 }
357 }
358 else if (std::holds_alternative<std::shared_ptr<const Constant<T>>>(_g))
359 {
360 auto g = std::get<std::shared_ptr<const Constant<T>>>(_g);
361 std::vector<T> value = g->value;
362 int bs = _function_space->dofmap()->bs();
363 std::int32_t x_size = x.size();
364 std::for_each(_dofs0.cbegin(), _dofs0.cend(),
365 [x_size, bs, scale, &value, &x](auto dof)
366 {
367 if (dof < x_size)
368 x[dof] = scale * value[dof % bs];
369 });
370 }
371 }
372
377 void set(xtl::span<T> x, const xtl::span<const T>& x0,
378 double scale = 1.0) const
379 {
380 if (std::holds_alternative<std::shared_ptr<const Function<T>>>(_g))
381 {
382 auto g = std::get<std::shared_ptr<const Function<T>>>(_g);
383 assert(g);
384 xtl::span<const T> values = g->x()->array();
385 assert(x.size() <= x0.size());
386 auto dofs1_g = _dofs1_g.empty() ? xtl::span(_dofs0) : xtl::span(_dofs1_g);
387 std::int32_t x_size = x.size();
388 for (std::size_t i = 0; i < _dofs0.size(); ++i)
389 {
390 if (_dofs0[i] < x_size)
391 {
392 assert(dofs1_g[i] < (std::int32_t)values.size());
393 x[_dofs0[i]] = scale * (values[dofs1_g[i]] - x0[_dofs0[i]]);
394 }
395 }
396 }
397 else if (std::holds_alternative<std::shared_ptr<const Constant<T>>>(_g))
398 {
399 auto g = std::get<std::shared_ptr<const Constant<T>>>(_g);
400 const std::vector<T>& value = g->value;
401 std::int32_t bs = _function_space->dofmap()->bs();
402 std::for_each(_dofs0.cbegin(), _dofs0.cend(),
403 [&x, &x0, &value, scale, bs](auto dof)
404 {
405 if (dof < (std::int32_t)x.size())
406 x[dof] = scale * (value[dof % bs] - x0[dof]);
407 });
408 }
409 }
410
419 void dof_values(xtl::span<T> values) const
420 {
421 if (std::holds_alternative<std::shared_ptr<const Function<T>>>(_g))
422 {
423 auto g = std::get<std::shared_ptr<const Function<T>>>(_g);
424 assert(g);
425 xtl::span<const T> g_values = g->x()->array();
426 auto dofs1_g = _dofs1_g.empty() ? xtl::span(_dofs0) : xtl::span(_dofs1_g);
427 for (std::size_t i = 0; i < dofs1_g.size(); ++i)
428 values[_dofs0[i]] = g_values[dofs1_g[i]];
429 }
430 else if (std::holds_alternative<std::shared_ptr<const Constant<T>>>(_g))
431 {
432 auto g = std::get<std::shared_ptr<const Constant<T>>>(_g);
433 assert(g);
434 const std::vector<T>& g_value = g->value;
435 const std::int32_t bs = _function_space->dofmap()->bs();
436 for (std::size_t i = 0; i < _dofs0.size(); ++i)
437 values[_dofs0[i]] = g_value[_dofs0[i] % bs];
438 }
439 }
440
447 void mark_dofs(const xtl::span<std::int8_t>& markers) const
448 {
449 for (std::size_t i = 0; i < _dofs0.size(); ++i)
450 {
451 assert(_dofs0[i] < (std::int32_t)markers.size());
452 markers[_dofs0[i]] = true;
453 }
454 }
455
456private:
457 // The function space (possibly a sub function space)
458 std::shared_ptr<const FunctionSpace> _function_space;
459
460 // The function
461 std::variant<std::shared_ptr<const Function<T>>,
462 std::shared_ptr<const Constant<T>>>
463 _g;
464
465 // Dof indices (_dofs0) in _function_space and (_dofs1_g) in the
466 // space of _g. _dofs1_g may be empty if _dofs0 can be re-used
467 std::vector<std::int32_t> _dofs0, _dofs1_g;
468
469 // The first _owned_indices in _dofs are owned by this process
470 int _owned_indices0 = -1;
471 int _owned_indices1 = -1;
472};
473} // namespace dolfinx::fem
Constant value which can be attached to a Form. Constants may be scalar (rank 0), vector (rank 1),...
Definition: Constant.h:21
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:244
void set(xtl::span< T > x, double scale=1.0) const
Set bc entries in x to scale * x_bc
Definition: DirichletBC.h:340
DirichletBC & operator=(const DirichletBC &bc)=default
Assignment operator.
std::pair< xtl::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:324
DirichletBC(const DirichletBC &bc)=default
Copy constructor.
DirichletBC(DirichletBC &&bc)=default
Move constructor.
void dof_values(xtl::span< T > values) const
Definition: DirichletBC.h:419
std::variant< std::shared_ptr< const Function< T > >, std::shared_ptr< const Constant< T > > > value() const
Return boundary value function g.
Definition: DirichletBC.h:313
void mark_dofs(const xtl::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:447
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:212
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:269
~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 tensor constant.
Definition: DirichletBC.h:190
std::shared_ptr< const fem::FunctionSpace > function_space() const
The function space to which boundary conditions are applied.
Definition: DirichletBC.h:304
void set(xtl::span< T > x, const xtl::span< const T > &x0, double scale=1.0) const
Set bc entries in x to scale * (x0 - x_bc)
Definition: DirichletBC.h:377
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:24
std::vector< std::int32_t > locate_dofs_topological(const FunctionSpace &V, int dim, const xtl::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:205
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:437