Reference documentation for deal.II version 9.3.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 
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  : 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],
165  update_flags);
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);
199  }
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 
261 namespace 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,
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>
282  const QCollection<dim> & q_collection,
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
337  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
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>
373  const hp::QCollection<dim - 1> & q_collection,
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>
384  const hp::MappingCollection<dim, spacedim> & mapping,
386  const std::vector<hp::QCollection<dim - 1>> &q_collection,
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>
399  const hp::QCollection<dim - 1> & q_collection,
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>
411  const std::vector<hp::QCollection<dim - 1>> &q_collection,
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
471  const typename Triangulation<dim, spacedim>::face_iterator &face,
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
485  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
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
520  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
521  const typename Triangulation<dim, spacedim>::face_iterator &face,
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>
538  const hp::QCollection<dim - 1> & q_collection,
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
609  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
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 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
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
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
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:3030
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
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:423
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:3888
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:370
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
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:293
::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:175
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
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:4342
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
UpdateFlags update_flags
Definition: fe_values.h:3934