Reference documentation for deal.II version Git 1dc1051882 2021-04-22 23:57:03 +0200
\(\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_interface_values.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 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 #ifndef dealii_fe_interface_values_h
17 #define dealii_fe_interface_values_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <deal.II/fe/fe_values.h>
22 #include <deal.II/fe/mapping_q1.h>
23 
25 
27 
52 template <int dim, int spacedim = dim>
54 {
55 public:
59  const unsigned int n_quadrature_points;
60 
61 
69  const Quadrature<dim - 1> & quadrature,
70  const UpdateFlags update_flags);
71 
79  const hp::QCollection<dim - 1> & quadrature,
80  const UpdateFlags update_flags);
81 
89  const Quadrature<dim - 1> & quadrature,
90  const UpdateFlags update_flags);
91 
123  template <class CellIteratorType>
124  void
125  reinit(const CellIteratorType & cell,
126  const unsigned int face_no,
127  const unsigned int sub_face_no,
128  const typename identity<CellIteratorType>::type &cell_neighbor,
129  const unsigned int face_no_neighbor,
130  const unsigned int sub_face_no_neighbor);
131 
143  template <class CellIteratorType>
144  void
145  reinit(const CellIteratorType &cell, const unsigned int face_no);
146 
155  get_fe_face_values(const unsigned int cell_index) const;
156 
160  const Mapping<dim, spacedim> &
161  get_mapping() const;
162 
167  get_fe() const;
168 
172  const Quadrature<dim - 1> &
173  get_quadrature() const;
174 
179  get_update_flags() const;
180 
192  bool
193  at_boundary() const;
194 
206  double
207  JxW(const unsigned int quadrature_point) const;
208 
214  const std::vector<double> &
215  get_JxW_values() const;
216 
226  const std::vector<Tensor<1, spacedim>> &
227  get_normal_vectors() const;
228 
234  const std::vector<Point<spacedim>> &
235  get_quadrature_points() const;
236 
237 
248  unsigned
249  n_current_interface_dofs() const;
250 
261  std::vector<types::global_dof_index>
263 
275  std::array<unsigned int, 2>
276  interface_dof_to_dof_indices(const unsigned int interface_dof_index) const;
277 
287  normal(const unsigned int q_point_index) const;
288 
319  double
320  shape_value(const bool here_or_there,
321  const unsigned int interface_dof_index,
322  const unsigned int q_point,
323  const unsigned int component = 0) const;
324 
339  double
340  jump(const unsigned int interface_dof_index,
341  const unsigned int q_point,
342  const unsigned int component = 0) const;
343 
353  double
354  average(const unsigned int interface_dof_index,
355  const unsigned int q_point,
356  const unsigned int component = 0) const;
357 
368  average_gradient(const unsigned int interface_dof_index,
369  const unsigned int q_point,
370  const unsigned int component = 0) const;
371 
383  average_hessian(const unsigned int interface_dof_index,
384  const unsigned int q_point,
385  const unsigned int component = 0) const;
386 
397  jump_gradient(const unsigned int interface_dof_index,
398  const unsigned int q_point,
399  const unsigned int component = 0) const;
400 
412  jump_hessian(const unsigned int interface_dof_index,
413  const unsigned int q_point,
414  const unsigned int component = 0) const;
415 
426  jump_3rd_derivative(const unsigned int interface_dof_index,
427  const unsigned int q_point,
428  const unsigned int component = 0) const;
429 
434 private:
438  std::vector<types::global_dof_index> interface_dof_indices;
439 
445  std::vector<std::array<unsigned int, 2>> dofmap;
446 
451 
456 
461 
466 
472 
479 };
480 
481 
482 
483 #ifndef DOXYGEN
484 
485 /*---------------------- Inline functions ---------------------*/
486 
487 template <int dim, int spacedim>
489  const Mapping<dim, spacedim> & mapping,
491  const Quadrature<dim - 1> & quadrature,
492  const UpdateFlags update_flags)
493  : n_quadrature_points(quadrature.size())
494  , internal_fe_face_values(mapping, fe, quadrature, update_flags)
495  , internal_fe_subface_values(mapping, fe, quadrature, update_flags)
496  , internal_fe_face_values_neighbor(mapping, fe, quadrature, update_flags)
497  , internal_fe_subface_values_neighbor(mapping, fe, quadrature, update_flags)
498  , fe_face_values(nullptr)
499  , fe_face_values_neighbor(nullptr)
500 {}
501 
502 template <int dim, int spacedim>
504  const Mapping<dim, spacedim> & mapping,
506  const hp::QCollection<dim - 1> & quadrature,
507  const UpdateFlags update_flags)
509  , internal_fe_face_values(mapping, fe, quadrature, update_flags)
510  , internal_fe_subface_values(mapping, fe, quadrature, update_flags)
511  , internal_fe_face_values_neighbor(mapping, fe, quadrature[0], update_flags)
513  fe,
514  quadrature[0],
515  update_flags)
516  , fe_face_values(nullptr)
517  , fe_face_values_neighbor(nullptr)
518 {}
519 
520 
521 
522 template <int dim, int spacedim>
525  const Quadrature<dim - 1> & quadrature,
526  const UpdateFlags update_flags)
527  : n_quadrature_points(quadrature.size())
529  fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
530  fe,
531  quadrature,
532  update_flags)
534  fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
535  fe,
536  quadrature,
537  update_flags)
539  fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
540  fe,
541  quadrature,
542  update_flags)
544  fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
545  fe,
546  quadrature,
547  update_flags)
548  , fe_face_values(nullptr)
549  , fe_face_values_neighbor(nullptr)
550 {}
551 
552 
553 
554 template <int dim, int spacedim>
555 template <class CellIteratorType>
556 void
558  const CellIteratorType & cell,
559  const unsigned int face_no,
560  const unsigned int sub_face_no,
561  const typename identity<CellIteratorType>::type &cell_neighbor,
562  const unsigned int face_no_neighbor,
563  const unsigned int sub_face_no_neighbor)
564 {
565  if (sub_face_no == numbers::invalid_unsigned_int)
566  {
567  internal_fe_face_values.reinit(cell, face_no);
569  }
570  else
571  {
572  internal_fe_subface_values.reinit(cell, face_no, sub_face_no);
574  }
575  if (sub_face_no_neighbor == numbers::invalid_unsigned_int)
576  {
577  internal_fe_face_values_neighbor.reinit(cell_neighbor, face_no_neighbor);
579  }
580  else
581  {
582  internal_fe_subface_values_neighbor.reinit(cell_neighbor,
583  face_no_neighbor,
584  sub_face_no_neighbor);
586  }
587 
588  AssertDimension(fe_face_values->n_quadrature_points,
589  fe_face_values_neighbor->n_quadrature_points);
590 
591  const_cast<unsigned int &>(this->n_quadrature_points) =
592  fe_face_values->n_quadrature_points;
593 
594  // Set up dof mapping and remove duplicates (for continuous elements).
595  {
596  // Get dof indices first:
597  std::vector<types::global_dof_index> v(
598  fe_face_values->get_fe().n_dofs_per_cell());
599  cell->get_active_or_mg_dof_indices(v);
600  std::vector<types::global_dof_index> v2(
601  fe_face_values_neighbor->get_fe().n_dofs_per_cell());
602  cell_neighbor->get_active_or_mg_dof_indices(v2);
603 
604 
605 
606  // Fill a map from the global dof index to the left and right
607  // local index.
608  std::map<types::global_dof_index, std::pair<unsigned int, unsigned int>>
609  tempmap;
610  std::pair<unsigned int, unsigned int> invalid_entry(
612 
613  for (unsigned int i = 0; i < v.size(); ++i)
614  {
615  // If not already existing, add an invalid entry:
616  auto result = tempmap.insert(std::make_pair(v[i], invalid_entry));
617  result.first->second.first = i;
618  }
619 
620  for (unsigned int i = 0; i < v2.size(); ++i)
621  {
622  // If not already existing, add an invalid entry:
623  auto result = tempmap.insert(std::make_pair(v2[i], invalid_entry));
624  result.first->second.second = i;
625  }
626 
627  // Transfer from the map to the sorted std::vectors.
628  dofmap.resize(tempmap.size());
629  interface_dof_indices.resize(tempmap.size());
630  unsigned int idx = 0;
631  for (auto &x : tempmap)
632  {
633  interface_dof_indices[idx] = x.first;
634  dofmap[idx] = {{x.second.first, x.second.second}};
635  ++idx;
636  }
637  }
638 }
639 
640 
641 
642 template <int dim, int spacedim>
643 template <class CellIteratorType>
644 void
645 FEInterfaceValues<dim, spacedim>::reinit(const CellIteratorType &cell,
646  const unsigned int face_no)
647 {
648  internal_fe_face_values.reinit(cell, face_no);
650  fe_face_values_neighbor = nullptr;
651 
652  interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell());
653  cell->get_active_or_mg_dof_indices(interface_dof_indices);
654 
655 
656  dofmap.resize(interface_dof_indices.size());
657 
658  for (unsigned int i = 0; i < interface_dof_indices.size(); ++i)
659  {
661  }
662 }
663 
664 
665 
666 template <int dim, int spacedim>
667 inline double
668 FEInterfaceValues<dim, spacedim>::JxW(const unsigned int q) const
669 {
670  Assert(fe_face_values != nullptr,
671  ExcMessage("This call requires a call to reinit() first."));
672  return fe_face_values->JxW(q);
673 }
674 
675 
676 
677 template <int dim, int spacedim>
678 const std::vector<double> &
680 {
681  Assert(fe_face_values != nullptr,
682  ExcMessage("This call requires a call to reinit() first."));
683  return fe_face_values->get_JxW_values();
684 }
685 
686 
687 
688 template <int dim, int spacedim>
689 const std::vector<Tensor<1, spacedim>> &
691 {
692  Assert(fe_face_values != nullptr,
693  ExcMessage("This call requires a call to reinit() first."));
694  return fe_face_values->get_normal_vectors();
695 }
696 
697 
698 
699 template <int dim, int spacedim>
702 {
703  return internal_fe_face_values.get_mapping();
704 }
705 
706 
707 
708 template <int dim, int spacedim>
711 {
712  return internal_fe_face_values.get_fe();
713 }
714 
715 
716 
717 template <int dim, int spacedim>
718 const Quadrature<dim - 1> &
720 {
721  return internal_fe_face_values.get_quadrature();
722 }
723 
724 
725 
726 template <int dim, int spacedim>
727 const std::vector<Point<spacedim>> &
729 {
730  Assert(fe_face_values != nullptr,
731  ExcMessage("This call requires a call to reinit() first."));
732  return fe_face_values->get_quadrature_points();
733 }
734 
735 
736 
737 template <int dim, int spacedim>
740 {
741  return internal_fe_face_values.get_update_flags();
742 }
743 
744 
745 
746 template <int dim, int spacedim>
747 unsigned
749 {
750  Assert(
751  interface_dof_indices.size() > 0,
752  ExcMessage(
753  "n_current_interface_dofs() is only available after a call to reinit()."));
754  return interface_dof_indices.size();
755 }
756 
757 
758 
759 template <int dim, int spacedim>
760 bool
762 {
763  return fe_face_values_neighbor == nullptr;
764 }
765 
766 
767 
768 template <int dim, int spacedim>
769 std::vector<types::global_dof_index>
771 {
772  return interface_dof_indices;
773 }
774 
775 
776 
777 template <int dim, int spacedim>
778 std::array<unsigned int, 2>
780  const unsigned int interface_dof_index) const
781 {
782  AssertIndexRange(interface_dof_index, n_current_interface_dofs());
783  return dofmap[interface_dof_index];
784 }
785 
786 
787 
788 template <int dim, int spacedim>
791  const unsigned int cell_index) const
792 {
793  AssertIndexRange(cell_index, 2);
794  Assert(
795  cell_index == 0 || !at_boundary(),
796  ExcMessage(
797  "You are on a boundary, so you can only ask for the first FEFaceValues object."));
798 
799  return (cell_index == 0) ? *fe_face_values : *fe_face_values_neighbor;
800 }
801 
802 
803 
804 template <int dim, int spacedim>
806 FEInterfaceValues<dim, spacedim>::normal(const unsigned int q_point_index) const
807 {
808  return fe_face_values->normal_vector(q_point_index);
809 }
810 
811 
812 
813 template <int dim, int spacedim>
814 double
816  const bool here_or_there,
817  const unsigned int interface_dof_index,
818  const unsigned int q_point,
819  const unsigned int component) const
820 {
821  const auto dof_pair = dofmap[interface_dof_index];
822 
823  if (here_or_there && dof_pair[0] != numbers::invalid_unsigned_int)
824  return get_fe_face_values(0).shape_value_component(dof_pair[0],
825  q_point,
826  component);
827  if (!here_or_there && dof_pair[1] != numbers::invalid_unsigned_int)
828  return get_fe_face_values(1).shape_value_component(dof_pair[1],
829  q_point,
830  component);
831 
832  return 0.0;
833 }
834 
835 
836 
837 template <int dim, int spacedim>
838 double
839 FEInterfaceValues<dim, spacedim>::jump(const unsigned int interface_dof_index,
840  const unsigned int q_point,
841  const unsigned int component) const
842 {
843  const auto dof_pair = dofmap[interface_dof_index];
844 
845  double value = 0.0;
846 
847  if (dof_pair[0] != numbers::invalid_unsigned_int)
848  value += get_fe_face_values(0).shape_value_component(dof_pair[0],
849  q_point,
850  component);
851  if (dof_pair[1] != numbers::invalid_unsigned_int)
852  value -= get_fe_face_values(1).shape_value_component(dof_pair[1],
853  q_point,
854  component);
855  return value;
856 }
857 
858 
859 
860 template <int dim, int spacedim>
861 double
863  const unsigned int interface_dof_index,
864  const unsigned int q_point,
865  const unsigned int component) const
866 {
867  const auto dof_pair = dofmap[interface_dof_index];
868 
869  if (at_boundary())
870  return 1.0 * get_fe_face_values(0).shape_value_component(dof_pair[0],
871  q_point,
872  component);
873 
874  double value = 0.0;
875 
876  if (dof_pair[0] != numbers::invalid_unsigned_int)
877  value += 0.5 * get_fe_face_values(0).shape_value_component(dof_pair[0],
878  q_point,
879  component);
880  if (dof_pair[1] != numbers::invalid_unsigned_int)
881  value += 0.5 * get_fe_face_values(1).shape_value_component(dof_pair[1],
882  q_point,
883  component);
884 
885  return value;
886 }
887 
888 
889 
890 template <int dim, int spacedim>
893  const unsigned int interface_dof_index,
894  const unsigned int q_point,
895  const unsigned int component) const
896 {
897  const auto dof_pair = dofmap[interface_dof_index];
898 
899  if (at_boundary())
900  return get_fe_face_values(0).shape_grad_component(dof_pair[0],
901  q_point,
902  component);
903 
904  Tensor<1, spacedim> value;
905 
906  if (dof_pair[0] != numbers::invalid_unsigned_int)
907  value += 0.5 * get_fe_face_values(0).shape_grad_component(dof_pair[0],
908  q_point,
909  component);
910  if (dof_pair[1] != numbers::invalid_unsigned_int)
911  value += 0.5 * get_fe_face_values(1).shape_grad_component(dof_pair[1],
912  q_point,
913  component);
914 
915  return value;
916 }
917 
918 
919 
920 template <int dim, int spacedim>
923  const unsigned int interface_dof_index,
924  const unsigned int q_point,
925  const unsigned int component) const
926 {
927  const auto dof_pair = dofmap[interface_dof_index];
928 
929  if (at_boundary())
930  return get_fe_face_values(0).shape_hessian_component(dof_pair[0],
931  q_point,
932  component);
933 
934  Tensor<2, spacedim> value;
935 
936  if (dof_pair[0] != numbers::invalid_unsigned_int)
937  value += 0.5 * get_fe_face_values(0).shape_hessian_component(dof_pair[0],
938  q_point,
939  component);
940  if (dof_pair[1] != numbers::invalid_unsigned_int)
941  value += 0.5 * get_fe_face_values(1).shape_hessian_component(dof_pair[1],
942  q_point,
943  component);
944 
945  return value;
946 }
947 
948 
949 
950 template <int dim, int spacedim>
953  const unsigned int interface_dof_index,
954  const unsigned int q_point,
955  const unsigned int component) const
956 {
957  const auto dof_pair = dofmap[interface_dof_index];
958 
959  if (at_boundary())
960  return get_fe_face_values(0).shape_grad_component(dof_pair[0],
961  q_point,
962  component);
963 
964  Tensor<1, spacedim> value;
965 
966  if (dof_pair[0] != numbers::invalid_unsigned_int)
967  value += 1.0 * get_fe_face_values(0).shape_grad_component(dof_pair[0],
968  q_point,
969  component);
970  if (dof_pair[1] != numbers::invalid_unsigned_int)
971  value += -1.0 * get_fe_face_values(1).shape_grad_component(dof_pair[1],
972  q_point,
973  component);
974 
975  return value;
976 }
977 
978 
979 
980 template <int dim, int spacedim>
983  const unsigned int interface_dof_index,
984  const unsigned int q_point,
985  const unsigned int component) const
986 {
987  const auto dof_pair = dofmap[interface_dof_index];
988 
989  if (at_boundary())
990  return get_fe_face_values(0).shape_hessian_component(dof_pair[0],
991  q_point,
992  component);
993 
994  Tensor<2, spacedim> value;
995 
996  if (dof_pair[0] != numbers::invalid_unsigned_int)
997  value += 1.0 * get_fe_face_values(0).shape_hessian_component(dof_pair[0],
998  q_point,
999  component);
1000  if (dof_pair[1] != numbers::invalid_unsigned_int)
1001  value += -1.0 * get_fe_face_values(1).shape_hessian_component(dof_pair[1],
1002  q_point,
1003  component);
1004 
1005  return value;
1006 }
1007 
1008 
1009 template <int dim, int spacedim>
1012  const unsigned int interface_dof_index,
1013  const unsigned int q_point,
1014  const unsigned int component) const
1015 {
1016  const auto dof_pair = dofmap[interface_dof_index];
1017 
1018  if (at_boundary())
1019  return get_fe_face_values(0).shape_3rd_derivative_component(dof_pair[0],
1020  q_point,
1021  component);
1022 
1023  Tensor<3, spacedim> value;
1024 
1025  if (dof_pair[0] != numbers::invalid_unsigned_int)
1026  value +=
1027  1.0 * get_fe_face_values(0).shape_3rd_derivative_component(dof_pair[0],
1028  q_point,
1029  component);
1030  if (dof_pair[1] != numbers::invalid_unsigned_int)
1031  value +=
1032  -1.0 * get_fe_face_values(1).shape_3rd_derivative_component(dof_pair[1],
1033  q_point,
1034  component);
1035 
1036  return value;
1037 }
1038 
1039 
1040 #endif // DOXYGEN
1041 
1043 
1044 #endif
Tensor< 2, spacedim > jump_hessian(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component=0) const
static const unsigned int invalid_unsigned_int
Definition: types.h:196
std::vector< std::array< unsigned int, 2 > > dofmap
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1630
double average(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component=0) const
double JxW(const unsigned int quadrature_point) const
void reinit(const CellIteratorType &cell, const unsigned int face_no, const unsigned int sub_face_no, const typename identity< CellIteratorType >::type &cell_neighbor, const unsigned int face_no_neighbor, const unsigned int sub_face_no_neighbor)
std::array< unsigned int, 2 > interface_dof_to_dof_indices(const unsigned int interface_dof_index) const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1698
Tensor< 1, spacedim > average_gradient(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component=0) const
FEFaceValues< dim, spacedim > internal_fe_face_values
Tensor< 1, spacedim > normal(const unsigned int q_point_index) const
Tensor< 1, spacedim > jump_gradient(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component=0) const
FEFaceValues< dim, spacedim > internal_fe_face_values_neighbor
Tensor< 2, spacedim > average_hessian(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component=0) const
std::vector< types::global_dof_index > interface_dof_indices
const Mapping< dim, spacedim > & get_mapping() const
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int max_n_quadrature_points() const
Definition: q_collection.h:181
#define Assert(cond, exc)
Definition: exceptions.h:1473
double shape_value(const bool here_or_there, const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component=0) const
UpdateFlags
unsigned n_current_interface_dofs() const
Abstract base class for mapping classes.
Definition: mapping.h:303
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
FEFaceValuesBase< dim, spacedim > * fe_face_values_neighbor
const std::vector< Tensor< 1, spacedim > > & get_normal_vectors() const
const Quadrature< dim - 1 > & get_quadrature() const
FEInterfaceValues(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const Quadrature< dim - 1 > &quadrature, const UpdateFlags update_flags)
unsigned int size() const
bool at_boundary() const
FESubfaceValues< dim, spacedim > internal_fe_subface_values
std::vector< types::global_dof_index > get_interface_dof_indices() const
const unsigned int n_quadrature_points
Tensor< 3, spacedim > jump_3rd_derivative(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component=0) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
UpdateFlags get_update_flags() const
const FEFaceValuesBase< dim, spacedim > & get_fe_face_values(const unsigned int cell_index) const
const std::vector< double > & get_JxW_values() const
const FiniteElement< dim, spacedim > & get_fe() const
const std::vector< Point< spacedim > > & get_quadrature_points() const
double jump(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component=0) const
FESubfaceValues< dim, spacedim > internal_fe_subface_values_neighbor
FEFaceValuesBase< dim, spacedim > * fe_face_values
ReferenceCell reference_cell() const