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