Reference documentation for deal.II version Git 953c9590e9 2020-10-28 17:25:14 -0400
\(\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  // -------------------------- FEValuesBase -------------------------
29 
30  template <int dim, int q_dim, class FEValuesType>
33  & mapping_collection,
35  const QCollection<q_dim> & q_collection,
36  const UpdateFlags update_flags)
37  : fe_collection(&fe_collection)
38  , mapping_collection(&mapping_collection)
39  , q_collection(q_collection)
40  , fe_values_table(fe_collection.size(),
41  mapping_collection.size(),
42  q_collection.size())
43  , present_fe_values_index(numbers::invalid_unsigned_int,
46  , update_flags(update_flags)
47  {}
48 
49 
50  template <int dim, int q_dim, class FEValuesType>
53  const QCollection<q_dim> & q_collection,
55  : fe_collection(&fe_collection)
56  , mapping_collection(
57  &::hp::StaticMappingQ1<dim, FEValuesType::space_dimension>::
58  mapping_collection)
59  , q_collection(q_collection)
60  , fe_values_table(fe_collection.size(), 1, q_collection.size())
61  , present_fe_values_index(numbers::invalid_unsigned_int,
64  , update_flags(update_flags)
65  {}
66 
67 
68 
69  template <int dim, int q_dim, class FEValuesType>
72  : fe_collection(other.fe_collection)
73  , mapping_collection(other.mapping_collection)
74  , q_collection(other.q_collection)
75  , fe_values_table(fe_collection->size(),
76  mapping_collection->size(),
77  q_collection.size())
78  , present_fe_values_index(other.present_fe_values_index)
79  , update_flags(other.update_flags)
80  {
81  // We've already resized the `fe_values_table` correctly above, but right
82  // now it just contains nullptrs. Create copies of the objects that
83  // `other.fe_values_table` stores
84  Threads::TaskGroup<> task_group;
85  for (unsigned int fe_index = 0; fe_index < fe_collection->size();
86  ++fe_index)
87  for (unsigned int m_index = 0; m_index < mapping_collection->size();
88  ++m_index)
89  for (unsigned int q_index = 0; q_index < q_collection.size(); ++q_index)
90  if (other.fe_values_table[fe_index][m_index][q_index].get() !=
91  nullptr)
92  task_group += Threads::new_task([&, fe_index, m_index, q_index]() {
93  fe_values_table[fe_index][m_index][q_index] =
94  std::make_unique<FEValuesType>((*mapping_collection)[m_index],
95  (*fe_collection)[fe_index],
96  q_collection[q_index],
97  update_flags);
98  });
99 
100  task_group.join_all();
101  }
102 
103 
104 
105  template <int dim, int q_dim, class FEValuesType>
106  FEValuesType &
108  const unsigned int fe_index,
109  const unsigned int mapping_index,
110  const unsigned int q_index)
111  {
112  AssertIndexRange(fe_index, fe_collection->size());
113  AssertIndexRange(mapping_index, mapping_collection->size());
114  AssertIndexRange(q_index, q_collection.size());
115 
116 
117  // set the triple of indices
118  // that we want to work with
119  present_fe_values_index = TableIndices<3>(fe_index, mapping_index, q_index);
120 
121  // first check whether we
122  // already have an object for
123  // this particular combination
124  // of indices
125  if (fe_values_table(present_fe_values_index).get() == nullptr)
126  fe_values_table(present_fe_values_index) =
127  std::make_unique<FEValuesType>((*mapping_collection)[mapping_index],
128  (*fe_collection)[fe_index],
129  q_collection[q_index],
130  update_flags);
131 
132  // now there definitely is one!
133  return *fe_values_table(present_fe_values_index);
134  }
135 
136 
137 
138  template <int dim, int q_dim, class FEValuesType>
139  void
141  const std::vector<unsigned int> &fe_indices,
142  const std::vector<unsigned int> &mapping_indices,
143  const std::vector<unsigned int> &q_indices)
144  {
145  AssertDimension(fe_indices.size(), mapping_indices.size());
146  AssertDimension(fe_indices.size(), q_indices.size());
147 
148  Threads::TaskGroup<> task_group;
149  for (unsigned int i = 0; i < fe_indices.size(); ++i)
150  {
151  const unsigned int fe_index = fe_indices[i],
152  mapping_index = mapping_indices[i],
153  q_index = q_indices[i];
154 
155  AssertIndexRange(fe_index, fe_collection->size());
156  AssertIndexRange(mapping_index, mapping_collection->size());
157  AssertIndexRange(q_index, q_collection.size());
158 
159  task_group +=
160  Threads::new_task([&, fe_index, mapping_index, q_index]() {
161  fe_values_table(TableIndices<3>(fe_index, mapping_index, q_index)) =
162  std::make_unique<FEValuesType>(
163  (*mapping_collection)[mapping_index],
164  (*fe_collection)[fe_index],
165  q_collection[q_index],
166  update_flags);
167  });
168  }
169 
170  task_group.join_all();
171  }
172 
173 
174 
175  template <int dim, int q_dim, class FEValuesType>
176  void
178  {
179  const unsigned int size = fe_collection->size();
180  std::vector<unsigned int> indices(size);
181  std::iota(indices.begin(), indices.end(), 0);
182 
183  precalculate_fe_values(/*fe_indices=*/indices,
184  /*mapping_indices=*/
185  (mapping_collection->size() > 1) ?
186  indices :
187  std::vector<unsigned int>(size, 0),
188  /*q_indices=*/
189  (q_collection.size() > 1) ?
190  indices :
191  std::vector<unsigned int>(size, 0));
192  }
193 } // namespace hp
194 
195 
196 namespace hp
197 {
198  // -------------------------- FEValues -------------------------
199 
200 
201  template <int dim, int spacedim>
204  const FECollection<dim, spacedim> & fe_collection,
205  const QCollection<dim> & q_collection,
207  : hp::FEValuesBase<dim, dim, ::FEValues<dim, spacedim>>(mapping,
208  fe_collection,
209  q_collection,
210  update_flags)
211  {}
212 
213 
214  template <int dim, int spacedim>
219  : hp::FEValuesBase<dim, dim, ::FEValues<dim, spacedim>>(fe_collection,
220  q_collection,
221  update_flags)
222  {}
223 
224 
225  template <int dim, int spacedim>
226  template <bool lda>
227  void
230  const unsigned int q_index,
231  const unsigned int mapping_index,
232  const unsigned int fe_index)
233  {
234  // determine which indices we
235  // should actually use
236  unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
237  real_fe_index = fe_index;
238 
239  if (real_q_index == numbers::invalid_unsigned_int)
240  {
241  if (this->q_collection.size() > 1)
242  real_q_index = cell->active_fe_index();
243  else
244  real_q_index = 0;
245  }
246 
247  if (real_mapping_index == numbers::invalid_unsigned_int)
248  {
249  if (this->mapping_collection->size() > 1)
250  real_mapping_index = cell->active_fe_index();
251  else
252  real_mapping_index = 0;
253  }
254 
255  if (real_fe_index == numbers::invalid_unsigned_int)
256  real_fe_index = cell->active_fe_index();
257 
258  // some checks
259  AssertIndexRange(real_q_index, this->q_collection.size());
260  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
261  AssertIndexRange(real_fe_index, this->fe_collection->size());
262 
263  // now finally actually get the
264  // corresponding object and
265  // initialize it
266  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
267  .reinit(cell);
268  }
269 
270 
271 
272  template <int dim, int spacedim>
273  void
275  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
276  const unsigned int q_index,
277  const unsigned int mapping_index,
278  const unsigned int fe_index)
279  {
280  // determine which indices we
281  // should actually use
282  unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
283  real_fe_index = fe_index;
284 
285  if (real_q_index == numbers::invalid_unsigned_int)
286  real_q_index = 0;
287 
288  if (real_mapping_index == numbers::invalid_unsigned_int)
289  real_mapping_index = 0;
290 
291  if (real_fe_index == numbers::invalid_unsigned_int)
292  real_fe_index = 0;
293 
294  // some checks
295  AssertIndexRange(real_q_index, this->q_collection.size());
296  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
297  AssertIndexRange(real_fe_index, this->fe_collection->size());
298 
299  // now finally actually get the
300  // corresponding object and
301  // initialize it
302  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
303  .reinit(cell);
304  }
305 
306 
307  // -------------------------- FEFaceValues -------------------------
308 
309 
310  template <int dim, int spacedim>
316  : hp::FEValuesBase<dim, dim - 1, ::FEFaceValues<dim, spacedim>>(
317  mapping,
318  fe_collection,
319  q_collection,
320  update_flags)
321  {}
322 
323 
324  template <int dim, int spacedim>
329  : hp::FEValuesBase<dim, dim - 1, ::FEFaceValues<dim, spacedim>>(
330  fe_collection,
331  q_collection,
332  update_flags)
333  {}
334 
335 
336  template <int dim, int spacedim>
337  template <bool lda>
338  void
341  const unsigned int face_no,
342  const unsigned int q_index,
343  const unsigned int mapping_index,
344  const unsigned int fe_index)
345  {
346  // determine which indices we
347  // should actually use
348  unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
349  real_fe_index = fe_index;
350 
351  if (real_q_index == numbers::invalid_unsigned_int)
352  {
353  if (this->q_collection.size() > 1)
354  real_q_index = cell->active_fe_index();
355  else
356  real_q_index = 0;
357  }
358 
359  if (real_mapping_index == numbers::invalid_unsigned_int)
360  {
361  if (this->mapping_collection->size() > 1)
362  real_mapping_index = cell->active_fe_index();
363  else
364  real_mapping_index = 0;
365  }
366 
367  if (real_fe_index == numbers::invalid_unsigned_int)
368  real_fe_index = cell->active_fe_index();
369 
370  // some checks
371  AssertIndexRange(real_q_index, this->q_collection.size());
372  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
373  AssertIndexRange(real_fe_index, this->fe_collection->size());
374 
375  // now finally actually get the
376  // corresponding object and
377  // initialize it
378  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
379  .reinit(cell, face_no);
380  }
381 
382 
383 
384  template <int dim, int spacedim>
385  template <bool lda>
386  void
389  const typename Triangulation<dim, spacedim>::face_iterator &face,
390  const unsigned int q_index,
391  const unsigned int mapping_index,
392  const unsigned int fe_index)
393  {
394  const auto face_n = cell->face_iterator_to_index(face);
395  reinit(cell, face_n, q_index, mapping_index, fe_index);
396  }
397 
398 
399 
400  template <int dim, int spacedim>
401  void
403  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
404  const unsigned int face_no,
405  const unsigned int q_index,
406  const unsigned int mapping_index,
407  const unsigned int fe_index)
408  {
409  // determine which indices we
410  // should actually use
411  unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
412  real_fe_index = fe_index;
413 
414  if (real_q_index == numbers::invalid_unsigned_int)
415  real_q_index = 0;
416 
417  if (real_mapping_index == numbers::invalid_unsigned_int)
418  real_mapping_index = 0;
419 
420  if (real_fe_index == numbers::invalid_unsigned_int)
421  real_fe_index = 0;
422 
423  // some checks
424  AssertIndexRange(real_q_index, this->q_collection.size());
425  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
426  AssertIndexRange(real_fe_index, this->fe_collection->size());
427 
428  // now finally actually get the
429  // corresponding object and
430  // initialize it
431  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
432  .reinit(cell, face_no);
433  }
434 
435 
436 
437  template <int dim, int spacedim>
438  void
440  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
441  const typename Triangulation<dim, spacedim>::face_iterator &face,
442  const unsigned int q_index,
443  const unsigned int mapping_index,
444  const unsigned int fe_index)
445  {
446  const auto face_n = cell->face_iterator_to_index(face);
447  reinit(cell, face_n, q_index, mapping_index, fe_index);
448  }
449 
450 
451  // -------------------------- FESubfaceValues -------------------------
452 
453 
454  template <int dim, int spacedim>
460  : hp::FEValuesBase<dim, dim - 1, ::FESubfaceValues<dim, spacedim>>(
461  mapping,
462  fe_collection,
463  q_collection,
464  update_flags)
465  {}
466 
467 
468  template <int dim, int spacedim>
470  const hp::FECollection<dim, spacedim> &fe_collection,
471  const hp::QCollection<dim - 1> & q_collection,
472  const UpdateFlags update_flags)
473  : hp::FEValuesBase<dim, dim - 1, ::FESubfaceValues<dim, spacedim>>(
474  fe_collection,
475  q_collection,
476  update_flags)
477  {}
478 
479 
480  template <int dim, int spacedim>
481  template <bool lda>
482  void
485  const unsigned int face_no,
486  const unsigned int subface_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  {
498  if (this->q_collection.size() > 1)
499  real_q_index = cell->active_fe_index();
500  else
501  real_q_index = 0;
502  }
503 
504  if (real_mapping_index == numbers::invalid_unsigned_int)
505  {
506  if (this->mapping_collection->size() > 1)
507  real_mapping_index = cell->active_fe_index();
508  else
509  real_mapping_index = 0;
510  }
511 
512  if (real_fe_index == numbers::invalid_unsigned_int)
513  real_fe_index = cell->active_fe_index();
514 
515  // some checks
516  AssertIndexRange(real_q_index, this->q_collection.size());
517  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
518  AssertIndexRange(real_fe_index, this->fe_collection->size());
519 
520  // now finally actually get the
521  // corresponding object and
522  // initialize it
523  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
524  .reinit(cell, face_no, subface_no);
525  }
526 
527 
528 
529  template <int dim, int spacedim>
530  void
532  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
533  const unsigned int face_no,
534  const unsigned int subface_no,
535  const unsigned int q_index,
536  const unsigned int mapping_index,
537  const unsigned int fe_index)
538  {
539  // determine which indices we
540  // should actually use
541  unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
542  real_fe_index = fe_index;
543 
544  if (real_q_index == numbers::invalid_unsigned_int)
545  real_q_index = 0;
546 
547  if (real_mapping_index == numbers::invalid_unsigned_int)
548  real_mapping_index = 0;
549 
550  if (real_fe_index == numbers::invalid_unsigned_int)
551  real_fe_index = 0;
552 
553  // some checks
554  AssertIndexRange(real_q_index, this->q_collection.size());
555  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
556  AssertIndexRange(real_fe_index, this->fe_collection->size());
557 
558  // now finally actually get the
559  // corresponding object and
560  // initialize it
561  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
562  .reinit(cell, face_no, subface_no);
563  }
564 } // namespace hp
565 
566 
567 // explicit instantiations
568 #include "fe_values.inst"
569 
570 
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:483
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1580
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:202
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:3072
#define AssertIndexRange(index, range)
Definition: exceptions.h:1648
const SmartPointer< const MappingCollection< dim, ::FEValues< dim, spacedim > ::space_dimension >, FEValuesBase< dim, q_dim, ::FEValues< dim, spacedim > > > mapping_collection
Definition: fe_values.h:195
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:339
static const unsigned int space_dimension
Definition: fe_values.h:2085
const SmartPointer< const Mapping< dim, spacedim >, FEValuesBase< dim, spacedim > > mapping
Definition: fe_values.h:3497
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:311
unsigned int size() const
Definition: q_collection.h:174
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:369
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:187
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:228
::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:107
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:368
Table< 3, std::unique_ptr< FEValuesType > > fe_values_table
Definition: fe_values.h:214
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:31
FEValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim > &quadrature, const UpdateFlags update_flags)
Definition: fe_values.cc:4386
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:455
UpdateFlags update_flags
Definition: fe_values.h:3543