Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3057-gb5b616642a 2025-04-10 20:30:00+00:00
\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}} \newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=} \newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]} \newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
dof_accessor.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1998 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_dof_accessor_h
16#define dealii_dof_accessor_h
17
18
19#include <deal.II/base/config.h>
20
21#include <deal.II/base/types.h>
22
26
28
32
33#include <boost/container/small_vector.hpp>
34
35#include <limits>
36#include <set>
37#include <type_traits>
38#include <vector>
39
40
42
43// Forward declarations
44#ifndef DOXYGEN
45template <typename number>
46class FullMatrix;
47template <typename number>
48class Vector;
49template <typename number>
51
52template <typename Accessor>
53class TriaRawIterator;
54
55template <int, int>
56class FiniteElement;
57
58namespace internal
59{
60 namespace DoFCellAccessorImplementation
61 {
62 struct Implementation;
63 }
64
65 namespace DoFHandlerImplementation
66 {
67 struct Implementation;
68 namespace Policy
69 {
70 struct Implementation;
71 }
72 } // namespace DoFHandlerImplementation
73
74 namespace hp
75 {
76 namespace DoFHandlerImplementation
77 {
78 struct Implementation;
79 }
80 } // namespace hp
81} // namespace internal
82#endif
83
84
85namespace internal
86{
87 namespace DoFAccessorImplementation
88 {
105 template <int structdim, int dim, int spacedim>
114
115
121 template <int dim, int spacedim>
122 struct Inheritance<dim, dim, spacedim>
123 {
129 };
130
131 struct Implementation;
132 } // namespace DoFAccessorImplementation
133} // namespace internal
134
135
136/* -------------------------------------------------------------------------- */
137
138
139
212template <int structdim, int dim, int spacedim, bool level_dof_access>
213class DoFAccessor : public ::internal::DoFAccessorImplementation::
214 Inheritance<structdim, dim, spacedim>::BaseClass
215{
216public:
221 static constexpr unsigned int dimension = dim;
222
227 static constexpr unsigned int space_dimension = spacedim;
228
233 using BaseClass = typename ::internal::DoFAccessorImplementation::
234 Inheritance<structdim, dimension, space_dimension>::BaseClass;
235
240
252
269 const int level,
270 const int index,
272
277 default;
278
282 DoFAccessor( // NOLINT
284 default; // NOLINT
285
289 ~DoFAccessor() = default;
290
303 template <int structdim2, int dim2, int spacedim2>
305
310 template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
313
317 template <bool level_dof_access2>
319
331 delete;
332
337 operator=( // NOLINT
339 default; // NOLINT
340
350
354 template <bool level_dof_access2>
355 void
357
362 void
364
369 static bool
371
383 child(const unsigned int c) const;
384
390 typename ::internal::DoFHandlerImplementation::
391 Iterators<dim, spacedim, level_dof_access>::line_iterator
392 line(const unsigned int i) const;
393
399 typename ::internal::DoFHandlerImplementation::
400 Iterators<dim, spacedim, level_dof_access>::quad_iterator
401 quad(const unsigned int i) const;
402
448 void
450 std::vector<types::global_dof_index> &dof_indices,
451 const types::fe_index fe_index = numbers::invalid_fe_index) const;
452
459 void
461 const int level,
462 std::vector<types::global_dof_index> &dof_indices,
463 const types::fe_index fe_index = numbers::invalid_fe_index) const;
464
468 void
470 const int level,
471 const std::vector<types::global_dof_index> &dof_indices,
473
497 const unsigned int vertex,
498 const unsigned int i,
499 const types::fe_index fe_index = numbers::invalid_fe_index) const;
500
508 const int level,
509 const unsigned int vertex,
510 const unsigned int i,
511 const types::fe_index fe_index = numbers::invalid_fe_index) const;
512
541 dof_index(const unsigned int i,
542 const types::fe_index fe_index = numbers::invalid_fe_index) const;
543
548 mg_dof_index(const int level, const unsigned int i) const;
549
571 unsigned int
573
582 nth_active_fe_index(const unsigned int n) const;
583
590 std::set<types::fe_index>
592
602 bool
603 fe_index_is_active(const types::fe_index fe_index) const;
604
611 get_fe(const types::fe_index fe_index) const;
612
623 "This accessor object has not been "
624 "associated with any DoFHandler object.");
657
658protected:
663
664public:
678 template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
679 bool
682
686 template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
687 bool
690
691protected:
695 void
697
716 void
718 const unsigned int i,
719 const types::global_dof_index index,
720 const types::fe_index fe_index = numbers::invalid_fe_index) const;
721
722 void
724 const unsigned int i,
725 const types::global_dof_index index) const;
726
727 void
729 const int level,
730 const unsigned int vertex,
731 const unsigned int i,
732 const types::global_dof_index index,
733 const types::fe_index fe_index = numbers::invalid_fe_index) const;
734
735 // Iterator classes need to be friends because they need to access
736 // operator== and operator!=.
737 template <typename>
738 friend class TriaRawIterator;
739 template <int, int, int, bool>
740 friend class DoFAccessor;
741
742private:
743 // Make the DoFHandler class a friend so that it can call the set_xxx()
744 // functions.
745 friend class DoFHandler<dim, spacedim>;
746
747 friend struct ::internal::DoFHandlerImplementation::Policy::
748 Implementation;
749 friend struct ::internal::DoFHandlerImplementation::Implementation;
750 friend struct ::internal::hp::DoFHandlerImplementation::Implementation;
751 friend struct ::internal::DoFCellAccessorImplementation::Implementation;
752 friend struct ::internal::DoFAccessorImplementation::Implementation;
753};
754
755
756
764template <int spacedim, bool level_dof_access>
765class DoFAccessor<0, 1, spacedim, level_dof_access>
766 : public TriaAccessor<0, 1, spacedim>
767{
768public:
773 static constexpr unsigned int dimension = 1;
774
779 static constexpr unsigned int space_dimension = spacedim;
780
786
791
802 DoFAccessor();
803
821 const Triangulation<1, spacedim> *tria,
822 const typename TriaAccessor<0, 1, spacedim>::VertexKind vertex_kind,
823 const unsigned int vertex_index,
825
833 const int = 0,
834 const int = 0,
836
849 template <int structdim2, int dim2, int spacedim2>
851
856 template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
859
864
868 // NOLINTNEXTLINE OSX does not compile with noexcept
870
874 ~DoFAccessor() = default;
875
887
893 default;
894
902 const DoFHandler<1, spacedim> &
903 get_dof_handler() const;
904
908 template <bool level_dof_access2>
909 void
910 copy_from(const DoFAccessor<0, 1, spacedim, level_dof_access2> &a);
911
916 void
917 copy_from(const TriaAccessorBase<0, 1, spacedim> &da);
918
931 TriaIterator<DoFAccessor<0, 1, spacedim, level_dof_access>>
932 child(const unsigned int c) const;
933
940 typename ::internal::DoFHandlerImplementation::
941 Iterators<1, spacedim, level_dof_access>::line_iterator
942 line(const unsigned int i) const;
943
950 typename ::internal::DoFHandlerImplementation::
951 Iterators<1, spacedim, level_dof_access>::quad_iterator
952 quad(const unsigned int i) const;
953
997 void
999 std::vector<types::global_dof_index> &dof_indices,
1000 const types::fe_index fe_index = numbers::invalid_fe_index) const;
1001
1008 void
1010 const int level,
1011 std::vector<types::global_dof_index> &dof_indices,
1012 const types::fe_index fe_index = numbers::invalid_fe_index) const;
1013
1032 types::global_dof_index
1034 const unsigned int vertex,
1035 const unsigned int i,
1036 const types::fe_index fe_index = numbers::invalid_fe_index) const;
1037
1055 types::global_dof_index
1056 dof_index(const unsigned int i,
1057 const types::fe_index fe_index = numbers::invalid_fe_index) const;
1058
1077 unsigned int
1078 n_active_fe_indices() const;
1079
1087 types::fe_index
1088 nth_active_fe_index(const unsigned int n) const;
1089
1098 bool
1099 fe_index_is_active(const types::fe_index fe_index) const;
1100
1106 const FiniteElement<1, spacedim> &
1107 get_fe(const types::fe_index fe_index) const;
1108
1119 "This accessor object has not been "
1120 "associated with any DoFHandler object.");
1153
1154protected:
1159
1163 template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
1164 bool
1165 operator==(
1166 const DoFAccessor<structdim2, dim2, spacedim2, level_dof_access2> &) const;
1167
1171 template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
1172 bool
1173 operator!=(
1174 const DoFAccessor<structdim2, dim2, spacedim2, level_dof_access2> &) const;
1175
1179 void
1180 set_dof_handler(DoFHandler<1, spacedim> *dh);
1181
1200 void
1202 const unsigned int i,
1203 const types::global_dof_index index,
1204 const types::fe_index fe_index = numbers::invalid_fe_index) const;
1205
1206 // Iterator classes need to be friends because they need to access
1207 // operator== and operator!=.
1208 template <typename>
1209 friend class TriaRawIterator;
1210
1211
1212 // Make the DoFHandler class a friend so that it can call the set_xxx()
1213 // functions.
1214 friend class DoFHandler<1, spacedim>;
1215
1216 friend struct ::internal::DoFHandlerImplementation::Policy::
1218 friend struct ::internal::DoFHandlerImplementation::Implementation;
1219 friend struct ::internal::hp::DoFHandlerImplementation::Implementation;
1220 friend struct ::internal::DoFCellAccessorImplementation::Implementation;
1221};
1222
1223
1224
1225/* -------------------------------------------------------------------------- */
1226
1227
1248template <int structdim, int dim, int spacedim = dim>
1249class DoFInvalidAccessor : public InvalidAccessor<structdim, dim, spacedim>
1250{
1251public:
1257
1265 DoFInvalidAccessor(const void *parent = nullptr,
1266 const int level = -1,
1267 const int index = -1,
1268 const AccessorData *local_data = nullptr);
1269
1278
1283 template <typename OtherAccessor>
1284 DoFInvalidAccessor(const OtherAccessor &);
1285
1292 dof_index(const unsigned int i,
1293 const types::fe_index fe_index =
1295
1301 void
1303 const unsigned int i,
1304 const types::global_dof_index index,
1305 const types::fe_index fe_index = numbers::invalid_fe_index) const;
1306};
1307
1308
1309
1310/* -------------------------------------------------------------------------- */
1311
1312
1324template <int dimension_, int space_dimension_, bool level_dof_access>
1325class DoFCellAccessor : public DoFAccessor<dimension_,
1326 dimension_,
1327 space_dimension_,
1328 level_dof_access>
1329{
1330public:
1334 static const unsigned int dim = dimension_;
1335
1339 static const unsigned int spacedim = space_dimension_;
1340
1341
1346
1353
1358
1364 dimension_,
1365 space_dimension_,
1366 level_dof_access>>;
1367
1379 const int level,
1380 const int index,
1381 const AccessorData *local_data);
1382
1395 template <int structdim2, int dim2, int spacedim2>
1397
1402 template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
1403 explicit DoFCellAccessor(
1405
1411 default;
1412
1418 &&) = default; // NOLINT
1419
1423 ~DoFCellAccessor() = default;
1424
1437 delete;
1438
1443 operator=( // NOLINT
1445 &&) = default; // NOLINT
1446
1461 parent() const;
1462
1476 neighbor(const unsigned int i) const;
1477
1484 periodic_neighbor(const unsigned int i) const;
1485
1492 neighbor_or_periodic_neighbor(const unsigned int i) const;
1493
1500 child(const unsigned int i) const;
1501
1505 boost::container::small_vector<
1509 child_iterators() const;
1510
1518 face(const unsigned int i) const;
1519
1523 boost::container::small_vector<face_iterator,
1525 face_iterators() const;
1526
1534 neighbor_child_on_subface(const unsigned int face_no,
1535 const unsigned int subface_no) const;
1536
1544 periodic_neighbor_child_on_subface(const unsigned int face_no,
1545 const unsigned int subface_no) const;
1546
1575 template <class InputVector, typename number>
1576 void
1577 get_dof_values(const InputVector &values, Vector<number> &local_values) const;
1578
1596 template <typename Number, typename ForwardIterator>
1597 void
1598 get_dof_values(const ReadVector<Number> &values,
1599 ForwardIterator local_values_begin,
1600 ForwardIterator local_values_end) const;
1601
1622 template <class InputVector, typename ForwardIterator>
1623 void
1624 get_dof_values(
1626 const InputVector &values,
1627 ForwardIterator local_values_begin,
1628 ForwardIterator local_values_end) const;
1629
1654 template <class OutputVector, typename number>
1655 void
1656 set_dof_values(const Vector<number> &local_values,
1657 OutputVector &values) const;
1658
1690 template <typename Number>
1691 void
1692 get_interpolated_dof_values(
1693 const ReadVector<Number> &values,
1694 Vector<Number> &interpolated_values,
1695 const types::fe_index fe_index = numbers::invalid_fe_index) const;
1696
1758 template <class OutputVector, typename number>
1759 void
1760 set_dof_values_by_interpolation(
1761 const Vector<number> &local_values,
1762 OutputVector &values,
1764 const bool perform_check = false) const;
1765
1775 template <class OutputVector, typename number>
1776 void
1777 distribute_local_to_global_by_interpolation(
1778 const Vector<number> &local_values,
1779 OutputVector &values,
1780 const types::fe_index fe_index = numbers::invalid_fe_index) const;
1781
1796 template <typename number, typename OutputVector>
1797 void
1798 distribute_local_to_global(const Vector<number> &local_source,
1799 OutputVector &global_destination) const;
1800
1815 template <typename ForwardIterator, typename OutputVector>
1816 void
1817 distribute_local_to_global(ForwardIterator local_source_begin,
1818 ForwardIterator local_source_end,
1819 OutputVector &global_destination) const;
1820
1834 template <typename ForwardIterator, typename OutputVector>
1835 void
1836 distribute_local_to_global(
1838 ForwardIterator local_source_begin,
1839 ForwardIterator local_source_end,
1840 OutputVector &global_destination) const;
1841
1848 template <typename number, typename OutputMatrix>
1849 void
1850 distribute_local_to_global(const FullMatrix<number> &local_source,
1851 OutputMatrix &global_destination) const;
1852
1857 template <typename number, typename OutputMatrix, typename OutputVector>
1858 void
1859 distribute_local_to_global(const FullMatrix<number> &local_matrix,
1860 const Vector<number> &local_vector,
1861 OutputMatrix &global_matrix,
1862 OutputVector &global_vector) const;
1863
1888 void
1889 get_active_or_mg_dof_indices(
1890 std::vector<types::global_dof_index> &dof_indices) const;
1891
1923 void
1924 get_dof_indices(std::vector<types::global_dof_index> &dof_indices) const;
1925
1930 void
1931 get_mg_dof_indices(std::vector<types::global_dof_index> &dof_indices) const;
1932
1958 get_fe() const;
1959
1986 active_fe_index() const;
1987
2014 void
2015 set_active_fe_index(const types::fe_index i) const;
2024 void
2025 set_dof_indices(const std::vector<types::global_dof_index> &dof_indices);
2026
2030 void
2031 set_mg_dof_indices(const std::vector<types::global_dof_index> &dof_indices);
2032
2059 get_future_fe() const;
2060
2080 future_fe_index() const;
2081
2090 void
2091 set_future_fe_index(const types::fe_index i) const;
2092
2099 bool
2100 future_fe_index_set() const;
2101
2110 void
2111 clear_future_fe_index() const;
2116private:
2117 friend struct ::internal::DoFCellAccessorImplementation::Implementation;
2118};
2119
2120
2121template <int structdim, int dim, int spacedim, bool level_dof_access>
2122inline bool
2127
2128
2129
2130template <int structdim, int dim, int spacedim>
2131template <typename OtherAccessor>
2133 const OtherAccessor &)
2134{
2135 Assert(false,
2136 ExcMessage("You are attempting an illegal conversion between "
2137 "iterator/accessor types. The constructor you call "
2138 "only exists to make certain template constructs "
2139 "easier to write as dimension independent code but "
2140 "the conversion is not valid in the current context."));
2141}
2142
2143
2144
2145/*------------------------- Functions: DoFAccessor ---------------------------*/
2146
2147
2148template <int structdim, int dim, int spacedim, bool level_dof_access>
2153
2154
2155
2156template <int structdim, int dim, int spacedim, bool level_dof_access>
2158 const Triangulation<dim, spacedim> *tria,
2159 const int level,
2160 const int index,
2162 : ::internal::DoFAccessorImplementation::
2163 Inheritance<structdim, dim, spacedim>::BaseClass(tria, level, index)
2164 , dof_handler(const_cast<DoFHandler<dim, spacedim> *>(dof_handler))
2165{
2166 Assert(
2167 tria == nullptr || &dof_handler->get_triangulation() == tria,
2168 ExcMessage(
2169 "You can't create a DoF accessor in which the DoFHandler object "
2170 "uses a different triangulation than the one you pass as argument."));
2171}
2172
2173
2174
2175template <int structdim, int dim, int spacedim, bool level_dof_access>
2176template <int structdim2, int dim2, int spacedim2>
2182
2183
2184
2185template <int structdim, int dim, int spacedim, bool level_dof_access>
2186template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
2189 : BaseClass(other)
2190 , dof_handler(nullptr)
2191{
2192 Assert(false,
2193 ExcMessage(
2194 "You are trying to assign iterators that are incompatible. "
2195 "The reason for incompatibility is that they refer to objects of "
2196 "different dimensionality (e.g., assigning a line iterator "
2197 "to a quad iterator)."));
2198}
2199
2200
2201
2202template <int structdim, int dim, int spacedim, bool level_dof_access>
2203template <bool level_dof_access2>
2206 : BaseClass(other)
2207 , dof_handler(const_cast<DoFHandler<dim, spacedim> *>(other.dof_handler))
2208{}
2209
2210
2211
2212template <int structdim, int dim, int spacedim, bool level_dof_access>
2213inline void
2216{
2217 Assert(dh != nullptr, ExcInvalidObject());
2218 this->dof_handler = dh;
2219}
2220
2221
2222
2223template <int structdim, int dim, int spacedim, bool level_dof_access>
2224inline const DoFHandler<dim, spacedim> &
2226{
2227 Assert(this->dof_handler != nullptr, ExcInvalidObject());
2228 return *this->dof_handler;
2229}
2230
2231
2232
2233template <int structdim, int dim, int spacedim, bool level_dof_access>
2234inline void
2237{
2238 Assert(this->dof_handler != nullptr, ExcInvalidObject());
2239 BaseClass::copy_from(da);
2240}
2241
2242
2243
2244template <int structdim, int dim, int spacedim, bool level_dof_access>
2245template <bool level_dof_access2>
2246inline void
2253
2254
2255
2256template <int structdim, int dim, int spacedim, bool level_dof_access>
2257template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
2258inline bool
2261{
2262 Assert(structdim == structdim2, ExcCantCompareIterators());
2263 Assert(this->dof_handler == a.dof_handler, ExcCantCompareIterators());
2264 return (BaseClass::operator==(a));
2265}
2266
2267
2268
2269template <int structdim, int dim, int spacedim, bool level_dof_access>
2270template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
2271inline bool
2274{
2275 Assert(structdim == structdim2, ExcCantCompareIterators());
2276 Assert(this->dof_handler == a.dof_handler, ExcCantCompareIterators());
2277 return (BaseClass::operator!=(a));
2278}
2279
2280
2281
2282template <int structdim, int dim, int spacedim, bool level_dof_access>
2285 const unsigned int i) const
2286{
2287 Assert(static_cast<unsigned int>(this->present_level) <
2288 this->dof_handler->object_dof_indices.size(),
2289 ExcMessage("DoFHandler not initialized"));
2290
2293
2295 *t, this->dof_handler);
2296 return q;
2297}
2298
2299
2300namespace internal
2301{
2302 namespace DoFAccessorImplementation
2303 {
2308 template <int structdim, int dim, int spacedim, bool level_dof_access>
2312 const types::fe_index fe_index)
2313 {
2314 if (cell.get_dof_handler().has_hp_capabilities() == false)
2315 {
2316 // No hp enabled, and the argument is at its default value -> we
2317 // can translate to the default active fe index
2318 Assert(
2319 (fe_index == numbers::invalid_fe_index) ||
2321 ExcMessage(
2322 "It is not possible to specify a FE index if no hp support is used!"));
2323
2325 }
2326 else
2327 {
2328 // Otherwise: If anything other than the default is provided by
2329 // the caller, then we should take just that. As an exception, if
2330 // we are on a cell (rather than a face/edge/vertex), then we know
2331 // that there is only one active fe index on this cell and we can
2332 // use that:
2333 if ((dim == structdim) && (fe_index == numbers::invalid_fe_index))
2334 {
2336
2337 return cell.nth_active_fe_index(0);
2338 }
2339
2340 Assert((fe_index != numbers::invalid_fe_index),
2341 ExcMessage(
2342 "You need to specify a FE index if hp support is used!"));
2343
2344 return fe_index;
2345 }
2346 }
2347
2353 {
2363 boost::container::small_vector<::types::global_dof_index, 27>;
2364
2375 template <int dim,
2376 int spacedim,
2377 int structdim,
2378 typename GlobalIndexType,
2379 typename DoFPProcessor>
2380 static void
2382 const unsigned int obj_level,
2383 const unsigned int obj_index,
2384 const types::fe_index fe_index,
2385 const unsigned int local_index,
2386 const std::integral_constant<int, structdim> &,
2387 GlobalIndexType &global_index,
2388 const DoFPProcessor &process)
2389 {
2390 Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
2391
2392 // 1) no hp used -> fe_index == 0
2393 if (dof_handler.hp_capability_enabled == false)
2394 {
2395 AssertDimension(fe_index,
2397
2398 process(
2399 dof_handler.object_dof_indices
2400 [obj_level][structdim]
2401 [dof_handler.object_dof_ptr[obj_level][structdim][obj_index] +
2402 local_index],
2403 global_index);
2404
2405 return;
2406 }
2407
2408 // 2) cell and hp is used -> there is only one fe_index
2409 if (structdim == dim)
2410 {
2411 process(
2412 dof_handler.object_dof_indices
2413 [obj_level][structdim]
2414 [dof_handler.object_dof_ptr[obj_level][structdim][obj_index] +
2415 local_index],
2416 global_index);
2417 return;
2418 }
2419
2420 // 3) general entity and hp is used
2421 AssertIndexRange(obj_level, dof_handler.object_dof_indices.size());
2422 AssertIndexRange(structdim,
2423 dof_handler.object_dof_indices[obj_level].size());
2424
2426
2427 AssertIndexRange(structdim, dof_handler.hp_object_fe_ptr.size());
2428 AssertIndexRange(obj_index,
2429 dof_handler.hp_object_fe_ptr[structdim].size());
2430
2431 const auto ptr =
2432 std::find(dof_handler.hp_object_fe_indices[structdim].begin() +
2433 dof_handler.hp_object_fe_ptr[structdim][obj_index],
2434 dof_handler.hp_object_fe_indices[structdim].begin() +
2435 dof_handler.hp_object_fe_ptr[structdim][obj_index + 1],
2436 fe_index);
2437
2438 Assert(ptr != dof_handler.hp_object_fe_indices[structdim].begin() +
2439 dof_handler.hp_object_fe_ptr[structdim][obj_index + 1],
2440 ExcMessage(
2441 "You are requesting an active FE index that is not assigned "
2442 "to any of the cells connected to this entity."));
2443
2444 const types::fe_index fe_index_ =
2445 std::distance(dof_handler.hp_object_fe_indices[structdim].begin() +
2446 dof_handler.hp_object_fe_ptr[structdim][obj_index],
2447 ptr);
2448
2450 dof_handler.hp_capability_enabled ?
2451 (dof_handler.hp_object_fe_ptr[structdim][obj_index] + fe_index_) :
2452 obj_index,
2453 dof_handler.object_dof_ptr[obj_level][structdim].size());
2454
2456 dof_handler.object_dof_ptr
2457 [obj_level][structdim]
2458 [dof_handler.hp_capability_enabled ?
2459 (dof_handler.hp_object_fe_ptr[structdim][obj_index] +
2460 fe_index_) :
2461 obj_index] +
2462 local_index,
2463 dof_handler.object_dof_indices[obj_level][structdim].size());
2464
2465 process(dof_handler.object_dof_indices
2466 [obj_level][structdim]
2467 [dof_handler.object_dof_ptr
2468 [obj_level][structdim]
2469 [dof_handler.hp_capability_enabled ?
2470 (dof_handler.hp_object_fe_ptr[structdim][obj_index] +
2471 fe_index_) :
2472 obj_index] +
2473 local_index],
2474 global_index);
2475 }
2476
2486 template <int dim, int spacedim, int structdim>
2487 static std::pair<unsigned int, unsigned int>
2489 const unsigned int obj_level,
2490 const unsigned int obj_index,
2491 const types::fe_index fe_index,
2492 const std::integral_constant<int, structdim> &)
2493 {
2494 Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
2495
2496 // determine range of dofs in global data structure
2497 // 1) cell
2498 if (structdim == dim)
2499 {
2500 const unsigned int ptr_0 =
2501 dof_handler.object_dof_ptr[obj_level][structdim][obj_index];
2502 const unsigned int length =
2503 dof_handler.get_fe(fe_index).template n_dofs_per_object<dim>(0);
2504
2505 return {ptr_0, length};
2506 }
2507
2508 // 2) hp is not used -> fe_index == 0
2509 if (dof_handler.hp_capability_enabled == false)
2510 {
2511 AssertDimension(fe_index,
2513
2514 const unsigned int ptr_0 =
2515 dof_handler.object_dof_ptr[obj_level][structdim][obj_index];
2516 const unsigned int length =
2517 dof_handler.object_dof_ptr[obj_level][structdim][obj_index + 1] -
2518 ptr_0;
2519
2520 return {ptr_0, length};
2521 }
2522
2523 // 3) hp is used
2524 AssertIndexRange(obj_level, dof_handler.object_dof_indices.size());
2525 AssertIndexRange(structdim,
2526 dof_handler.object_dof_indices[obj_level].size());
2527
2528 AssertIndexRange(structdim, dof_handler.hp_object_fe_ptr.size());
2529 AssertIndexRange(obj_index,
2530 dof_handler.hp_object_fe_ptr[structdim].size());
2531
2532 const auto fe_index_local_ptr =
2533 std::find(dof_handler.hp_object_fe_indices[structdim].begin() +
2534 dof_handler.hp_object_fe_ptr[structdim][obj_index],
2535 dof_handler.hp_object_fe_indices[structdim].begin() +
2536 dof_handler.hp_object_fe_ptr[structdim][obj_index + 1],
2537 fe_index);
2538
2539 Assert(fe_index_local_ptr !=
2540 dof_handler.hp_object_fe_indices[structdim].begin() +
2541 dof_handler.hp_object_fe_ptr[structdim][obj_index + 1],
2542 ExcMessage(
2543 "You tried to call a function accessing DoF indices, but "
2544 "they appear not be available (yet) or inconsistent. "
2545 "Did you call distribute_dofs() first? Alternatively, if "
2546 "you are using different elements on different cells (i.e., "
2547 "you are using the hp capabilities of deal.II), did you "
2548 "change the active_fe_index of a cell since you last "
2549 "called distribute_dofs()?"));
2550
2551 const types::fe_index fe_index_local =
2552 std::distance(dof_handler.hp_object_fe_indices[structdim].begin() +
2553 dof_handler.hp_object_fe_ptr[structdim][obj_index],
2554 fe_index_local_ptr);
2555
2557 dof_handler.hp_object_fe_ptr[structdim][obj_index] + fe_index_local,
2558 dof_handler.object_dof_ptr[obj_level][structdim].size());
2559
2560 const unsigned int ptr_0 =
2561 dof_handler
2562 .object_dof_ptr[obj_level][structdim]
2563 [dof_handler.hp_object_fe_ptr[structdim][obj_index] +
2564 fe_index_local];
2565 const unsigned int ptr_1 =
2566 dof_handler
2567 .object_dof_ptr[obj_level][structdim]
2568 [dof_handler.hp_object_fe_ptr[structdim][obj_index] +
2569 fe_index_local + 1];
2570
2571 return {ptr_0, ptr_1 - ptr_0};
2572 }
2573
2574 template <int dim, int spacedim, int structdim, bool level_dof_access>
2575 static std::pair<unsigned int, unsigned int>
2577 const ::DoFAccessor<structdim, dim, spacedim, level_dof_access>
2578 accessor,
2579 const types::fe_index fe_index)
2580 {
2581 return process_object_range(accessor.get_dof_handler(),
2582 accessor.level(),
2583 accessor.index(),
2584 fe_index,
2585 std::integral_constant<int, structdim>());
2586 }
2587
2588 template <int dim, int spacedim, int structdim>
2589 static std::pair<unsigned int, unsigned int>
2591 const unsigned int)
2592 {
2594
2595 return {0, 0};
2596 }
2597
2598
2599
2608 template <int dim,
2609 int spacedim,
2610 int structdim,
2611 typename DoFProcessor,
2612 typename DoFMapping>
2613 static DEAL_II_ALWAYS_INLINE void
2615 const unsigned int obj_level,
2616 const unsigned int obj_index,
2617 const types::fe_index fe_index,
2618 const DoFMapping &mapping,
2619 const std::integral_constant<int, structdim> &dd,
2620 types::global_dof_index *&dof_indices_ptr,
2621 const DoFProcessor &process)
2622 {
2623 Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
2624
2625 // determine range of dofs in global data structure
2626 const auto range =
2627 process_object_range(dof_handler, obj_level, obj_index, fe_index, dd);
2628 if (range.second == 0)
2629 return;
2630
2631 std::vector<types::global_dof_index> &object_dof_indices =
2632 dof_handler
2633 .object_dof_indices[structdim < dim ? 0 : obj_level][structdim];
2634 AssertIndexRange(range.first, object_dof_indices.size());
2636 object_dof_indices.data() + range.first;
2637
2638 // process dofs
2639 for (unsigned int i = 0; i < range.second; ++i)
2640 {
2641 process(
2642 stored_indices[(structdim == 0 || structdim == dim) ? i :
2643 mapping(i)],
2644 dof_indices_ptr);
2645 if (dof_indices_ptr != nullptr)
2646 ++dof_indices_ptr;
2647 }
2648 }
2649
2650
2651
2662 template <int dim, int spacedim, int structdim>
2663 static void
2665 const unsigned int obj_level,
2666 const unsigned int obj_index,
2667 const types::fe_index fe_index,
2668 const unsigned int local_index,
2669 const std::integral_constant<int, structdim> &dd,
2670 const types::global_dof_index global_index)
2671 {
2672 process_dof_index(dof_handler,
2673 obj_level,
2674 obj_index,
2675 fe_index,
2676 local_index,
2677 dd,
2678 global_index,
2679 [](auto &ptr, const auto &value) { ptr = value; });
2680 }
2681
2682
2693 template <int dim, int spacedim, int structdim>
2696 const unsigned int obj_level,
2697 const unsigned int obj_index,
2698 const types::fe_index fe_index,
2699 const unsigned int local_index,
2700 const std::integral_constant<int, structdim> &dd)
2701 {
2702 types::global_dof_index global_index;
2703 process_dof_index(dof_handler,
2704 obj_level,
2705 obj_index,
2706 fe_index,
2707 local_index,
2708 dd,
2709 global_index,
2710 [](const auto &ptr, auto &value) { value = ptr; });
2711 return global_index;
2712 }
2713
2714
2715 template <int dim, int spacedim>
2718 const int level,
2719 const unsigned int vertex_index,
2720 const unsigned int i)
2721 {
2722 Assert(dof_handler.hp_capability_enabled == false,
2723 ExcMessage(
2724 "DoFHandler in hp-mode does not implement multilevel DoFs."));
2725
2726 return dof_handler.mg_vertex_dofs[vertex_index].access_index(
2727 level, i, dof_handler.get_fe().n_dofs_per_vertex());
2728 }
2729
2730
2731
2741 template <int dim, int spacedim, int structdim>
2742 static unsigned int
2744 const unsigned int obj_level,
2745 const unsigned int obj_index,
2746 const std::integral_constant<int, structdim> &)
2747 {
2748 (void)obj_level;
2749
2750 Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
2751
2752 // 1) no hp used -> fe_index == 0
2753 if (dof_handler.hp_capability_enabled == false)
2754 return 1;
2755
2756 // 2) cell and hp is used -> there is only one fe_index
2757 if (structdim == dim)
2758 return 1;
2759
2760 // 3) general entity and hp is used
2761 AssertIndexRange(structdim, dof_handler.hp_object_fe_ptr.size());
2762 AssertIndexRange(obj_index + 1,
2763 dof_handler.hp_object_fe_ptr[structdim].size());
2764
2765 return dof_handler.hp_object_fe_ptr[structdim][obj_index + 1] -
2766 dof_handler.hp_object_fe_ptr[structdim][obj_index];
2767 }
2768
2769
2770
2780 template <int dim, int spacedim, int structdim>
2781 static types::fe_index
2783 const unsigned int obj_level,
2784 const unsigned int obj_index,
2785 const unsigned int local_index,
2786 const std::integral_constant<int, structdim> &)
2787 {
2788 Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
2789
2790 // for cells only one active FE index available
2791 Assert(((structdim == dim) &&
2793 false,
2795
2796 // 1) no hp used -> fe_index == 0
2797 if (dof_handler.hp_capability_enabled == false)
2799
2800 // 2) cell and hp is used -> there is only one fe_index
2801 if (structdim == dim)
2802 return dof_handler.hp_cell_active_fe_indices[obj_level][obj_index];
2803
2804 // 3) general entity and hp is used
2805 AssertIndexRange(structdim, dof_handler.hp_object_fe_indices.size());
2806 AssertIndexRange(structdim, dof_handler.hp_object_fe_ptr.size());
2807 AssertIndexRange(obj_index,
2808 dof_handler.hp_object_fe_ptr[structdim].size());
2809 AssertIndexRange(dof_handler.hp_object_fe_ptr[structdim][obj_index] +
2810 local_index,
2811 dof_handler.hp_object_fe_indices[structdim].size());
2812
2813 return dof_handler.hp_object_fe_indices
2814 [structdim]
2815 [dof_handler.hp_object_fe_ptr[structdim][obj_index] + local_index];
2816 }
2817
2818
2819
2832 template <int dim, int spacedim, int structdim>
2833 static std::set<types::fe_index>
2835 const unsigned int obj_level,
2836 const unsigned int obj_index,
2837 const std::integral_constant<int, structdim> &t)
2838 {
2839 Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
2840
2841 // 1) no hp used -> fe_index == 0
2842 if (dof_handler.hp_capability_enabled == false)
2844
2845 // 2) cell and hp is used -> there is only one fe_index
2846 if (structdim == dim)
2847 return {dof_handler.hp_cell_active_fe_indices[obj_level][obj_index]};
2848
2849 // 3) general entity and hp is used
2850 std::set<types::fe_index> active_fe_indices;
2851 for (unsigned int i = 0;
2852 i < n_active_fe_indices(dof_handler, obj_level, obj_index, t);
2853 ++i)
2854 active_fe_indices.insert(
2855 nth_active_fe_index(dof_handler, obj_level, obj_index, i, t));
2856 return active_fe_indices;
2857 }
2858
2859
2860
2861 template <int dim, int spacedim, int structdim>
2862 static bool
2864 const unsigned int obj_level,
2865 const unsigned int obj_index,
2866 const types::fe_index fe_index,
2867 const std::integral_constant<int, structdim> &)
2868 {
2869 Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
2870
2871 // 1) no hp used -> fe_index == 0
2872 if (dof_handler.hp_capability_enabled == false)
2874
2875 // 2) cell and hp is used -> there is only one fe_index
2876 if (structdim == dim)
2877 return dof_handler.hp_cell_active_fe_indices[obj_level][obj_index] ==
2878 fe_index;
2879
2880 // 3) general entity and hp is used
2881 return std::find(
2882 dof_handler.hp_object_fe_indices[structdim].begin() +
2883 dof_handler.hp_object_fe_ptr[structdim][obj_index],
2884 dof_handler.hp_object_fe_indices[structdim].begin() +
2885 dof_handler.hp_object_fe_ptr[structdim][obj_index + 1],
2886 fe_index) !=
2887 (dof_handler.hp_object_fe_indices[structdim].begin() +
2888 dof_handler.hp_object_fe_ptr[structdim][obj_index + 1]);
2889 }
2890
2891
2892
2893 template <typename InputVector, typename ForwardIterator>
2894 static void
2895 extract_subvector_to(const InputVector &values,
2896 const types::global_dof_index *cache,
2897 const types::global_dof_index *cache_end,
2898 ForwardIterator local_values_begin)
2899 {
2900 values.extract_subvector_to(cache, cache_end, local_values_begin);
2901 }
2902
2903
2904
2905#ifdef DEAL_II_WITH_TRILINOS
2906 static std::vector<unsigned int>
2908 const types::global_dof_index *v_end)
2909 {
2910 // initialize original index locations
2911 std::vector<unsigned int> idx(v_end - v_begin);
2912 std::iota(idx.begin(), idx.end(), 0u);
2913
2914 // sort indices based on comparing values in v
2915 std::sort(idx.begin(),
2916 idx.end(),
2917 [&v_begin](unsigned int i1, unsigned int i2) {
2918 return *(v_begin + i1) < *(v_begin + i2);
2919 });
2920
2921 return idx;
2922 }
2923
2924
2925
2926# ifdef DEAL_II_TRILINOS_WITH_TPETRA
2927 template <typename ForwardIterator, typename Number, typename MemorySpace>
2928 static void
2931 &values,
2932 const types::global_dof_index *cache_begin,
2933 const types::global_dof_index *cache_end,
2934 ForwardIterator local_values_begin)
2935 {
2936 std::vector<unsigned int> sorted_indices_pos =
2937 sort_indices(cache_begin, cache_end);
2938 const unsigned int cache_size = cache_end - cache_begin;
2939 std::vector<types::global_dof_index> cache_indices(cache_size);
2940 for (unsigned int i = 0; i < cache_size; ++i)
2941 cache_indices[i] = *(cache_begin + sorted_indices_pos[i]);
2942
2943 IndexSet index_set(cache_indices.back() + 1);
2944 index_set.add_indices(cache_indices.begin(), cache_indices.end());
2945 index_set.compress();
2946 LinearAlgebra::ReadWriteVector<Number> read_write_vector(index_set);
2947 read_write_vector.import_elements(values, VectorOperation::insert);
2948
2949 // Copy the elements from read_write_vector and reorder them.
2950 for (unsigned int i = 0; i < cache_size; ++i, ++local_values_begin)
2951 *local_values_begin = read_write_vector[sorted_indices_pos[i]];
2952 }
2953# endif
2954
2955
2956
2957 template <typename ForwardIterator>
2958 static void
2960 const types::global_dof_index *cache_begin,
2961 const types::global_dof_index *cache_end,
2962 ForwardIterator local_values_begin)
2963 {
2964 std::vector<unsigned int> sorted_indices_pos =
2965 sort_indices(cache_begin, cache_end);
2966 const unsigned int cache_size = cache_end - cache_begin;
2967 std::vector<types::global_dof_index> cache_indices(cache_size);
2968 for (unsigned int i = 0; i < cache_size; ++i)
2969 cache_indices[i] = *(cache_begin + sorted_indices_pos[i]);
2970
2971 IndexSet index_set(cache_indices.back() + 1);
2972 index_set.add_indices(cache_indices.begin(), cache_indices.end());
2973 index_set.compress();
2974 LinearAlgebra::ReadWriteVector<double> read_write_vector(index_set);
2975 read_write_vector.import_elements(values, VectorOperation::insert);
2976
2977 // Copy the elements from read_write_vector and reorder them.
2978 for (unsigned int i = 0; i < cache_size; ++i, ++local_values_begin)
2979 *local_values_begin = read_write_vector[sorted_indices_pos[i]];
2980 }
2981#endif
2982
2987 template <int dim, int spacedim, bool level_dof_access, int structdim>
2988 static unsigned int
2990 const ::DoFAccessor<structdim, dim, spacedim, level_dof_access>
2991 &accessor,
2992 const types::fe_index fe_index_,
2993 const bool count_level_dofs)
2994 {
2995 // note: we cannot rely on the template parameter level_dof_access here,
2996 // since the function get_mg_dof_indices()/set_mg_dof_indices() can be
2997 // called even if level_dof_access==false.
2998 if (count_level_dofs)
2999 {
3000 const auto &fe = accessor.get_fe(fe_index_);
3001
3002 const unsigned int //
3003 dofs_per_vertex = fe.n_dofs_per_vertex(), //
3004 dofs_per_line = fe.n_dofs_per_line(), //
3005 dofs_per_quad = fe.n_dofs_per_quad(0 /*dummy*/), //
3006 dofs_per_hex = fe.n_dofs_per_hex(); //
3007
3008 unsigned int index = 0;
3009
3010 // 1) VERTEX dofs
3011 index += dofs_per_vertex * accessor.n_vertices();
3012
3013 // 2) LINE dofs
3014 if (structdim == 2 || structdim == 3)
3015 index += dofs_per_line * accessor.n_lines();
3016
3017 // 3) FACE dofs
3018 if (structdim == 3)
3019 index += dofs_per_quad * accessor.n_faces();
3020
3021 // 4) INNER dofs
3022 const unsigned int interior_dofs =
3023 structdim == 1 ? dofs_per_line :
3024 (structdim == 2 ? dofs_per_quad : dofs_per_hex);
3025
3026 index += interior_dofs;
3027
3028 return index;
3029 }
3030 else
3031 {
3032 const auto fe_index =
3034 accessor, fe_index_);
3035
3036 unsigned int index = 0;
3037
3038 // 1) VERTEX dofs
3039 for (const auto vertex : accessor.vertex_indices())
3040 index += process_object_range(accessor.get_dof_handler(),
3041 0,
3042 accessor.vertex_index(vertex),
3043 fe_index,
3044 std::integral_constant<int, 0>())
3045 .second;
3046
3047 // 2) LINE dofs
3048 if (structdim == 2 || structdim == 3)
3049 for (const auto line : accessor.line_indices())
3050 index +=
3051 process_object_range(*accessor.line(line), fe_index).second;
3052
3053 // 3) FACE dofs
3054 if (structdim == 3)
3055 for (const auto face : accessor.face_indices())
3056 index +=
3057 process_object_range(*accessor.quad(face), fe_index).second;
3058
3059 // 4) INNER dofs
3060 index += process_object_range(accessor, fe_index).second;
3061
3062 return index;
3063 }
3064 }
3065
3066
3067
3068 // The next few internal helper functions are needed to support various
3069 // DoFIndicesType kinds, e.g. actual vectors of DoFIndices or empty
3070 // types that we use when we only want to work on the internally stored
3071 // DoFs and never extract any number.
3072 template <typename ArrayType>
3073 static unsigned int
3074 get_array_length(const ArrayType &array)
3075 {
3076 return array.size();
3077 }
3078
3079 static unsigned int
3080 get_array_length(const std::tuple<> &)
3081 {
3082 return 0;
3083 }
3084
3085 template <typename ArrayType>
3087 get_array_ptr(const ArrayType &array)
3088 {
3089 return const_cast<types::global_dof_index *>(array.data());
3090 }
3091
3093 get_array_ptr(const std::tuple<> &)
3094 {
3095 return nullptr;
3096 }
3097
3098
3099
3105 template <int dim,
3106 int spacedim,
3107 bool level_dof_access,
3108 int structdim,
3109 typename DoFIndicesType,
3110 typename DoFOperation,
3111 typename DoFProcessor>
3112 static void
3114 const ::DoFAccessor<structdim, dim, spacedim, level_dof_access>
3115 &accessor,
3116 const DoFIndicesType &const_dof_indices,
3117 const types::fe_index fe_index_,
3118 const DoFOperation &dof_operation,
3119 const DoFProcessor &dof_processor,
3120 const bool count_level_dofs)
3121 {
3122 const types::fe_index fe_index =
3124 accessor, fe_index_);
3125
3126 // we cannot rely on the template parameter level_dof_access here, since
3127 // the function get_mg_dof_indices()/set_mg_dof_indices() can be called
3128 // even if level_dof_access==false.
3129 (void)count_level_dofs;
3130
3131 const auto &fe = accessor.get_fe(fe_index);
3132
3133 // we want to pass in rvalue 'std::tuple<>' types as `DoFIndicesType`,
3134 // but we need non-const references for std::vector<> types, so get in
3135 // a const reference here and immediately cast the constness away -
3136 // note that any use of the dereferenced invalid type will result in a
3137 // segfault
3138 types::global_dof_index *dof_indices_ptr =
3139 get_array_ptr(const_dof_indices);
3140 types::global_dof_index *end_dof_indices =
3141 dof_indices_ptr + get_array_length(const_dof_indices);
3142
3143 // 1) VERTEX dofs, only step into the functions if we actually have
3144 // DoFs on them
3145 if (fe.n_dofs_per_vertex() > 0)
3146 for (const auto vertex : accessor.vertex_indices())
3147 dof_operation.process_vertex_dofs(*accessor.dof_handler,
3148 accessor.vertex_index(vertex),
3149 fe_index,
3150 dof_indices_ptr,
3151 dof_processor);
3152
3153 // 2) copy dof numbers from the LINE. for lines with the wrong
3154 // orientation (which might occur in 3d), we have already made sure that
3155 // we're ok by picking the correct vertices (this happens automatically
3156 // in the vertex() function). however, if the line is in wrong
3157 // orientation, we look at it in flipped orientation and we will have to
3158 // adjust the shape function indices that we see to correspond to the
3159 // correct (face/cell-local) ordering.
3160 if ((structdim == 2 || structdim == 3) && fe.n_dofs_per_line() > 0)
3161 {
3162 const auto line_indices = internal::TriaAccessorImplementation::
3163 Implementation::get_line_indices_of_cell(accessor);
3164 const auto line_orientations =
3165 internal::TriaAccessorImplementation::Implementation::
3166 get_line_orientations_of_cell(accessor);
3167
3168 for (const auto line : accessor.line_indices())
3169 {
3170 const auto line_orientation = line_orientations[line];
3171 if (line_orientation == numbers::default_geometric_orientation)
3172 dof_operation.process_dofs(
3173 accessor.get_dof_handler(),
3174 0,
3175 line_indices[line],
3176 fe_index,
3177 [](const auto d) { return d; },
3178 std::integral_constant<int, 1>(),
3179 dof_indices_ptr,
3180 dof_processor);
3181 else
3182 dof_operation.process_dofs(
3183 accessor.get_dof_handler(),
3184 0,
3185 line_indices[line],
3186 fe_index,
3187 [&fe, line_orientation](const auto d) {
3188 return fe.adjust_line_dof_index_for_line_orientation(
3189 d, line_orientation);
3190 },
3191 std::integral_constant<int, 1>(),
3192 dof_indices_ptr,
3193 dof_processor);
3194 }
3195 }
3196
3197 // 3) copy dof numbers from the FACE. for faces with the wrong
3198 // orientation, we have already made sure that we're ok by picking the
3199 // correct lines and vertices (this happens automatically in the line()
3200 // and vertex() functions). however, if the face is in wrong
3201 // orientation, we look at it in flipped orientation and we will have to
3202 // adjust the shape function indices that we see to correspond to the
3203 // correct (cell-local) ordering. The same applies, if the face_rotation
3204 // or face_orientation is non-standard
3205 if (structdim == 3 && fe.max_dofs_per_quad() > 0)
3206 for (const auto face_no : accessor.face_indices())
3207 {
3208 const auto combined_orientation =
3209 accessor.combined_face_orientation(face_no);
3210 const unsigned int quad_index = accessor.quad_index(face_no);
3211 if (combined_orientation ==
3213 dof_operation.process_dofs(
3214 accessor.get_dof_handler(),
3215 0,
3216 quad_index,
3217 fe_index,
3218 [](const auto d) { return d; },
3219 std::integral_constant<int, 2>(),
3220 dof_indices_ptr,
3221 dof_processor);
3222 else
3223 dof_operation.process_dofs(
3224 accessor.get_dof_handler(),
3225 0,
3226 quad_index,
3227 fe_index,
3228 [&](const auto d) {
3229 return fe.adjust_quad_dof_index_for_face_orientation(
3230 d, face_no, combined_orientation);
3231 },
3232 std::integral_constant<int, 2>(),
3233 dof_indices_ptr,
3234 dof_processor);
3235 }
3236
3237 // 4) INNER dofs - here we need to make sure that the shortcut to not
3238 // run the function does not miss the faces of wedge and pyramid
3239 // elements where n_dofs_per_object might not return the largest
3240 // possible value
3241 if (((dim == 3 && structdim == 2) ?
3242 fe.max_dofs_per_quad() :
3243 fe.template n_dofs_per_object<structdim>()) > 0)
3244 dof_operation.process_dofs(
3245 accessor.get_dof_handler(),
3246 accessor.level(),
3247 accessor.index(),
3248 fe_index,
3249 [&](const auto d) { return d; },
3250 std::integral_constant<int, structdim>(),
3251 dof_indices_ptr,
3252 dof_processor);
3253
3254 if (dof_indices_ptr != nullptr)
3255 {
3256 AssertDimension(n_dof_indices(accessor, fe_index, count_level_dofs),
3257 dof_indices_ptr - get_array_ptr(const_dof_indices));
3258 }
3259
3260 // PM: This is a part that should not be reached since it indicates that
3261 // an object (and/or its subobjects) is not active. However,
3262 // unfortunately this function is called by
3263 // DoFTools::set_periodicity_constraints() indirectly by
3264 // get_dof_indices() also for artificial faces to determine if a face
3265 // is artificial.
3267 for (; dof_indices_ptr < end_dof_indices; ++dof_indices_ptr)
3268 dof_processor(invalid_index, dof_indices_ptr);
3269 }
3270
3271
3272
3277 template <int dim, int spacedim>
3279 {
3283 template <typename DoFProcessor>
3286 const unsigned int vertex_index,
3287 const types::fe_index fe_index,
3288 types::global_dof_index *&dof_indices_ptr,
3289 const DoFProcessor &dof_processor) const
3290 {
3292 dof_handler,
3293 0,
3294 vertex_index,
3295 fe_index,
3296 [](const auto d) {
3298 return d;
3299 },
3300 std::integral_constant<int, 0>(),
3301 dof_indices_ptr,
3302 dof_processor);
3303 }
3304
3308 template <int structdim, typename DoFMapping, typename DoFProcessor>
3311 const unsigned int obj_level,
3312 const unsigned int obj_index,
3313 const types::fe_index fe_index,
3314 const DoFMapping &mapping,
3315 const std::integral_constant<int, structdim>,
3316 types::global_dof_index *&dof_indices_ptr,
3317 const DoFProcessor &dof_processor) const
3318 {
3320 dof_handler,
3321 obj_level,
3322 obj_index,
3323 fe_index,
3324 mapping,
3325 std::integral_constant<int, std::min(structdim, dim)>(),
3326 dof_indices_ptr,
3327 dof_processor);
3328 }
3329 };
3330
3331
3332
3337 template <int dim, int spacedim>
3339 {
3343 MGDoFIndexProcessor(const unsigned int level)
3344 : level(level)
3345 {}
3346
3350 template <typename DoFProcessor>
3353 const unsigned int vertex_index,
3354 const types::fe_index,
3355 types::global_dof_index *&dof_indices_ptr,
3356 const DoFProcessor &dof_processor) const
3357 {
3358 const unsigned int n_indices =
3359 dof_handler.get_fe(0).template n_dofs_per_object<0>();
3360 types::global_dof_index *stored_indices =
3361 &dof_handler.mg_vertex_dofs[vertex_index].access_index(level,
3362 0,
3363 n_indices);
3364 for (unsigned int d = 0; d < n_indices; ++d, ++dof_indices_ptr)
3365 dof_processor(stored_indices[d], dof_indices_ptr);
3366 }
3367
3371 template <int structdim, typename DoFMapping, typename DoFProcessor>
3374 const unsigned int,
3375 const unsigned int obj_index,
3376 const types::fe_index fe_index,
3377 const DoFMapping &mapping,
3378 const std::integral_constant<int, structdim>,
3379 types::global_dof_index *&dof_indices_ptr,
3380 const DoFProcessor &dof_processor) const
3381 {
3382 const unsigned int n_indices =
3383 dof_handler.get_fe(0).template n_dofs_per_object<structdim>();
3384 types::global_dof_index *stored_indices = &get_mg_dof_index(
3385 dof_handler,
3386 dof_handler.mg_levels[level],
3387 dof_handler.mg_faces,
3388 obj_index,
3389 fe_index,
3390 0,
3391 std::integral_constant<int, std::min(structdim, dim)>());
3392 for (unsigned int d = 0; d < n_indices; ++d, ++dof_indices_ptr)
3393 dof_processor(stored_indices[structdim < dim ? mapping(d) : d],
3394 dof_indices_ptr);
3395 }
3396
3397 private:
3398 const unsigned int level;
3399 };
3400
3401
3402
3403 template <int dim, int spacedim, bool level_dof_access, int structdim>
3404 static void
3406 const ::DoFAccessor<structdim, dim, spacedim, level_dof_access>
3407 &accessor,
3408 std::vector<types::global_dof_index> &dof_indices,
3409 const types::fe_index fe_index)
3410 {
3412 accessor,
3413 dof_indices,
3414 fe_index,
3416 [](auto stored_index, auto dof_ptr) { *dof_ptr = stored_index; },
3417 false);
3418 }
3419
3420
3421
3422 template <int dim, int spacedim, bool level_dof_access, int structdim>
3423 static void
3425 const ::DoFAccessor<structdim, dim, spacedim, level_dof_access>
3426 &accessor,
3427 const std::vector<types::global_dof_index> &dof_indices,
3428 const types::fe_index fe_index)
3429 {
3430 // Note: this function is as general as `get_dof_indices()`. This
3431 // assert is placed here since it is currently only used by the
3432 // function DoFCellAccessor::set_dof_indices(), which is called by
3433 // internal::DoFHandlerImplementation::Policy::Implementation::distribute_dofs().
3434 // In the case of new use cases, this assert can be removed.
3435 Assert(
3436 dim == structdim,
3437 ExcMessage(
3438 "This function is intended to be used for DoFCellAccessor, i.e., "
3439 "dimension == structdim."));
3440
3442 accessor,
3443 dof_indices,
3444 fe_index,
3446 [](auto &stored_index, auto dof_ptr) { stored_index = *dof_ptr; },
3447 false);
3448 }
3449
3450
3451
3452 template <int dim, int spacedim, bool level_dof_access, int structdim>
3453 static void
3455 const ::DoFAccessor<structdim, dim, spacedim, level_dof_access>
3456 &accessor,
3457 const int level,
3458 std::vector<types::global_dof_index> &dof_indices,
3459 const types::fe_index fe_index)
3460 {
3462 ExcMessage("MG DoF indices cannot be queried in hp case"));
3464 accessor,
3465 dof_indices,
3466 fe_index,
3468 [](auto stored_index, auto dof_ptr) { *dof_ptr = stored_index; },
3469 true);
3470 }
3471
3472
3473
3474 template <int dim, int spacedim, bool level_dof_access, int structdim>
3475 static void
3477 const ::DoFAccessor<structdim, dim, spacedim, level_dof_access>
3478 &accessor,
3479 const int level,
3480 const std::vector<types::global_dof_index> &dof_indices,
3481 const types::fe_index fe_index)
3482 {
3484 ExcMessage("MG DoF indices cannot be queried in hp case"));
3485
3486 // Note: this function is as general as `get_mg_dof_indices()`. This
3487 // assert is placed here since it is currently only used by the
3488 // function DoFCellAccessor::set_mg_dof_indices(), which is called by
3489 // internal::DoFHandlerImplementation::Policy::Implementation::distribute_mg_dofs().
3490 // In the case of new use cases, this assert can be removed.
3491 Assert(dim == structdim,
3492 ExcMessage("This function is intended to be used for "
3493 "DoFCellAccessor, i.e., dimension == structdim."));
3494
3496 accessor,
3497 dof_indices,
3498 fe_index,
3500 [](auto &stored_index, auto dof_ptr) { stored_index = *dof_ptr; },
3501 true);
3502 }
3503
3504
3505
3506 template <int dim, int spacedim>
3509 const DoFHandler<dim, spacedim> &dof_handler,
3511 &mg_level,
3513 &,
3514 const unsigned int obj_index,
3515 const types::fe_index fe_index,
3516 const unsigned int local_index,
3517 const std::integral_constant<int, dim>)
3518 {
3519 Assert(dof_handler.hp_capability_enabled == false,
3521
3522 return mg_level->dof_object.access_dof_index(
3523 static_cast<const DoFHandler<dim, spacedim> &>(dof_handler),
3524 obj_index,
3525 fe_index,
3526 local_index);
3527 }
3528
3529
3530
3531 template <int dim, int spacedim, std::enable_if_t<(dim > 1), int> = 0>
3534 const DoFHandler<dim, spacedim> &dof_handler,
3536 &,
3538 &mg_faces,
3539 const unsigned int obj_index,
3540 const types::fe_index fe_index,
3541 const unsigned int local_index,
3542 const std::integral_constant<int, 1>)
3543 {
3544 return mg_faces->lines.access_dof_index(
3545 static_cast<const DoFHandler<dim, spacedim> &>(dof_handler),
3546 obj_index,
3547 fe_index,
3548 local_index);
3549 }
3550
3551
3552
3553 template <int spacedim>
3556 const DoFHandler<3, spacedim> &dof_handler,
3558 &,
3560 &mg_faces,
3561 const unsigned int obj_index,
3562 const types::fe_index fe_index,
3563 const unsigned int local_index,
3564 const std::integral_constant<int, 2>)
3565 {
3566 Assert(dof_handler.hp_capability_enabled == false,
3568 return mg_faces->quads.access_dof_index(
3569 static_cast<const DoFHandler<3, spacedim> &>(dof_handler),
3570 obj_index,
3571 fe_index,
3572 local_index);
3573 }
3574 };
3575
3576
3577
3578 template <int dim, int spacedim, bool level_dof_access>
3579 void
3581 const ::DoFCellAccessor<dim, spacedim, level_dof_access> &accessor,
3583 const unsigned int fe_index);
3584 } // namespace DoFAccessorImplementation
3585} // namespace internal
3586
3587
3588
3589template <int structdim, int dim, int spacedim, bool level_dof_access>
3592 const unsigned int i,
3593 const types::fe_index fe_index_) const
3594{
3595 const auto fe_index =
3597 fe_index_);
3598
3599 // access the respective DoF
3600 return ::internal::DoFAccessorImplementation::Implementation::
3601 get_dof_index(*this->dof_handler,
3602 this->level(),
3603 this->index(),
3604 fe_index,
3605 i,
3606 std::integral_constant<int, structdim>());
3607}
3608
3609
3610template <int structdim, int dim, int spacedim, bool level_dof_access>
3613 const int level,
3614 const unsigned int i) const
3615{
3617 *this->dof_handler,
3618 this->dof_handler->mg_levels[level],
3619 this->dof_handler->mg_faces,
3620 this->index(),
3621 0,
3622 i,
3623 std::integral_constant<int, structdim>());
3624}
3625
3626
3627template <int structdim, int dim, int spacedim, bool level_dof_access>
3628inline void
3630 const unsigned int i,
3631 const types::global_dof_index index,
3632 const types::fe_index fe_index_) const
3633{
3634 const auto fe_index =
3636 fe_index_);
3637
3638 // access the respective DoF
3640 *this->dof_handler,
3641 this->level(),
3642 this->index(),
3643 fe_index,
3644 i,
3645 std::integral_constant<int, structdim>(),
3646 index);
3647}
3648
3649
3650
3651template <int structdim, int dim, int spacedim, bool level_dof_access>
3652inline void
3654 const int level,
3655 const unsigned int i,
3656 const types::global_dof_index index) const
3657{
3659 *this->dof_handler,
3660 this->dof_handler->mg_levels[level],
3661 this->dof_handler->mg_faces,
3662 this->index(),
3663 0,
3664 i,
3665 std::integral_constant<int, structdim>()) = index;
3666}
3667
3668
3669
3670template <int structdim, int dim, int spacedim, bool level_dof_access>
3671inline unsigned int
3673 const
3674{
3675 // access the respective DoF
3676 return ::internal::DoFAccessorImplementation::Implementation::
3677 n_active_fe_indices(*this->dof_handler,
3678 this->level(),
3679 this->index(),
3680 std::integral_constant<int, structdim>());
3681}
3682
3683
3684
3685template <int structdim, int dim, int spacedim, bool level_dof_access>
3686inline types::fe_index
3688 const unsigned int n) const
3689{
3690 // access the respective DoF
3691 return ::internal::DoFAccessorImplementation::Implementation::
3692 nth_active_fe_index(*this->dof_handler,
3693 this->level(),
3694 this->index(),
3695 n,
3696 std::integral_constant<int, structdim>());
3697}
3698
3699
3700
3701template <int structdim, int dim, int spacedim, bool level_dof_access>
3702inline std::set<types::fe_index>
3704 const
3705{
3706 std::set<types::fe_index> active_fe_indices;
3707 for (unsigned int i = 0; i < n_active_fe_indices(); ++i)
3708 active_fe_indices.insert(nth_active_fe_index(i));
3709 return active_fe_indices;
3710}
3711
3712
3713
3714template <int structdim, int dim, int spacedim, bool level_dof_access>
3715inline bool
3717 const types::fe_index fe_index) const
3718{
3719 // access the respective DoF
3720 return ::internal::DoFAccessorImplementation::Implementation::
3721 fe_index_is_active(*this->dof_handler,
3722 this->level(),
3723 this->index(),
3724 fe_index,
3725 std::integral_constant<int, structdim>());
3726}
3727
3728
3729
3730template <int structdim, int dim, int spacedim, bool level_dof_access>
3733 const unsigned int vertex,
3734 const unsigned int i,
3735 const types::fe_index fe_index_) const
3736{
3737 const types::fe_index fe_index =
3738 (((this->dof_handler->hp_capability_enabled == false) &&
3739 (fe_index_ == numbers::invalid_fe_index)) ?
3740 // No hp enabled, and the argument is at its default value -> we
3741 // can translate to the default active fe index
3743 // Otherwise: If anything other than the default is provided by
3744 // the caller, then we should take just that. As an exception, if
3745 // we are on a cell (rather than a face/edge/vertex), then we know
3746 // that there is only one active fe index on this cell and we can
3747 // use that:
3748 ((dim == structdim) && (fe_index_ == numbers::invalid_fe_index) ?
3749 this->nth_active_fe_index(0) :
3750 fe_index_));
3751
3752 return ::internal::DoFAccessorImplementation::Implementation::
3753 get_dof_index(*this->dof_handler,
3754 0,
3755 this->vertex_index(vertex),
3756 fe_index,
3757 i,
3758 std::integral_constant<int, 0>());
3759}
3760
3761
3762template <int structdim, int dim, int spacedim, bool level_dof_access>
3765 const int level,
3766 const unsigned int vertex,
3767 const unsigned int i,
3768 const types::fe_index fe_index_) const
3769{
3770 const auto fe_index =
3772 fe_index_);
3773 (void)fe_index;
3774 Assert(this->dof_handler != nullptr, ExcInvalidObject());
3775 Assert(this->dof_handler->mg_vertex_dofs.size() > 0,
3776 ExcMessage("Multigrid DoF indices can only be accessed after "
3777 "DoFHandler::distribute_mg_dofs() has been called!"));
3778 AssertIndexRange(vertex, this->n_vertices());
3779 AssertIndexRange(i, this->dof_handler->get_fe(fe_index).n_dofs_per_vertex());
3780
3781 Assert(dof_handler->hp_capability_enabled == false,
3782 ExcMessage(
3783 "DoFHandler in hp-mode does not implement multilevel DoFs."));
3784
3785 return this->dof_handler->mg_vertex_dofs[this->vertex_index(vertex)]
3786 .access_index(level, i, this->dof_handler->get_fe().n_dofs_per_vertex());
3787}
3788
3789
3790
3791template <int structdim, int dim, int spacedim, bool level_dof_access>
3792inline void
3795 const unsigned int vertex,
3796 const unsigned int i,
3797 const types::global_dof_index index,
3798 const types::fe_index fe_index_) const
3799{
3800 const auto fe_index =
3802 fe_index_);
3803 (void)fe_index;
3804 Assert(this->dof_handler != nullptr, ExcInvalidObject());
3805 AssertIndexRange(vertex, this->n_vertices());
3806 AssertIndexRange(i, this->dof_handler->get_fe(fe_index).n_dofs_per_vertex());
3807
3808 Assert(dof_handler->hp_capability_enabled == false,
3809 ExcMessage(
3810 "DoFHandler in hp-mode does not implement multilevel DoFs."));
3811
3812 this->dof_handler->mg_vertex_dofs[this->vertex_index(vertex)].access_index(
3813 level, i, this->dof_handler->get_fe().n_dofs_per_vertex()) = index;
3814}
3815
3816
3817
3818template <int structdim, int dim, int spacedim, bool level_dof_access>
3819inline const FiniteElement<dim, spacedim> &
3821 const types::fe_index fe_index) const
3822{
3823 Assert(fe_index_is_active(fe_index) == true,
3824 ExcMessage("This function can only be called for active FE indices"));
3825
3826 return this->dof_handler->get_fe(fe_index);
3827}
3828
3829
3830
3831template <int structdim, int dim, int spacedim, bool level_dof_access>
3832inline typename ::internal::DoFHandlerImplementation::
3833 Iterators<dim, spacedim, level_dof_access>::line_iterator
3835 const unsigned int i) const
3836{
3837 // if we are asking for a particular line and this object refers to
3838 // a line, then the only valid index is i==0 and we should return
3839 // *this
3840 if (structdim == 1)
3841 {
3842 Assert(i == 0,
3843 ExcMessage("You can only ask for line zero if the "
3844 "current object is a line itself."));
3845 return typename ::internal::DoFHandlerImplementation::
3846 Iterators<dim, spacedim, level_dof_access>::cell_iterator(
3847 &this->get_triangulation(),
3848 this->level(),
3849 this->index(),
3850 &this->get_dof_handler());
3851 }
3852
3853 // otherwise we need to be in structdim>=2
3854 Assert(structdim > 1, ExcImpossibleInDim(structdim));
3855 Assert(dim > 1, ExcImpossibleInDim(dim));
3856
3857 // checking of 'i' happens in line_index(i)
3858 return typename ::internal::DoFHandlerImplementation::
3859 Iterators<dim, spacedim, level_dof_access>::line_iterator(
3860 this->tria,
3861 0, // only sub-objects are allowed, which have no level
3862 this->line_index(i),
3863 this->dof_handler);
3864}
3865
3866
3867template <int structdim, int dim, int spacedim, bool level_dof_access>
3868inline typename ::internal::DoFHandlerImplementation::
3869 Iterators<dim, spacedim, level_dof_access>::quad_iterator
3871 const unsigned int i) const
3872{
3873 // if we are asking for a
3874 // particular quad and this object
3875 // refers to a quad, then the only
3876 // valid index is i==0 and we
3877 // should return *this
3878 if (structdim == 2)
3879 {
3880 Assert(i == 0,
3881 ExcMessage("You can only ask for quad zero if the "
3882 "current object is a quad itself."));
3883 return typename ::internal::DoFHandlerImplementation::
3884 Iterators<dim, spacedim>::cell_iterator(&this->get_triangulation(),
3885 this->level(),
3886 this->index(),
3887 &this->get_dof_handler());
3888 }
3889
3890 // otherwise we need to be in structdim>=3
3891 Assert(structdim > 2, ExcImpossibleInDim(structdim));
3892 Assert(dim > 2, ExcImpossibleInDim(dim));
3893
3894 // checking of 'i' happens in quad_index(i)
3895 return typename ::internal::DoFHandlerImplementation::
3896 Iterators<dim, spacedim, level_dof_access>::quad_iterator(
3897 this->tria,
3898 0, // only sub-objects are allowed, which have no level
3899 this->quad_index(i),
3900 this->dof_handler);
3901}
3902
3903
3904/*----------------- Functions: DoFAccessor<0,1,spacedim> --------------------*/
3905
3906
3907template <int spacedim, bool level_dof_access>
3909{
3910 Assert(false, ExcInvalidObject());
3911}
3912
3913
3914
3915template <int spacedim, bool level_dof_access>
3917 const Triangulation<1, spacedim> *tria,
3918 const typename TriaAccessor<0, 1, spacedim>::VertexKind vertex_kind,
3919 const unsigned int vertex_index,
3920 const DoFHandler<1, spacedim> *dof_handler)
3921 : BaseClass(tria, vertex_kind, vertex_index)
3922 , dof_handler(const_cast<DoFHandler<1, spacedim> *>(dof_handler))
3923{}
3924
3925
3926
3927template <int spacedim, bool level_dof_access>
3929 const Triangulation<1, spacedim> *tria,
3930 const int level,
3931 const int index,
3932 const DoFHandler<1, spacedim> *dof_handler)
3933 // This is the constructor signature for "ordinary" (non-vertex)
3934 // accessors and we shouldn't be calling it altogether. But it is also
3935 // the constructor that the default-constructor of TriaRawIterator
3936 // calls when default-constructing an iterator object. If so, this
3937 // happens with level==-2 and index==-2, and this is the only case we
3938 // would like to support. We do this by just forwarding to the
3939 // other constructor of this class, and then asserting the condition
3940 // on level and index.
3941 : DoFAccessor<0, 1, spacedim, level_dof_access>(
3942 tria,
3943 TriaAccessor<0, 1, spacedim>::interior_vertex,
3944 0U,
3945 dof_handler)
3946{
3947 (void)level;
3948 (void)index;
3949 Assert((tria == nullptr) && (level == -2) && (index == -2) &&
3950 (dof_handler == nullptr),
3951 ExcMessage(
3952 "This constructor can not be called for face iterators in 1d, "
3953 "except to default-construct iterator objects."));
3954}
3955
3956
3957
3958template <int spacedim, bool level_dof_access>
3959template <int structdim2, int dim2, int spacedim2>
3965
3966
3967
3968template <int spacedim, bool level_dof_access>
3969template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
3975
3976
3977
3978template <int spacedim, bool level_dof_access>
3979inline void
3982{
3983 Assert(dh != nullptr, ExcInvalidObject());
3984 this->dof_handler = dh;
3985}
3986
3987
3988
3989template <int spacedim, bool level_dof_access>
3990inline void
3992 const unsigned int /*i*/,
3993 const types::global_dof_index /*index*/,
3994 const types::fe_index /*fe_index*/) const
3995{
3997}
3998
3999
4000
4001template <int spacedim, bool level_dof_access>
4002inline const DoFHandler<1, spacedim> &
4004{
4005 return *this->dof_handler;
4006}
4007
4008
4009
4010template <int spacedim, bool level_dof_access>
4011inline void
4013 std::vector<types::global_dof_index> &dof_indices,
4014 const types::fe_index fe_index_) const
4015{
4016 const auto fe_index =
4018 fe_index_);
4019
4020 for (unsigned int i = 0; i < dof_indices.size(); ++i)
4021 dof_indices[i] = ::internal::DoFAccessorImplementation::
4022 Implementation::get_dof_index(*dof_handler,
4023 0,
4024 this->global_vertex_index,
4025 fe_index,
4026 i,
4027 std::integral_constant<int, 0>());
4028}
4029
4030
4031
4032template <int spacedim, bool level_dof_access>
4033inline void
4035 const int level,
4036 std::vector<types::global_dof_index> &dof_indices,
4037 const types::fe_index fe_index_) const
4038{
4039 const auto fe_index =
4041 fe_index_);
4042 (void)fe_index;
4044
4045 for (unsigned int i = 0; i < dof_indices.size(); ++i)
4046 dof_indices[i] =
4048 mg_vertex_dof_index(*dof_handler, level, this->global_vertex_index, i);
4049}
4050
4051
4052
4053template <int spacedim, bool level_dof_access>
4056 const unsigned int vertex,
4057 const unsigned int i,
4058 const types::fe_index fe_index_) const
4059{
4060 const auto fe_index =
4062 fe_index_);
4063
4064 (void)vertex;
4065 AssertIndexRange(vertex, 1);
4066 return ::internal::DoFAccessorImplementation::Implementation::
4067 get_dof_index(*dof_handler,
4068 0,
4069 this->global_vertex_index,
4070 fe_index,
4071 i,
4072 std::integral_constant<int, 0>());
4073}
4074
4075
4076
4077template <int spacedim, bool level_dof_access>
4080 const unsigned int i,
4081 const types::fe_index fe_index_) const
4082{
4083 const auto fe_index =
4085 fe_index_);
4086
4087 return ::internal::DoFAccessorImplementation::Implementation::
4088 get_dof_index(*this->dof_handler,
4089 0,
4090 this->vertex_index(0),
4091 fe_index,
4092 i,
4093 std::integral_constant<int, 0>());
4094}
4095
4096
4097
4098template <int spacedim, bool level_dof_access>
4099inline unsigned int
4104
4105
4106
4107template <int spacedim, bool level_dof_access>
4108inline types::fe_index
4110 const unsigned int /*n*/) const
4111{
4112 return 0;
4113}
4114
4115
4116
4117template <int spacedim, bool level_dof_access>
4118inline bool
4120 const types::fe_index /*fe_index*/) const
4121{
4123 return false;
4124}
4125
4126
4127
4128template <int spacedim, bool level_dof_access>
4129inline const FiniteElement<1, spacedim> &
4131 const types::fe_index fe_index) const
4132{
4133 Assert(this->dof_handler != nullptr, ExcInvalidObject());
4134 return dof_handler->get_fe(fe_index);
4135}
4136
4137
4138
4139template <int spacedim, bool level_dof_access>
4140inline void
4143{
4144 Assert(this->dof_handler != nullptr, ExcInvalidObject());
4145 BaseClass::copy_from(da);
4146}
4147
4148
4149
4150template <int spacedim, bool level_dof_access>
4151template <bool level_dof_access2>
4152inline void
4155{
4156 BaseClass::copy_from(a);
4157 set_dof_handler(a.dof_handler);
4158}
4159
4160
4161
4162template <int spacedim, bool level_dof_access>
4169
4170
4171
4172template <int spacedim, bool level_dof_access>
4173inline typename ::internal::DoFHandlerImplementation::
4174 Iterators<1, spacedim, level_dof_access>::line_iterator
4176 const unsigned int /*c*/) const
4177{
4179 return typename ::internal::DoFHandlerImplementation::
4180 Iterators<1, spacedim, level_dof_access>::line_iterator();
4181}
4182
4183
4184
4185template <int spacedim, bool level_dof_access>
4186inline typename ::internal::DoFHandlerImplementation::
4187 Iterators<1, spacedim, level_dof_access>::quad_iterator
4189 const unsigned int /*c*/) const
4190{
4192 return typename ::internal::DoFHandlerImplementation::
4193 Iterators<1, spacedim, level_dof_access>::quad_iterator();
4194}
4195
4196
4197
4198template <int spacedim, bool level_dof_access>
4199template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
4200inline bool
4203{
4204 Assert(structdim2 == 0, ExcCantCompareIterators());
4205 Assert(this->dof_handler == a.dof_handler, ExcCantCompareIterators());
4206 return (BaseClass::operator==(a));
4207}
4208
4209
4210
4211template <int spacedim, bool level_dof_access>
4212template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
4213inline bool
4216{
4217 Assert(structdim2 == 0, ExcCantCompareIterators());
4218 Assert(this->dof_handler == a.dof_handler, ExcCantCompareIterators());
4219 return (BaseClass::operator!=(a));
4220}
4221
4222
4223
4224/*------------------------- Functions: DoFCellAccessor -----------------------*/
4225
4226
4227namespace internal
4228{
4229 namespace DoFCellAccessorImplementation
4230 {
4236 {
4241 template <int dim, int spacedim, bool level_dof_access>
4242 static types::fe_index
4245 {
4246 if (accessor.dof_handler->hp_capability_enabled == false)
4248
4249 Assert(accessor.dof_handler != nullptr,
4250 (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
4251 Assert(static_cast<unsigned int>(accessor.level()) <
4252 accessor.dof_handler->hp_cell_future_fe_indices.size(),
4253 ExcMessage("DoFHandler not initialized"));
4254
4255 return accessor.dof_handler
4256 ->hp_cell_active_fe_indices[accessor.level()][accessor.present_index];
4257 }
4258
4259
4260
4265 template <int dim, int spacedim, bool level_dof_access>
4266 static void
4269 const types::fe_index i)
4270 {
4271 if (accessor.dof_handler->hp_capability_enabled == false)
4272 {
4274 return;
4275 }
4276
4277 Assert(accessor.dof_handler != nullptr,
4278 (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
4279 Assert(static_cast<unsigned int>(accessor.level()) <
4280 accessor.dof_handler->hp_cell_future_fe_indices.size(),
4281 ExcMessage("DoFHandler not initialized"));
4283 ExcMessage("Invalid finite element index."));
4284
4285 accessor.dof_handler
4286 ->hp_cell_active_fe_indices[accessor.level()]
4287 [accessor.present_index] = i;
4288 }
4289
4290
4291
4296 template <int dim, int spacedim, bool level_dof_access>
4297 static types::fe_index
4300 {
4301 if (accessor.dof_handler->hp_capability_enabled == false)
4303
4304 Assert(accessor.dof_handler != nullptr,
4305 (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
4306 Assert(static_cast<unsigned int>(accessor.level()) <
4307 accessor.dof_handler->hp_cell_future_fe_indices.size(),
4308 ExcMessage("DoFHandler not initialized"));
4309
4310 if (future_fe_index_set(accessor))
4311 return accessor.dof_handler
4312 ->hp_cell_future_fe_indices[accessor.level()]
4313 [accessor.present_index];
4314 else
4315 return accessor.dof_handler
4316 ->hp_cell_active_fe_indices[accessor.level()]
4317 [accessor.present_index];
4318 }
4319
4320
4325 template <int dim, int spacedim, bool level_dof_access>
4326 static void
4329 const types::fe_index i)
4330 {
4331 if (accessor.dof_handler->hp_capability_enabled == false)
4332 {
4334 return;
4335 }
4336
4337 Assert(accessor.dof_handler != nullptr,
4338 (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
4339 Assert(static_cast<unsigned int>(accessor.level()) <
4340 accessor.dof_handler->hp_cell_future_fe_indices.size(),
4341 ExcMessage("DoFHandler not initialized"));
4343 ExcMessage("Invalid finite element index."));
4344
4345 accessor.dof_handler
4346 ->hp_cell_future_fe_indices[accessor.level()]
4347 [accessor.present_index] = i;
4348 }
4349
4350
4351
4356 template <int dim, int spacedim, bool level_dof_access>
4357 static bool
4360 {
4361 if (accessor.dof_handler->hp_capability_enabled == false)
4362 return false;
4363
4364 Assert(accessor.dof_handler != nullptr,
4365 (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
4366 Assert(static_cast<unsigned int>(accessor.level()) <
4367 accessor.dof_handler->hp_cell_future_fe_indices.size(),
4368 ExcMessage("DoFHandler not initialized"));
4369
4370 return accessor.dof_handler
4371 ->hp_cell_future_fe_indices[accessor.level()]
4372 [accessor.present_index] !=
4374 }
4375
4376
4377
4382 template <int dim, int spacedim, bool level_dof_access>
4383 static void
4386 {
4387 if (accessor.dof_handler->hp_capability_enabled == false)
4388 return;
4389
4390 Assert(accessor.dof_handler != nullptr,
4391 (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
4392 Assert(static_cast<unsigned int>(accessor.level()) <
4393 accessor.dof_handler->hp_cell_future_fe_indices.size(),
4394 ExcMessage("DoFHandler not initialized"));
4395
4396 accessor.dof_handler
4397 ->hp_cell_future_fe_indices[accessor.level()]
4398 [accessor.present_index] =
4400 }
4401 };
4402 } // namespace DoFCellAccessorImplementation
4403} // namespace internal
4404
4405
4406
4407template <int dimension_, int space_dimension_, bool level_dof_access>
4410 const int level,
4411 const int index,
4412 const AccessorData *local_data)
4413 : DoFAccessor<dimension_, dimension_, space_dimension_, level_dof_access>(
4414 tria,
4415 level,
4416 index,
4417 local_data)
4418{}
4419
4420
4421
4422template <int dimension_, int space_dimension_, bool level_dof_access>
4423template <int structdim2, int dim2, int spacedim2>
4429
4430
4431
4432template <int dimension_, int space_dimension_, bool level_dof_access>
4433template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
4439
4440
4441
4442template <int dimension_, int space_dimension_, bool level_dof_access>
4443inline TriaIterator<
4446 const unsigned int i) const
4447{
4449 q(this->tria,
4450 this->neighbor_level(i),
4451 this->neighbor_index(i),
4452 this->dof_handler);
4453
4454 if constexpr (running_in_debug_mode())
4455 {
4457 Assert(q->used(), ExcInternalError());
4458 }
4459 return q;
4460}
4461
4462
4463
4464template <int dimension_, int space_dimension_, bool level_dof_access>
4465inline TriaIterator<
4468 const unsigned int i) const
4469{
4471 q(this->tria, this->level() + 1, this->child_index(i), this->dof_handler);
4472
4473 if constexpr (running_in_debug_mode())
4474 {
4476 Assert(q->used(), ExcInternalError());
4477 }
4478 return q;
4479}
4480
4481
4482
4483template <int dimension_, int space_dimension_, bool level_dof_access>
4484inline boost::container::small_vector<
4488 child_iterators() const
4489{
4490 boost::container::small_vector<
4494 child_iterators(this->n_children());
4495
4496 for (unsigned int i = 0; i < this->n_children(); ++i)
4497 child_iterators[i] = this->child(i);
4498
4499 return child_iterators;
4500}
4501
4502
4503
4504template <int dimension_, int space_dimension_, bool level_dof_access>
4505inline TriaIterator<
4508{
4510 q(this->tria, this->level() - 1, this->parent_index(), this->dof_handler);
4511
4512 return q;
4513}
4514
4515
4516
4517namespace internal
4518{
4519 namespace DoFCellAccessorImplementation
4520 {
4521 template <int dim, int spacedim, bool level_dof_access>
4522 inline TriaIterator<
4523 ::DoFAccessor<dim - 1, dim, spacedim, level_dof_access>>
4525 const ::DoFCellAccessor<dim, spacedim, level_dof_access> &cell,
4526 const unsigned int i,
4527 const std::integral_constant<int, 1>)
4528 {
4530 &cell.get_triangulation(),
4531 ((i == 0) && cell.at_boundary(0) ?
4533 ((i == 1) && cell.at_boundary(1) ?
4536 cell.vertex_index(i),
4537 &cell.get_dof_handler());
4538 return ::TriaIterator<
4540 }
4541
4542
4543 template <int dim, int spacedim, bool level_dof_access>
4544 inline TriaIterator<
4545 ::DoFAccessor<dim - 1, dim, spacedim, level_dof_access>>
4547 const ::DoFCellAccessor<dim, spacedim, level_dof_access> &cell,
4548 const unsigned int i,
4549 const std::integral_constant<int, 2>)
4550 {
4551 return cell.line(i);
4552 }
4553
4554
4555 template <int dim, int spacedim, bool level_dof_access>
4556 inline TriaIterator<
4557 ::DoFAccessor<dim - 1, dim, spacedim, level_dof_access>>
4559 const ::DoFCellAccessor<dim, spacedim, level_dof_access> &cell,
4560 const unsigned int i,
4561 const std::integral_constant<int, 3>)
4562 {
4563 return cell.quad(i);
4564 }
4565 } // namespace DoFCellAccessorImplementation
4566} // namespace internal
4567
4568
4569
4570template <int dimension_, int space_dimension_, bool level_dof_access>
4571inline typename DoFCellAccessor<dimension_,
4572 space_dimension_,
4573 level_dof_access>::face_iterator
4575 const unsigned int i) const
4576{
4577 AssertIndexRange(i, this->n_faces());
4578
4579 return ::internal::DoFCellAccessorImplementation::get_face(
4580 *this, i, std::integral_constant<int, dimension_>());
4581}
4582
4583
4584
4585template <int dimension_, int space_dimension_, bool level_dof_access>
4586inline boost::container::small_vector<
4588 face_iterator,
4591 face_iterators() const
4592{
4593 boost::container::small_vector<
4595 face_iterator,
4597 face_iterators(this->n_faces());
4598
4599 for (const unsigned int i : this->face_indices())
4600 face_iterators[i] =
4602 *this, i, std::integral_constant<int, dimension_>());
4603
4604 return face_iterators;
4605}
4606
4607
4608
4609template <int dimension_, int space_dimension_, bool level_dof_access>
4610inline void
4613 std::vector<types::global_dof_index> &dof_indices) const
4614{
4615 if (level_dof_access)
4616 get_mg_dof_indices(dof_indices);
4617 else
4618 get_dof_indices(dof_indices);
4619}
4620
4621
4622
4623template <int dimension_, int space_dimension_, bool level_dof_access>
4624template <class InputVector, typename number>
4625inline void
4627 const InputVector &values,
4628 Vector<number> &local_values) const
4629{
4630 get_dof_values(values, local_values.begin(), local_values.end());
4631}
4632
4633
4634
4635template <int dimension_, int space_dimension_, bool level_dof_access>
4636template <typename Number, typename ForwardIterator>
4637inline void
4639 const ReadVector<Number> &values,
4640 ForwardIterator local_values_begin,
4641 ForwardIterator local_values_end) const
4642{
4643 (void)local_values_end;
4644 Assert(this->is_artificial() == false,
4645 ExcMessage("Can't ask for DoF indices on artificial cells."));
4646 Assert(this->is_active(), ExcMessage("Cell must be active."));
4647 Assert(this->dof_handler != nullptr, typename BaseClass::ExcInvalidObject());
4648
4649 Assert(static_cast<unsigned int>(local_values_end - local_values_begin) ==
4650 this->get_fe().n_dofs_per_cell(),
4652 Assert(values.size() == this->get_dof_handler().n_dofs(),
4654
4656 dof_indices(this->get_fe().n_dofs_per_cell());
4658 *this, dof_indices, this->active_fe_index());
4659
4660 boost::container::small_vector<Number, 27> values_temp(local_values_end -
4661 local_values_begin);
4662 auto view = make_array_view(values_temp.begin(), values_temp.end());
4663 values.extract_subvector_to(make_array_view(dof_indices.begin(),
4664 dof_indices.end()),
4665 view);
4666 using view_type = std::remove_reference_t<decltype(*local_values_begin)>;
4667 ArrayView<view_type> values_view2(&*local_values_begin,
4668 local_values_end - local_values_begin);
4669 std::copy(values_temp.begin(), values_temp.end(), values_view2.begin());
4670}
4671
4672
4673
4674template <int dimension_, int space_dimension_, bool level_dof_access>
4675template <class InputVector, typename ForwardIterator>
4676inline void
4679 const InputVector &values,
4680 ForwardIterator local_values_begin,
4681 ForwardIterator local_values_end) const
4682{
4683 Assert(this->is_artificial() == false,
4684 ExcMessage("Can't ask for DoF indices on artificial cells."));
4685 Assert(this->is_active(), ExcMessage("Cell must be active."));
4686
4687 Assert(static_cast<unsigned int>(local_values_end - local_values_begin) ==
4688 this->get_fe().n_dofs_per_cell(),
4690 Assert(values.size() == this->get_dof_handler().n_dofs(),
4692
4693
4695 dof_indices(this->get_fe().n_dofs_per_cell());
4697 *this, dof_indices, this->active_fe_index());
4698
4699 constraints.get_dof_values(values,
4700 dof_indices.data(),
4701 local_values_begin,
4702 local_values_end);
4703}
4704
4705
4706
4707template <int dimension_, int space_dimension_, bool level_dof_access>
4708template <class OutputVector, typename number>
4709inline void
4711 const Vector<number> &local_values,
4712 OutputVector &values) const
4713{
4714 Assert(this->is_artificial() == false,
4715 ExcMessage("Can't ask for DoF indices on artificial cells."));
4716 Assert(this->is_active(), ExcMessage("Cell must be active."));
4717
4718 Assert(static_cast<unsigned int>(local_values.size()) ==
4719 this->get_fe().n_dofs_per_cell(),
4721 Assert(values.size() == this->get_dof_handler().n_dofs(),
4723
4724
4725 Assert(this->dof_handler != nullptr, typename BaseClass::ExcInvalidObject());
4727 dof_indices(this->get_fe().n_dofs_per_cell());
4729 *this, dof_indices, this->active_fe_index());
4730
4731 for (unsigned int i = 0; i < this->get_fe().n_dofs_per_cell(); ++i)
4733 dof_indices[i],
4734 values);
4735}
4736
4737
4738
4739template <int dimension_, int space_dimension_, bool level_dof_access>
4742{
4743 Assert(this->dof_handler != nullptr, typename BaseClass::ExcInvalidObject());
4744 Assert((this->dof_handler->hp_capability_enabled == false) ||
4745 this->is_active(),
4746 ExcMessage(
4747 "For DoFHandler objects in hp-mode, finite elements are only "
4748 "associated with active cells. Consequently, you can not ask "
4749 "for the active finite element on cells with children."));
4750
4751 const auto &fe = this->dof_handler->get_fe(active_fe_index());
4752
4753 Assert(this->reference_cell() == fe.reference_cell(),
4754 internal::ExcNonMatchingReferenceCellTypes(this->reference_cell(),
4755 fe.reference_cell()));
4756
4757 return fe;
4758}
4759
4760
4761
4762template <int dimension_, int space_dimension_, bool level_dof_access>
4763inline types::fe_index
4765 active_fe_index() const
4766{
4767 Assert((this->dof_handler->hp_capability_enabled == false) ||
4768 this->is_active(),
4769 ExcMessage(
4770 "You can not ask for the active FE index on a cell that has "
4771 "children because no degrees of freedom are assigned "
4772 "to this cell and, consequently, no finite element "
4773 "is associated with it."));
4774 Assert((this->dof_handler->hp_capability_enabled == false) ||
4775 (this->is_locally_owned() || this->is_ghost()),
4776 ExcMessage("You can only query active FE index information on cells "
4777 "that are either locally owned or (after distributing "
4778 "degrees of freedom) are ghost cells."));
4779
4780 return ::internal::DoFCellAccessorImplementation::Implementation::
4781 active_fe_index(*this);
4782}
4783
4784
4785
4786template <int dimension_, int space_dimension_, bool level_dof_access>
4787inline void
4790{
4791 Assert((this->dof_handler->hp_capability_enabled == false) ||
4792 this->is_active(),
4793 ExcMessage("You can not set the active FE index on a cell that has "
4794 "children because no degrees of freedom will be assigned "
4795 "to this cell."));
4796
4797 Assert((this->dof_handler->hp_capability_enabled == false) ||
4798 this->is_locally_owned(),
4799 ExcMessage("You can only set active FE index information on cells "
4800 "that are locally owned. On ghost cells, this information "
4801 "will automatically be propagated from the owning process "
4802 "of that cell, and there is no information at all on "
4803 "artificial cells."));
4804
4806 set_active_fe_index(*this, i);
4807}
4808
4809
4810
4811template <int dimension_, int space_dimension_, bool level_dof_access>
4814 const
4815{
4816 Assert(this->dof_handler != nullptr, typename BaseClass::ExcInvalidObject());
4817 Assert((this->dof_handler->hp_capability_enabled == false) ||
4818 this->is_active(),
4819 ExcMessage(
4820 "For DoFHandler objects in hp-mode, finite elements are only "
4821 "associated with active cells. Consequently, you can not ask "
4822 "for the future finite element on cells with children."));
4823
4824 return this->dof_handler->get_fe(future_fe_index());
4825}
4826
4827
4828
4829template <int dimension_, int space_dimension_, bool level_dof_access>
4830inline types::fe_index
4832 future_fe_index() const
4833{
4834 Assert((this->dof_handler->hp_capability_enabled == false) ||
4835 (this->has_children() == false),
4836 ExcMessage(
4837 "You can not ask for the future FE index on a cell that has "
4838 "children because no degrees of freedom are assigned "
4839 "to this cell and, consequently, no finite element "
4840 "is associated with it."));
4841 Assert((this->dof_handler->hp_capability_enabled == false) ||
4842 (this->is_locally_owned()),
4843 ExcMessage("You can only query future FE index information on cells "
4844 "that are locally owned."));
4845
4846 return ::internal::DoFCellAccessorImplementation::Implementation::
4847 future_fe_index(*this);
4848}
4849
4850
4851
4852template <int dimension_, int space_dimension_, bool level_dof_access>
4853inline void
4856{
4857 Assert((this->dof_handler->hp_capability_enabled == false) ||
4858 (this->has_children() == false),
4859 ExcMessage("You can not set the future FE index on a cell that has "
4860 "children because no degrees of freedom will be assigned "
4861 "to this cell."));
4862
4863 Assert((this->dof_handler->hp_capability_enabled == false) ||
4864 this->is_locally_owned(),
4865 ExcMessage("You can only set future FE index information on cells "
4866 "that are locally owned."));
4867
4869 set_future_fe_index(*this, i);
4870}
4871
4872
4873
4874template <int dimension_, int space_dimension_, bool level_dof_access>
4875inline bool
4878{
4879 Assert((this->dof_handler->hp_capability_enabled == false) ||
4880 (this->has_children() == false),
4881 ExcMessage(
4882 "You can not ask for the future FE index on a cell that has "
4883 "children because no degrees of freedom are assigned "
4884 "to this cell and, consequently, no finite element "
4885 "is associated with it."));
4886 Assert((this->dof_handler->hp_capability_enabled == false) ||
4887 (this->is_locally_owned()),
4888 ExcMessage("You can only query future FE index information on cells "
4889 "that are locally owned."));
4890
4891 return ::internal::DoFCellAccessorImplementation::Implementation::
4892 future_fe_index_set(*this);
4893}
4894
4895
4896
4897template <int dimension_, int space_dimension_, bool level_dof_access>
4898inline void
4901{
4902 Assert((this->dof_handler->hp_capability_enabled == false) ||
4903 (this->has_children() == false),
4904 ExcMessage(
4905 "You can not ask for the future FE index on a cell that has "
4906 "children because no degrees of freedom are assigned "
4907 "to this cell and, consequently, no finite element "
4908 "is associated with it."));
4909 Assert((this->dof_handler->hp_capability_enabled == false) ||
4910 (this->is_locally_owned()),
4911 ExcMessage("You can only query future FE index information on cells "
4912 "that are locally owned."));
4913
4916}
4917
4918
4919
4920template <int dimension_, int space_dimension_, bool level_dof_access>
4921template <typename number, typename OutputVector>
4922inline void
4925 OutputVector &global_destination) const
4926{
4927 this->distribute_local_to_global(local_source.begin(),
4928 local_source.end(),
4929 global_destination);
4930}
4931
4932
4933
4934template <int dimension_, int space_dimension_, bool level_dof_access>
4935template <typename ForwardIterator, typename OutputVector>
4936inline void
4938 distribute_local_to_global(ForwardIterator local_source_begin,
4939 ForwardIterator local_source_end,
4940 OutputVector &global_destination) const
4941{
4942 Assert(this->dof_handler != nullptr,
4943 (typename std::decay_t<decltype(*this)>::ExcInvalidObject()));
4944 Assert(static_cast<unsigned int>(local_source_end - local_source_begin) ==
4945 this->get_fe().n_dofs_per_cell(),
4946 (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
4947 Assert(this->dof_handler->n_dofs() == global_destination.size(),
4948 (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
4949
4950 Assert(!this->has_children(), ExcMessage("Cell must be active"));
4951
4952 Assert(
4953 internal::ArrayViewHelper::is_contiguous(local_source_begin,
4954 local_source_end),
4955 ExcMessage(
4956 "This function can not be called with iterator types that do not point to contiguous memory."));
4957
4958 const unsigned int n_dofs = local_source_end - local_source_begin;
4959
4961 dof_indices(n_dofs);
4963 *this, dof_indices, this->active_fe_index());
4964
4965 // distribute cell vector
4966 global_destination.add(n_dofs, dof_indices.data(), &(*local_source_begin));
4967}
4968
4969
4970
4971template <int dimension_, int space_dimension_, bool level_dof_access>
4972template <typename ForwardIterator, typename OutputVector>
4973inline void
4977 ForwardIterator local_source_begin,
4978 ForwardIterator local_source_end,
4979 OutputVector &global_destination) const
4980{
4981 Assert(this->dof_handler != nullptr,
4982 (typename std::decay_t<decltype(*this)>::ExcInvalidObject()));
4983 Assert(local_source_end - local_source_begin ==
4984 this->get_fe().n_dofs_per_cell(),
4985 (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
4986 Assert(this->dof_handler->n_dofs() == global_destination.size(),
4987 (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
4988
4989 Assert(!this->has_children(), ExcMessage("Cell must be active."));
4990
4992 dof_indices(this->get_fe().n_dofs_per_cell());
4994 *this, dof_indices, this->active_fe_index());
4995
4996 // distribute cell vector
4997 constraints.distribute_local_to_global(local_source_begin,
4998 local_source_end,
4999 dof_indices.data(),
5000 global_destination);
5001}
5002
5003
5004
5005template <int dimension_, int space_dimension_, bool level_dof_access>
5006template <typename number, typename OutputMatrix>
5007inline void
5010 OutputMatrix &global_destination) const
5011{
5012 Assert(this->dof_handler != nullptr,
5013 (typename std::decay_t<decltype(*this)>::ExcInvalidObject()));
5014 Assert(local_source.m() == this->get_fe().n_dofs_per_cell(),
5015 (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
5016 Assert(local_source.n() == this->get_fe().n_dofs_per_cell(),
5017 (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
5018 Assert(this->dof_handler->n_dofs() == global_destination.m(),
5019 (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
5020 Assert(this->dof_handler->n_dofs() == global_destination.n(),
5021 (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
5022
5023 Assert(!this->has_children(), ExcMessage("Cell must be active."));
5024
5025 const unsigned int n_dofs = local_source.m();
5026
5028 dof_indices(n_dofs);
5030 *this, dof_indices, this->active_fe_index());
5031
5032 // distribute cell matrix
5033 for (unsigned int i = 0; i < n_dofs; ++i)
5034 global_destination.add(dof_indices[i],
5035 n_dofs,
5036 dof_indices.data(),
5037 &local_source(i, 0));
5038}
5039
5040
5041
5042template <int dimension_, int space_dimension_, bool level_dof_access>
5043template <typename number, typename OutputMatrix, typename OutputVector>
5044inline void
5047 const Vector<number> &local_vector,
5048 OutputMatrix &global_matrix,
5049 OutputVector &global_vector) const
5050{
5051 Assert(this->dof_handler != nullptr,
5052 (typename std::decay_t<decltype(*this)>::ExcInvalidObject()));
5053 Assert(local_matrix.m() == this->get_fe().n_dofs_per_cell(),
5054 (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
5055 Assert(local_matrix.n() == this->get_fe().n_dofs_per_cell(),
5056 (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
5057 Assert(this->dof_handler->n_dofs() == global_matrix.m(),
5058 (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
5059 Assert(this->dof_handler->n_dofs() == global_matrix.n(),
5060 (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
5061 Assert(local_vector.size() == this->get_fe().n_dofs_per_cell(),
5062 (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
5063 Assert(this->dof_handler->n_dofs() == global_vector.size(),
5064 (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
5065
5066 Assert(!this->has_children(), ExcMessage("Cell must be active."));
5067
5068 const unsigned int n_dofs = this->get_fe().n_dofs_per_cell();
5070 dof_indices(n_dofs);
5072 *this, dof_indices, this->active_fe_index());
5073
5074 // distribute cell matrices
5075 for (unsigned int i = 0; i < n_dofs; ++i)
5076 {
5077 global_matrix.add(dof_indices[i],
5078 n_dofs,
5079 dof_indices.data(),
5080 &local_matrix(i, 0));
5081 global_vector(dof_indices[i]) += local_vector(i);
5082 }
5083}
5084
5085
5086
5088
5089
5090#endif
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:954
void get_dof_values(const VectorType &global_vector, ForwardIteratorInd local_indices_begin, ForwardIteratorVec local_vector_begin, ForwardIteratorVec local_vector_end) const
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
iterator begin() const
Definition array_view.h:707
DoFAccessor< 0, 1, spacedim, level_dof_access > & operator=(const DoFAccessor< 0, 1, spacedim, level_dof_access > &da)=delete
DoFAccessor< 0, 1, spacedim, level_dof_access > & operator=(DoFAccessor< 0, 1, spacedim, level_dof_access > &&) noexcept=default
DoFAccessor(DoFAccessor< 0, 1, spacedim, level_dof_access > &&)=default
DoFAccessor(const DoFAccessor< 0, 1, spacedim, level_dof_access > &)=default
typename::internal::DoFHandlerImplementation::Iterators< dim, spacedim, level_dof_access >::line_iterator line(const unsigned int i) const
friend struct ::internal::DoFHandlerImplementation::Policy::Implementation
static constexpr unsigned int space_dimension
const DoFHandler< dim, spacedim > & get_dof_handler() const
void set_mg_dof_index(const int level, const unsigned int i, const types::global_dof_index index) const
friend class DoFAccessor
DoFAccessor< structdim, dim, spacedim, level_dof_access > & operator=(const DoFAccessor< structdim, dim, spacedim, level_dof_access > &da)=delete
types::fe_index nth_active_fe_index(const unsigned int n) const
DoFAccessor(DoFAccessor< structdim, dim, spacedim, level_dof_access > &&)
DoFAccessor(const DoFAccessor< structdim, dim, spacedim, level_dof_access2 > &)
TriaIterator< DoFAccessor< structdim, dim, spacedim, level_dof_access > > child(const unsigned int c) const
void set_mg_vertex_dof_index(const int level, const unsigned int vertex, const unsigned int i, const types::global_dof_index index, const types::fe_index fe_index=numbers::invalid_fe_index) const
~DoFAccessor()=default
static constexpr unsigned int dimension
void set_dof_handler(DoFHandler< dim, spacedim > *dh)
DoFAccessor(const DoFAccessor< structdim2, dim2, spacedim2, level_dof_access2 > &)
DoFAccessor< structdim, dim, spacedim, level_dof_access > & operator=(DoFAccessor< structdim, dim, spacedim, level_dof_access > &&)
std::set< types::fe_index > get_active_fe_indices() const
DoFAccessor(const DoFAccessor< structdim, dim, spacedim, level_dof_access > &)=default
void set_mg_dof_indices(const int level, const std::vector< types::global_dof_index > &dof_indices, const types::fe_index fe_index=numbers::invalid_fe_index)
void get_mg_dof_indices(const int level, std::vector< types::global_dof_index > &dof_indices, const types::fe_index fe_index=numbers::invalid_fe_index) const
types::global_dof_index mg_vertex_dof_index(const int level, const unsigned int vertex, const unsigned int i, const types::fe_index fe_index=numbers::invalid_fe_index) const
DoFAccessor(const Triangulation< dim, spacedim > *tria, const int level, const int index, const DoFHandler< dim, spacedim > *dof_handler)
bool operator!=(const DoFAccessor< structdim2, dim2, spacedim2, level_dof_access2 > &) const
void copy_from(const TriaAccessorBase< structdim, dim, spacedim > &da)
types::global_dof_index dof_index(const unsigned int i, const types::fe_index fe_index=numbers::invalid_fe_index) const
typename::internal::DoFHandlerImplementation::Iterators< dim, spacedim, level_dof_access >::quad_iterator quad(const unsigned int i) const
void copy_from(const DoFAccessor< structdim, dim, spacedim, level_dof_access2 > &a)
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index fe_index) const
void get_dof_indices(std::vector< types::global_dof_index > &dof_indices, const types::fe_index fe_index=numbers::invalid_fe_index) const
bool operator==(const DoFAccessor< structdim2, dim2, spacedim2, level_dof_access2 > &) const
types::global_dof_index mg_dof_index(const int level, const unsigned int i) const
void set_dof_index(const unsigned int i, const types::global_dof_index index, const types::fe_index fe_index=numbers::invalid_fe_index) const
unsigned int n_active_fe_indices() const
static bool is_level_cell()
typename ::internal::DoFAccessorImplementation::Inheritance< structdim, dimension, space_dimension >::BaseClass BaseClass
types::global_dof_index vertex_dof_index(const unsigned int vertex, const unsigned int i, const types::fe_index fe_index=numbers::invalid_fe_index) const
DoFHandler< dim, spacedim > * dof_handler
DoFAccessor(const InvalidAccessor< structdim2, dim2, spacedim2 > &)
bool fe_index_is_active(const types::fe_index fe_index) const
boost::container::small_vector< TriaIterator< DoFCellAccessor< dimension_, space_dimension_, level_dof_access > >, GeometryInfo< dimension_ >::max_children_per_cell > child_iterators() const
DoFCellAccessor(DoFCellAccessor< dimension_, space_dimension_, level_dof_access > &&)=default
void get_active_or_mg_dof_indices(std::vector< types::global_dof_index > &dof_indices) const
void get_dof_values(const InputVector &values, Vector< number > &local_values) const
boost::container::small_vector< face_iterator, GeometryInfo< dimension_ >::faces_per_cell > face_iterators() const
types::fe_index future_fe_index() const
DoFCellAccessor< dimension_, space_dimension_, level_dof_access > & operator=(const DoFCellAccessor< dimension_, space_dimension_, level_dof_access > &da)=delete
const FiniteElement< dimension_, space_dimension_ > & get_fe() const
TriaIterator< DoFCellAccessor< dimension_, space_dimension_, level_dof_access > > child(const unsigned int i) const
DoFCellAccessor(const Triangulation< dimension_, space_dimension_ > *tria, const int level, const int index, const AccessorData *local_data)
void set_future_fe_index(const types::fe_index i) const
void distribute_local_to_global(const Vector< number > &local_source, OutputVector &global_destination) const
const FiniteElement< dimension_, space_dimension_ > & get_future_fe() const
TriaIterator< DoFCellAccessor< dimension_, space_dimension_, level_dof_access > > neighbor(const unsigned int i) const
void set_active_fe_index(const types::fe_index i) const
face_iterator face(const unsigned int i) const
void clear_future_fe_index() const
void set_dof_values(const Vector< number > &local_values, OutputVector &values) const
DoFCellAccessor(const DoFCellAccessor< dimension_, space_dimension_, level_dof_access > &)=default
~DoFCellAccessor()=default
bool future_fe_index_set() const
DoFCellAccessor< dimension_, space_dimension_, level_dof_access > & operator=(DoFCellAccessor< dimension_, space_dimension_, level_dof_access > &&)=default
types::fe_index active_fe_index() const
TriaIterator< DoFCellAccessor< dimension_, space_dimension_, level_dof_access > > parent() const
std::vector< std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel< dim > > > mg_levels
std::vector< std::vector< types::fe_index > > hp_cell_future_fe_indices
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
std::vector< std::vector< types::fe_index > > hp_cell_active_fe_indices
std::vector< MGVertexDoFs > mg_vertex_dofs
std::vector< std::array< std::vector< types::global_dof_index >, dim+1 > > object_dof_indices
const Triangulation< dim, spacedim > & get_triangulation() const
std::vector< std::array< std::vector< offset_type >, dim+1 > > object_dof_ptr
std::unique_ptr<::internal::DoFHandlerImplementation::DoFFaces< dim > > mg_faces
bool hp_capability_enabled
bool has_hp_capabilities() const
std::array< std::vector< offset_type >, dim+1 > hp_object_fe_ptr
std::array< std::vector< types::fe_index >, dim+1 > hp_object_fe_indices
DoFInvalidAccessor(const void *parent=nullptr, const int level=-1, const int index=-1, const AccessorData *local_data=nullptr)
typename InvalidAccessor< structdim, dim, spacedim >::AccessorData AccessorData
unsigned int n_dofs_per_vertex() const
size_type n() const
size_type m() const
void compress() const
Definition index_set.h:1795
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1842
void import_elements(const ::Vector< Number > &vec, VectorOperation::values operation, const std::shared_ptr< const Utilities::MPI::CommunicationPatternBase > &communication_pattern={})
TriaIterator< TriaAccessor< structdim, dim, spacedim > > child(const unsigned int i) const
IteratorState::IteratorStates state() const
virtual size_type size() const override
iterator end()
iterator begin()
#define DEAL_II_ALWAYS_INLINE
Definition config.h:161
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_RESTRICT
Definition config.h:162
constexpr bool running_in_debug_mode()
Definition config.h:73
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
Point< 2 > second
Definition grid_out.cc:4633
unsigned int level
Definition grid_out.cc:4635
static ::ExceptionBase & ExcVectorNotEmpty()
#define DeclException0(Exception0)
static ::ExceptionBase & ExcNonMatchingReferenceCellTypes(ReferenceCell arg1, ReferenceCell arg2)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
static ::ExceptionBase & ExcNotActive()
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcMatrixDoesNotMatch()
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
static ::ExceptionBase & ExcInvalidObject()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcCantCompareIterators()
static ::ExceptionBase & ExcVectorDoesNotMatch()
static ::ExceptionBase & ExcMessage(std::string arg1)
@ past_the_end
Iterator reached end of container.
Definition hp.h:117
types::fe_index get_fe_index_or_default(const DoFAccessor< structdim, dim, spacedim, level_dof_access > &cell, const types::fe_index fe_index)
void get_cell_dof_indices(const ::DoFCellAccessor< dim, spacedim, level_dof_access > &accessor, Implementation::dof_index_vector_type &dof_indices, const unsigned int fe_index)
TriaIterator< ::DoFAccessor< dim - 1, dim, spacedim, level_dof_access > > get_face(const ::DoFCellAccessor< dim, spacedim, level_dof_access > &cell, const unsigned int i, const std::integral_constant< int, 1 >)
constexpr types::global_dof_index invalid_dof_index
Definition types.h:269
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:352
constexpr types::fe_index invalid_fe_index
Definition types.h:260
STL namespace.
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition types.h:32
unsigned short int fe_index
Definition types.h:72
void process_vertex_dofs(DoFHandler< dim, spacedim > &dof_handler, const unsigned int vertex_index, const types::fe_index fe_index, types::global_dof_index *&dof_indices_ptr, const DoFProcessor &dof_processor) const
void process_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int obj_level, const unsigned int obj_index, const types::fe_index fe_index, const DoFMapping &mapping, const std::integral_constant< int, structdim >, types::global_dof_index *&dof_indices_ptr, const DoFProcessor &dof_processor) const
void process_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int, const unsigned int obj_index, const types::fe_index fe_index, const DoFMapping &mapping, const std::integral_constant< int, structdim >, types::global_dof_index *&dof_indices_ptr, const DoFProcessor &dof_processor) const
void process_vertex_dofs(DoFHandler< dim, spacedim > &dof_handler, const unsigned int vertex_index, const types::fe_index, types::global_dof_index *&dof_indices_ptr, const DoFProcessor &dof_processor) const
static types::global_dof_index * get_array_ptr(const std::tuple<> &)
static void process_dof_indices(const ::DoFAccessor< structdim, dim, spacedim, level_dof_access > &accessor, const DoFIndicesType &const_dof_indices, const types::fe_index fe_index_, const DoFOperation &dof_operation, const DoFProcessor &dof_processor, const bool count_level_dofs)
static void get_dof_indices(const ::DoFAccessor< structdim, dim, spacedim, level_dof_access > &accessor, std::vector< types::global_dof_index > &dof_indices, const types::fe_index fe_index)
static void extract_subvector_to(const LinearAlgebra::TpetraWrappers::Vector< Number, MemorySpace > &values, const types::global_dof_index *cache_begin, const types::global_dof_index *cache_end, ForwardIterator local_values_begin)
static std::set< types::fe_index > get_active_fe_indices(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int obj_level, const unsigned int obj_index, const std::integral_constant< int, structdim > &t)
static types::global_dof_index * get_array_ptr(const ArrayType &array)
static unsigned int get_array_length(const std::tuple<> &)
static unsigned int n_active_fe_indices(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int obj_level, const unsigned int obj_index, const std::integral_constant< int, structdim > &)
static void set_dof_index(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int obj_level, const unsigned int obj_index, const types::fe_index fe_index, const unsigned int local_index, const std::integral_constant< int, structdim > &dd, const types::global_dof_index global_index)
static void extract_subvector_to(const LinearAlgebra::EpetraWrappers::Vector &values, const types::global_dof_index *cache_begin, const types::global_dof_index *cache_end, ForwardIterator local_values_begin)
static void get_mg_dof_indices(const ::DoFAccessor< structdim, dim, spacedim, level_dof_access > &accessor, const int level, std::vector< types::global_dof_index > &dof_indices, const types::fe_index fe_index)
static bool fe_index_is_active(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int obj_level, const unsigned int obj_index, const types::fe_index fe_index, const std::integral_constant< int, structdim > &)
static types::global_dof_index & get_mg_dof_index(const DoFHandler< dim, spacedim > &dof_handler, const std::unique_ptr< internal::DoFHandlerImplementation::DoFLevel< dim > > &, const std::unique_ptr< internal::DoFHandlerImplementation::DoFFaces< dim > > &mg_faces, const unsigned int obj_index, const types::fe_index fe_index, const unsigned int local_index, const std::integral_constant< int, 1 >)
static types::global_dof_index & get_mg_dof_index(const DoFHandler< dim, spacedim > &dof_handler, const std::unique_ptr< internal::DoFHandlerImplementation::DoFLevel< dim > > &mg_level, const std::unique_ptr< internal::DoFHandlerImplementation::DoFFaces< dim > > &, const unsigned int obj_index, const types::fe_index fe_index, const unsigned int local_index, const std::integral_constant< int, dim >)
static std::pair< unsigned int, unsigned int > process_object_range(const ::DoFAccessor< structdim, dim, spacedim, level_dof_access > accessor, const types::fe_index fe_index)
static types::global_dof_index & mg_vertex_dof_index(DoFHandler< dim, spacedim > &dof_handler, const int level, const unsigned int vertex_index, const unsigned int i)
static void set_mg_dof_indices(const ::DoFAccessor< structdim, dim, spacedim, level_dof_access > &accessor, const int level, const std::vector< types::global_dof_index > &dof_indices, const types::fe_index fe_index)
boost::container::small_vector<::types::global_dof_index, 27 > dof_index_vector_type
static void set_dof_indices(const ::DoFAccessor< structdim, dim, spacedim, level_dof_access > &accessor, const std::vector< types::global_dof_index > &dof_indices, const types::fe_index fe_index)
static void process_dof_index(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int obj_level, const unsigned int obj_index, const types::fe_index fe_index, const unsigned int local_index, const std::integral_constant< int, structdim > &, GlobalIndexType &global_index, const DoFPProcessor &process)
static types::global_dof_index & get_mg_dof_index(const DoFHandler< 3, spacedim > &dof_handler, const std::unique_ptr< internal::DoFHandlerImplementation::DoFLevel< 3 > > &, const std::unique_ptr< internal::DoFHandlerImplementation::DoFFaces< 3 > > &mg_faces, const unsigned int obj_index, const types::fe_index fe_index, const unsigned int local_index, const std::integral_constant< int, 2 >)
static std::pair< unsigned int, unsigned int > process_object_range(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int obj_level, const unsigned int obj_index, const types::fe_index fe_index, const std::integral_constant< int, structdim > &)
static void extract_subvector_to(const InputVector &values, const types::global_dof_index *cache, const types::global_dof_index *cache_end, ForwardIterator local_values_begin)
static std::vector< unsigned int > sort_indices(const types::global_dof_index *v_begin, const types::global_dof_index *v_end)
static types::fe_index nth_active_fe_index(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int obj_level, const unsigned int obj_index, const unsigned int local_index, const std::integral_constant< int, structdim > &)
static void process_object(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int obj_level, const unsigned int obj_index, const types::fe_index fe_index, const DoFMapping &mapping, const std::integral_constant< int, structdim > &dd, types::global_dof_index *&dof_indices_ptr, const DoFProcessor &process)
static types::global_dof_index get_dof_index(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int obj_level, const unsigned int obj_index, const types::fe_index fe_index, const unsigned int local_index, const std::integral_constant< int, structdim > &dd)
static unsigned int n_dof_indices(const ::DoFAccessor< structdim, dim, spacedim, level_dof_access > &accessor, const types::fe_index fe_index_, const bool count_level_dofs)
static std::pair< unsigned int, unsigned int > process_object_range(::DoFInvalidAccessor< structdim, dim, spacedim >, const unsigned int)
static unsigned int get_array_length(const ArrayType &array)
static types::fe_index future_fe_index(const DoFCellAccessor< dim, spacedim, level_dof_access > &accessor)
static void clear_future_fe_index(const DoFCellAccessor< dim, spacedim, level_dof_access > &accessor)
static bool future_fe_index_set(const DoFCellAccessor< dim, spacedim, level_dof_access > &accessor)
static void set_active_fe_index(const DoFCellAccessor< dim, spacedim, level_dof_access > &accessor, const types::fe_index i)
static void set_future_fe_index(const DoFCellAccessor< dim, spacedim, level_dof_access > &accessor, const types::fe_index i)
static types::fe_index active_fe_index(const DoFCellAccessor< dim, spacedim, level_dof_access > &accessor)