version 0.4.4
utils/dirichletvalues.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers ikarus@ibb.uni-stuttgart.de
2// SPDX-License-Identifier: LGPL-3.0-or-later
3
10#pragma once
11#include <algorithm>
12#include <concepts>
13#include <cstddef>
14#include <functional>
15#include <type_traits>
16#include <utility>
17#include <vector>
18
19#include <dune/common/float_cmp.hh>
20#include <dune/common/hybridutilities.hh>
21#include <dune/common/indices.hh>
22#include <dune/functions/backends/istlvectorbackend.hh>
23#include <dune/functions/functionspacebases/boundarydofs.hh>
24#include <dune/functions/functionspacebases/flatmultiindex.hh>
25#include <dune/functions/functionspacebases/interpolate.hh>
26#include <dune/functions/functionspacebases/subspacebasis.hh>
27
28#include <Eigen/Core>
29
30#include <autodiff/forward/real/real.hpp>
31
33
34namespace Dune {
35template <class K, int N>
37}
38
39namespace Ikarus {
40
44template <class>
46
47template <Concepts::EigenVector T>
49{
50 using SizeType = Eigen::Index;
51};
52
53template <>
54struct DeriveSizeType<std::vector<bool>>
55{
56 using SizeType = std::vector<bool>::size_type;
57};
58
59namespace Impl {
60 template <typename Tree>
61 struct PreBasisInfo
62 {
63 };
64
65 template <typename Tree>
66 requires(Tree::isLeaf)
67 struct PreBasisInfo<Tree>
68 {
69 static constexpr std::size_t size = 0;
70 using NodalSolutionType = double;
71 };
72
73 template <typename Tree>
74 requires(Tree::isPower)
75 struct PreBasisInfo<Tree>
76 {
77 static constexpr std::size_t size = Tree::degree();
78 using NodalSolutionType = Eigen::Vector<double, size>;
79 };
80
81 template <typename Tree>
82 requires(Tree::isComposite)
83 struct PreBasisInfo<Tree>
84 {
85 using ChildTreeType = Tree::template Child<0>::Type;
86 static_assert(not ChildTreeType::isComposite,
87 "DirichletValues is not implemented to handle a composite basis within a composite basis.");
88
89 static constexpr std::size_t size = PreBasisInfo<ChildTreeType>::size;
90 using NodalSolutionType = PreBasisInfo<ChildTreeType>::NodalSolutionType;
91 };
92} // namespace Impl
93
106template <typename B, typename FC = std::vector<bool>>
108{
109public:
110 using Basis = std::remove_cvref_t<B>;
111 using FlagsType = FC;
112 using BackendType = decltype(Dune::Functions::istlVectorBackend(std::declval<FlagsType&>()));
114 static constexpr int worldDimension = Basis::GridView::dimensionworld;
115 using LocalView = Basis::LocalView;
116 using Tree = LocalView::Tree;
117 using NodalSolutionType = typename Impl::PreBasisInfo<Tree>::NodalSolutionType;
118 static constexpr std::size_t numberOfChildrenAtNode = Impl::PreBasisInfo<Tree>::size;
119 explicit DirichletValues(const B& basis)
120 : basis_{basis},
121 dirichletFlagsBackend_{dirichletFlags_} {
122 dirichletFlagsBackend_.resize(basis_);
123 std::fill(dirichletFlags_.begin(), dirichletFlags_.end(), false);
124 if constexpr (Tree::isComposite) {
125 Dune::Hybrid::forEach(
126 Dune::Hybrid::integralRange(Dune::index_constant<std::tuple_size_v<typename Basis::PreBasis::SubPreBases>>()),
127 [&](const auto i) {
128 static_assert(not Tree::template Child<i>::Type::isComposite,
129 "DirichletValues is not implemented to handle a composite basis within a composite basis.");
130 });
131 }
132 }
133
143 template <typename F, typename TreePath = Dune::TypeTree::HybridTreePath<>>
144 void fixBoundaryDOFs(F&& f, TreePath&& treePath = {}) {
145 using namespace Dune::Functions;
146 using SubSpaceLocalView =
147 typename std::remove_cvref_t<decltype(subspaceBasis(basis_, std::forward<TreePath>(treePath)))>::LocalView;
148
150 auto lambda = [&](auto&& indexGlobal) { f(dirichletFlagsBackend_, indexGlobal); };
151 Dune::Functions::forEachBoundaryDOF(subspaceBasis(basis_, std::forward<TreePath>(treePath)), lambda);
152 } else if constexpr (Concepts::IsFunctorWithArgs<F, BackendType, int, SubSpaceLocalView>) {
153 auto lambda = [&](auto&& localIndex, auto&& localView) { f(dirichletFlagsBackend_, localIndex, localView); };
154 Dune::Functions::forEachBoundaryDOF(subspaceBasis(basis_, std::forward<TreePath>(treePath)), lambda);
155 } else if constexpr (Concepts::IsFunctorWithArgs<F, BackendType, int, SubSpaceLocalView,
156 typename Basis::GridView::Intersection>) {
157 auto lambda = [&](auto&& localIndex, auto&& localView, auto&& intersection) {
158 f(dirichletFlagsBackend_, localIndex, localView, intersection);
159 };
160 Dune::Functions::forEachBoundaryDOF(subspaceBasis(basis_, std::forward<TreePath>(treePath)), lambda);
161 } else {
162 static_assert(Dune::AlwaysFalse<F>(), "fixBoundaryDOFs: A function with this signature is not supported");
163 }
164 }
165
174 template <typename F>
175 void fixDOFs(F&& f) {
176 f(basis_, dirichletFlagsBackend_);
177 }
178
185 template <typename MultiIndex>
186 requires(not std::integral<MultiIndex>)
187 void setSingleDOF(const MultiIndex i, bool flag) {
188 dirichletFlagsBackend_[i] = flag;
189 }
190
198 void setSingleDOF(std::size_t i, bool flag)
199 requires(std::same_as<typename Basis::MultiIndex, Dune::Functions::FlatMultiIndex<size_t>>)
200 {
201 dirichletFlags_[i] = flag;
202 }
203
207 void reset() { std::fill(dirichletFlags_.begin(), dirichletFlags_.end(), false); }
208
209 /* \brief Returns the local basis object */
210 const auto& basis() const { return basis_; }
211
212 /* \brief Returns a boolean values, if the give multiIndex is constrained */
213 template <typename MultiIndex>
214 requires(not std::integral<MultiIndex>)
215 [[nodiscard]] bool isConstrained(const MultiIndex& multiIndex) const {
216 return dirichletFlagsBackend_[multiIndex];
217 }
218
219 /* \brief Returns a boolean values, if the i-th degree of freedom is constrained */
220 [[nodiscard]] bool isConstrained(std::size_t i) const
221 requires(std::same_as<typename Basis::MultiIndex, Dune::Functions::FlatMultiIndex<size_t>>)
222 {
223 return dirichletFlags_[i];
224 }
225
226 /* \brief Returns how many degrees of freedoms are fixed */
227 auto fixedDOFsize() const { return std::ranges::count(dirichletFlags_, true); }
228
229 /* \brief Returns how many degrees of freedom there are */
230 auto size() const { return dirichletFlags_.size(); }
231
232 /* \brief Returns the underlying dof flag container */
233 auto& container() const { return dirichletFlags_; }
234
246 template <typename F>
247 void storeInhomogeneousBoundaryCondition(F&& f, double lambda = 1.0) {
248 auto derivativeLambda = [&](const auto& globalCoord, const double& lambda) {
249 autodiff::real lambdaDual = lambda;
250 lambdaDual[1] = 1; // Setting the derivative in lambda direction to 1
251 return derivative(f(globalCoord, lambdaDual));
252 };
253 dirichletFunctions_.push_back({f, derivativeLambda});
254 setInhomogeneousBoundaryConditionFlag(lambda);
255 }
256
262 void setZeroAtConstrainedDofs(Eigen::VectorXd& xIh) const {
263 for (Eigen::Index i = 0; i < xIh.size(); ++i)
264 if (this->isConstrained(i))
265 xIh[i] = 0.0;
266 }
267
277 void evaluateInhomogeneousBoundaryCondition(Eigen::VectorXd& xIh, const double& lambda) const {
278 Eigen::VectorXd inhomogeneousBoundaryVectorDummy;
279 inhomogeneousBoundaryVectorDummy.setZero(this->size());
280 xIh.resizeLike(inhomogeneousBoundaryVectorDummy);
281 xIh.setZero();
282 auto runInterpolateImpl = [&](const auto& basis) {
283 for (auto& f : dirichletFunctions_) {
284 interpolate(basis, inhomogeneousBoundaryVectorDummy,
285 [&](const auto& globalCoord) { return f.value(globalCoord, lambda); });
286 xIh += inhomogeneousBoundaryVectorDummy;
287 }
288 };
289 runInterpolate(runInterpolateImpl);
290 }
291
301 void evaluateInhomogeneousBoundaryConditionDerivative(Eigen::VectorXd& xIh, const double& lambda) const {
302 Eigen::VectorXd inhomogeneousBoundaryVectorDummy;
303 inhomogeneousBoundaryVectorDummy.setZero(this->size());
304 xIh.resizeLike(inhomogeneousBoundaryVectorDummy);
305 xIh.setZero();
306 auto runInterpolateImpl = [&](const auto& basis) {
307 for (auto& f : dirichletFunctions_) {
308 interpolate(basis, inhomogeneousBoundaryVectorDummy,
309 [&](const auto& globalCoord) { return f.derivative(globalCoord, lambda); });
310 xIh += inhomogeneousBoundaryVectorDummy;
311 }
312 };
313 runInterpolate(runInterpolateImpl);
314 }
315
316private:
317 Basis basis_;
318 FlagsType dirichletFlags_;
319 BackendType dirichletFlagsBackend_;
320 struct DirichletFunctions
321 {
322 using Signature = std::function<NodalSolutionType(const Dune::FieldVector<double, worldDimension>&, const double&)>;
323 Signature value;
324 Signature derivative;
325 };
326 std::vector<DirichletFunctions> dirichletFunctions_;
327
328 void setInhomogeneousBoundaryConditionFlag(double lambda) {
329 Eigen::VectorXd inhomogeneousBoundaryVectorDummy;
330 inhomogeneousBoundaryVectorDummy.setZero(this->size());
331 this->evaluateInhomogeneousBoundaryCondition(inhomogeneousBoundaryVectorDummy, lambda);
332 for (const std::size_t i : Dune::range(this->size()))
333 if (Dune::FloatCmp::ne(inhomogeneousBoundaryVectorDummy[i], 0.0))
334 this->setSingleDOF(i, true);
335 }
336
337 template <typename F>
338 void runInterpolate(F&& f) const {
339 if constexpr (Tree::isComposite)
340 f(subspaceBasis(basis_, Dune::Indices::_0));
341 else
342 f(basis_);
343 }
344};
345
346} // namespace Ikarus
Definition: assemblermanipulatorbuildingblocks.hh:22
Definition: utils/dirichletvalues.hh:34
Definition: utils/dirichletvalues.hh:36
A helper struct to derive the SizeType of the underlying container.
Definition: utils/dirichletvalues.hh:45
Eigen::Index SizeType
Definition: utils/dirichletvalues.hh:50
std::vector< bool >::size_type SizeType
Definition: utils/dirichletvalues.hh:56
Class for handling Dirichlet boundary conditions in Ikarus.
Definition: utils/dirichletvalues.hh:108
auto size() const
Definition: utils/dirichletvalues.hh:230
void setZeroAtConstrainedDofs(Eigen::VectorXd &xIh) const
Function to write zeros at constrained Dirichlet entries.
Definition: utils/dirichletvalues.hh:262
void setSingleDOF(const MultiIndex i, bool flag)
Fixes and unfixes (set boolean value to true or false) a specific degree of freedom.
Definition: utils/dirichletvalues.hh:187
bool isConstrained(const MultiIndex &multiIndex) const
Definition: utils/dirichletvalues.hh:215
bool isConstrained(std::size_t i) const
Definition: utils/dirichletvalues.hh:220
void setSingleDOF(std::size_t i, bool flag)
Fixes or unfixes (set boolean value to true or false) a specific degree of freedom.
Definition: utils/dirichletvalues.hh:198
std::remove_cvref_t< B > Basis
Definition: utils/dirichletvalues.hh:110
void fixBoundaryDOFs(F &&f, TreePath &&treePath={})
Function to fix (set boolean values to true or false) degrees of freedom on the boundary.
Definition: utils/dirichletvalues.hh:144
decltype(Dune::Functions::istlVectorBackend(std::declval< FlagsType & >())) BackendType
Definition: utils/dirichletvalues.hh:112
void evaluateInhomogeneousBoundaryConditionDerivative(Eigen::VectorXd &xIh, const double &lambda) const
Function to evaluate all stored inhomogeneous Dirichlet boundary derivative functions.
Definition: utils/dirichletvalues.hh:301
void storeInhomogeneousBoundaryCondition(F &&f, double lambda=1.0)
Function to insert a function for inhomogeneous Dirichlet boundary conditions.
Definition: utils/dirichletvalues.hh:247
LocalView::Tree Tree
Definition: utils/dirichletvalues.hh:116
static constexpr int worldDimension
Definition: utils/dirichletvalues.hh:114
void fixDOFs(F &&f)
Function to fix (set boolean values to true or false) degrees of freedom.
Definition: utils/dirichletvalues.hh:175
FC FlagsType
Definition: utils/dirichletvalues.hh:111
void evaluateInhomogeneousBoundaryCondition(Eigen::VectorXd &xIh, const double &lambda) const
Function to evaluate all stored inhomogeneous Dirichlet boundary functions.
Definition: utils/dirichletvalues.hh:277
Basis::LocalView LocalView
Definition: utils/dirichletvalues.hh:115
void reset()
Resets all degrees of freedom.
Definition: utils/dirichletvalues.hh:207
typename DeriveSizeType< FlagsType >::SizeType SizeType
Definition: utils/dirichletvalues.hh:113
const auto & basis() const
Definition: utils/dirichletvalues.hh:210
auto fixedDOFsize() const
Definition: utils/dirichletvalues.hh:227
DirichletValues(const B &basis)
Definition: utils/dirichletvalues.hh:119
static constexpr std::size_t numberOfChildrenAtNode
Definition: utils/dirichletvalues.hh:118
auto & container() const
Definition: utils/dirichletvalues.hh:233
typename Impl::PreBasisInfo< Tree >::NodalSolutionType NodalSolutionType
Definition: utils/dirichletvalues.hh:117
Concept defining the requirements for functors with arguments.
Definition: utils/concepts.hh:357
Several concepts.