Reference documentation for deal.II version Git adf7f9fdfe 2021-01-28 13:54:50 +0100
\(\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())
530  .template get_default_linear_mapping<dim, spacedim>(),
531  fe,
532  quadrature,
533  update_flags)
536  .template get_default_linear_mapping<dim, spacedim>(),
537  fe,
538  quadrature,
539  update_flags)
542  .template get_default_linear_mapping<dim, spacedim>(),
543  fe,
544  quadrature,
545  update_flags)
548  .template get_default_linear_mapping<dim, spacedim>(),
549  fe,
550  quadrature,
551  update_flags)
552  , fe_face_values(nullptr)
553  , fe_face_values_neighbor(nullptr)
554 {}
555 
556 
557 
558 template <int dim, int spacedim>
559 template <class CellIteratorType>
560 void
562  const CellIteratorType & cell,
563  const unsigned int face_no,
564  const unsigned int sub_face_no,
565  const typename identity<CellIteratorType>::type &cell_neighbor,
566  const unsigned int face_no_neighbor,
567  const unsigned int sub_face_no_neighbor)
568 {
569  if (sub_face_no == numbers::invalid_unsigned_int)
570  {
571  internal_fe_face_values.reinit(cell, face_no);
573  }
574  else
575  {
576  internal_fe_subface_values.reinit(cell, face_no, sub_face_no);
578  }
579  if (sub_face_no_neighbor == numbers::invalid_unsigned_int)
580  {
581  internal_fe_face_values_neighbor.reinit(cell_neighbor, face_no_neighbor);
583  }
584  else
585  {
586  internal_fe_subface_values_neighbor.reinit(cell_neighbor,
587  face_no_neighbor,
588  sub_face_no_neighbor);
590  }
591 
592  AssertDimension(fe_face_values->n_quadrature_points,
593  fe_face_values_neighbor->n_quadrature_points);
594 
595  const_cast<unsigned int &>(this->n_quadrature_points) =
596  fe_face_values->n_quadrature_points;
597 
598  // Set up dof mapping and remove duplicates (for continuous elements).
599  {
600  // Get dof indices first:
601  std::vector<types::global_dof_index> v(
602  fe_face_values->get_fe().n_dofs_per_cell());
603  cell->get_active_or_mg_dof_indices(v);
604  std::vector<types::global_dof_index> v2(
605  fe_face_values_neighbor->get_fe().n_dofs_per_cell());
606  cell_neighbor->get_active_or_mg_dof_indices(v2);
607 
608 
609 
610  // Fill a map from the global dof index to the left and right
611  // local index.
612  std::map<types::global_dof_index, std::pair<unsigned int, unsigned int>>
613  tempmap;
614  std::pair<unsigned int, unsigned int> invalid_entry(
616 
617  for (unsigned int i = 0; i < v.size(); ++i)
618  {
619  // If not already existing, add an invalid entry:
620  auto result = tempmap.insert(std::make_pair(v[i], invalid_entry));
621  result.first->second.first = i;
622  }
623 
624  for (unsigned int i = 0; i < v2.size(); ++i)
625  {
626  // If not already existing, add an invalid entry:
627  auto result = tempmap.insert(std::make_pair(v2[i], invalid_entry));
628  result.first->second.second = i;
629  }
630 
631  // Transfer from the map to the sorted std::vectors.
632  dofmap.resize(tempmap.size());
633  interface_dof_indices.resize(tempmap.size());
634  unsigned int idx = 0;
635  for (auto &x : tempmap)
636  {
637  interface_dof_indices[idx] = x.first;
638  dofmap[idx] = {{x.second.first, x.second.second}};
639  ++idx;
640  }
641  }
642 }
643 
644 
645 
646 template <int dim, int spacedim>
647 template <class CellIteratorType>
648 void
649 FEInterfaceValues<dim, spacedim>::reinit(const CellIteratorType &cell,
650  const unsigned int face_no)
651 {
652  internal_fe_face_values.reinit(cell, face_no);
654  fe_face_values_neighbor = nullptr;
655 
656  interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell());
657  cell->get_active_or_mg_dof_indices(interface_dof_indices);
658 
659 
660  dofmap.resize(interface_dof_indices.size());
661 
662  for (unsigned int i = 0; i < interface_dof_indices.size(); ++i)
663  {
665  }
666 }
667 
668 
669 
670 template <int dim, int spacedim>
671 inline double
672 FEInterfaceValues<dim, spacedim>::JxW(const unsigned int q) const
673 {
674  Assert(fe_face_values != nullptr,
675  ExcMessage("This call requires a call to reinit() first."));
676  return fe_face_values->JxW(q);
677 }
678 
679 
680 
681 template <int dim, int spacedim>
682 const std::vector<double> &
684 {
685  Assert(fe_face_values != nullptr,
686  ExcMessage("This call requires a call to reinit() first."));
687  return fe_face_values->get_JxW_values();
688 }
689 
690 
691 
692 template <int dim, int spacedim>
693 const std::vector<Tensor<1, spacedim>> &
695 {
696  Assert(fe_face_values != nullptr,
697  ExcMessage("This call requires a call to reinit() first."));
698  return fe_face_values->get_normal_vectors();
699 }
700 
701 
702 
703 template <int dim, int spacedim>
706 {
707  return internal_fe_face_values.get_mapping();
708 }
709 
710 
711 
712 template <int dim, int spacedim>
715 {
716  return internal_fe_face_values.get_fe();
717 }
718 
719 
720 
721 template <int dim, int spacedim>
722 const Quadrature<dim - 1> &
724 {
725  return internal_fe_face_values.get_quadrature();
726 }
727 
728 
729 
730 template <int dim, int spacedim>
731 const std::vector<Point<spacedim>> &
733 {
734  Assert(fe_face_values != nullptr,
735  ExcMessage("This call requires a call to reinit() first."));
736  return fe_face_values->get_quadrature_points();
737 }
738 
739 
740 
741 template <int dim, int spacedim>
744 {
745  return internal_fe_face_values.get_update_flags();
746 }
747 
748 
749 
750 template <int dim, int spacedim>
751 unsigned
753 {
754  Assert(
755  interface_dof_indices.size() > 0,
756  ExcMessage(
757  "n_current_interface_dofs() is only available after a call to reinit()."));
758  return interface_dof_indices.size();
759 }
760 
761 
762 
763 template <int dim, int spacedim>
764 bool
766 {
767  return fe_face_values_neighbor == nullptr;
768 }
769 
770 
771 
772 template <int dim, int spacedim>
773 std::vector<types::global_dof_index>
775 {
776  return interface_dof_indices;
777 }
778 
779 
780 
781 template <int dim, int spacedim>
782 std::array<unsigned int, 2>
784  const unsigned int interface_dof_index) const
785 {
786  AssertIndexRange(interface_dof_index, n_current_interface_dofs());
787  return dofmap[interface_dof_index];
788 }
789 
790 
791 
792 template <int dim, int spacedim>
795  const unsigned int cell_index) const
796 {
797  AssertIndexRange(cell_index, 2);
798  Assert(
799  cell_index == 0 || !at_boundary(),
800  ExcMessage(
801  "You are on a boundary, so you can only ask for the first FEFaceValues object."));
802 
803  return (cell_index == 0) ? *fe_face_values : *fe_face_values_neighbor;
804 }
805 
806 
807 
808 template <int dim, int spacedim>
810 FEInterfaceValues<dim, spacedim>::normal(const unsigned int q_point_index) const
811 {
812  return fe_face_values->normal_vector(q_point_index);
813 }
814 
815 
816 
817 template <int dim, int spacedim>
818 double
820  const bool here_or_there,
821  const unsigned int interface_dof_index,
822  const unsigned int q_point,
823  const unsigned int component) const
824 {
825  const auto dof_pair = dofmap[interface_dof_index];
826 
827  if (here_or_there && dof_pair[0] != numbers::invalid_unsigned_int)
828  return get_fe_face_values(0).shape_value_component(dof_pair[0],
829  q_point,
830  component);
831  if (!here_or_there && dof_pair[1] != numbers::invalid_unsigned_int)
832  return get_fe_face_values(1).shape_value_component(dof_pair[1],
833  q_point,
834  component);
835 
836  return 0.0;
837 }
838 
839 
840 
841 template <int dim, int spacedim>
842 double
843 FEInterfaceValues<dim, spacedim>::jump(const unsigned int interface_dof_index,
844  const unsigned int q_point,
845  const unsigned int component) const
846 {
847  const auto dof_pair = dofmap[interface_dof_index];
848 
849  double value = 0.0;
850 
851  if (dof_pair[0] != numbers::invalid_unsigned_int)
852  value += get_fe_face_values(0).shape_value_component(dof_pair[0],
853  q_point,
854  component);
855  if (dof_pair[1] != numbers::invalid_unsigned_int)
856  value -= get_fe_face_values(1).shape_value_component(dof_pair[1],
857  q_point,
858  component);
859  return value;
860 }
861 
862 
863 
864 template <int dim, int spacedim>
865 double
867  const unsigned int interface_dof_index,
868  const unsigned int q_point,
869  const unsigned int component) const
870 {
871  const auto dof_pair = dofmap[interface_dof_index];
872 
873  if (at_boundary())
874  return 1.0 * get_fe_face_values(0).shape_value_component(dof_pair[0],
875  q_point,
876  component);
877 
878  double value = 0.0;
879 
880  if (dof_pair[0] != numbers::invalid_unsigned_int)
881  value += 0.5 * get_fe_face_values(0).shape_value_component(dof_pair[0],
882  q_point,
883  component);
884  if (dof_pair[1] != numbers::invalid_unsigned_int)
885  value += 0.5 * get_fe_face_values(1).shape_value_component(dof_pair[1],
886  q_point,
887  component);
888 
889  return value;
890 }
891 
892 
893 
894 template <int dim, int spacedim>
897  const unsigned int interface_dof_index,
898  const unsigned int q_point,
899  const unsigned int component) const
900 {
901  const auto dof_pair = dofmap[interface_dof_index];
902 
903  if (at_boundary())
904  return get_fe_face_values(0).shape_grad_component(dof_pair[0],
905  q_point,
906  component);
907 
908  Tensor<1, spacedim> value;
909 
910  if (dof_pair[0] != numbers::invalid_unsigned_int)
911  value += 0.5 * get_fe_face_values(0).shape_grad_component(dof_pair[0],
912  q_point,
913  component);
914  if (dof_pair[1] != numbers::invalid_unsigned_int)
915  value += 0.5 * get_fe_face_values(1).shape_grad_component(dof_pair[1],
916  q_point,
917  component);
918 
919  return value;
920 }
921 
922 
923 
924 template <int dim, int spacedim>
927  const unsigned int interface_dof_index,
928  const unsigned int q_point,
929  const unsigned int component) const
930 {
931  const auto dof_pair = dofmap[interface_dof_index];
932 
933  if (at_boundary())
934  return get_fe_face_values(0).shape_hessian_component(dof_pair[0],
935  q_point,
936  component);
937 
938  Tensor<2, spacedim> value;
939 
940  if (dof_pair[0] != numbers::invalid_unsigned_int)
941  value += 0.5 * get_fe_face_values(0).shape_hessian_component(dof_pair[0],
942  q_point,
943  component);
944  if (dof_pair[1] != numbers::invalid_unsigned_int)
945  value += 0.5 * get_fe_face_values(1).shape_hessian_component(dof_pair[1],
946  q_point,
947  component);
948 
949  return value;
950 }
951 
952 
953 
954 template <int dim, int spacedim>
957  const unsigned int interface_dof_index,
958  const unsigned int q_point,
959  const unsigned int component) const
960 {
961  const auto dof_pair = dofmap[interface_dof_index];
962 
963  if (at_boundary())
964  return get_fe_face_values(0).shape_grad_component(dof_pair[0],
965  q_point,
966  component);
967 
968  Tensor<1, spacedim> value;
969 
970  if (dof_pair[0] != numbers::invalid_unsigned_int)
971  value += 1.0 * get_fe_face_values(0).shape_grad_component(dof_pair[0],
972  q_point,
973  component);
974  if (dof_pair[1] != numbers::invalid_unsigned_int)
975  value += -1.0 * get_fe_face_values(1).shape_grad_component(dof_pair[1],
976  q_point,
977  component);
978 
979  return value;
980 }
981 
982 
983 
984 template <int dim, int spacedim>
987  const unsigned int interface_dof_index,
988  const unsigned int q_point,
989  const unsigned int component) const
990 {
991  const auto dof_pair = dofmap[interface_dof_index];
992 
993  if (at_boundary())
994  return get_fe_face_values(0).shape_hessian_component(dof_pair[0],
995  q_point,
996  component);
997 
998  Tensor<2, spacedim> value;
999 
1000  if (dof_pair[0] != numbers::invalid_unsigned_int)
1001  value += 1.0 * get_fe_face_values(0).shape_hessian_component(dof_pair[0],
1002  q_point,
1003  component);
1004  if (dof_pair[1] != numbers::invalid_unsigned_int)
1005  value += -1.0 * get_fe_face_values(1).shape_hessian_component(dof_pair[1],
1006  q_point,
1007  component);
1008 
1009  return value;
1010 }
1011 
1012 
1013 template <int dim, int spacedim>
1016  const unsigned int interface_dof_index,
1017  const unsigned int q_point,
1018  const unsigned int component) const
1019 {
1020  const auto dof_pair = dofmap[interface_dof_index];
1021 
1022  if (at_boundary())
1023  return get_fe_face_values(0).shape_3rd_derivative_component(dof_pair[0],
1024  q_point,
1025  component);
1026 
1027  Tensor<3, spacedim> value;
1028 
1029  if (dof_pair[0] != numbers::invalid_unsigned_int)
1030  value +=
1031  1.0 * get_fe_face_values(0).shape_3rd_derivative_component(dof_pair[0],
1032  q_point,
1033  component);
1034  if (dof_pair[1] != numbers::invalid_unsigned_int)
1035  value +=
1036  -1.0 * get_fe_face_values(1).shape_3rd_derivative_component(dof_pair[1],
1037  q_point,
1038  component);
1039 
1040  return value;
1041 }
1042 
1043 
1044 #endif // DOXYGEN
1045 
1047 
1048 #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:1623
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:1691
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
ReferenceCell::Type reference_cell_type() const
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int max_n_quadrature_points() const
Definition: q_collection.h:209
#define Assert(cond, exc)
Definition: exceptions.h:1466
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:380
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:379
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