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