Reference documentation for deal.II version 9.4.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\}}\)
fe_values.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2003 - 2021 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
17
19
21
22#include <memory>
23
25
26namespace hp
27{
28 namespace
29 {
34 template <int q_dim>
35 std::vector<QCollection<q_dim>>
36 translate(const QCollection<q_dim> &q_collection)
37 {
38 std::vector<QCollection<q_dim>> q_collections;
39
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, class FEValuesType>
68 & mapping_collection,
70 const QCollection<q_dim> & q_collection,
71 const UpdateFlags update_flags)
72 : fe_collection(&fe_collection)
73 , mapping_collection(&mapping_collection)
74 , q_collection(q_collection)
75 , q_collections(translate(q_collection))
76 , fe_values_table(fe_collection.size(),
77 mapping_collection.size(),
78 q_collection.size())
79 , present_fe_values_index(numbers::invalid_unsigned_int,
82 , update_flags(update_flags)
83 {}
84
85 template <int dim, int q_dim, class FEValuesType>
88 & mapping_collection,
90 const std::vector<QCollection<q_dim>> & q_collections,
91 const UpdateFlags update_flags)
92 : fe_collection(&fe_collection)
93 , mapping_collection(&mapping_collection)
94 , q_collection(translate(q_collections))
95 , q_collections(q_collections)
96 , fe_values_table(fe_collection.size(),
97 mapping_collection.size(),
98 q_collection.size())
99 , present_fe_values_index(numbers::invalid_unsigned_int,
102 , update_flags(update_flags)
103 {}
104
105
106 template <int dim, int q_dim, class FEValuesType>
109 const QCollection<q_dim> & q_collection,
110 const UpdateFlags update_flags)
111 : FEValuesBase(
112 ::hp::StaticMappingQ1<dim, FEValuesType::space_dimension>::
113 mapping_collection,
114 fe_collection,
115 q_collection,
116 update_flags)
117 {}
119
120 template <int dim, int q_dim, class FEValuesType>
123 const std::vector<QCollection<q_dim>> & q_collection,
124 const UpdateFlags update_flags)
125 : FEValuesBase(
126 ::hp::StaticMappingQ1<dim, FEValuesType::space_dimension>::
127 mapping_collection,
128 fe_collection,
129 q_collection,
130 update_flags)
131 {}
132
133
134
135 template <int dim, int q_dim, class FEValuesType>
138 : Subscriptor(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],
167 });
168
169 task_group.join_all();
170 }
171
172
173
174 template <int dim, int q_dim, class FEValuesType>
175 FEValuesType &
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, class FEValuesType>
205 void
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, class FEValuesType>
242 void
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>
270 const FECollection<dim, spacedim> & fe_collection,
271 const QCollection<dim> & q_collection,
272 const UpdateFlags update_flags)
273 : hp::FEValuesBase<dim, dim, ::FEValues<dim, spacedim>>(mapping,
274 fe_collection,
275 q_collection,
276 update_flags)
277 {}
278
279
280 template <int dim, int spacedim>
282 const FECollection<dim, spacedim> &fe_collection,
283 const QCollection<dim> & q_collection,
284 const UpdateFlags update_flags)
285 : hp::FEValuesBase<dim, dim, ::FEValues<dim, spacedim>>(fe_collection,
286 q_collection,
287 update_flags)
288 {}
289
290
291 template <int dim, int spacedim>
292 template <bool lda>
293 void
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
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>
373 const hp::FECollection<dim, spacedim> & fe_collection,
374 const hp::QCollection<dim - 1> & q_collection,
375 const UpdateFlags update_flags)
376 : hp::FEValuesBase<dim, dim - 1, ::FEFaceValues<dim, spacedim>>(
377 mapping,
378 fe_collection,
379 q_collection,
380 update_flags)
381 {}
382
383 template <int dim, int spacedim>
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>
399 const hp::FECollection<dim, spacedim> &fe_collection,
400 const hp::QCollection<dim - 1> & q_collection,
401 const UpdateFlags update_flags)
402 : hp::FEValuesBase<dim, dim - 1, ::FEFaceValues<dim, spacedim>>(
403 fe_collection,
404 q_collection,
405 update_flags)
406 {}
407
408
409 template <int dim, int spacedim>
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
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
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
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
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>
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>
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
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
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 "fe_values.inst"
644
645
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)
Definition: fe_values.cc:2863
UpdateFlags update_flags
Definition: fe_values.h:3985
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
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
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
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
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
Table< 3, std::unique_ptr< FEValuesType > > fe_values_table
Definition: fe_values.h:246
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
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
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
UpdateFlags
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
Task< RT > new_task(const std::function< RT()> &function)
Definition: hp.h:118
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
static const unsigned int invalid_unsigned_int
Definition: types.h:201