Reference documentation for deal.II version GIT 0b65fff18a 2023-09-27 19:30:02+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 
155 template <int dim, int spacedim>
157 
158 
159 
160 template <int dim, int spacedim>
163  const Quadrature<dim> &q,
164  const UpdateFlags update_flags)
165  : FEValuesBase<dim, spacedim>(q.size(),
166  fe.n_dofs_per_cell(),
168  mapping,
169  fe)
170  , quadrature(q)
171 {
173 }
174 
175 
176 
177 template <int dim, int spacedim>
180  const hp::QCollection<dim> &q,
181  const UpdateFlags update_flags)
182  : FEValues(mapping, fe, q[0], update_flags)
183 {
184  AssertDimension(q.size(), 1);
185 }
186 
187 
188 
189 template <int dim, int spacedim>
191  const Quadrature<dim> &q,
192  const UpdateFlags update_flags)
193  : FEValuesBase<dim, spacedim>(
194  q.size(),
195  fe.n_dofs_per_cell(),
197  fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
198  fe)
199  , quadrature(q)
200 {
202 }
203 
204 
205 
206 template <int dim, int spacedim>
208  const hp::QCollection<dim> &q,
209  const UpdateFlags update_flags)
210  : FEValues(fe, q[0], update_flags)
211 {
212  AssertDimension(q.size(), 1);
213 }
214 
215 
216 
217 template <int dim, int spacedim>
218 void
220 {
221  // You can compute normal vectors to the cells only in the
222  // codimension one case.
223  if (dim != spacedim - 1)
224  Assert((update_flags & update_normal_vectors) == false,
225  ExcMessage("You can only pass the 'update_normal_vectors' "
226  "flag to FEFaceValues or FESubfaceValues objects, "
227  "but not to an FEValues object unless the "
228  "triangulation it refers to is embedded in a higher "
229  "dimensional space."));
230 
231  const UpdateFlags flags = this->compute_update_flags(update_flags);
232 
233  // initialize the base classes
234  if (flags & update_mapping)
235  this->mapping_output.initialize(this->max_n_quadrature_points, flags);
236  this->finite_element_output.initialize(this->max_n_quadrature_points,
237  *this->fe,
238  flags);
239 
240  // then get objects into which the FE and the Mapping can store
241  // intermediate data used across calls to reinit. we can do this in parallel
243  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>>
244  fe_get_data = Threads::new_task([&]() {
245  return this->fe->get_data(flags,
246  *this->mapping,
247  quadrature,
248  this->finite_element_output);
249  });
250 
252  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
253  mapping_get_data;
254  if (flags & update_mapping)
255  mapping_get_data = Threads::new_task(
256  [&]() { return this->mapping->get_data(flags, quadrature); });
257 
258  this->update_flags = flags;
259 
260  // then collect answers from the two task above
261  this->fe_data = std::move(fe_get_data.return_value());
262  if (flags & update_mapping)
263  this->mapping_data = std::move(mapping_get_data.return_value());
264  else
265  this->mapping_data =
266  std::make_unique<typename Mapping<dim, spacedim>::InternalDataBase>();
267 }
268 
269 
270 
271 template <int dim, int spacedim>
272 void
274  const typename Triangulation<dim, spacedim>::cell_iterator &cell)
275 {
276  // Check that mapping and reference cell type are compatible:
277  Assert(this->get_mapping().is_compatible_with(cell->reference_cell()),
278  ExcMessage(
279  "You are trying to call FEValues::reinit() with a cell of type " +
280  cell->reference_cell().to_string() +
281  " with a Mapping that is not compatible with it."));
282 
283  // no FE in this cell, so no assertion
284  // necessary here
285  this->maybe_invalidate_previous_present_cell(cell);
286  this->check_cell_similarity(cell);
287 
288  this->present_cell = {cell};
289 
290  // this was the part of the work that is dependent on the actual
291  // data type of the iterator. now pass on to the function doing
292  // the real work.
293  do_reinit();
294 }
295 
296 
297 
298 template <int dim, int spacedim>
299 template <bool lda>
300 void
303 {
304  // assert that the finite elements passed to the constructor and
305  // used by the DoFHandler used by this cell, are the same
306  Assert(static_cast<const FiniteElementData<dim> &>(*this->fe) ==
307  static_cast<const FiniteElementData<dim> &>(cell->get_fe()),
309 
310  // Check that mapping and reference cell type are compatible:
311  Assert(this->get_mapping().is_compatible_with(cell->reference_cell()),
312  ExcMessage(
313  "You are trying to call FEValues::reinit() with a cell of type " +
314  cell->reference_cell().to_string() +
315  " with a Mapping that is not compatible with it."));
316 
317  this->maybe_invalidate_previous_present_cell(cell);
318  this->check_cell_similarity(cell);
319 
320  this->present_cell = {cell};
321 
322  // this was the part of the work that is dependent on the actual
323  // data type of the iterator. now pass on to the function doing
324  // the real work.
325  do_reinit();
326 }
327 
328 
329 
330 template <int dim, int spacedim>
331 void
333 {
334  // first call the mapping and let it generate the data
335  // specific to the mapping. also let it inspect the
336  // cell similarity flag and, if necessary, update
337  // it
338  if (this->update_flags & update_mapping)
339  {
340  this->cell_similarity =
341  this->get_mapping().fill_fe_values(this->present_cell,
342  this->cell_similarity,
343  quadrature,
344  *this->mapping_data,
345  this->mapping_output);
346  }
347 
348  // then call the finite element and, with the data
349  // already filled by the mapping, let it compute the
350  // data for the mapped shape function values, gradients,
351  // etc.
352  this->get_fe().fill_fe_values(this->present_cell,
353  this->cell_similarity,
354  this->quadrature,
355  this->get_mapping(),
356  *this->mapping_data,
357  this->mapping_output,
358  *this->fe_data,
359  this->finite_element_output);
360 }
361 
362 
363 
364 template <int dim, int spacedim>
365 std::size_t
367 {
370 }
371 
372 
373 /*------------------------------- FEFaceValuesBase --------------------------*/
374 
375 
376 template <int dim, int spacedim>
378  const unsigned int dofs_per_cell,
379  const UpdateFlags flags,
380  const Mapping<dim, spacedim> &mapping,
382  const Quadrature<dim - 1> &quadrature)
383  : FEFaceValuesBase<dim, spacedim>(dofs_per_cell,
384  flags,
385  mapping,
386  fe,
387  hp::QCollection<dim - 1>(quadrature))
388 {}
389 
390 
391 
392 template <int dim, int spacedim>
394  const unsigned int dofs_per_cell,
395  const UpdateFlags,
396  const Mapping<dim, spacedim> &mapping,
398  const hp::QCollection<dim - 1> &quadrature)
399  : FEValuesBase<dim, spacedim>(quadrature.max_n_quadrature_points(),
400  dofs_per_cell,
402  mapping,
403  fe)
404  , present_face_index(numbers::invalid_unsigned_int)
405  , quadrature(quadrature)
406 {
407  Assert(quadrature.size() == 1 ||
408  quadrature.size() == fe.reference_cell().n_faces(),
409  ExcInternalError());
410 }
411 
412 
413 
414 template <int dim, int spacedim>
415 const std::vector<Tensor<1, spacedim>> &
417 {
418  Assert(this->update_flags & update_boundary_forms,
420  "update_boundary_forms")));
421  return this->mapping_output.boundary_forms;
422 }
423 
424 
425 
426 template <int dim, int spacedim>
427 std::size_t
429 {
432 }
433 
434 
435 /*------------------------------- FEFaceValues -------------------------------*/
436 
437 template <int dim, int spacedim>
438 const unsigned int FEFaceValues<dim, spacedim>::dimension;
439 
440 
441 
442 template <int dim, int spacedim>
444 
445 
446 
447 template <int dim, int spacedim>
449  const Mapping<dim, spacedim> &mapping,
451  const Quadrature<dim - 1> &quadrature,
452  const UpdateFlags update_flags)
453  : FEFaceValues<dim, spacedim>(mapping,
454  fe,
455  hp::QCollection<dim - 1>(quadrature),
456  update_flags)
457 {}
458 
459 
460 
461 template <int dim, int spacedim>
463  const Mapping<dim, spacedim> &mapping,
465  const hp::QCollection<dim - 1> &quadrature,
466  const UpdateFlags update_flags)
467  : FEFaceValuesBase<dim, spacedim>(fe.n_dofs_per_cell(),
468  update_flags,
469  mapping,
470  fe,
471  quadrature)
472 {
474 }
475 
476 
477 
478 template <int dim, int spacedim>
481  const Quadrature<dim - 1> &quadrature,
482  const UpdateFlags update_flags)
483  : FEFaceValues<dim, spacedim>(fe,
484  hp::QCollection<dim - 1>(quadrature),
485  update_flags)
486 {}
487 
488 
489 
490 template <int dim, int spacedim>
493  const hp::QCollection<dim - 1> &quadrature,
494  const UpdateFlags update_flags)
495  : FEFaceValuesBase<dim, spacedim>(
496  fe.n_dofs_per_cell(),
497  update_flags,
498  fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
499  fe,
500  quadrature)
501 {
503 }
504 
505 
506 
507 template <int dim, int spacedim>
508 void
510 {
511  const UpdateFlags flags = this->compute_update_flags(update_flags);
512 
513  // initialize the base classes
514  if (flags & update_mapping)
515  this->mapping_output.initialize(this->max_n_quadrature_points, flags);
516  this->finite_element_output.initialize(this->max_n_quadrature_points,
517  *this->fe,
518  flags);
519 
520  // then get objects into which the FE and the Mapping can store
521  // intermediate data used across calls to reinit. this can be done in parallel
522 
523  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase> (
524  FiniteElement<dim, spacedim>::*finite_element_get_face_data)(
525  const UpdateFlags,
526  const Mapping<dim, spacedim> &,
527  const hp::QCollection<dim - 1> &,
529  spacedim>
531 
532  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> (
533  Mapping<dim, spacedim>::*mapping_get_face_data)(
534  const UpdateFlags, const hp::QCollection<dim - 1> &) const =
536 
537 
539  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>>
540  fe_get_data = Threads::new_task(finite_element_get_face_data,
541  *this->fe,
542  flags,
543  *this->mapping,
544  this->quadrature,
545  this->finite_element_output);
547  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
548  mapping_get_data;
549  if (flags & update_mapping)
550  mapping_get_data = Threads::new_task(mapping_get_face_data,
551  *this->mapping,
552  flags,
553  this->quadrature);
554 
555  this->update_flags = flags;
556 
557  // then collect answers from the two task above
558  this->fe_data = std::move(fe_get_data.return_value());
559  if (flags & update_mapping)
560  this->mapping_data = std::move(mapping_get_data.return_value());
561  else
562  this->mapping_data =
563  std::make_unique<typename Mapping<dim, spacedim>::InternalDataBase>();
564 }
565 
566 
567 
568 template <int dim, int spacedim>
569 template <bool lda>
570 void
573  const unsigned int face_no)
574 {
575  // assert that the finite elements passed to the constructor and
576  // used by the DoFHandler used by this cell, are the same
577  Assert(static_cast<const FiniteElementData<dim> &>(*this->fe) ==
578  static_cast<const FiniteElementData<dim> &>(
579  cell->get_dof_handler().get_fe(cell->active_fe_index())),
581 
583 
584  this->maybe_invalidate_previous_present_cell(cell);
585  this->present_cell = {cell};
586 
587  // this was the part of the work that is dependent on the actual
588  // data type of the iterator. now pass on to the function doing
589  // the real work.
590  do_reinit(face_no);
591 }
592 
593 
594 
595 template <int dim, int spacedim>
596 template <bool lda>
597 void
600  const typename Triangulation<dim, spacedim>::face_iterator &face)
601 {
602  const auto face_n = cell->face_iterator_to_index(face);
603  reinit(cell, face_n);
604 }
605 
606 
607 
608 template <int dim, int spacedim>
609 void
611  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
612  const unsigned int face_no)
613 {
615 
616  this->maybe_invalidate_previous_present_cell(cell);
617  this->present_cell = {cell};
618 
619  // this was the part of the work that is dependent on the actual
620  // data type of the iterator. now pass on to the function doing
621  // the real work.
622  do_reinit(face_no);
623 }
624 
625 
626 
627 template <int dim, int spacedim>
628 void
630  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
631  const typename Triangulation<dim, spacedim>::face_iterator &face)
632 {
633  const auto face_n = cell->face_iterator_to_index(face);
634  reinit(cell, face_n);
635 }
636 
637 
638 
639 template <int dim, int spacedim>
640 void
641 FEFaceValues<dim, spacedim>::do_reinit(const unsigned int face_no)
642 {
643  this->present_face_no = face_no;
644 
645  // first of all, set the present_face_index (if available)
646  const typename Triangulation<dim, spacedim>::cell_iterator cell =
647  this->present_cell;
648  this->present_face_index = cell->face_index(face_no);
649 
650  if (this->update_flags & update_mapping)
651  {
652  this->get_mapping().fill_fe_face_values(this->present_cell,
653  face_no,
654  this->quadrature,
655  *this->mapping_data,
656  this->mapping_output);
657  }
658 
659  this->get_fe().fill_fe_face_values(this->present_cell,
660  face_no,
661  this->quadrature,
662  this->get_mapping(),
663  *this->mapping_data,
664  this->mapping_output,
665  *this->fe_data,
666  this->finite_element_output);
667 
668  const_cast<unsigned int &>(this->n_quadrature_points) =
669  this->quadrature[this->quadrature.size() == 1 ? 0 : face_no].size();
670 }
671 
672 
673 /* ---------------------------- FESubFaceValues ---------------------------- */
674 
675 
676 template <int dim, int spacedim>
678 
679 
680 
681 template <int dim, int spacedim>
683 
684 
685 
686 template <int dim, int spacedim>
688  const Mapping<dim, spacedim> &mapping,
690  const Quadrature<dim - 1> &quadrature,
691  const UpdateFlags update_flags)
692  : FEFaceValuesBase<dim, spacedim>(fe.n_dofs_per_cell(),
693  update_flags,
694  mapping,
695  fe,
696  quadrature)
697 {
699 }
700 
701 
702 
703 template <int dim, int spacedim>
705  const Mapping<dim, spacedim> &mapping,
707  const hp::QCollection<dim - 1> &quadrature,
708  const UpdateFlags update_flags)
709  : FESubfaceValues(mapping, fe, quadrature[0], update_flags)
710 {
712 }
713 
714 
715 
716 template <int dim, int spacedim>
719  const Quadrature<dim - 1> &quadrature,
720  const UpdateFlags update_flags)
721  : FEFaceValuesBase<dim, spacedim>(
722  fe.n_dofs_per_cell(),
723  update_flags,
724  fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
725  fe,
726  quadrature)
727 {
729 }
730 
731 
732 
733 template <int dim, int spacedim>
736  const hp::QCollection<dim - 1> &quadrature,
737  const UpdateFlags update_flags)
738  : FESubfaceValues(fe, quadrature[0], update_flags)
739 {
741 }
742 
743 
744 
745 template <int dim, int spacedim>
746 void
748 {
749  const UpdateFlags flags = this->compute_update_flags(update_flags);
750 
751  // initialize the base classes
752  if (flags & update_mapping)
753  this->mapping_output.initialize(this->max_n_quadrature_points, flags);
754  this->finite_element_output.initialize(this->max_n_quadrature_points,
755  *this->fe,
756  flags);
757 
758  // then get objects into which the FE and the Mapping can store
759  // intermediate data used across calls to reinit. this can be done
760  // in parallel
762  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>>
763  fe_get_data =
765  *this->fe,
766  flags,
767  *this->mapping,
768  this->quadrature[0],
769  this->finite_element_output);
771  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
772  mapping_get_data;
773  if (flags & update_mapping)
774  mapping_get_data =
776  *this->mapping,
777  flags,
778  this->quadrature[0]);
779 
780  this->update_flags = flags;
781 
782  // then collect answers from the two task above
783  this->fe_data = std::move(fe_get_data.return_value());
784  if (flags & update_mapping)
785  this->mapping_data = std::move(mapping_get_data.return_value());
786  else
787  this->mapping_data =
788  std::make_unique<typename Mapping<dim, spacedim>::InternalDataBase>();
789 }
790 
791 
792 
793 template <int dim, int spacedim>
794 template <bool lda>
795 void
798  const unsigned int face_no,
799  const unsigned int subface_no)
800 {
801  // assert that the finite elements passed to the constructor and
802  // used by the DoFHandler used by this cell, are the same
803  Assert(static_cast<const FiniteElementData<dim> &>(*this->fe) ==
804  static_cast<const FiniteElementData<dim> &>(
805  cell->get_dof_handler().get_fe(cell->active_fe_index())),
808  // We would like to check for subface_no < cell->face(face_no)->n_children(),
809  // but unfortunately the current function is also called for
810  // faces without children (see tests/fe/mapping.cc). Therefore,
811  // we must use following workaround of two separate assertions
812  Assert(cell->face(face_no)->has_children() ||
814  ExcIndexRange(subface_no,
815  0,
817  Assert(!cell->face(face_no)->has_children() ||
818  subface_no < cell->face(face_no)->n_active_descendants(),
819  ExcIndexRange(subface_no,
820  0,
821  cell->face(face_no)->n_active_descendants()));
822  Assert(cell->has_children() == false,
823  ExcMessage("You can't use subface data for cells that are "
824  "already refined. Iterate over their children "
825  "instead in these cases."));
826 
827  this->maybe_invalidate_previous_present_cell(cell);
828  this->present_cell = {cell};
829 
830  // this was the part of the work that is dependent on the actual
831  // data type of the iterator. now pass on to the function doing
832  // the real work.
833  do_reinit(face_no, subface_no);
834 }
835 
836 
837 
838 template <int dim, int spacedim>
839 template <bool lda>
840 void
843  const typename Triangulation<dim, spacedim>::face_iterator &face,
844  const typename Triangulation<dim, spacedim>::face_iterator &subface)
845 {
846  reinit(cell,
847  cell->face_iterator_to_index(face),
848  face->child_iterator_to_index(subface));
849 }
850 
851 
852 
853 template <int dim, int spacedim>
854 void
856  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
857  const unsigned int face_no,
858  const unsigned int subface_no)
859 {
861  // We would like to check for subface_no < cell->face(face_no)->n_children(),
862  // but unfortunately the current function is also called for
863  // faces without children for periodic faces, which have hanging nodes on
864  // the other side (see include/deal.II/matrix_free/mapping_info.templates.h).
865  AssertIndexRange(subface_no,
866  (cell->has_periodic_neighbor(face_no) ?
867  cell->periodic_neighbor(face_no)
868  ->face(cell->periodic_neighbor_face_no(face_no))
869  ->n_children() :
870  cell->face(face_no)->n_children()));
871 
872  this->maybe_invalidate_previous_present_cell(cell);
873  this->present_cell = {cell};
874 
875  // this was the part of the work that is dependent on the actual
876  // data type of the iterator. now pass on to the function doing
877  // the real work.
878  do_reinit(face_no, subface_no);
879 }
880 
881 
882 
883 template <int dim, int spacedim>
884 void
886  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
887  const typename Triangulation<dim, spacedim>::face_iterator &face,
888  const typename Triangulation<dim, spacedim>::face_iterator &subface)
889 {
890  reinit(cell,
891  cell->face_iterator_to_index(face),
892  face->child_iterator_to_index(subface));
893 }
894 
895 
896 
897 template <int dim, int spacedim>
898 void
899 FESubfaceValues<dim, spacedim>::do_reinit(const unsigned int face_no,
900  const unsigned int subface_no)
901 {
902  this->present_face_no = face_no;
903 
904  // first of all, set the present_face_index (if available)
905  const typename Triangulation<dim, spacedim>::cell_iterator cell =
906  this->present_cell;
907 
908  if (!cell->face(face_no)->has_children())
909  // no subfaces at all, so set present_face_index to this face rather
910  // than any subface
911  this->present_face_index = cell->face_index(face_no);
912  else if (dim != 3)
913  this->present_face_index = cell->face(face_no)->child_index(subface_no);
914  else
915  {
916  // this is the same logic we use in cell->neighbor_child_on_subface(). See
917  // there for an explanation of the different cases
918  unsigned int subface_index = numbers::invalid_unsigned_int;
919  switch (cell->subface_case(face_no))
920  {
924  subface_index = cell->face(face_no)->child_index(subface_no);
925  break;
928  subface_index = cell->face(face_no)
929  ->child(subface_no / 2)
930  ->child_index(subface_no % 2);
931  break;
934  switch (subface_no)
935  {
936  case 0:
937  case 1:
938  subface_index =
939  cell->face(face_no)->child(0)->child_index(subface_no);
940  break;
941  case 2:
942  subface_index = cell->face(face_no)->child_index(1);
943  break;
944  default:
945  Assert(false, ExcInternalError());
946  }
947  break;
950  switch (subface_no)
951  {
952  case 0:
953  subface_index = cell->face(face_no)->child_index(0);
954  break;
955  case 1:
956  case 2:
957  subface_index =
958  cell->face(face_no)->child(1)->child_index(subface_no - 1);
959  break;
960  default:
961  Assert(false, ExcInternalError());
962  }
963  break;
964  default:
965  Assert(false, ExcInternalError());
966  break;
967  }
968  Assert(subface_index != numbers::invalid_unsigned_int,
969  ExcInternalError());
970  this->present_face_index = subface_index;
971  }
972 
973  // now ask the mapping and the finite element to do the actual work
974  if (this->update_flags & update_mapping)
975  {
976  this->get_mapping().fill_fe_subface_values(this->present_cell,
977  face_no,
978  subface_no,
979  this->quadrature[0],
980  *this->mapping_data,
981  this->mapping_output);
982  }
983 
984  this->get_fe().fill_fe_subface_values(this->present_cell,
985  face_no,
986  subface_no,
987  this->quadrature[0],
988  this->get_mapping(),
989  *this->mapping_data,
990  this->mapping_output,
991  *this->fe_data,
992  this->finite_element_output);
993 }
994 
995 
996 /*------------------------- Explicit Instantiations --------------------------*/
997 
998 #include "fe_values.inst"
999 
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)
Definition: fe_values.cc:377
std::size_t memory_consumption() const
Definition: fe_values.cc:428
const std::vector< Tensor< 1, spacedim > > & get_boundary_forms() const
Definition: fe_values.cc:416
const hp::QCollection< dim - 1 > quadrature
Definition: fe_values.h:303
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access >> &cell, const unsigned int face_no)
void initialize(const UpdateFlags update_flags)
Definition: fe_values.cc:509
void do_reinit(const unsigned int face_no)
Definition: fe_values.cc:641
FEFaceValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim - 1 > &quadrature, const UpdateFlags update_flags)
Definition: fe_values.cc:448
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:687
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)
Definition: fe_values.cc:747
void do_reinit(const unsigned int face_no, const unsigned int subface_no)
Definition: fe_values.cc:899
const SmartPointer< const FiniteElement< dim, spacedim >, FEValuesBase< dim, spacedim > > fe
FEValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim > &quadrature, const UpdateFlags update_flags)
Definition: fe_values.cc:161
void do_reinit()
Definition: fe_values.cc:332
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access >> &cell)
void initialize(const UpdateFlags update_flags)
Definition: fe_values.cc:219
std::size_t memory_consumption() const
Definition: fe_values.cc:366
unsigned int n_dofs_per_cell() const
unsigned int n_components() const
const ComponentMask & get_nonzero_components(const unsigned int i) const
virtual std::unique_ptr< InternalDataBase > get_data(const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const Quadrature< dim > &quadrature, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const =0
unsigned int n_nonzero_components(const unsigned int i) 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:1616
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
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:213
T signaling_nan()