Reference documentation for deal.II version 9.3.3
\(\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 : 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, class 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);
200
201
202
203 template <int dim, int q_dim, class 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, class 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
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:3030
UpdateFlags update_flags
Definition: fe_values.h:3934
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
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: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:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
UpdateFlags
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
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:196