Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
fe_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 
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 Quadrature<dim - 1> & quadrature,
80  const UpdateFlags update_flags);
81 
113  template <class CellIteratorType>
114  void
115  reinit(const CellIteratorType & cell,
116  const unsigned int face_no,
117  const unsigned int sub_face_no,
118  const typename identity<CellIteratorType>::type &cell_neighbor,
119  const unsigned int face_no_neighbor,
120  const unsigned int sub_face_no_neighbor);
121 
133  template <class CellIteratorType>
134  void
135  reinit(const CellIteratorType &cell, const unsigned int face_no);
136 
145  get_fe_face_values(const unsigned int cell_index) const;
146 
150  const Quadrature<dim - 1> &
151  get_quadrature() const;
152 
157  get_update_flags() const;
158 
170  bool
171  at_boundary() const;
172 
184  double
185  JxW(const unsigned int quadrature_point) const;
186 
192  const std::vector<double> &
193  get_JxW_values() const;
194 
204  const std::vector<Tensor<1, spacedim>> &
205  get_normal_vectors() const;
206 
212  const std::vector<Point<spacedim>> &
213  get_quadrature_points() const;
214 
215 
226  unsigned
227  n_current_interface_dofs() const;
228 
239  std::vector<types::global_dof_index>
241 
253  std::array<unsigned int, 2>
254  interface_dof_to_dof_indices(const unsigned int interface_dof_index) const;
255 
265  normal(const unsigned int q_point_index) const;
266 
297  double
298  shape_value(const bool here_or_there,
299  const unsigned int interface_dof_index,
300  const unsigned int q_point,
301  const unsigned int component = 0) const;
302 
317  double
318  jump(const unsigned int interface_dof_index,
319  const unsigned int q_point,
320  const unsigned int component = 0) const;
321 
331  double
332  average(const unsigned int interface_dof_index,
333  const unsigned int q_point,
334  const unsigned int component = 0) const;
335 
346  average_gradient(const unsigned int interface_dof_index,
347  const unsigned int q_point,
348  const unsigned int component = 0) const;
349 
361  average_hessian(const unsigned int interface_dof_index,
362  const unsigned int q_point,
363  const unsigned int component = 0) const;
364 
375  jump_gradient(const unsigned int interface_dof_index,
376  const unsigned int q_point,
377  const unsigned int component = 0) const;
378 
390  jump_hessian(const unsigned int interface_dof_index,
391  const unsigned int q_point,
392  const unsigned int component = 0) const;
393 
404  jump_3rd_derivative(const unsigned int interface_dof_index,
405  const unsigned int q_point,
406  const unsigned int component = 0) const;
407 
412 private:
416  std::vector<types::global_dof_index> interface_dof_indices;
417 
423  std::vector<std::array<unsigned int, 2>> dofmap;
424 
429 
434 
439 
444 
450 
457 };
458 
459 
460 
461 #ifndef DOXYGEN
462 
463 /*---------------------- Inline functions ---------------------*/
464 
465 template <int dim, int spacedim>
467  const Mapping<dim, spacedim> & mapping,
469  const Quadrature<dim - 1> & quadrature,
470  const UpdateFlags update_flags)
471  : n_quadrature_points(quadrature.size())
472  , internal_fe_face_values(mapping, fe, quadrature, update_flags)
473  , internal_fe_subface_values(mapping, fe, quadrature, update_flags)
474  , internal_fe_face_values_neighbor(mapping, fe, quadrature, update_flags)
475  , internal_fe_subface_values_neighbor(mapping, fe, quadrature, update_flags)
476  , fe_face_values(nullptr)
477  , fe_face_values_neighbor(nullptr)
478 {}
479 
480 
481 
482 template <int dim, int spacedim>
485  const Quadrature<dim - 1> & quadrature,
486  const UpdateFlags update_flags)
487  : n_quadrature_points(quadrature.size())
489  fe,
490  quadrature,
491  update_flags)
493  fe,
494  quadrature,
495  update_flags)
497  fe,
498  quadrature,
499  update_flags)
501  fe,
502  quadrature,
503  update_flags)
504  , fe_face_values(nullptr)
505  , fe_face_values_neighbor(nullptr)
506 {}
507 
508 
509 
510 template <int dim, int spacedim>
511 template <class CellIteratorType>
512 void
514  const CellIteratorType & cell,
515  const unsigned int face_no,
516  const unsigned int sub_face_no,
517  const typename identity<CellIteratorType>::type &cell_neighbor,
518  const unsigned int face_no_neighbor,
519  const unsigned int sub_face_no_neighbor)
520 {
521  if (sub_face_no == numbers::invalid_unsigned_int)
522  {
523  internal_fe_face_values.reinit(cell, face_no);
525  }
526  else
527  {
528  internal_fe_subface_values.reinit(cell, face_no, sub_face_no);
530  }
531  if (sub_face_no_neighbor == numbers::invalid_unsigned_int)
532  {
533  internal_fe_face_values_neighbor.reinit(cell_neighbor, face_no_neighbor);
535  }
536  else
537  {
539  face_no_neighbor,
540  sub_face_no_neighbor);
542  }
543 
544  // Set up dof mapping and remove duplicates (for continuous elements).
545  {
546  // Get dof indices first:
547  std::vector<types::global_dof_index> v(
548  fe_face_values->get_fe().n_dofs_per_cell());
549  cell->get_active_or_mg_dof_indices(v);
550  std::vector<types::global_dof_index> v2(
551  fe_face_values_neighbor->get_fe().n_dofs_per_cell());
552  cell_neighbor->get_active_or_mg_dof_indices(v2);
553 
554 
555 
556  // Fill a map from the global dof index to the left and right
557  // local index.
558  std::map<types::global_dof_index, std::pair<unsigned int, unsigned int>>
559  tempmap;
560  std::pair<unsigned int, unsigned int> invalid_entry(
562 
563  for (unsigned int i = 0; i < v.size(); ++i)
564  {
565  // If not already existing, add an invalid entry:
566  auto result = tempmap.insert(std::make_pair(v[i], invalid_entry));
567  result.first->second.first = i;
568  }
569 
570  for (unsigned int i = 0; i < v2.size(); ++i)
571  {
572  // If not already existing, add an invalid entry:
573  auto result = tempmap.insert(std::make_pair(v2[i], invalid_entry));
574  result.first->second.second = i;
575  }
576 
577  // Transfer from the map to the sorted std::vectors.
578  dofmap.resize(tempmap.size());
579  interface_dof_indices.resize(tempmap.size());
580  unsigned int idx = 0;
581  for (auto &x : tempmap)
582  {
583  interface_dof_indices[idx] = x.first;
584  dofmap[idx] = {{x.second.first, x.second.second}};
585  ++idx;
586  }
587  }
588 }
589 
590 
591 
592 template <int dim, int spacedim>
593 template <class CellIteratorType>
594 void
595 FEInterfaceValues<dim, spacedim>::reinit(const CellIteratorType &cell,
596  const unsigned int face_no)
597 {
598  internal_fe_face_values.reinit(cell, face_no);
600  fe_face_values_neighbor = nullptr;
601 
602  interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell());
603  cell->get_active_or_mg_dof_indices(interface_dof_indices);
604 
605 
606  dofmap.resize(interface_dof_indices.size());
607 
608  for (unsigned int i = 0; i < interface_dof_indices.size(); ++i)
609  {
611  }
612 }
613 
614 
615 
616 template <int dim, int spacedim>
617 inline double
618 FEInterfaceValues<dim, spacedim>::JxW(const unsigned int q) const
619 {
620  Assert(fe_face_values != nullptr,
621  ExcMessage("This call requires a call to reinit() first."));
622  return fe_face_values->JxW(q);
623 }
624 
625 
626 
627 template <int dim, int spacedim>
628 const std::vector<double> &
630 {
631  Assert(fe_face_values != nullptr,
632  ExcMessage("This call requires a call to reinit() first."));
633  return fe_face_values->get_JxW_values();
634 }
635 
636 
637 
638 template <int dim, int spacedim>
639 const std::vector<Tensor<1, spacedim>> &
641 {
642  Assert(fe_face_values != nullptr,
643  ExcMessage("This call requires a call to reinit() first."));
645 }
646 
647 
648 
649 template <int dim, int spacedim>
650 const Quadrature<dim - 1> &
652 {
654 }
655 
656 
657 
658 template <int dim, int spacedim>
659 const std::vector<Point<spacedim>> &
661 {
662  Assert(fe_face_values != nullptr,
663  ExcMessage("This call requires a call to reinit() first."));
665 }
666 
667 
668 
669 template <int dim, int spacedim>
672 {
674 }
675 
676 
677 
678 template <int dim, int spacedim>
679 unsigned
681 {
682  Assert(
683  interface_dof_indices.size() > 0,
684  ExcMessage(
685  "n_current_interface_dofs() is only available after a call to reinit()."));
686  return interface_dof_indices.size();
687 }
688 
689 
690 
691 template <int dim, int spacedim>
692 bool
694 {
695  return fe_face_values_neighbor == nullptr;
696 }
697 
698 
699 
700 template <int dim, int spacedim>
701 std::vector<types::global_dof_index>
703 {
704  return interface_dof_indices;
705 }
706 
707 
708 
709 template <int dim, int spacedim>
710 std::array<unsigned int, 2>
712  const unsigned int interface_dof_index) const
713 {
714  AssertIndexRange(interface_dof_index, n_current_interface_dofs());
715  return dofmap[interface_dof_index];
716 }
717 
718 
719 
720 template <int dim, int spacedim>
723  const unsigned int cell_index) const
724 {
725  AssertIndexRange(cell_index, 2);
726  Assert(
727  cell_index == 0 || !at_boundary(),
728  ExcMessage(
729  "You are on a boundary, so you can only ask for the first FEFaceValues object."));
730 
731  return (cell_index == 0) ? *fe_face_values : *fe_face_values_neighbor;
732 }
733 
734 
735 
736 template <int dim, int spacedim>
738 FEInterfaceValues<dim, spacedim>::normal(const unsigned int q_point_index) const
739 {
740  return fe_face_values->normal_vector(q_point_index);
741 }
742 
743 
744 
745 template <int dim, int spacedim>
746 double
748  const bool here_or_there,
749  const unsigned int interface_dof_index,
750  const unsigned int q_point,
751  const unsigned int component) const
752 {
753  const auto dof_pair = dofmap[interface_dof_index];
754 
755  if (here_or_there && dof_pair[0] != numbers::invalid_unsigned_int)
756  return get_fe_face_values(0).shape_value_component(dof_pair[0],
757  q_point,
758  component);
759  if (!here_or_there && dof_pair[1] != numbers::invalid_unsigned_int)
760  return get_fe_face_values(1).shape_value_component(dof_pair[1],
761  q_point,
762  component);
763 
764  return 0.0;
765 }
766 
767 
768 
769 template <int dim, int spacedim>
770 double
771 FEInterfaceValues<dim, spacedim>::jump(const unsigned int interface_dof_index,
772  const unsigned int q_point,
773  const unsigned int component) const
774 {
775  const auto dof_pair = dofmap[interface_dof_index];
776 
777  double value = 0.0;
778 
779  if (dof_pair[0] != numbers::invalid_unsigned_int)
780  value += get_fe_face_values(0).shape_value_component(dof_pair[0],
781  q_point,
782  component);
783  if (dof_pair[1] != numbers::invalid_unsigned_int)
784  value -= get_fe_face_values(1).shape_value_component(dof_pair[1],
785  q_point,
786  component);
787  return value;
788 }
789 
790 
791 
792 template <int dim, int spacedim>
793 double
795  const unsigned int interface_dof_index,
796  const unsigned int q_point,
797  const unsigned int component) const
798 {
799  const auto dof_pair = dofmap[interface_dof_index];
800 
801  if (at_boundary())
802  return 1.0 * get_fe_face_values(0).shape_value_component(dof_pair[0],
803  q_point,
804  component);
805 
806  double value = 0.0;
807 
808  if (dof_pair[0] != numbers::invalid_unsigned_int)
809  value += 0.5 * get_fe_face_values(0).shape_value_component(dof_pair[0],
810  q_point,
811  component);
812  if (dof_pair[1] != numbers::invalid_unsigned_int)
813  value += 0.5 * get_fe_face_values(1).shape_value_component(dof_pair[1],
814  q_point,
815  component);
816 
817  return value;
818 }
819 
820 
821 
822 template <int dim, int spacedim>
825  const unsigned int interface_dof_index,
826  const unsigned int q_point,
827  const unsigned int component) const
828 {
829  const auto dof_pair = dofmap[interface_dof_index];
830 
831  if (at_boundary())
832  return get_fe_face_values(0).shape_grad_component(dof_pair[0],
833  q_point,
834  component);
835 
837 
838  if (dof_pair[0] != numbers::invalid_unsigned_int)
839  value += 0.5 * get_fe_face_values(0).shape_grad_component(dof_pair[0],
840  q_point,
841  component);
842  if (dof_pair[1] != numbers::invalid_unsigned_int)
843  value += 0.5 * get_fe_face_values(1).shape_grad_component(dof_pair[1],
844  q_point,
845  component);
846 
847  return value;
848 }
849 
850 
851 
852 template <int dim, int spacedim>
855  const unsigned int interface_dof_index,
856  const unsigned int q_point,
857  const unsigned int component) const
858 {
859  const auto dof_pair = dofmap[interface_dof_index];
860 
861  if (at_boundary())
862  return get_fe_face_values(0).shape_hessian_component(dof_pair[0],
863  q_point,
864  component);
865 
867 
868  if (dof_pair[0] != numbers::invalid_unsigned_int)
869  value += 0.5 * get_fe_face_values(0).shape_hessian_component(dof_pair[0],
870  q_point,
871  component);
872  if (dof_pair[1] != numbers::invalid_unsigned_int)
873  value += 0.5 * get_fe_face_values(1).shape_hessian_component(dof_pair[1],
874  q_point,
875  component);
876 
877  return value;
878 }
879 
880 
881 
882 template <int dim, int spacedim>
885  const unsigned int interface_dof_index,
886  const unsigned int q_point,
887  const unsigned int component) const
888 {
889  const auto dof_pair = dofmap[interface_dof_index];
890 
891  if (at_boundary())
892  return get_fe_face_values(0).shape_grad_component(dof_pair[0],
893  q_point,
894  component);
895 
897 
898  if (dof_pair[0] != numbers::invalid_unsigned_int)
899  value += 1.0 * get_fe_face_values(0).shape_grad_component(dof_pair[0],
900  q_point,
901  component);
902  if (dof_pair[1] != numbers::invalid_unsigned_int)
903  value += -1.0 * get_fe_face_values(1).shape_grad_component(dof_pair[1],
904  q_point,
905  component);
906 
907  return value;
908 }
909 
910 
911 
912 template <int dim, int spacedim>
915  const unsigned int interface_dof_index,
916  const unsigned int q_point,
917  const unsigned int component) const
918 {
919  const auto dof_pair = dofmap[interface_dof_index];
920 
921  if (at_boundary())
922  return get_fe_face_values(0).shape_hessian_component(dof_pair[0],
923  q_point,
924  component);
925 
927 
928  if (dof_pair[0] != numbers::invalid_unsigned_int)
929  value += 1.0 * get_fe_face_values(0).shape_hessian_component(dof_pair[0],
930  q_point,
931  component);
932  if (dof_pair[1] != numbers::invalid_unsigned_int)
933  value += -1.0 * get_fe_face_values(1).shape_hessian_component(dof_pair[1],
934  q_point,
935  component);
936 
937  return value;
938 }
939 
940 
941 template <int dim, int spacedim>
944  const unsigned int interface_dof_index,
945  const unsigned int q_point,
946  const unsigned int component) const
947 {
948  const auto dof_pair = dofmap[interface_dof_index];
949 
950  if (at_boundary())
951  return get_fe_face_values(0).shape_3rd_derivative_component(dof_pair[0],
952  q_point,
953  component);
954 
956 
957  if (dof_pair[0] != numbers::invalid_unsigned_int)
958  value +=
959  1.0 * get_fe_face_values(0).shape_3rd_derivative_component(dof_pair[0],
960  q_point,
961  component);
962  if (dof_pair[1] != numbers::invalid_unsigned_int)
963  value +=
964  -1.0 * get_fe_face_values(1).shape_3rd_derivative_component(dof_pair[1],
965  q_point,
966  component);
967 
968  return value;
969 }
970 
971 
972 #endif // DOXYGEN
973 
975 
976 #endif
void reinit(const TriaIterator< DoFCellAccessor< DoFHandlerType< dim, spacedim >, level_dof_access >> &cell, const unsigned int face_no, const unsigned int subface_no)
Definition: fe_values.cc:4941
static const unsigned int invalid_unsigned_int
Definition: types.h:191
Tensor< 2, dim > average_hessian(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component=0) const
std::vector< std::array< unsigned int, 2 > > dofmap
FESubfaceValues< dim > internal_fe_subface_values
double average(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component=0) const
FEFaceValues< dim > internal_fe_face_values_neighbor
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
Tensor< 1, dim > average_gradient(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component=0) const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
FEFaceValuesBase< dim > * fe_face_values_neighbor
Tensor< 3, dim > jump_3rd_derivative(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component=0) const
const FiniteElement< dim, spacedim > & get_fe() const
Tensor< 1, spacedim > normal(const unsigned int q_point_index) const
UpdateFlags get_update_flags() const
std::vector< types::global_dof_index > interface_dof_indices
const Quadrature< dim - 1 > & get_quadrature() const
const std::vector< Point< spacedim > > & get_quadrature_points() const
FEFaceValues< dim > internal_fe_face_values
static ::ExceptionBase & ExcMessage(std::string arg1)
FESubfaceValues< dim > internal_fe_subface_values_neighbor
#define Assert(cond, exc)
Definition: exceptions.h:1419
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:302
FEFaceValuesBase< dim > * fe_face_values
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
const std::vector< Tensor< 1, spacedim > > & get_normal_vectors() const
Tensor< 1, dim > jump_gradient(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component=0) 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
std::vector< types::global_dof_index > get_interface_dof_indices() const
const unsigned int n_quadrature_points
const std::vector< double > & get_JxW_values() const
double JxW(const unsigned int quadrature_point) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
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
void reinit(const TriaIterator< DoFCellAccessor< DoFHandlerType< dim, spacedim >, level_dof_access >> &cell, const unsigned int face_no)
Definition: fe_values.cc:4740
const std::vector< Tensor< 1, spacedim > > & get_normal_vectors() const
Definition: fe_values.cc:4199
static const bool value
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
Tensor< 2, dim > jump_hessian(const unsigned int interface_dof_index, const unsigned int q_point, const unsigned int component=0) const
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const