Reference documentation for deal.II version GIT 29f9da0a34 2023-12-07 10:00:01+00:00
\(\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) 1998 - 2023 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 
19 #include <deal.II/base/numbers.h>
23 
25 
27 
28 #include <deal.II/fe/fe.h>
29 #include <deal.II/fe/fe_values.h>
30 #include <deal.II/fe/mapping.h>
31 
34 
35 #include <deal.II/lac/vector.h>
36 
37 #include <boost/container/small_vector.hpp>
38 
39 #include <iomanip>
40 #include <memory>
41 #include <type_traits>
42 
44 
45 
46 namespace internal
47 {
48  template <int dim, int spacedim>
49  inline std::vector<unsigned int>
51  {
52  std::vector<unsigned int> shape_function_to_row_table(
54  unsigned int row = 0;
55  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
56  {
57  // loop over all components that are nonzero for this particular
58  // shape function. if a component is zero then we leave the
59  // value in the table unchanged (at the invalid value)
60  // otherwise it is mapped to the next free entry
61  unsigned int nth_nonzero_component = 0;
62  for (unsigned int c = 0; c < fe.n_components(); ++c)
63  if (fe.get_nonzero_components(i)[c] == true)
64  {
65  shape_function_to_row_table[i * fe.n_components() + c] =
66  row + nth_nonzero_component;
67  ++nth_nonzero_component;
68  }
69  row += fe.n_nonzero_components(i);
70  }
71 
72  return shape_function_to_row_table;
73  }
74 } // namespace internal
75 
76 
77 
78 namespace internal
79 {
80  namespace FEValuesImplementation
81  {
82  template <int dim, int spacedim>
83  void
85  const unsigned int n_quadrature_points,
87  const UpdateFlags flags)
88  {
89  // initialize the table mapping from shape function number to
90  // the rows in the tables storing the data by shape function and
91  // nonzero component
92  this->shape_function_to_row_table =
94 
95  // count the total number of non-zero components accumulated
96  // over all shape functions
97  unsigned int n_nonzero_shape_components = 0;
98  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
99  n_nonzero_shape_components += fe.n_nonzero_components(i);
100  Assert(n_nonzero_shape_components >= fe.n_dofs_per_cell(),
101  ExcInternalError());
102 
103  // with the number of rows now known, initialize those fields
104  // that we will need to their correct size
105  if (flags & update_values)
106  {
107  this->shape_values.reinit(n_nonzero_shape_components,
108  n_quadrature_points);
109  this->shape_values.fill(numbers::signaling_nan<double>());
110  }
111 
112  if (flags & update_gradients)
113  {
114  this->shape_gradients.reinit(n_nonzero_shape_components,
115  n_quadrature_points);
116  this->shape_gradients.fill(
118  }
119 
120  if (flags & update_hessians)
121  {
122  this->shape_hessians.reinit(n_nonzero_shape_components,
123  n_quadrature_points);
124  this->shape_hessians.fill(
126  }
127 
128  if (flags & update_3rd_derivatives)
129  {
130  this->shape_3rd_derivatives.reinit(n_nonzero_shape_components,
131  n_quadrature_points);
132  this->shape_3rd_derivatives.fill(
134  }
135  }
136 
137 
138 
139  template <int dim, int spacedim>
140  std::size_t
142  {
143  return (
145  MemoryConsumption::memory_consumption(shape_gradients) +
146  MemoryConsumption::memory_consumption(shape_hessians) +
147  MemoryConsumption::memory_consumption(shape_3rd_derivatives) +
148  MemoryConsumption::memory_consumption(shape_function_to_row_table));
149  }
150  } // namespace FEValuesImplementation
151 } // namespace internal
152 
153 /*------------------------------- FEValues -------------------------------*/
154 #ifndef DOXYGEN
155 
156 template <int dim, int spacedim>
158 
159 
160 
161 template <int dim, int spacedim>
164  const Quadrature<dim> &q,
165  const UpdateFlags update_flags)
166  : FEValuesBase<dim, spacedim>(q.size(),
167  fe.n_dofs_per_cell(),
169  mapping,
170  fe)
171  , quadrature(q)
172 {
173  initialize(update_flags);
174 }
175 
176 
177 
178 template <int dim, int spacedim>
181  const hp::QCollection<dim> &q,
182  const UpdateFlags update_flags)
183  : FEValues(mapping, fe, q[0], update_flags)
184 {
185  AssertDimension(q.size(), 1);
186 }
187 
188 
189 
190 template <int dim, int spacedim>
192  const Quadrature<dim> &q,
193  const UpdateFlags update_flags)
194  : FEValuesBase<dim, spacedim>(
195  q.size(),
196  fe.n_dofs_per_cell(),
198  fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
199  fe)
200  , quadrature(q)
201 {
202  initialize(update_flags);
203 }
204 
205 
206 
207 template <int dim, int spacedim>
209  const hp::QCollection<dim> &q,
210  const UpdateFlags update_flags)
211  : FEValues(fe, q[0], update_flags)
212 {
213  AssertDimension(q.size(), 1);
214 }
215 
216 
217 
218 template <int dim, int spacedim>
219 void
221 {
222  // You can compute normal vectors to the cells only in the
223  // codimension one case.
224  if (dim != spacedim - 1)
225  Assert((update_flags & update_normal_vectors) == false,
226  ExcMessage("You can only pass the 'update_normal_vectors' "
227  "flag to FEFaceValues or FESubfaceValues objects, "
228  "but not to an FEValues object unless the "
229  "triangulation it refers to is embedded in a higher "
230  "dimensional space."));
231 
232  const UpdateFlags flags = this->compute_update_flags(update_flags);
233 
234  // initialize the base classes
235  if (flags & update_mapping)
236  this->mapping_output.initialize(this->max_n_quadrature_points, flags);
237  this->finite_element_output.initialize(this->max_n_quadrature_points,
238  *this->fe,
239  flags);
240 
241  // then get objects into which the FE and the Mapping can store
242  // intermediate data used across calls to reinit. we can do this in parallel
244  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>>
245  fe_get_data = Threads::new_task([&]() {
246  return this->fe->get_data(flags,
247  *this->mapping,
248  quadrature,
249  this->finite_element_output);
250  });
251 
253  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
254  mapping_get_data;
255  if (flags & update_mapping)
256  mapping_get_data = Threads::new_task(
257  [&]() { return this->mapping->get_data(flags, quadrature); });
258 
259  this->update_flags = flags;
260 
261  // then collect answers from the two task above
262  this->fe_data = std::move(fe_get_data.return_value());
263  if (flags & update_mapping)
264  this->mapping_data = std::move(mapping_get_data.return_value());
265  else
266  this->mapping_data =
267  std::make_unique<typename Mapping<dim, spacedim>::InternalDataBase>();
268 }
269 
270 
271 
272 template <int dim, int spacedim>
273 void
275  const typename Triangulation<dim, spacedim>::cell_iterator &cell)
276 {
277  // Check that mapping and reference cell type are compatible:
278  Assert(this->get_mapping().is_compatible_with(cell->reference_cell()),
279  ExcMessage(
280  "You are trying to call FEValues::reinit() with a cell of type " +
281  cell->reference_cell().to_string() +
282  " with a Mapping that is not compatible with it."));
283 
284  // no FE in this cell, so no assertion
285  // necessary here
286  this->maybe_invalidate_previous_present_cell(cell);
287  this->check_cell_similarity(cell);
288 
289  this->present_cell = {cell};
290 
291  // this was the part of the work that is dependent on the actual
292  // data type of the iterator. now pass on to the function doing
293  // the real work.
294  do_reinit();
295 }
296 
297 
298 
299 template <int dim, int spacedim>
300 template <bool lda>
301 void
304 {
305  // assert that the finite elements passed to the constructor and
306  // used by the DoFHandler used by this cell, are the same
307  Assert(static_cast<const FiniteElementData<dim> &>(*this->fe) ==
308  static_cast<const FiniteElementData<dim> &>(cell->get_fe()),
310 
311  // Check that mapping and reference cell type are compatible:
312  Assert(this->get_mapping().is_compatible_with(cell->reference_cell()),
313  ExcMessage(
314  "You are trying to call FEValues::reinit() with a cell of type " +
315  cell->reference_cell().to_string() +
316  " with a Mapping that is not compatible with it."));
317 
318  this->maybe_invalidate_previous_present_cell(cell);
319  this->check_cell_similarity(cell);
320 
321  this->present_cell = {cell};
322 
323  // this was the part of the work that is dependent on the actual
324  // data type of the iterator. now pass on to the function doing
325  // the real work.
326  do_reinit();
327 }
328 
329 
330 
331 template <int dim, int spacedim>
332 void
334 {
335  // first call the mapping and let it generate the data
336  // specific to the mapping. also let it inspect the
337  // cell similarity flag and, if necessary, update
338  // it
339  if (this->update_flags & update_mapping)
340  {
341  this->cell_similarity =
342  this->get_mapping().fill_fe_values(this->present_cell,
343  this->cell_similarity,
344  quadrature,
345  *this->mapping_data,
346  this->mapping_output);
347  }
348 
349  // then call the finite element and, with the data
350  // already filled by the mapping, let it compute the
351  // data for the mapped shape function values, gradients,
352  // etc.
353  this->get_fe().fill_fe_values(this->present_cell,
354  this->cell_similarity,
355  this->quadrature,
356  this->get_mapping(),
357  *this->mapping_data,
358  this->mapping_output,
359  *this->fe_data,
360  this->finite_element_output);
361 }
362 
363 
364 
365 template <int dim, int spacedim>
366 std::size_t
368 {
371 }
372 
373 #endif
374 /*------------------------------- FEFaceValuesBase --------------------------*/
375 #ifndef DOXYGEN
376 
377 template <int dim, int spacedim>
379  const unsigned int dofs_per_cell,
380  const UpdateFlags flags,
381  const Mapping<dim, spacedim> &mapping,
383  const Quadrature<dim - 1> &quadrature)
384  : FEFaceValuesBase<dim, spacedim>(dofs_per_cell,
385  flags,
386  mapping,
387  fe,
388  hp::QCollection<dim - 1>(quadrature))
389 {}
390 
391 
392 
393 template <int dim, int spacedim>
395  const unsigned int dofs_per_cell,
396  const UpdateFlags,
397  const Mapping<dim, spacedim> &mapping,
399  const hp::QCollection<dim - 1> &quadrature)
400  : FEValuesBase<dim, spacedim>(quadrature.max_n_quadrature_points(),
401  dofs_per_cell,
403  mapping,
404  fe)
405  , present_face_index(numbers::invalid_unsigned_int)
406  , quadrature(quadrature)
407 {
408  Assert(quadrature.size() == 1 ||
409  quadrature.size() == fe.reference_cell().n_faces(),
410  ExcInternalError());
411 }
412 
413 
414 
415 template <int dim, int spacedim>
416 const std::vector<Tensor<1, spacedim>> &
418 {
419  Assert(this->update_flags & update_boundary_forms,
421  "update_boundary_forms")));
422  return this->mapping_output.boundary_forms;
423 }
424 
425 
426 
427 template <int dim, int spacedim>
428 std::size_t
430 {
433 }
434 #endif
435 
436 /*------------------------------- FEFaceValues -------------------------------*/
437 
438 #ifndef DOXYGEN
439 
440 template <int dim, int spacedim>
441 const unsigned int FEFaceValues<dim, spacedim>::dimension;
442 
443 
444 
445 template <int dim, int spacedim>
447 
448 
449 
450 template <int dim, int spacedim>
452  const Mapping<dim, spacedim> &mapping,
454  const Quadrature<dim - 1> &quadrature,
455  const UpdateFlags update_flags)
456  : FEFaceValues<dim, spacedim>(mapping,
457  fe,
458  hp::QCollection<dim - 1>(quadrature),
459  update_flags)
460 {}
461 
462 
463 
464 template <int dim, int spacedim>
466  const Mapping<dim, spacedim> &mapping,
468  const hp::QCollection<dim - 1> &quadrature,
469  const UpdateFlags update_flags)
470  : FEFaceValuesBase<dim, spacedim>(fe.n_dofs_per_cell(),
471  update_flags,
472  mapping,
473  fe,
474  quadrature)
475 {
476  initialize(update_flags);
477 }
478 
479 
480 
481 template <int dim, int spacedim>
484  const Quadrature<dim - 1> &quadrature,
485  const UpdateFlags update_flags)
486  : FEFaceValues<dim, spacedim>(fe,
487  hp::QCollection<dim - 1>(quadrature),
488  update_flags)
489 {}
490 
491 
492 
493 template <int dim, int spacedim>
496  const hp::QCollection<dim - 1> &quadrature,
497  const UpdateFlags update_flags)
498  : FEFaceValuesBase<dim, spacedim>(
499  fe.n_dofs_per_cell(),
500  update_flags,
501  fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
502  fe,
503  quadrature)
504 {
505  initialize(update_flags);
506 }
507 
508 
509 
510 template <int dim, int spacedim>
511 void
513 {
514  const UpdateFlags flags = this->compute_update_flags(update_flags);
515 
516  // initialize the base classes
517  if (flags & update_mapping)
518  this->mapping_output.initialize(this->max_n_quadrature_points, flags);
519  this->finite_element_output.initialize(this->max_n_quadrature_points,
520  *this->fe,
521  flags);
522 
523  // then get objects into which the FE and the Mapping can store
524  // intermediate data used across calls to reinit. this can be done in parallel
525 
526  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase> (
527  FiniteElement<dim, spacedim>::*finite_element_get_face_data)(
528  const UpdateFlags,
529  const Mapping<dim, spacedim> &,
530  const hp::QCollection<dim - 1> &,
532  spacedim>
534 
535  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> (
536  Mapping<dim, spacedim>::*mapping_get_face_data)(
537  const UpdateFlags, const hp::QCollection<dim - 1> &) const =
539 
540 
542  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>>
543  fe_get_data = Threads::new_task(finite_element_get_face_data,
544  *this->fe,
545  flags,
546  *this->mapping,
547  this->quadrature,
548  this->finite_element_output);
550  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
551  mapping_get_data;
552  if (flags & update_mapping)
553  mapping_get_data = Threads::new_task(mapping_get_face_data,
554  *this->mapping,
555  flags,
556  this->quadrature);
557 
558  this->update_flags = flags;
559 
560  // then collect answers from the two task above
561  this->fe_data = std::move(fe_get_data.return_value());
562  if (flags & update_mapping)
563  this->mapping_data = std::move(mapping_get_data.return_value());
564  else
565  this->mapping_data =
566  std::make_unique<typename Mapping<dim, spacedim>::InternalDataBase>();
567 }
568 
569 
570 
571 template <int dim, int spacedim>
572 template <bool lda>
573 void
576  const unsigned int face_no)
577 {
578  // assert that the finite elements passed to the constructor and
579  // used by the DoFHandler used by this cell, are the same
580  Assert(static_cast<const FiniteElementData<dim> &>(*this->fe) ==
581  static_cast<const FiniteElementData<dim> &>(
582  cell->get_dof_handler().get_fe(cell->active_fe_index())),
584 
586 
587  this->maybe_invalidate_previous_present_cell(cell);
588  this->present_cell = {cell};
589 
590  // this was the part of the work that is dependent on the actual
591  // data type of the iterator. now pass on to the function doing
592  // the real work.
593  do_reinit(face_no);
594 }
595 
596 
597 
598 template <int dim, int spacedim>
599 template <bool lda>
600 void
603  const typename Triangulation<dim, spacedim>::face_iterator &face)
604 {
605  const auto face_n = cell->face_iterator_to_index(face);
606  reinit(cell, face_n);
607 }
608 
609 
610 
611 template <int dim, int spacedim>
612 void
614  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
615  const unsigned int face_no)
616 {
618 
619  this->maybe_invalidate_previous_present_cell(cell);
620  this->present_cell = {cell};
621 
622  // this was the part of the work that is dependent on the actual
623  // data type of the iterator. now pass on to the function doing
624  // the real work.
625  do_reinit(face_no);
626 }
627 
628 
629 
630 template <int dim, int spacedim>
631 void
633  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
634  const typename Triangulation<dim, spacedim>::face_iterator &face)
635 {
636  const auto face_n = cell->face_iterator_to_index(face);
637  reinit(cell, face_n);
638 }
639 
640 
641 
642 template <int dim, int spacedim>
643 void
644 FEFaceValues<dim, spacedim>::do_reinit(const unsigned int face_no)
645 {
646  this->present_face_no = face_no;
647 
648  // first of all, set the present_face_index (if available)
649  const typename Triangulation<dim, spacedim>::cell_iterator cell =
650  this->present_cell;
651  this->present_face_index = cell->face_index(face_no);
652 
653  if (this->update_flags & update_mapping)
654  {
655  this->get_mapping().fill_fe_face_values(this->present_cell,
656  face_no,
657  this->quadrature,
658  *this->mapping_data,
659  this->mapping_output);
660  }
661 
662  this->get_fe().fill_fe_face_values(this->present_cell,
663  face_no,
664  this->quadrature,
665  this->get_mapping(),
666  *this->mapping_data,
667  this->mapping_output,
668  *this->fe_data,
669  this->finite_element_output);
670 
671  const_cast<unsigned int &>(this->n_quadrature_points) =
672  this->quadrature[this->quadrature.size() == 1 ? 0 : face_no].size();
673 }
674 
675 
676 /* ---------------------------- FESubFaceValues ---------------------------- */
677 
678 
679 template <int dim, int spacedim>
681 
682 
683 
684 template <int dim, int spacedim>
686 
687 
688 
689 template <int dim, int spacedim>
691  const Mapping<dim, spacedim> &mapping,
693  const Quadrature<dim - 1> &quadrature,
694  const UpdateFlags update_flags)
695  : FEFaceValuesBase<dim, spacedim>(fe.n_dofs_per_cell(),
696  update_flags,
697  mapping,
698  fe,
699  quadrature)
700 {
701  initialize(update_flags);
702 }
703 
704 
705 
706 template <int dim, int spacedim>
708  const Mapping<dim, spacedim> &mapping,
710  const hp::QCollection<dim - 1> &quadrature,
711  const UpdateFlags update_flags)
712  : FESubfaceValues(mapping, fe, quadrature[0], update_flags)
713 {
714  AssertDimension(quadrature.size(), 1);
715 }
716 
717 
718 
719 template <int dim, int spacedim>
722  const Quadrature<dim - 1> &quadrature,
723  const UpdateFlags update_flags)
724  : FEFaceValuesBase<dim, spacedim>(
725  fe.n_dofs_per_cell(),
726  update_flags,
727  fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
728  fe,
729  quadrature)
730 {
731  initialize(update_flags);
732 }
733 
734 
735 
736 template <int dim, int spacedim>
739  const hp::QCollection<dim - 1> &quadrature,
740  const UpdateFlags update_flags)
741  : FESubfaceValues(fe, quadrature[0], update_flags)
742 {
743  AssertDimension(quadrature.size(), 1);
744 }
745 
746 
747 
748 template <int dim, int spacedim>
749 void
751 {
752  const UpdateFlags flags = this->compute_update_flags(update_flags);
753 
754  // initialize the base classes
755  if (flags & update_mapping)
756  this->mapping_output.initialize(this->max_n_quadrature_points, flags);
757  this->finite_element_output.initialize(this->max_n_quadrature_points,
758  *this->fe,
759  flags);
760 
761  // then get objects into which the FE and the Mapping can store
762  // intermediate data used across calls to reinit. this can be done
763  // in parallel
765  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>>
766  fe_get_data =
768  *this->fe,
769  flags,
770  *this->mapping,
771  this->quadrature[0],
772  this->finite_element_output);
774  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
775  mapping_get_data;
776  if (flags & update_mapping)
777  mapping_get_data =
779  *this->mapping,
780  flags,
781  this->quadrature[0]);
782 
783  this->update_flags = flags;
784 
785  // then collect answers from the two task above
786  this->fe_data = std::move(fe_get_data.return_value());
787  if (flags & update_mapping)
788  this->mapping_data = std::move(mapping_get_data.return_value());
789  else
790  this->mapping_data =
791  std::make_unique<typename Mapping<dim, spacedim>::InternalDataBase>();
792 }
793 
794 
795 
796 template <int dim, int spacedim>
797 template <bool lda>
798 void
801  const unsigned int face_no,
802  const unsigned int subface_no)
803 {
804  // assert that the finite elements passed to the constructor and
805  // used by the DoFHandler used by this cell, are the same
806  Assert(static_cast<const FiniteElementData<dim> &>(*this->fe) ==
807  static_cast<const FiniteElementData<dim> &>(
808  cell->get_dof_handler().get_fe(cell->active_fe_index())),
811  // We would like to check for subface_no < cell->face(face_no)->n_children(),
812  // but unfortunately the current function is also called for
813  // faces without children (see tests/fe/mapping.cc). Therefore,
814  // we must use following workaround of two separate assertions
815  Assert(cell->face(face_no)->has_children() ||
817  ExcIndexRange(subface_no,
818  0,
820  Assert(!cell->face(face_no)->has_children() ||
821  subface_no < cell->face(face_no)->n_active_descendants(),
822  ExcIndexRange(subface_no,
823  0,
824  cell->face(face_no)->n_active_descendants()));
825  Assert(cell->has_children() == false,
826  ExcMessage("You can't use subface data for cells that are "
827  "already refined. Iterate over their children "
828  "instead in these cases."));
829 
830  this->maybe_invalidate_previous_present_cell(cell);
831  this->present_cell = {cell};
832 
833  // this was the part of the work that is dependent on the actual
834  // data type of the iterator. now pass on to the function doing
835  // the real work.
836  do_reinit(face_no, subface_no);
837 }
838 
839 
840 
841 template <int dim, int spacedim>
842 template <bool lda>
843 void
846  const typename Triangulation<dim, spacedim>::face_iterator &face,
847  const typename Triangulation<dim, spacedim>::face_iterator &subface)
848 {
849  reinit(cell,
850  cell->face_iterator_to_index(face),
851  face->child_iterator_to_index(subface));
852 }
853 
854 
855 
856 template <int dim, int spacedim>
857 void
859  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
860  const unsigned int face_no,
861  const unsigned int subface_no)
862 {
864  // We would like to check for subface_no < cell->face(face_no)->n_children(),
865  // but unfortunately the current function is also called for
866  // faces without children for periodic faces, which have hanging nodes on
867  // the other side (see include/deal.II/matrix_free/mapping_info.templates.h).
868  AssertIndexRange(subface_no,
869  (cell->has_periodic_neighbor(face_no) ?
870  cell->periodic_neighbor(face_no)
871  ->face(cell->periodic_neighbor_face_no(face_no))
872  ->n_children() :
873  cell->face(face_no)->n_children()));
874 
875  this->maybe_invalidate_previous_present_cell(cell);
876  this->present_cell = {cell};
877 
878  // this was the part of the work that is dependent on the actual
879  // data type of the iterator. now pass on to the function doing
880  // the real work.
881  do_reinit(face_no, subface_no);
882 }
883 
884 
885 
886 template <int dim, int spacedim>
887 void
889  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
890  const typename Triangulation<dim, spacedim>::face_iterator &face,
891  const typename Triangulation<dim, spacedim>::face_iterator &subface)
892 {
893  reinit(cell,
894  cell->face_iterator_to_index(face),
895  face->child_iterator_to_index(subface));
896 }
897 
898 
899 
900 template <int dim, int spacedim>
901 void
902 FESubfaceValues<dim, spacedim>::do_reinit(const unsigned int face_no,
903  const unsigned int subface_no)
904 {
905  this->present_face_no = face_no;
906 
907  // first of all, set the present_face_index (if available)
908  const typename Triangulation<dim, spacedim>::cell_iterator cell =
909  this->present_cell;
910 
911  if (!cell->face(face_no)->has_children())
912  // no subfaces at all, so set present_face_index to this face rather
913  // than any subface
914  this->present_face_index = cell->face_index(face_no);
915  else if (dim != 3)
916  this->present_face_index = cell->face(face_no)->child_index(subface_no);
917  else
918  {
919  // this is the same logic we use in cell->neighbor_child_on_subface(). See
920  // there for an explanation of the different cases
921  unsigned int subface_index = numbers::invalid_unsigned_int;
922  switch (cell->subface_case(face_no))
923  {
927  subface_index = cell->face(face_no)->child_index(subface_no);
928  break;
931  subface_index = cell->face(face_no)
932  ->child(subface_no / 2)
933  ->child_index(subface_no % 2);
934  break;
937  switch (subface_no)
938  {
939  case 0:
940  case 1:
941  subface_index =
942  cell->face(face_no)->child(0)->child_index(subface_no);
943  break;
944  case 2:
945  subface_index = cell->face(face_no)->child_index(1);
946  break;
947  default:
948  Assert(false, ExcInternalError());
949  }
950  break;
953  switch (subface_no)
954  {
955  case 0:
956  subface_index = cell->face(face_no)->child_index(0);
957  break;
958  case 1:
959  case 2:
960  subface_index =
961  cell->face(face_no)->child(1)->child_index(subface_no - 1);
962  break;
963  default:
964  Assert(false, ExcInternalError());
965  }
966  break;
967  default:
968  Assert(false, ExcInternalError());
969  break;
970  }
971  Assert(subface_index != numbers::invalid_unsigned_int,
972  ExcInternalError());
973  this->present_face_index = subface_index;
974  }
975 
976  // now ask the mapping and the finite element to do the actual work
977  if (this->update_flags & update_mapping)
978  {
979  this->get_mapping().fill_fe_subface_values(this->present_cell,
980  face_no,
981  subface_no,
982  this->quadrature[0],
983  *this->mapping_data,
984  this->mapping_output);
985  }
986 
987  this->get_fe().fill_fe_subface_values(this->present_cell,
988  face_no,
989  subface_no,
990  this->quadrature[0],
991  this->get_mapping(),
992  *this->mapping_data,
993  this->mapping_output,
994  *this->fe_data,
995  this->finite_element_output);
996 }
997 #endif
998 
999 /*------------------------- Explicit Instantiations --------------------------*/
1000 
1001 #include "fe_values.inst"
1002 
const std::vector< Tensor< 1, spacedim > > & get_boundary_forms() const
FEFaceValuesBase(const unsigned int dofs_per_cell, const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim - 1 > &quadrature)
std::size_t memory_consumption() const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access >> &cell, const unsigned int face_no)
void initialize(const UpdateFlags update_flags)
void do_reinit(const unsigned int face_no)
FEFaceValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim - 1 > &quadrature, const UpdateFlags update_flags)
FESubfaceValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim - 1 > &face_quadrature, const UpdateFlags update_flags)
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access >> &cell, const unsigned int face_no, const unsigned int subface_no)
void initialize(const UpdateFlags update_flags)
void do_reinit(const unsigned int face_no, const unsigned int subface_no)
FEValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim > &quadrature, const UpdateFlags update_flags)
void do_reinit()
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access >> &cell)
void initialize(const UpdateFlags update_flags)
std::size_t memory_consumption() const
unsigned int n_dofs_per_cell() const
unsigned int n_components() const
ReferenceCell reference_cell() const
const ComponentMask & get_nonzero_components(const unsigned int i) const
unsigned int n_nonzero_components(const unsigned int i) const
unsigned int n_faces() const
unsigned int size() const
Definition: collection.h:265
void initialize(const unsigned int n_quadrature_points, const FiniteElement< dim, spacedim > &fe, const UpdateFlags flags)
Definition: fe_values.cc:84
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1631
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcMessage(std::string arg1)
UpdateFlags
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_3rd_derivatives
Third derivatives of shape functions.
@ update_mapping
@ update_gradients
Shape function gradients.
@ update_default
No update.
@ update_boundary_forms
Outer normal vector, not normalized.
Task< RT > new_task(const std::function< RT()> &function)
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Definition: mapping.cc:285
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
Definition: hp.h:118
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
std::vector< unsigned int > make_shape_function_to_row_table(const FiniteElement< dim, spacedim > &fe)
Definition: fe_values.cc:50
static const unsigned int invalid_unsigned_int
Definition: types.h:221
T signaling_nan()