Reference documentation for deal.II version 9.2.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 - 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 
16 #include <deal.II/base/std_cxx14/memory.h>
17 #include <deal.II/base/thread_management.h>
18 
19 #include <deal.II/fe/mapping_q1.h>
20 
21 #include <deal.II/hp/fe_values.h>
22 
23 DEAL_II_NAMESPACE_OPEN
24 
25 namespace hp
26 {
27  // -------------------------- FEValuesBase -------------------------
28 
29  template <int dim, int q_dim, class FEValuesType>
30  FEValuesBase<dim, q_dim, FEValuesType>::FEValuesBase(
31  const MappingCollection<dim, FEValuesType::space_dimension>
32  & mapping_collection,
33  const FECollection<dim, FEValuesType::space_dimension> &fe_collection,
34  const QCollection<q_dim> & q_collection,
35  const UpdateFlags update_flags)
36  : fe_collection(&fe_collection)
37  , mapping_collection(&mapping_collection)
38  , q_collection(q_collection)
39  , fe_values_table(fe_collection.size(),
40  mapping_collection.size(),
41  q_collection.size())
42  , present_fe_values_index(numbers::invalid_unsigned_int,
43  numbers::invalid_unsigned_int,
44  numbers::invalid_unsigned_int)
45  , update_flags(update_flags)
46  {}
47 
48 
49  template <int dim, int q_dim, class FEValuesType>
50  FEValuesBase<dim, q_dim, FEValuesType>::FEValuesBase(
51  const FECollection<dim, FEValuesType::space_dimension> &fe_collection,
52  const QCollection<q_dim> & q_collection,
53  const UpdateFlags update_flags)
54  : fe_collection(&fe_collection)
55  , mapping_collection(
56  &::hp::StaticMappingQ1<dim, FEValuesType::space_dimension>::
57  mapping_collection)
58  , q_collection(q_collection)
59  , fe_values_table(fe_collection.size(), 1, q_collection.size())
60  , present_fe_values_index(numbers::invalid_unsigned_int,
61  numbers::invalid_unsigned_int,
62  numbers::invalid_unsigned_int)
63  , update_flags(update_flags)
64  {}
65 
66 
67 
68  template <int dim, int q_dim, class FEValuesType>
69  FEValuesBase<dim, q_dim, FEValuesType>::FEValuesBase(
70  const FEValuesBase<dim, q_dim, FEValuesType> &other)
71  : fe_collection(other.fe_collection)
72  , mapping_collection(other.mapping_collection)
73  , q_collection(other.q_collection)
74  , fe_values_table(fe_collection->size(),
75  mapping_collection->size(),
76  q_collection.size())
77  , present_fe_values_index(other.present_fe_values_index)
78  , update_flags(other.update_flags)
79  {
80  // We've already resized the `fe_values_table` correctly above, but right
81  // now it just contains nullptrs. Create copies of the objects that
82  // `other.fe_values_table` stores
83  Threads::TaskGroup<> task_group;
84  for (unsigned int fe_index = 0; fe_index < fe_collection->size();
85  ++fe_index)
86  for (unsigned int m_index = 0; m_index < mapping_collection->size();
87  ++m_index)
88  for (unsigned int q_index = 0; q_index < q_collection.size(); ++q_index)
89  if (other.fe_values_table[fe_index][m_index][q_index].get() !=
90  nullptr)
91  task_group += Threads::new_task([&, fe_index, m_index, q_index]() {
92  fe_values_table[fe_index][m_index][q_index] =
93  std_cxx14::make_unique<FEValuesType>(
94  (*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 &
107  FEValuesBase<dim, q_dim, FEValuesType>::select_fe_values(
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_cxx14::make_unique<FEValuesType>(
128  (*mapping_collection)[mapping_index],
129  (*fe_collection)[fe_index],
130  q_collection[q_index],
131  update_flags);
132 
133  // now there definitely is one!
134  return *fe_values_table(present_fe_values_index);
135  }
136 
137 
138 
139  template <int dim, int q_dim, class FEValuesType>
140  void
141  FEValuesBase<dim, q_dim, FEValuesType>::precalculate_fe_values(
142  const std::vector<unsigned int> &fe_indices,
143  const std::vector<unsigned int> &mapping_indices,
144  const std::vector<unsigned int> &q_indices)
145  {
146  AssertDimension(fe_indices.size(), mapping_indices.size());
147  AssertDimension(fe_indices.size(), q_indices.size());
148 
149  Threads::TaskGroup<> task_group;
150  for (unsigned int i = 0; i < fe_indices.size(); ++i)
151  {
152  const unsigned int fe_index = fe_indices[i],
153  mapping_index = mapping_indices[i],
154  q_index = q_indices[i];
155 
156  AssertIndexRange(fe_index, fe_collection->size());
157  AssertIndexRange(mapping_index, mapping_collection->size());
158  AssertIndexRange(q_index, q_collection.size());
159 
160  task_group +=
161  Threads::new_task([&, fe_index, mapping_index, q_index]() {
162  fe_values_table(TableIndices<3>(fe_index, mapping_index, q_index)) =
163  std_cxx14::make_unique<FEValuesType>(
164  (*mapping_collection)[mapping_index],
165  (*fe_collection)[fe_index],
166  q_collection[q_index],
167  update_flags);
168  });
169  }
170 
171  task_group.join_all();
172  }
173 
174 
175 
176  template <int dim, int q_dim, class FEValuesType>
177  void
178  FEValuesBase<dim, q_dim, FEValuesType>::precalculate_fe_values()
179  {
180  const unsigned int size = fe_collection->size();
181  std::vector<unsigned int> indices(size);
182  std::iota(indices.begin(), indices.end(), 0);
183 
184  precalculate_fe_values(/*fe_indices=*/indices,
185  /*mapping_indices=*/
186  (mapping_collection->size() > 1) ?
187  indices :
188  std::vector<unsigned int>(size, 0),
189  /*q_indices=*/
190  (q_collection.size() > 1) ?
191  indices :
192  std::vector<unsigned int>(size, 0));
193  }
194 } // namespace hp
195 
196 
197 namespace hp
198 {
199  // -------------------------- FEValues -------------------------
200 
201 
202  template <int dim, int spacedim>
203  FEValues<dim, spacedim>::FEValues(
204  const MappingCollection<dim, spacedim> &mapping,
205  const FECollection<dim, spacedim> & fe_collection,
206  const QCollection<dim> & q_collection,
207  const UpdateFlags update_flags)
208  : hp::FEValuesBase<dim, dim, ::FEValues<dim, spacedim>>(mapping,
209  fe_collection,
210  q_collection,
211  update_flags)
212  {}
213 
214 
215  template <int dim, int spacedim>
216  FEValues<dim, spacedim>::FEValues(
217  const FECollection<dim, spacedim> &fe_collection,
218  const QCollection<dim> & q_collection,
219  const UpdateFlags update_flags)
220  : hp::FEValuesBase<dim, dim, ::FEValues<dim, spacedim>>(fe_collection,
221  q_collection,
222  update_flags)
223  {}
224 
225 
226  template <int dim, int spacedim>
227  template <typename DoFHandlerType, bool lda>
228  void
229  FEValues<dim, spacedim>::reinit(
230  const TriaIterator<DoFCellAccessor<DoFHandlerType, lda>> cell,
231  const unsigned int q_index,
232  const unsigned int mapping_index,
233  const unsigned int fe_index)
234  {
235  // determine which indices we
236  // should actually use
237  unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
238  real_fe_index = fe_index;
239 
240  if (real_q_index == numbers::invalid_unsigned_int)
241  {
242  if (this->q_collection.size() > 1)
243  real_q_index = cell->active_fe_index();
244  else
245  real_q_index = 0;
246  }
247 
248  if (real_mapping_index == numbers::invalid_unsigned_int)
249  {
250  if (this->mapping_collection->size() > 1)
251  real_mapping_index = cell->active_fe_index();
252  else
253  real_mapping_index = 0;
254  }
255 
256  if (real_fe_index == numbers::invalid_unsigned_int)
257  real_fe_index = cell->active_fe_index();
258 
259  // some checks
260  AssertIndexRange(real_q_index, this->q_collection.size());
261  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
262  AssertIndexRange(real_fe_index, this->fe_collection->size());
263 
264  // now finally actually get the
265  // corresponding object and
266  // initialize it
267  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
268  .reinit(cell);
269  }
270 
271 
272 
273  template <int dim, int spacedim>
274  void
275  FEValues<dim, spacedim>::reinit(
276  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
277  const unsigned int q_index,
278  const unsigned int mapping_index,
279  const unsigned int fe_index)
280  {
281  // determine which indices we
282  // should actually use
283  unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
284  real_fe_index = fe_index;
285 
286  if (real_q_index == numbers::invalid_unsigned_int)
287  real_q_index = 0;
288 
289  if (real_mapping_index == numbers::invalid_unsigned_int)
290  real_mapping_index = 0;
291 
292  if (real_fe_index == numbers::invalid_unsigned_int)
293  real_fe_index = 0;
294 
295  // some checks
296  AssertIndexRange(real_q_index, this->q_collection.size());
297  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
298  AssertIndexRange(real_fe_index, this->fe_collection->size());
299 
300  // now finally actually get the
301  // corresponding object and
302  // initialize it
303  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
304  .reinit(cell);
305  }
306 
307 
308  // -------------------------- FEFaceValues -------------------------
309 
310 
311  template <int dim, int spacedim>
312  FEFaceValues<dim, spacedim>::FEFaceValues(
313  const hp::MappingCollection<dim, spacedim> &mapping,
314  const hp::FECollection<dim, spacedim> & fe_collection,
315  const hp::QCollection<dim - 1> & q_collection,
316  const UpdateFlags update_flags)
317  : hp::FEValuesBase<dim, dim - 1, ::FEFaceValues<dim, spacedim>>(
318  mapping,
319  fe_collection,
320  q_collection,
321  update_flags)
322  {}
323 
324 
325  template <int dim, int spacedim>
326  FEFaceValues<dim, spacedim>::FEFaceValues(
327  const hp::FECollection<dim, spacedim> &fe_collection,
328  const hp::QCollection<dim - 1> & q_collection,
329  const UpdateFlags update_flags)
330  : hp::FEValuesBase<dim, dim - 1, ::FEFaceValues<dim, spacedim>>(
331  fe_collection,
332  q_collection,
333  update_flags)
334  {}
335 
336 
337  template <int dim, int spacedim>
338  template <typename DoFHandlerType, bool lda>
339  void
340  FEFaceValues<dim, spacedim>::reinit(
341  const TriaIterator<DoFCellAccessor<DoFHandlerType, lda>> cell,
342  const unsigned int face_no,
343  const unsigned int q_index,
344  const unsigned int mapping_index,
345  const unsigned int fe_index)
346  {
347  // determine which indices we
348  // should actually use
349  unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
350  real_fe_index = fe_index;
351 
352  if (real_q_index == numbers::invalid_unsigned_int)
353  {
354  if (this->q_collection.size() > 1)
355  real_q_index = cell->active_fe_index();
356  else
357  real_q_index = 0;
358  }
359 
360  if (real_mapping_index == numbers::invalid_unsigned_int)
361  {
362  if (this->mapping_collection->size() > 1)
363  real_mapping_index = cell->active_fe_index();
364  else
365  real_mapping_index = 0;
366  }
367 
368  if (real_fe_index == numbers::invalid_unsigned_int)
369  real_fe_index = cell->active_fe_index();
370 
371  // some checks
372  AssertIndexRange(real_q_index, this->q_collection.size());
373  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
374  AssertIndexRange(real_fe_index, this->fe_collection->size());
375 
376  // now finally actually get the
377  // corresponding object and
378  // initialize it
379  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
380  .reinit(cell, face_no);
381  }
382 
383 
384 
385  template <int dim, int spacedim>
386  void
387  FEFaceValues<dim, spacedim>::reinit(
388  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
389  const unsigned int face_no,
390  const unsigned int q_index,
391  const unsigned int mapping_index,
392  const unsigned int fe_index)
393  {
394  // determine which indices we
395  // should actually use
396  unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
397  real_fe_index = fe_index;
398 
399  if (real_q_index == numbers::invalid_unsigned_int)
400  real_q_index = 0;
401 
402  if (real_mapping_index == numbers::invalid_unsigned_int)
403  real_mapping_index = 0;
404 
405  if (real_fe_index == numbers::invalid_unsigned_int)
406  real_fe_index = 0;
407 
408  // some checks
409  AssertIndexRange(real_q_index, this->q_collection.size());
410  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
411  AssertIndexRange(real_fe_index, this->fe_collection->size());
412 
413  // now finally actually get the
414  // corresponding object and
415  // initialize it
416  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
417  .reinit(cell, face_no);
418  }
419 
420 
421  // -------------------------- FESubfaceValues -------------------------
422 
423 
424  template <int dim, int spacedim>
425  FESubfaceValues<dim, spacedim>::FESubfaceValues(
426  const hp::MappingCollection<dim, spacedim> &mapping,
427  const hp::FECollection<dim, spacedim> & fe_collection,
428  const hp::QCollection<dim - 1> & q_collection,
429  const UpdateFlags update_flags)
430  : hp::FEValuesBase<dim, dim - 1, ::FESubfaceValues<dim, spacedim>>(
431  mapping,
432  fe_collection,
433  q_collection,
434  update_flags)
435  {}
436 
437 
438  template <int dim, int spacedim>
439  FESubfaceValues<dim, spacedim>::FESubfaceValues(
440  const hp::FECollection<dim, spacedim> &fe_collection,
441  const hp::QCollection<dim - 1> & q_collection,
442  const UpdateFlags update_flags)
443  : hp::FEValuesBase<dim, dim - 1, ::FESubfaceValues<dim, spacedim>>(
444  fe_collection,
445  q_collection,
446  update_flags)
447  {}
448 
449 
450  template <int dim, int spacedim>
451  template <typename DoFHandlerType, bool lda>
452  void
453  FESubfaceValues<dim, spacedim>::reinit(
454  const TriaIterator<DoFCellAccessor<DoFHandlerType, lda>> cell,
455  const unsigned int face_no,
456  const unsigned int subface_no,
457  const unsigned int q_index,
458  const unsigned int mapping_index,
459  const unsigned int fe_index)
460  {
461  // determine which indices we
462  // should actually use
463  unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
464  real_fe_index = fe_index;
465 
466  if (real_q_index == numbers::invalid_unsigned_int)
467  {
468  if (this->q_collection.size() > 1)
469  real_q_index = cell->active_fe_index();
470  else
471  real_q_index = 0;
472  }
473 
474  if (real_mapping_index == numbers::invalid_unsigned_int)
475  {
476  if (this->mapping_collection->size() > 1)
477  real_mapping_index = cell->active_fe_index();
478  else
479  real_mapping_index = 0;
480  }
481 
482  if (real_fe_index == numbers::invalid_unsigned_int)
483  real_fe_index = cell->active_fe_index();
484 
485  // some checks
486  AssertIndexRange(real_q_index, this->q_collection.size());
487  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
488  AssertIndexRange(real_fe_index, this->fe_collection->size());
489 
490  // now finally actually get the
491  // corresponding object and
492  // initialize it
493  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
494  .reinit(cell, face_no, subface_no);
495  }
496 
497 
498 
499  template <int dim, int spacedim>
500  void
501  FESubfaceValues<dim, spacedim>::reinit(
502  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
503  const unsigned int face_no,
504  const unsigned int subface_no,
505  const unsigned int q_index,
506  const unsigned int mapping_index,
507  const unsigned int fe_index)
508  {
509  // determine which indices we
510  // should actually use
511  unsigned int real_q_index = q_index, real_mapping_index = mapping_index,
512  real_fe_index = fe_index;
513 
514  if (real_q_index == numbers::invalid_unsigned_int)
515  real_q_index = 0;
516 
517  if (real_mapping_index == numbers::invalid_unsigned_int)
518  real_mapping_index = 0;
519 
520  if (real_fe_index == numbers::invalid_unsigned_int)
521  real_fe_index = 0;
522 
523  // some checks
524  AssertIndexRange(real_q_index, this->q_collection.size());
525  AssertIndexRange(real_mapping_index, this->mapping_collection->size());
526  AssertIndexRange(real_fe_index, this->fe_collection->size());
527 
528  // now finally actually get the
529  // corresponding object and
530  // initialize it
531  this->select_fe_values(real_fe_index, real_mapping_index, real_q_index)
532  .reinit(cell, face_no, subface_no);
533  }
534 } // namespace hp
535 
536 
537 // explicit instantiations
538 #include "fe_values.inst"
539 
540 
541 DEAL_II_NAMESPACE_CLOSE
hp::FEValues
Definition: fe_values.h:281
hp::FEValues::FEValues
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:203
TableIndices< 3 >
FEValuesBase::update_flags
UpdateFlags update_flags
Definition: fe_values.h:3556
mapping_q1.h
FEValues::FEValues
FEValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim > &quadrature, const UpdateFlags update_flags)
Definition: fe_values.cc:4390
hp::FEValuesBase
Definition: fe_values.h:66
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
hp::QCollection< q_dim >
Threads::new_task
Task< RT > new_task(const std::function< RT()> &function)
Definition: thread_management.h:1647
thread_management.h
hp::StaticMappingQ1
Definition: mapping_collection.h:142
hp::FESubfaceValues::FESubfaceValues
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:425
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
FEValues::reinit
void reinit(const TriaIterator< DoFCellAccessor< DoFHandlerType< dim, spacedim >, level_dof_access >> &cell)
DoFCellAccessor
Definition: dof_accessor.h:1321
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
hp::FEValuesBase::fe_values_table
Table< 3, std::unique_ptr< FEValuesType > > fe_values_table
Definition: fe_values.h:211
numbers
Definition: numbers.h:207
hp::MappingCollection< dim, FEValuesType::space_dimension >
UpdateFlags
UpdateFlags
Definition: fe_update_flags.h:66
FEFaceValues::reinit
void reinit(const TriaIterator< DoFCellAccessor< DoFHandlerType< dim, spacedim >, level_dof_access >> &cell, const unsigned int face_no)
Threads::TaskGroup::join_all
void join_all() const
Definition: thread_management.h:1833
hp::FEValuesBase::FEValuesBase
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:30
hp::FESubfaceValues
Definition: fe_values.h:515
Threads::TaskGroup
Definition: thread_management.h:1798
hp
Definition: hp.h:117
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:191
hp::FEFaceValues
Definition: fe_values.h:408
FESubfaceValues::FESubfaceValues
FESubfaceValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim - 1 > &face_quadrature, const UpdateFlags update_flags)
Definition: fe_values.cc:4856
FEFaceValues::FEFaceValues
FEFaceValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim - 1 > &quadrature, const UpdateFlags update_flags)
Definition: fe_values.cc:4657
memory.h
fe_values.h
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
hp::FESubfaceValues::reinit
void reinit(const TriaIterator< DoFCellAccessor< DoFHandlerType, 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:453
TriaIterator
Definition: tria_iterator.h:578
FEValuesBase::FEValuesBase
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:3076
hp::FECollection< dim, FEValuesType::space_dimension >
Generated by   doxygen 1.8.17