deal.II version 9.8.0
\(\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\}}\)
Loading...
Searching...
No Matches
fe_values.cc
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2005 - 2025 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13#include <deal.II/base/thread_management.h>
14
15#include <deal.II/fe/mapping_q1.h>
16
17#include <deal.II/hp/fe_values.h>
18
19#include <memory>
20#include <numeric>
21
22
23DEAL_II_NAMESPACE_OPEN
24
25namespace hp
26{
27 namespace
28 {
33 template <int q_dim>
34 std::vector<QCollection<q_dim>>
35 translate(const QCollection<q_dim> &q_collection)
36 {
37 std::vector<QCollection<q_dim>> q_collections;
38
39 q_collections.reserve(q_collection.size());
40 for (unsigned int q = 0; q < q_collection.size(); ++q)
41 q_collections.emplace_back(q_collection[q]);
42
43 return q_collections;
44 }
45
50 template <int q_dim>
51 QCollection<q_dim>
52 translate(const std::vector<QCollection<q_dim>> &q_collections)
53 {
54 QCollection<q_dim> result;
55
56 for (unsigned int q = 0; q < q_collections.size(); ++q)
57 result.push_back(q_collections[q][0]);
58
59 return result;
60 }
61 } // namespace
62
63 // -------------------------- FEValuesBase -------------------------
64
65 template <int dim, int q_dim, typename FEValuesType>
84
85 template <int dim, int q_dim, typename FEValuesType>
104
105
106 template <int dim, int q_dim, typename FEValuesType>
118
119
120 template <int dim, int q_dim, typename FEValuesType>
132
133
134
135 template <int dim, int q_dim, typename FEValuesType>
136 FEValuesBase<dim, q_dim, FEValuesType>::FEValuesBase(
137 const FEValuesBase<dim, q_dim, FEValuesType> &other)
138 : EnableObserverPointer(other)
139 , fe_collection(other.fe_collection)
140 , mapping_collection(other.mapping_collection)
141 , q_collection(other.q_collection)
142 , q_collections(other.q_collections)
143 , fe_values_table(other.fe_values_table.size(0),
144 other.fe_values_table.size(1),
145 other.fe_values_table.size(2))
146 , present_fe_values_index(other.present_fe_values_index)
147 , update_flags(other.update_flags)
148 {
149 // We've already resized the `fe_values_table` correctly above, but right
150 // now it just contains nullptrs. Create copies of the objects that
151 // `other.fe_values_table` stores
152 Threads::TaskGroup<> task_group;
153 for (unsigned int fe_index = 0; fe_index < other.fe_values_table.size(0);
154 ++fe_index)
155 for (unsigned int m_index = 0; m_index < other.fe_values_table.size(1);
156 ++m_index)
157 for (unsigned int q_index = 0; q_index < other.fe_values_table.size(2);
158 ++q_index)
159 if (other.fe_values_table[fe_index][m_index][q_index].get() !=
160 nullptr)
161 task_group += Threads::new_task([&, fe_index, m_index, q_index]() {
162 fe_values_table[fe_index][m_index][q_index] =
163 std::make_unique<FEValuesType>((*mapping_collection)[m_index],
164 (*fe_collection)[fe_index],
165 q_collections[q_index],
166 update_flags);
167 });
168
169 task_group.join_all();
170 }
171
172
173
174 template <int dim, int q_dim, typename FEValuesType>
175 FEValuesType &
176 FEValuesBase<dim, q_dim, FEValuesType>::select_fe_values(
177 const unsigned int fe_index,
178 const unsigned int mapping_index,
179 const unsigned int q_index)
180 {
181 AssertIndexRange(fe_index, fe_collection->size());
182 AssertIndexRange(mapping_index, mapping_collection->size());
183 AssertIndexRange(q_index, q_collections.size());
184
185
186 // set the triple of indices that we want to work with
187 present_fe_values_index = {fe_index, mapping_index, q_index};
188
189 // first check whether we already have an object for this particular
190 // combination of indices
191 if (fe_values_table(present_fe_values_index).get() == nullptr)
192 fe_values_table(present_fe_values_index) =
193 std::make_unique<FEValuesType>((*mapping_collection)[mapping_index],
194 (*fe_collection)[fe_index],
195 q_collections[q_index],
196 update_flags);
197
198 // now there definitely is one!
199 return *fe_values_table(present_fe_values_index);
200 }
201
202
203
204 template <int dim, int q_dim, typename FEValuesType>
205 void
206 FEValuesBase<dim, q_dim, FEValuesType>::precalculate_fe_values(
207 const std::vector<unsigned int> &fe_indices,
208 const std::vector<unsigned int> &mapping_indices,
209 const std::vector<unsigned int> &q_indices)
210 {
211 AssertDimension(fe_indices.size(), mapping_indices.size());
212 AssertDimension(fe_indices.size(), q_indices.size());
213
214 Threads::TaskGroup<> task_group;
215 for (unsigned int i = 0; i < fe_indices.size(); ++i)
216 {
217 const unsigned int fe_index = fe_indices[i],
218 mapping_index = mapping_indices[i],
219 q_index = q_indices[i];
220
221 AssertIndexRange(fe_index, fe_collection->size());
222 AssertIndexRange(mapping_index, mapping_collection->size());
223 AssertIndexRange(q_index, q_collections.size());
224
225 task_group +=
226 Threads::new_task([&, fe_index, mapping_index, q_index]() {
227 fe_values_table[fe_index][mapping_index][q_index] =
228 std::make_unique<FEValuesType>(
229 (*mapping_collection)[mapping_index],
230 (*fe_collection)[fe_index],
231 q_collections[q_index],
232 update_flags);
233 });
234 }
235
236 task_group.join_all();
237 }
238
239
240
241 template <int dim, int q_dim, typename FEValuesType>
242 void
243 FEValuesBase<dim, q_dim, FEValuesType>::precalculate_fe_values()
244 {
245 const unsigned int size = fe_collection->size();
246 std::vector<unsigned int> indices(size);
247 std::iota(indices.begin(), indices.end(), 0);
248
249 precalculate_fe_values(/*fe_indices=*/indices,
250 /*mapping_indices=*/
251 (mapping_collection->size() > 1) ?
252 indices :
253 std::vector<unsigned int>(size, 0),
254 /*q_indices=*/
255 (q_collections.size() > 1) ?
256 indices :
257 std::vector<unsigned int>(size, 0));
258 }
259} // namespace hp
260
261
262namespace hp
263{
264 // -------------------------- FEValues -------------------------
265
266
267 template <int dim, int spacedim>
278
279
280 template <int dim, int spacedim>
289
290
291 template <int dim, int spacedim>
292 template <bool lda>
293 void
294 FEValues<dim, spacedim>::reinit(
295 const TriaIterator<DoFCellAccessor<dim, spacedim, lda>> &cell,
296 const unsigned int q_index,
297 const unsigned int mapping_index,
298 const unsigned int fe_index)
299 {
300 // determine which indices we should actually use
301 unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
302 real_fe_index = fe_index;
303
304 if (real_q_index == numbers::invalid_unsigned_int)
305 {
306 if (this->q_collections.size() > 1)
307 real_q_index = cell->active_fe_index();
308 else
309 real_q_index = 0;
310 }
311
312 if (real_mapping_index == numbers::invalid_unsigned_int)
313 {
314 if (this->mapping_collection->size() > 1)
315 real_mapping_index = cell->active_fe_index();
316 else
317 real_mapping_index = 0;
318 }
319
320 if (real_fe_index == numbers::invalid_unsigned_int)
321 real_fe_index = cell->active_fe_index();
322
323 // some checks
324 AssertIndexRange(real_q_index, this->q_collections.size());
325 AssertIndexRange(real_mapping_index, this->mapping_collection->size());
326 AssertIndexRange(real_fe_index, this->fe_collection->size());
327
328 // now finally actually get the corresponding object and initialize it
329 this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
330 .reinit(cell);
331 }
332
333
334
335 template <int dim, int spacedim>
336 void
337 FEValues<dim, spacedim>::reinit(
338 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
339 const unsigned int q_index,
340 const unsigned int mapping_index,
341 const unsigned int fe_index)
342 {
343 // determine which indices we should actually use
344 unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
345 real_fe_index = fe_index;
346
347 if (real_q_index == numbers::invalid_unsigned_int)
348 real_q_index = 0;
349
350 if (real_mapping_index == numbers::invalid_unsigned_int)
351 real_mapping_index = 0;
352
353 if (real_fe_index == numbers::invalid_unsigned_int)
354 real_fe_index = 0;
355
356 // some checks
357 AssertIndexRange(real_q_index, this->q_collections.size());
358 AssertIndexRange(real_mapping_index, this->mapping_collection->size());
359 AssertIndexRange(real_fe_index, this->fe_collection->size());
360
361 // now finally actually get the corresponding object and initialize it
362 this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
363 .reinit(cell);
364 }
365
366
367 // -------------------------- FEFaceValues -------------------------
368
369
370 template <int dim, int spacedim>
382
383 template <int dim, int spacedim>
384 FEFaceValues<dim, spacedim>::FEFaceValues(
385 const hp::MappingCollection<dim, spacedim> &mapping,
386 const hp::FECollection<dim, spacedim> &fe_collection,
387 const std::vector<hp::QCollection<dim - 1>> &q_collection,
388 const UpdateFlags update_flags)
389 : hp::FEValuesBase<dim, dim - 1, ::FEFaceValues<dim, spacedim>>(
390 mapping,
391 fe_collection,
392 q_collection,
393 update_flags)
394 {}
395
396
397 template <int dim, int spacedim>
407
408
409 template <int dim, int spacedim>
410 FEFaceValues<dim, spacedim>::FEFaceValues(
411 const hp::FECollection<dim, spacedim> &fe_collection,
412 const std::vector<hp::QCollection<dim - 1>> &q_collection,
413 const UpdateFlags update_flags)
414 : hp::FEValuesBase<dim, dim - 1, ::FEFaceValues<dim, spacedim>>(
415 fe_collection,
416 q_collection,
417 update_flags)
418 {}
419
420
421 template <int dim, int spacedim>
422 template <bool lda>
423 void
424 FEFaceValues<dim, spacedim>::reinit(
425 const TriaIterator<DoFCellAccessor<dim, spacedim, lda>> &cell,
426 const unsigned int face_no,
427 const unsigned int q_index,
428 const unsigned int mapping_index,
429 const unsigned int fe_index)
430 {
431 // determine which indices we
432 // should actually use
433 unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
434 real_fe_index = fe_index;
435
436 if (real_q_index == numbers::invalid_unsigned_int)
437 {
438 if (this->q_collections.size() > 1)
439 real_q_index = cell->active_fe_index();
440 else
441 real_q_index = 0;
442 }
443
444 if (real_mapping_index == numbers::invalid_unsigned_int)
445 {
446 if (this->mapping_collection->size() > 1)
447 real_mapping_index = cell->active_fe_index();
448 else
449 real_mapping_index = 0;
450 }
451
452 if (real_fe_index == numbers::invalid_unsigned_int)
453 real_fe_index = cell->active_fe_index();
454
455 // some checks
456 AssertIndexRange(real_q_index, this->q_collections.size());
457 AssertIndexRange(real_mapping_index, this->mapping_collection->size());
458 AssertIndexRange(real_fe_index, this->fe_collection->size());
459
460 // now finally actually get the corresponding object and initialize it
461 this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
462 .reinit(cell, face_no);
463 }
464
465
466
467 template <int dim, int spacedim>
468 template <bool lda>
469 void
470 FEFaceValues<dim, spacedim>::reinit(
471 const TriaIterator<DoFCellAccessor<dim, spacedim, lda>> &cell,
472 const typename Triangulation<dim, spacedim>::face_iterator &face,
473 const unsigned int q_index,
474 const unsigned int mapping_index,
475 const unsigned int fe_index)
476 {
477 const auto face_n = cell->face_iterator_to_index(face);
478 reinit(cell, face_n, q_index, mapping_index, fe_index);
479 }
480
481
482
483 template <int dim, int spacedim>
484 void
485 FEFaceValues<dim, spacedim>::reinit(
486 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
487 const unsigned int face_no,
488 const unsigned int q_index,
489 const unsigned int mapping_index,
490 const unsigned int fe_index)
491 {
492 // determine which indices we
493 // should actually use
494 unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
495 real_fe_index = fe_index;
496
497 if (real_q_index == numbers::invalid_unsigned_int)
498 real_q_index = 0;
499
500 if (real_mapping_index == numbers::invalid_unsigned_int)
501 real_mapping_index = 0;
502
503 if (real_fe_index == numbers::invalid_unsigned_int)
504 real_fe_index = 0;
505
506 // some checks
507 AssertIndexRange(real_q_index, this->q_collections.size());
508 AssertIndexRange(real_mapping_index, this->mapping_collection->size());
509 AssertIndexRange(real_fe_index, this->fe_collection->size());
510
511 // now finally actually get the corresponding object and initialize it
512 this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
513 .reinit(cell, face_no);
514 }
515
516
517
518 template <int dim, int spacedim>
519 void
520 FEFaceValues<dim, spacedim>::reinit(
521 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
522 const typename Triangulation<dim, spacedim>::face_iterator &face,
523 const unsigned int q_index,
524 const unsigned int mapping_index,
525 const unsigned int fe_index)
526 {
527 const auto face_n = cell->face_iterator_to_index(face);
528 reinit(cell, face_n, q_index, mapping_index, fe_index);
529 }
530
531
532 // -------------------------- FESubfaceValues -------------------------
533
534
535 template <int dim, int spacedim>
536 FESubfaceValues<dim, spacedim>::FESubfaceValues(
537 const hp::MappingCollection<dim, spacedim> &mapping,
538 const hp::FECollection<dim, spacedim> &fe_collection,
539 const hp::QCollection<dim - 1> &q_collection,
540 const UpdateFlags update_flags)
541 : hp::FEValuesBase<dim, dim - 1, ::FESubfaceValues<dim, spacedim>>(
542 mapping,
543 fe_collection,
544 q_collection,
545 update_flags)
546 {}
547
548
549 template <int dim, int spacedim>
550 FESubfaceValues<dim, spacedim>::FESubfaceValues(
551 const hp::FECollection<dim, spacedim> &fe_collection,
552 const hp::QCollection<dim - 1> &q_collection,
553 const UpdateFlags update_flags)
554 : hp::FEValuesBase<dim, dim - 1, ::FESubfaceValues<dim, spacedim>>(
555 fe_collection,
556 q_collection,
557 update_flags)
558 {}
559
560
561 template <int dim, int spacedim>
562 template <bool lda>
563 void
564 FESubfaceValues<dim, spacedim>::reinit(
565 const TriaIterator<DoFCellAccessor<dim, spacedim, lda>> &cell,
566 const unsigned int face_no,
567 const unsigned int subface_no,
568 const unsigned int q_index,
569 const unsigned int mapping_index,
570 const unsigned int fe_index)
571 {
572 // determine which indices we should actually use
573 unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
574 real_fe_index = fe_index;
575
576 if (real_q_index == numbers::invalid_unsigned_int)
577 {
578 if (this->q_collections.size() > 1)
579 real_q_index = cell->active_fe_index();
580 else
581 real_q_index = 0;
582 }
583
584 if (real_mapping_index == numbers::invalid_unsigned_int)
585 {
586 if (this->mapping_collection->size() > 1)
587 real_mapping_index = cell->active_fe_index();
588 else
589 real_mapping_index = 0;
590 }
591
592 if (real_fe_index == numbers::invalid_unsigned_int)
593 real_fe_index = cell->active_fe_index();
594
595 // some checks
596 AssertIndexRange(real_q_index, this->q_collections.size());
597 AssertIndexRange(real_mapping_index, this->mapping_collection->size());
598 AssertIndexRange(real_fe_index, this->fe_collection->size());
599
600 // now finally actually get the corresponding object and initialize it
601 this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
602 .reinit(cell, face_no, subface_no);
603 }
604
605
606
607 template <int dim, int spacedim>
608 void
609 FESubfaceValues<dim, spacedim>::reinit(
610 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
611 const unsigned int face_no,
612 const unsigned int subface_no,
613 const unsigned int q_index,
614 const unsigned int mapping_index,
615 const unsigned int fe_index)
616 {
617 // determine which indices we should actually use
618 unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
619 real_fe_index = fe_index;
620
621 if (real_q_index == numbers::invalid_unsigned_int)
622 real_q_index = 0;
623
624 if (real_mapping_index == numbers::invalid_unsigned_int)
625 real_mapping_index = 0;
626
627 if (real_fe_index == numbers::invalid_unsigned_int)
628 real_fe_index = 0;
629
630 // some checks
631 AssertIndexRange(real_q_index, this->q_collections.size());
632 AssertIndexRange(real_mapping_index, this->mapping_collection->size());
633 AssertIndexRange(real_fe_index, this->fe_collection->size());
634
635 // now finally actually get the corresponding object and initialize it
636 this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
637 .reinit(cell, face_no, subface_no);
638 }
639} // namespace hp
640
641
642// explicit instantiations
643#include "hp/fe_values.inst"
644
645
646DEAL_II_NAMESPACE_CLOSE
CellAccessor::face_iterator_to_index
unsigned int face_iterator_to_index(const TriaIterator< TriaAccessor< dim - 1, dim, spacedim > > &face) const
hp::FEFaceValues::reinit
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:424
hp::FEFaceValues::FEFaceValues
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:371
hp::FESubfaceValues::reinit
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:564
hp::FESubfaceValues::FESubfaceValues
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:536
hp::FEValuesBase::fe_collection
const ObserverPointer< const FECollection< dim, FEValuesType::space_dimension >, FEValuesBase< dim, q_dim, FEValuesType > > fe_collection
Definition fe_values.h:206
hp::FEValuesBase::FEValuesBase
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:66
hp::FEValuesBase< dim, dim - 1, ::FEFaceValues< dim, spacedim > >::q_collections
const std::vector< QCollection< q_dim > > q_collections
Definition fe_values.h:229
hp::FEValuesBase::precalculate_fe_values
void precalculate_fe_values(const std::vector< unsigned int > &fe_indices, const std::vector< unsigned int > &mapping_indices, const std::vector< unsigned int > &q_indices)
Definition fe_values.cc:206
hp::FEValuesBase::select_fe_values
FEValuesType & select_fe_values(const unsigned int fe_index, const unsigned int mapping_index, const unsigned int q_index)
Definition fe_values.cc:176
hp::FEValuesBase::q_collection
const QCollection< q_dim > q_collection
Definition fe_values.h:219
hp::FEValuesBase::mapping_collection
const ObserverPointer< const MappingCollection< dim, FEValuesType::space_dimension >, FEValuesBase< dim, q_dim, FEValuesType > > mapping_collection
Definition fe_values.h:214
hp::FEValuesBase< dim, dim - 1, ::FEFaceValues< dim, spacedim > >::fe_values_table
Table< 3, std::unique_ptr< ::FEFaceValues< dim, spacedim > > > fe_values_table
Definition fe_values.h:243
hp::FEValuesBase::precalculate_fe_values
void precalculate_fe_values()
Definition fe_values.cc:243
hp::FEValuesBase::update_flags
const UpdateFlags update_flags
Definition fe_values.h:254
hp::FEValues::reinit
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:294
hp::FEValues::FEValues
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:268
hp::QCollection::push_back
void push_back(const Quadrature< dim_in > &new_quadrature)
Definition q_collection.h:213
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
space_dimension
const unsigned int DoFAccessor< structdim, dim, spacedim, level_dof_access >::space_dimension
Definition dof_accessor.cc:35
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exception_macros.h:829
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exception_macros.h:901
Triangulation::face_iterator
TriaIterator< TriaAccessor< dim - 1, dim, spacedim > > face_iterator
Definition tria.h:1663
Triangulation::cell_iterator
TriaIterator< CellAccessor< dim, spacedim > > cell_iterator
Definition tria.h:1624
Threads::new_task
Task< RT > new_task(const std::function< RT()> &function)
Definition thread_management.h:1173
fe_values.h
mapping_q1.h
mapping
MappingQ< dim, spacedim > StaticMappingQ1< dim, spacedim >::mapping
Definition mapping_q1.h:101
hp
Definition hp.h:115
hp::mapping_collection
MappingCollection< dim, spacedim > StaticMappingQ1< dim, spacedim >::mapping_collection
Definition mapping_collection.h:178
numbers::invalid_unsigned_int
constexpr unsigned int invalid_unsigned_int
Definition types.h:236
thread_management.h
Generated by doxygen 1.16.1