Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-2776-gf808a1a12e 2025-03-07 23:20: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
fe_values.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2005 - 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
16
18
20
21#include <memory>
22#include <numeric>
23
24
26
27namespace hp
28{
29 namespace
30 {
35 template <int q_dim>
36 std::vector<QCollection<q_dim>>
37 translate(const QCollection<q_dim> &q_collection)
38 {
39 std::vector<QCollection<q_dim>> q_collections;
40
41 for (unsigned int q = 0; q < q_collection.size(); ++q)
42 q_collections.emplace_back(q_collection[q]);
43
44 return q_collections;
45 }
46
51 template <int q_dim>
52 QCollection<q_dim>
53 translate(const std::vector<QCollection<q_dim>> &q_collections)
54 {
55 QCollection<q_dim> result;
56
57 for (unsigned int q = 0; q < q_collections.size(); ++q)
58 result.push_back(q_collections[q][0]);
59
60 return result;
61 }
62 } // namespace
63
64 // -------------------------- FEValuesBase -------------------------
65
66 template <int dim, int q_dim, typename FEValuesType>
69 &mapping_collection,
71 const QCollection<q_dim> &q_collection,
72 const UpdateFlags update_flags)
73 : fe_collection(&fe_collection)
74 , mapping_collection(&mapping_collection)
75 , q_collection(q_collection)
76 , q_collections(translate(q_collection))
77 , fe_values_table(fe_collection.size(),
78 mapping_collection.size(),
79 q_collection.size())
80 , present_fe_values_index(numbers::invalid_unsigned_int,
81 numbers::invalid_unsigned_int,
82 numbers::invalid_unsigned_int)
83 , update_flags(update_flags)
84 {}
85
86 template <int dim, int q_dim, typename FEValuesType>
89 &mapping_collection,
91 const std::vector<QCollection<q_dim>> &q_collections,
92 const UpdateFlags update_flags)
93 : fe_collection(&fe_collection)
94 , mapping_collection(&mapping_collection)
95 , q_collection(translate(q_collections))
96 , q_collections(q_collections)
97 , fe_values_table(fe_collection.size(),
98 mapping_collection.size(),
99 q_collection.size())
100 , present_fe_values_index(numbers::invalid_unsigned_int,
101 numbers::invalid_unsigned_int,
102 numbers::invalid_unsigned_int)
103 , update_flags(update_flags)
104 {}
105
106
107 template <int dim, int q_dim, typename FEValuesType>
110 const QCollection<q_dim> &q_collection,
111 const UpdateFlags update_flags)
112 : FEValuesBase(
113 ::hp::StaticMappingQ1<dim, FEValuesType::space_dimension>::
114 mapping_collection,
115 fe_collection,
116 q_collection,
117 update_flags)
118 {}
119
120
121 template <int dim, int q_dim, typename FEValuesType>
124 const std::vector<QCollection<q_dim>> &q_collection,
125 const UpdateFlags update_flags)
126 : FEValuesBase(
127 ::hp::StaticMappingQ1<dim, FEValuesType::space_dimension>::
128 mapping_collection,
129 fe_collection,
130 q_collection,
131 update_flags)
132 {}
133
134
136 template <int dim, int q_dim, typename FEValuesType>
139 : EnableObserverPointer(other)
140 , fe_collection(other.fe_collection)
141 , mapping_collection(other.mapping_collection)
142 , q_collection(other.q_collection)
143 , q_collections(other.q_collections)
144 , fe_values_table(other.fe_values_table.size(0),
145 other.fe_values_table.size(1),
146 other.fe_values_table.size(2))
147 , present_fe_values_index(other.present_fe_values_index)
148 , update_flags(other.update_flags)
149 {
150 // We've already resized the `fe_values_table` correctly above, but right
151 // now it just contains nullptrs. Create copies of the objects that
152 // `other.fe_values_table` stores
153 Threads::TaskGroup<> task_group;
154 for (unsigned int fe_index = 0; fe_index < other.fe_values_table.size(0);
155 ++fe_index)
156 for (unsigned int m_index = 0; m_index < other.fe_values_table.size(1);
157 ++m_index)
158 for (unsigned int q_index = 0; q_index < other.fe_values_table.size(2);
159 ++q_index)
160 if (other.fe_values_table[fe_index][m_index][q_index].get() !=
161 nullptr)
162 task_group += Threads::new_task([&, fe_index, m_index, q_index]() {
163 fe_values_table[fe_index][m_index][q_index] =
164 std::make_unique<FEValuesType>((*mapping_collection)[m_index],
165 (*fe_collection)[fe_index],
166 q_collections[q_index],
168 });
169
170 task_group.join_all();
171 }
172
173
174
175 template <int dim, int q_dim, typename FEValuesType>
176 FEValuesType &
178 const unsigned int fe_index,
179 const unsigned int mapping_index,
180 const unsigned int q_index)
181 {
182 AssertIndexRange(fe_index, fe_collection->size());
183 AssertIndexRange(mapping_index, mapping_collection->size());
184 AssertIndexRange(q_index, q_collections.size());
185
186
187 // set the triple of indices that we want to work with
188 present_fe_values_index = {fe_index, mapping_index, q_index};
189
190 // first check whether we already have an object for this particular
191 // combination of indices
192 if (fe_values_table(present_fe_values_index).get() == nullptr)
193 fe_values_table(present_fe_values_index) =
194 std::make_unique<FEValuesType>((*mapping_collection)[mapping_index],
195 (*fe_collection)[fe_index],
196 q_collections[q_index],
197 update_flags);
198
199 // now there definitely is one!
200 return *fe_values_table(present_fe_values_index);
201 }
202
203
204
205 template <int dim, int q_dim, typename FEValuesType>
206 void
208 const std::vector<unsigned int> &fe_indices,
209 const std::vector<unsigned int> &mapping_indices,
210 const std::vector<unsigned int> &q_indices)
211 {
212 AssertDimension(fe_indices.size(), mapping_indices.size());
213 AssertDimension(fe_indices.size(), q_indices.size());
214
215 Threads::TaskGroup<> task_group;
216 for (unsigned int i = 0; i < fe_indices.size(); ++i)
217 {
218 const unsigned int fe_index = fe_indices[i],
219 mapping_index = mapping_indices[i],
220 q_index = q_indices[i];
221
222 AssertIndexRange(fe_index, fe_collection->size());
223 AssertIndexRange(mapping_index, mapping_collection->size());
224 AssertIndexRange(q_index, q_collections.size());
225
226 task_group +=
227 Threads::new_task([&, fe_index, mapping_index, q_index]() {
228 fe_values_table[fe_index][mapping_index][q_index] =
229 std::make_unique<FEValuesType>(
230 (*mapping_collection)[mapping_index],
231 (*fe_collection)[fe_index],
232 q_collections[q_index],
233 update_flags);
234 });
235 }
236
237 task_group.join_all();
238 }
239
240
241
242 template <int dim, int q_dim, typename FEValuesType>
243 void
245 {
246 const unsigned int size = fe_collection->size();
247 std::vector<unsigned int> indices(size);
248 std::iota(indices.begin(), indices.end(), 0);
249
250 precalculate_fe_values(/*fe_indices=*/indices,
251 /*mapping_indices=*/
252 (mapping_collection->size() > 1) ?
253 indices :
254 std::vector<unsigned int>(size, 0),
255 /*q_indices=*/
256 (q_collections.size() > 1) ?
257 indices :
258 std::vector<unsigned int>(size, 0));
259 }
260} // namespace hp
261
262
263namespace hp
264{
265 // -------------------------- FEValues -------------------------
266
267
268 template <int dim, int spacedim>
271 const FECollection<dim, spacedim> &fe_collection,
272 const QCollection<dim> &q_collection,
273 const UpdateFlags update_flags)
274 : hp::FEValuesBase<dim, dim, ::FEValues<dim, spacedim>>(mapping,
275 fe_collection,
276 q_collection,
277 update_flags)
278 {}
279
280
281 template <int dim, int spacedim>
283 const FECollection<dim, spacedim> &fe_collection,
284 const QCollection<dim> &q_collection,
285 const UpdateFlags update_flags)
286 : hp::FEValuesBase<dim, dim, ::FEValues<dim, spacedim>>(fe_collection,
287 q_collection,
288 update_flags)
289 {}
290
291
292 template <int dim, int spacedim>
293 template <bool lda>
294 void
297 const unsigned int q_index,
298 const unsigned int mapping_index,
299 const unsigned int fe_index)
300 {
301 // determine which indices we should actually use
302 unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
303 real_fe_index = fe_index;
304
305 if (real_q_index == numbers::invalid_unsigned_int)
306 {
307 if (this->q_collections.size() > 1)
308 real_q_index = cell->active_fe_index();
309 else
310 real_q_index = 0;
311 }
312
313 if (real_mapping_index == numbers::invalid_unsigned_int)
314 {
315 if (this->mapping_collection->size() > 1)
316 real_mapping_index = cell->active_fe_index();
317 else
318 real_mapping_index = 0;
319 }
320
321 if (real_fe_index == numbers::invalid_unsigned_int)
322 real_fe_index = cell->active_fe_index();
323
324 // some checks
325 AssertIndexRange(real_q_index, this->q_collections.size());
326 AssertIndexRange(real_mapping_index, this->mapping_collection->size());
327 AssertIndexRange(real_fe_index, this->fe_collection->size());
328
329 // now finally actually get the corresponding object and initialize it
330 this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
331 .reinit(cell);
332 }
333
334
335
336 template <int dim, int spacedim>
337 void
340 const unsigned int q_index,
341 const unsigned int mapping_index,
342 const unsigned int fe_index)
343 {
344 // determine which indices we should actually use
345 unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
346 real_fe_index = fe_index;
347
348 if (real_q_index == numbers::invalid_unsigned_int)
349 real_q_index = 0;
350
351 if (real_mapping_index == numbers::invalid_unsigned_int)
352 real_mapping_index = 0;
353
354 if (real_fe_index == numbers::invalid_unsigned_int)
355 real_fe_index = 0;
356
357 // some checks
358 AssertIndexRange(real_q_index, this->q_collections.size());
359 AssertIndexRange(real_mapping_index, this->mapping_collection->size());
360 AssertIndexRange(real_fe_index, this->fe_collection->size());
361
362 // now finally actually get the corresponding object and initialize it
363 this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
364 .reinit(cell);
365 }
366
367
368 // -------------------------- FEFaceValues -------------------------
369
370
371 template <int dim, int spacedim>
374 const hp::FECollection<dim, spacedim> &fe_collection,
375 const hp::QCollection<dim - 1> &q_collection,
376 const UpdateFlags update_flags)
377 : hp::FEValuesBase<dim, dim - 1, ::FEFaceValues<dim, spacedim>>(
378 mapping,
379 fe_collection,
380 q_collection,
381 update_flags)
382 {}
383
384 template <int dim, int spacedim>
387 const hp::FECollection<dim, spacedim> &fe_collection,
388 const std::vector<hp::QCollection<dim - 1>> &q_collection,
389 const UpdateFlags update_flags)
390 : hp::FEValuesBase<dim, dim - 1, ::FEFaceValues<dim, spacedim>>(
391 mapping,
392 fe_collection,
393 q_collection,
394 update_flags)
395 {}
396
397
398 template <int dim, int spacedim>
400 const hp::FECollection<dim, spacedim> &fe_collection,
401 const hp::QCollection<dim - 1> &q_collection,
402 const UpdateFlags update_flags)
403 : hp::FEValuesBase<dim, dim - 1, ::FEFaceValues<dim, spacedim>>(
404 fe_collection,
405 q_collection,
406 update_flags)
407 {}
408
409
410 template <int dim, int spacedim>
412 const hp::FECollection<dim, spacedim> &fe_collection,
413 const std::vector<hp::QCollection<dim - 1>> &q_collection,
414 const UpdateFlags update_flags)
415 : hp::FEValuesBase<dim, dim - 1, ::FEFaceValues<dim, spacedim>>(
416 fe_collection,
417 q_collection,
418 update_flags)
419 {}
420
421
422 template <int dim, int spacedim>
423 template <bool lda>
424 void
427 const unsigned int face_no,
428 const unsigned int q_index,
429 const unsigned int mapping_index,
430 const unsigned int fe_index)
431 {
432 // determine which indices we
433 // should actually use
434 unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
435 real_fe_index = fe_index;
436
437 if (real_q_index == numbers::invalid_unsigned_int)
438 {
439 if (this->q_collections.size() > 1)
440 real_q_index = cell->active_fe_index();
441 else
442 real_q_index = 0;
443 }
444
445 if (real_mapping_index == numbers::invalid_unsigned_int)
446 {
447 if (this->mapping_collection->size() > 1)
448 real_mapping_index = cell->active_fe_index();
449 else
450 real_mapping_index = 0;
451 }
452
453 if (real_fe_index == numbers::invalid_unsigned_int)
454 real_fe_index = cell->active_fe_index();
455
456 // some checks
457 AssertIndexRange(real_q_index, this->q_collections.size());
458 AssertIndexRange(real_mapping_index, this->mapping_collection->size());
459 AssertIndexRange(real_fe_index, this->fe_collection->size());
460
461 // now finally actually get the corresponding object and initialize it
462 this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
463 .reinit(cell, face_no);
464 }
465
466
467
468 template <int dim, int spacedim>
469 template <bool lda>
470 void
474 const unsigned int q_index,
475 const unsigned int mapping_index,
476 const unsigned int fe_index)
477 {
478 const auto face_n = cell->face_iterator_to_index(face);
479 reinit(cell, face_n, q_index, mapping_index, fe_index);
480 }
481
482
483
484 template <int dim, int spacedim>
485 void
488 const unsigned int face_no,
489 const unsigned int q_index,
490 const unsigned int mapping_index,
491 const unsigned int fe_index)
492 {
493 // determine which indices we
494 // should actually use
495 unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
496 real_fe_index = fe_index;
497
498 if (real_q_index == numbers::invalid_unsigned_int)
499 real_q_index = 0;
500
501 if (real_mapping_index == numbers::invalid_unsigned_int)
502 real_mapping_index = 0;
503
504 if (real_fe_index == numbers::invalid_unsigned_int)
505 real_fe_index = 0;
506
507 // some checks
508 AssertIndexRange(real_q_index, this->q_collections.size());
509 AssertIndexRange(real_mapping_index, this->mapping_collection->size());
510 AssertIndexRange(real_fe_index, this->fe_collection->size());
511
512 // now finally actually get the corresponding object and initialize it
513 this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
514 .reinit(cell, face_no);
515 }
516
517
518
519 template <int dim, int spacedim>
520 void
524 const unsigned int q_index,
525 const unsigned int mapping_index,
526 const unsigned int fe_index)
527 {
528 const auto face_n = cell->face_iterator_to_index(face);
529 reinit(cell, face_n, q_index, mapping_index, fe_index);
530 }
531
532
533 // -------------------------- FESubfaceValues -------------------------
534
535
536 template <int dim, int spacedim>
539 const hp::FECollection<dim, spacedim> &fe_collection,
540 const hp::QCollection<dim - 1> &q_collection,
541 const UpdateFlags update_flags)
542 : hp::FEValuesBase<dim, dim - 1, ::FESubfaceValues<dim, spacedim>>(
543 mapping,
544 fe_collection,
545 q_collection,
546 update_flags)
547 {}
548
549
550 template <int dim, int spacedim>
552 const hp::FECollection<dim, spacedim> &fe_collection,
553 const hp::QCollection<dim - 1> &q_collection,
554 const UpdateFlags update_flags)
555 : hp::FEValuesBase<dim, dim - 1, ::FESubfaceValues<dim, spacedim>>(
556 fe_collection,
557 q_collection,
558 update_flags)
559 {}
560
561
562 template <int dim, int spacedim>
563 template <bool lda>
564 void
567 const unsigned int face_no,
568 const unsigned int subface_no,
569 const unsigned int q_index,
570 const unsigned int mapping_index,
571 const unsigned int fe_index)
572 {
573 // determine which indices we should actually use
574 unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
575 real_fe_index = fe_index;
576
577 if (real_q_index == numbers::invalid_unsigned_int)
578 {
579 if (this->q_collections.size() > 1)
580 real_q_index = cell->active_fe_index();
581 else
582 real_q_index = 0;
583 }
584
585 if (real_mapping_index == numbers::invalid_unsigned_int)
586 {
587 if (this->mapping_collection->size() > 1)
588 real_mapping_index = cell->active_fe_index();
589 else
590 real_mapping_index = 0;
591 }
592
593 if (real_fe_index == numbers::invalid_unsigned_int)
594 real_fe_index = cell->active_fe_index();
595
596 // some checks
597 AssertIndexRange(real_q_index, this->q_collections.size());
598 AssertIndexRange(real_mapping_index, this->mapping_collection->size());
599 AssertIndexRange(real_fe_index, this->fe_collection->size());
600
601 // now finally actually get the corresponding object and initialize it
602 this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
603 .reinit(cell, face_no, subface_no);
604 }
605
606
607
608 template <int dim, int spacedim>
609 void
612 const unsigned int face_no,
613 const unsigned int subface_no,
614 const unsigned int q_index,
615 const unsigned int mapping_index,
616 const unsigned int fe_index)
617 {
618 // determine which indices we should actually use
619 unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
620 real_fe_index = fe_index;
621
622 if (real_q_index == numbers::invalid_unsigned_int)
623 real_q_index = 0;
624
625 if (real_mapping_index == numbers::invalid_unsigned_int)
626 real_mapping_index = 0;
627
628 if (real_fe_index == numbers::invalid_unsigned_int)
629 real_fe_index = 0;
630
631 // some checks
632 AssertIndexRange(real_q_index, this->q_collections.size());
633 AssertIndexRange(real_mapping_index, this->mapping_collection->size());
634 AssertIndexRange(real_fe_index, this->fe_collection->size());
635
636 // now finally actually get the corresponding object and initialize it
637 this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
638 .reinit(cell, face_no, subface_no);
639 }
640} // namespace hp
641
642
643// explicit instantiations
644#include "hp/fe_values.inst"
645
646
FEValuesBase(const unsigned int n_q_points, const unsigned int dofs_per_cell, const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe)
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, lda > > &cell, const unsigned int face_no, const unsigned int q_index=numbers::invalid_unsigned_int, const unsigned int mapping_index=numbers::invalid_unsigned_int, const unsigned int fe_index=numbers::invalid_unsigned_int)
Definition fe_values.cc:425
FEFaceValues(const hp::MappingCollection< dim, spacedim > &mapping_collection, const hp::FECollection< dim, spacedim > &fe_collection, const hp::QCollection< dim - 1 > &q_collection, const UpdateFlags update_flags)
Definition fe_values.cc:372
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, lda > > &cell, const unsigned int face_no, const unsigned int subface_no, const unsigned int q_index=numbers::invalid_unsigned_int, const unsigned int mapping_index=numbers::invalid_unsigned_int, const unsigned int fe_index=numbers::invalid_unsigned_int)
Definition fe_values.cc:565
FESubfaceValues(const hp::MappingCollection< dim, spacedim > &mapping_collection, const hp::FECollection< dim, spacedim > &fe_collection, const hp::QCollection< dim - 1 > &q_collection, const UpdateFlags update_flags)
Definition fe_values.cc:537
const ObserverPointer< const FECollection< dim, FEValuesType::space_dimension >, FEValuesBase< dim, q_dim, FEValuesType > > fe_collection
Definition fe_values.h:208
FEValuesBase(const MappingCollection< dim, FEValuesType::space_dimension > &mapping_collection, const FECollection< dim, FEValuesType::space_dimension > &fe_collection, const QCollection< q_dim > &q_collection, const UpdateFlags update_flags)
Definition fe_values.cc:67
const std::vector< QCollection< q_dim > > q_collections
Definition fe_values.h:231
FEValuesType & select_fe_values(const unsigned int fe_index, const unsigned int mapping_index, const unsigned int q_index)
Definition fe_values.cc:177
const ObserverPointer< const MappingCollection< dim, FEValuesType::space_dimension >, FEValuesBase< dim, q_dim, FEValuesType > > mapping_collection
Definition fe_values.h:216
Table< 3, std::unique_ptr< FEValuesType > > fe_values_table
Definition fe_values.h:245
void precalculate_fe_values()
Definition fe_values.cc:244
const UpdateFlags update_flags
Definition fe_values.h:256
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, lda > > &cell, const unsigned int q_index=numbers::invalid_unsigned_int, const unsigned int mapping_index=numbers::invalid_unsigned_int, const unsigned int fe_index=numbers::invalid_unsigned_int)
Definition fe_values.cc:295
FEValues(const MappingCollection< dim, spacedim > &mapping_collection, const FECollection< dim, spacedim > &fe_collection, const QCollection< dim > &q_collection, const UpdateFlags update_flags)
Definition fe_values.cc:269
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
UpdateFlags
Task< RT > new_task(const std::function< RT()> &function)
std::size_t size
Definition mpi.cc:743
Definition hp.h:117
constexpr unsigned int invalid_unsigned_int
Definition types.h:232