Reference documentation for deal.II version Git 73c87d96ef 2021-11-30 22:54:44 +0100
\(\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 
18 #include <deal.II/fe/mapping_q1.h>
19 
20 #include <deal.II/hp/fe_values.h>
21 
22 #include <memory>
23 
25 
26 namespace 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,
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,
111  : FEValuesBase(
112  ::hp::StaticMappingQ1<dim, FEValuesType::space_dimension>::
113  mapping_collection,
114  fe_collection,
115  q_collection,
116  update_flags)
117  {}
118 
119 
120  template <int dim, int q_dim, class FEValuesType>
123  const std::vector<QCollection<q_dim>> & q_collection,
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],
166  update_flags);
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 
262 namespace 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,
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>
283  const QCollection<dim> & q_collection,
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
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>
374  const hp::QCollection<dim - 1> & q_collection,
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>
385  const hp::MappingCollection<dim, spacedim> & mapping,
387  const std::vector<hp::QCollection<dim - 1>> &q_collection,
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>
400  const hp::QCollection<dim - 1> & q_collection,
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>
412  const std::vector<hp::QCollection<dim - 1>> &q_collection,
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
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
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
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>
539  const hp::QCollection<dim - 1> & q_collection,
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
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 "fe_values.inst"
644 
645 
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
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1655
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
Task< RT > new_task(const std::function< RT()> &function)
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
#define AssertIndexRange(index, range)
Definition: exceptions.h:1720
const SmartPointer< const MappingCollection< dim, ::FEValues< dim, spacedim > ::space_dimension >, FEValuesBase< dim, q_dim, ::FEValues< dim, spacedim > > > mapping_collection
Definition: fe_values.h:217
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
static const unsigned int space_dimension
Definition: fe_values.h:2423
const SmartPointer< const Mapping< dim, spacedim >, FEValuesBase< dim, spacedim > > mapping
Definition: fe_values.h:3939
UpdateFlags
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
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
Definition: hp.h:117
const SmartPointer< const FECollection< dim, ::FEValues< dim, spacedim > ::space_dimension >, FEValuesBase< dim, q_dim, ::FEValues< dim, spacedim > > > fe_collection
Definition: fe_values.h:209
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< dim, spacedim > & select_fe_values(const unsigned int fe_index, const unsigned int mapping_index, const unsigned int q_index)
Definition: fe_values.cc:176
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
Table< 3, std::unique_ptr< FEValuesType > > fe_values_table
Definition: fe_values.h:246
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
FEValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim > &quadrature, const UpdateFlags update_flags)
Definition: fe_values.cc:4165
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
UpdateFlags update_flags
Definition: fe_values.h:3985