Reference documentation for deal.II version Git a6dcdbf211 2021-01-17 11:23:55 +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 - 2020 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
186  // that we want to work with
187  present_fe_values_index = TableIndices<3>(fe_index, mapping_index, q_index);
188 
189  // first check whether we
190  // already have an object for
191  // this particular combination
192  // of indices
193  if (fe_values_table(present_fe_values_index).get() == nullptr)
194  fe_values_table(present_fe_values_index) =
195  std::make_unique<FEValuesType>((*mapping_collection)[mapping_index],
196  (*fe_collection)[fe_index],
197  q_collections[q_index],
198  update_flags);
199 
200  // now there definitely is one!
201  return *fe_values_table(present_fe_values_index);
202  }
203 
204 
205 
206  template <int dim, int q_dim, class FEValuesType>
207  void
209  const std::vector<unsigned int> &fe_indices,
210  const std::vector<unsigned int> &mapping_indices,
211  const std::vector<unsigned int> &q_indices)
212  {
213  AssertDimension(fe_indices.size(), mapping_indices.size());
214  AssertDimension(fe_indices.size(), q_indices.size());
215 
216  Threads::TaskGroup<> task_group;
217  for (unsigned int i = 0; i < fe_indices.size(); ++i)
218  {
219  const unsigned int fe_index = fe_indices[i],
220  mapping_index = mapping_indices[i],
221  q_index = q_indices[i];
222 
223  AssertIndexRange(fe_index, fe_collection->size());
224  AssertIndexRange(mapping_index, mapping_collection->size());
225  AssertIndexRange(q_index, q_collections.size());
226 
227  task_group +=
228  Threads::new_task([&, fe_index, mapping_index, q_index]() {
229  fe_values_table(TableIndices<3>(fe_index, mapping_index, q_index)) =
230  std::make_unique<FEValuesType>(
231  (*mapping_collection)[mapping_index],
232  (*fe_collection)[fe_index],
233  q_collections[q_index],
234  update_flags);
235  });
236  }
237 
238  task_group.join_all();
239  }
240 
241 
242 
243  template <int dim, int q_dim, class FEValuesType>
244  void
246  {
247  const unsigned int size = fe_collection->size();
248  std::vector<unsigned int> indices(size);
249  std::iota(indices.begin(), indices.end(), 0);
250 
251  precalculate_fe_values(/*fe_indices=*/indices,
252  /*mapping_indices=*/
253  (mapping_collection->size() > 1) ?
254  indices :
255  std::vector<unsigned int>(size, 0),
256  /*q_indices=*/
257  (q_collections.size() > 1) ?
258  indices :
259  std::vector<unsigned int>(size, 0));
260  }
261 } // namespace hp
262 
263 
264 namespace hp
265 {
266  // -------------------------- FEValues -------------------------
267 
268 
269  template <int dim, int spacedim>
272  const FECollection<dim, spacedim> & fe_collection,
273  const QCollection<dim> & q_collection,
275  : hp::FEValuesBase<dim, dim, ::FEValues<dim, spacedim>>(mapping,
276  fe_collection,
277  q_collection,
278  update_flags)
279  {}
280 
281 
282  template <int dim, int spacedim>
285  const QCollection<dim> & q_collection,
287  : hp::FEValuesBase<dim, dim, ::FEValues<dim, spacedim>>(fe_collection,
288  q_collection,
289  update_flags)
290  {}
291 
292 
293  template <int dim, int spacedim>
294  template <bool lda>
295  void
298  const unsigned int q_index,
299  const unsigned int mapping_index,
300  const unsigned int fe_index)
301  {
302  // determine which indices we
303  // should actually use
304  unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
305  real_fe_index = fe_index;
306 
307  if (real_q_index == numbers::invalid_unsigned_int)
308  {
309  if (this->q_collections.size() > 1)
310  real_q_index = cell->active_fe_index();
311  else
312  real_q_index = 0;
313  }
314 
315  if (real_mapping_index == numbers::invalid_unsigned_int)
316  {
317  if (this->mapping_collection->size() > 1)
318  real_mapping_index = cell->active_fe_index();
319  else
320  real_mapping_index = 0;
321  }
322 
323  if (real_fe_index == numbers::invalid_unsigned_int)
324  real_fe_index = cell->active_fe_index();
325 
326  // some checks
327  AssertIndexRange(real_q_index, this->q_collections.size());
328  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
329  AssertIndexRange(real_fe_index, this->fe_collection->size());
330 
331  // now finally actually get the
332  // corresponding object and
333  // initialize it
334  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
335  .reinit(cell);
336  }
337 
338 
339 
340  template <int dim, int spacedim>
341  void
343  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
344  const unsigned int q_index,
345  const unsigned int mapping_index,
346  const unsigned int fe_index)
347  {
348  // determine which indices we
349  // should actually use
350  unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
351  real_fe_index = fe_index;
352 
353  if (real_q_index == numbers::invalid_unsigned_int)
354  real_q_index = 0;
355 
356  if (real_mapping_index == numbers::invalid_unsigned_int)
357  real_mapping_index = 0;
358 
359  if (real_fe_index == numbers::invalid_unsigned_int)
360  real_fe_index = 0;
361 
362  // some checks
363  AssertIndexRange(real_q_index, this->q_collections.size());
364  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
365  AssertIndexRange(real_fe_index, this->fe_collection->size());
366 
367  // now finally actually get the
368  // corresponding object and
369  // initialize it
370  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
371  .reinit(cell);
372  }
373 
374 
375  // -------------------------- FEFaceValues -------------------------
376 
377 
378  template <int dim, int spacedim>
382  const hp::QCollection<dim - 1> & q_collection,
384  : hp::FEValuesBase<dim, dim - 1, ::FEFaceValues<dim, spacedim>>(
385  mapping,
386  fe_collection,
387  q_collection,
388  update_flags)
389  {}
390 
391  template <int dim, int spacedim>
393  const hp::MappingCollection<dim, spacedim> & mapping,
395  const std::vector<hp::QCollection<dim - 1>> &q_collection,
397  : hp::FEValuesBase<dim, dim - 1, ::FEFaceValues<dim, spacedim>>(
398  mapping,
399  fe_collection,
400  q_collection,
401  update_flags)
402  {}
403 
404 
405  template <int dim, int spacedim>
408  const hp::QCollection<dim - 1> & q_collection,
410  : hp::FEValuesBase<dim, dim - 1, ::FEFaceValues<dim, spacedim>>(
411  fe_collection,
412  q_collection,
413  update_flags)
414  {}
415 
416 
417  template <int dim, int spacedim>
420  const std::vector<hp::QCollection<dim - 1>> &q_collection,
422  : hp::FEValuesBase<dim, dim - 1, ::FEFaceValues<dim, spacedim>>(
423  fe_collection,
424  q_collection,
425  update_flags)
426  {}
427 
428 
429  template <int dim, int spacedim>
430  template <bool lda>
431  void
434  const unsigned int face_no,
435  const unsigned int q_index,
436  const unsigned int mapping_index,
437  const unsigned int fe_index)
438  {
439  // determine which indices we
440  // should actually use
441  unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
442  real_fe_index = fe_index;
443 
444  if (real_q_index == numbers::invalid_unsigned_int)
445  {
446  if (this->q_collections.size() > 1)
447  real_q_index = cell->active_fe_index();
448  else
449  real_q_index = 0;
450  }
451 
452  if (real_mapping_index == numbers::invalid_unsigned_int)
453  {
454  if (this->mapping_collection->size() > 1)
455  real_mapping_index = cell->active_fe_index();
456  else
457  real_mapping_index = 0;
458  }
459 
460  if (real_fe_index == numbers::invalid_unsigned_int)
461  real_fe_index = cell->active_fe_index();
462 
463  // some checks
464  AssertIndexRange(real_q_index, this->q_collections.size());
465  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
466  AssertIndexRange(real_fe_index, this->fe_collection->size());
467 
468  // now finally actually get the
469  // corresponding object and
470  // initialize it
471  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
472  .reinit(cell, face_no);
473  }
474 
475 
476 
477  template <int dim, int spacedim>
478  template <bool lda>
479  void
482  const typename Triangulation<dim, spacedim>::face_iterator &face,
483  const unsigned int q_index,
484  const unsigned int mapping_index,
485  const unsigned int fe_index)
486  {
487  const auto face_n = cell->face_iterator_to_index(face);
488  reinit(cell, face_n, q_index, mapping_index, fe_index);
489  }
490 
491 
492 
493  template <int dim, int spacedim>
494  void
496  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
497  const unsigned int face_no,
498  const unsigned int q_index,
499  const unsigned int mapping_index,
500  const unsigned int fe_index)
501  {
502  // determine which indices we
503  // should actually use
504  unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
505  real_fe_index = fe_index;
506 
507  if (real_q_index == numbers::invalid_unsigned_int)
508  real_q_index = 0;
509 
510  if (real_mapping_index == numbers::invalid_unsigned_int)
511  real_mapping_index = 0;
512 
513  if (real_fe_index == numbers::invalid_unsigned_int)
514  real_fe_index = 0;
515 
516  // some checks
517  AssertIndexRange(real_q_index, this->q_collections.size());
518  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
519  AssertIndexRange(real_fe_index, this->fe_collection->size());
520 
521  // now finally actually get the
522  // corresponding object and
523  // initialize it
524  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
525  .reinit(cell, face_no);
526  }
527 
528 
529 
530  template <int dim, int spacedim>
531  void
533  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
534  const typename Triangulation<dim, spacedim>::face_iterator &face,
535  const unsigned int q_index,
536  const unsigned int mapping_index,
537  const unsigned int fe_index)
538  {
539  const auto face_n = cell->face_iterator_to_index(face);
540  reinit(cell, face_n, q_index, mapping_index, fe_index);
541  }
542 
543 
544  // -------------------------- FESubfaceValues -------------------------
545 
546 
547  template <int dim, int spacedim>
551  const hp::QCollection<dim - 1> & q_collection,
553  : hp::FEValuesBase<dim, dim - 1, ::FESubfaceValues<dim, spacedim>>(
554  mapping,
555  fe_collection,
556  q_collection,
557  update_flags)
558  {}
559 
560 
561  template <int dim, int spacedim>
563  const hp::FECollection<dim, spacedim> &fe_collection,
564  const hp::QCollection<dim - 1> & q_collection,
565  const UpdateFlags update_flags)
566  : hp::FEValuesBase<dim, dim - 1, ::FESubfaceValues<dim, spacedim>>(
567  fe_collection,
568  q_collection,
569  update_flags)
570  {}
571 
572 
573  template <int dim, int spacedim>
574  template <bool lda>
575  void
578  const unsigned int face_no,
579  const unsigned int subface_no,
580  const unsigned int q_index,
581  const unsigned int mapping_index,
582  const unsigned int fe_index)
583  {
584  // determine which indices we
585  // should actually use
586  unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
587  real_fe_index = fe_index;
588 
589  if (real_q_index == numbers::invalid_unsigned_int)
590  {
591  if (this->q_collections.size() > 1)
592  real_q_index = cell->active_fe_index();
593  else
594  real_q_index = 0;
595  }
596 
597  if (real_mapping_index == numbers::invalid_unsigned_int)
598  {
599  if (this->mapping_collection->size() > 1)
600  real_mapping_index = cell->active_fe_index();
601  else
602  real_mapping_index = 0;
603  }
604 
605  if (real_fe_index == numbers::invalid_unsigned_int)
606  real_fe_index = cell->active_fe_index();
607 
608  // some checks
609  AssertIndexRange(real_q_index, this->q_collections.size());
610  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
611  AssertIndexRange(real_fe_index, this->fe_collection->size());
612 
613  // now finally actually get the
614  // corresponding object and
615  // initialize it
616  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
617  .reinit(cell, face_no, subface_no);
618  }
619 
620 
621 
622  template <int dim, int spacedim>
623  void
625  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
626  const unsigned int face_no,
627  const unsigned int subface_no,
628  const unsigned int q_index,
629  const unsigned int mapping_index,
630  const unsigned int fe_index)
631  {
632  // determine which indices we
633  // should actually use
634  unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
635  real_fe_index = fe_index;
636 
637  if (real_q_index == numbers::invalid_unsigned_int)
638  real_q_index = 0;
639 
640  if (real_mapping_index == numbers::invalid_unsigned_int)
641  real_mapping_index = 0;
642 
643  if (real_fe_index == numbers::invalid_unsigned_int)
644  real_fe_index = 0;
645 
646  // some checks
647  AssertIndexRange(real_q_index, this->q_collections.size());
648  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
649  AssertIndexRange(real_fe_index, this->fe_collection->size());
650 
651  // now finally actually get the
652  // corresponding object and
653  // initialize it
654  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
655  .reinit(cell, face_no, subface_no);
656  }
657 } // namespace hp
658 
659 
660 // explicit instantiations
661 #include "fe_values.inst"
662 
663 
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:576
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1623
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:270
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:3043
#define AssertIndexRange(index, range)
Definition: exceptions.h:1691
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:432
static const unsigned int space_dimension
Definition: fe_values.h:2177
const SmartPointer< const Mapping< dim, spacedim >, FEValuesBase< dim, spacedim > > mapping
Definition: fe_values.h:3642
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:379
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:380
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:296
::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
size_type size(const unsigned int i) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:379
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:4355
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:548
UpdateFlags update_flags
Definition: fe_values.h:3688