Reference documentation for deal.II version GIT b6bf1e606d 2022-08-11 15:25: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_tools.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2022 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_tools_H
17 #define dealii_fe_tools_H
18 
19 
20 
21 #include <deal.II/base/config.h>
22 
27 #include <deal.II/base/tensor.h>
28 
30 
32 
34 
35 #include <memory>
36 #include <string>
37 #include <vector>
38 
39 
41 
42 // Forward declarations
43 #ifndef DOXYGEN
44 template <typename number>
45 class FullMatrix;
46 template <int dim>
47 class Quadrature;
48 template <int dim, int spacedim>
49 class FiniteElement;
50 template <int dim, int spacedim>
51 class DoFHandler;
52 template <int dim>
53 class FiniteElementData;
54 template <typename number>
55 class AffineConstraints;
56 #endif
57 
58 
77 namespace FETools
78 {
87  template <int dim, int spacedim = dim>
88  class FEFactoryBase : public Subscriptor
89  {
90  public:
94  virtual std::unique_ptr<FiniteElement<dim, spacedim>>
95  get(const unsigned int degree) const = 0;
96 
102  virtual std::unique_ptr<FiniteElement<dim, spacedim>>
103  get(const Quadrature<1> &quad) const = 0;
104 
108  virtual ~FEFactoryBase() override = default;
109  };
110 
120  template <class FE>
121  class FEFactory : public FEFactoryBase<FE::dimension, FE::space_dimension>
122  {
123  public:
127  virtual std::unique_ptr<FiniteElement<FE::dimension, FE::space_dimension>>
128  get(const unsigned int degree) const override;
129 
134  virtual std::unique_ptr<FiniteElement<FE::dimension, FE::space_dimension>>
135  get(const Quadrature<1> &quad) const override;
136  };
137 
153  template <int dim, int spacedim>
154  void
156  std::vector<unsigned int> & renumbering,
157  std::vector<std::vector<unsigned int>> &start_indices);
158 
174  template <int dim, int spacedim>
175  void
177  std::vector<types::global_dof_index> &renumbering,
178  std::vector<types::global_dof_index> &block_data,
179  bool return_start_indices = true);
180 
194  template <int dim, typename number, int spacedim>
195  void
197  const FiniteElement<dim, spacedim> &fe2,
198  FullMatrix<number> &interpolation_matrix);
199 
211  template <int dim, typename number, int spacedim>
212  void
214  const FiniteElement<dim, spacedim> &fe2,
215  FullMatrix<number> &interpolation_matrix);
216 
228  template <int dim, typename number, int spacedim>
229  void
231  const FiniteElement<dim, spacedim> &fe2,
232  FullMatrix<number> &difference_matrix);
233 
237  template <int dim, typename number, int spacedim>
238  void
240  const FiniteElement<dim, spacedim> &fe2,
242 
311  template <int dim, int spacedim>
314 
353  template <int dim, typename number, int spacedim>
354  void
356  const FiniteElement<dim, spacedim> & fe,
357  std::vector<std::vector<FullMatrix<number>>> &matrices,
358  const bool isotropic_only = false,
359  const double threshold = 1.e-12);
360 
382  template <int dim, typename number, int spacedim>
383  void
387  const unsigned int face_coarse,
388  const unsigned int face_fine,
389  const double threshold = 1.e-12);
390 
422  template <int dim, typename number, int spacedim>
423  void
425  const FiniteElement<dim, spacedim> & fe,
426  std::vector<std::vector<FullMatrix<number>>> &matrices,
427  const bool isotropic_only = false);
428 
514  template <int dim, int spacedim>
515  void
518  const Quadrature<dim> & lhs_quadrature,
519  const Quadrature<dim> & rhs_quadrature,
520  FullMatrix<double> & X);
521 
529  template <int dim, int spacedim>
530  void
533  const Quadrature<dim> & quadrature,
534  FullMatrix<double> & I_q);
535 
549  template <int dim>
550  void
552  const FullMatrix<double> & projection_matrix,
553  const std::vector<Tensor<1, dim>> &vector_of_tensors_at_qp,
554  std::vector<Tensor<1, dim>> & vector_of_tensors_at_nodes);
555 
556 
557 
561  template <int dim>
562  void
564  const FullMatrix<double> & projection_matrix,
565  const std::vector<SymmetricTensor<2, dim>> &vector_of_tensors_at_qp,
566  std::vector<SymmetricTensor<2, dim>> & vector_of_tensors_at_nodes);
567 
568 
569 
579  template <int dim, int spacedim>
580  void
583  const Quadrature<dim - 1> & lhs_quadrature,
584  const Quadrature<dim - 1> & rhs_quadrature,
586  const unsigned int face,
587  FullMatrix<double> & X);
588 
589 
590 
606  template <int dim, int spacedim, typename number>
607  void
609  const FiniteElement<dim, spacedim> &finite_element,
610  const std::vector<Vector<number>> & support_point_values,
611  std::vector<number> & dof_values);
612 
613 
614 
649  template <int dim, int spacedim, class InVector, class OutVector>
650  void
652  const InVector & u1,
653  const DoFHandler<dim, spacedim> &dof2,
654  OutVector & u2);
655 
672  template <int dim, int spacedim, class InVector, class OutVector>
673  void
675  const DoFHandler<dim, spacedim> & dof1,
676  const InVector & u1,
677  const DoFHandler<dim, spacedim> & dof2,
679  OutVector & u2);
680 
694  template <int dim, class InVector, class OutVector, int spacedim>
695  void
697  const InVector & u1,
698  const FiniteElement<dim, spacedim> &fe2,
699  OutVector & u1_interpolated);
700 
713  template <int dim, class InVector, class OutVector, int spacedim>
714  void
716  const DoFHandler<dim, spacedim> & dof1,
718  const InVector & u1,
719  const DoFHandler<dim, spacedim> & dof2,
721  OutVector & u1_interpolated);
722 
732  template <int dim, class InVector, class OutVector, int spacedim>
733  void
735  const InVector & z1,
736  const FiniteElement<dim, spacedim> &fe2,
737  OutVector & z1_difference);
738 
751  template <int dim, class InVector, class OutVector, int spacedim>
752  void
754  const DoFHandler<dim, spacedim> & dof1,
756  const InVector & z1,
757  const DoFHandler<dim, spacedim> & dof2,
759  OutVector & z1_difference);
760 
761 
762 
771  template <int dim, class InVector, class OutVector, int spacedim>
772  void
774  const InVector & u1,
775  const DoFHandler<dim, spacedim> &dof2,
776  OutVector & u2);
777 
838  template <int dim, class InVector, class OutVector, int spacedim>
839  void
841  const InVector & z1,
842  const DoFHandler<dim, spacedim> &dof2,
843  OutVector & z2);
844 
857  template <int dim, class InVector, class OutVector, int spacedim>
858  void
860  const DoFHandler<dim, spacedim> & dof1,
861  const InVector & z1,
862  const DoFHandler<dim, spacedim> & dof2,
864  OutVector & z2);
865 
881  template <int dim>
882  std::vector<unsigned int>
884 
891  template <int dim>
892  std::vector<unsigned int>
894 
970  namespace Compositing
971  {
988  template <int dim, int spacedim>
991  const std::vector<const FiniteElement<dim, spacedim> *> &fes,
992  const std::vector<unsigned int> & multiplicities,
993  const bool do_tensor_product = true);
994 
1000  template <int dim, int spacedim>
1003  const std::initializer_list<
1004  std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>>
1005  &fe_systems);
1006 
1010  template <int dim, int spacedim>
1013  const unsigned int N1,
1014  const FiniteElement<dim, spacedim> *fe2 = nullptr,
1015  const unsigned int N2 = 0,
1016  const FiniteElement<dim, spacedim> *fe3 = nullptr,
1017  const unsigned int N3 = 0,
1018  const FiniteElement<dim, spacedim> *fe4 = nullptr,
1019  const unsigned int N4 = 0,
1020  const FiniteElement<dim, spacedim> *fe5 = nullptr,
1021  const unsigned int N5 = 0);
1022 
1035  template <int dim, int spacedim>
1036  std::vector<bool>
1038  const std::vector<const FiniteElement<dim, spacedim> *> &fes,
1039  const std::vector<unsigned int> & multiplicities);
1040 
1046  template <int dim, int spacedim>
1047  std::vector<bool>
1049  const std::initializer_list<
1050  std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>>
1051  &fe_systems);
1052 
1067  template <int dim, int spacedim>
1068  std::vector<bool>
1070  const FiniteElement<dim, spacedim> *fe1,
1071  const unsigned int N1,
1072  const FiniteElement<dim, spacedim> *fe2 = nullptr,
1073  const unsigned int N2 = 0,
1074  const FiniteElement<dim, spacedim> *fe3 = nullptr,
1075  const unsigned int N3 = 0,
1076  const FiniteElement<dim, spacedim> *fe4 = nullptr,
1077  const unsigned int N4 = 0,
1078  const FiniteElement<dim, spacedim> *fe5 = nullptr,
1079  const unsigned int N5 = 0);
1080 
1081 
1098  template <int dim, int spacedim>
1099  std::vector<ComponentMask>
1101  const std::vector<const FiniteElement<dim, spacedim> *> &fes,
1102  const std::vector<unsigned int> & multiplicities,
1103  const bool do_tensor_product = true);
1104 
1110  template <int dim, int spacedim>
1111  std::vector<ComponentMask>
1113  const std::initializer_list<
1114  std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>>
1115  &fe_systems);
1116 
1135  template <int dim, int spacedim>
1136  std::vector<ComponentMask>
1138  const FiniteElement<dim, spacedim> *fe1,
1139  const unsigned int N1,
1140  const FiniteElement<dim, spacedim> *fe2 = nullptr,
1141  const unsigned int N2 = 0,
1142  const FiniteElement<dim, spacedim> *fe3 = nullptr,
1143  const unsigned int N3 = 0,
1144  const FiniteElement<dim, spacedim> *fe4 = nullptr,
1145  const unsigned int N4 = 0,
1146  const FiniteElement<dim, spacedim> *fe5 = nullptr,
1147  const unsigned int N5 = 0,
1148  const bool do_tensor_product = true);
1149 
1166  template <int dim, int spacedim>
1167  void
1169  std::vector<std::pair<std::pair<unsigned int, unsigned int>,
1170  unsigned int>> &system_to_base_table,
1171  std::vector<std::pair<unsigned int, unsigned int>>
1172  & system_to_component_table,
1173  std::vector<std::pair<std::pair<unsigned int, unsigned int>,
1174  unsigned int>> &component_to_base_table,
1175  const FiniteElement<dim, spacedim> & finite_element,
1176  const bool do_tensor_product = true);
1177 
1193  template <int dim, int spacedim>
1194  void
1196  std::vector<std::pair<std::pair<unsigned int, unsigned int>,
1197  unsigned int>> &face_system_to_base_table,
1198  std::vector<std::pair<unsigned int, unsigned int>>
1199  & face_system_to_component_table,
1200  const FiniteElement<dim, spacedim> &finite_element,
1201  const bool do_tensor_product = true,
1202  const unsigned int face_no = 0 /*TODO*/);
1203 
1204  } // namespace Compositing
1205 
1206 
1240  template <int dim, int spacedim = dim>
1241  std::unique_ptr<FiniteElement<dim, spacedim>>
1242  get_fe_by_name(const std::string &name);
1243 
1286  template <int dim, int spacedim>
1287  void
1288  add_fe_name(const std::string & name,
1289  const FEFactoryBase<dim, spacedim> *factory);
1290 
1301  std::string,
1302  << "Can't re-generate a finite element from the string '"
1303  << arg1 << "'.");
1304 
1317  char,
1318  int,
1319  << "The dimension " << arg1
1320  << " in the finite element string must match "
1321  << "the space dimension " << arg2 << '.');
1322 
1329 
1343 
1351  "You are using continuous elements on a grid with "
1352  "hanging nodes but without providing hanging node "
1353  "constraints. Use the respective function with "
1354  "additional AffineConstraints argument(s), instead.");
1367  int,
1368  int,
1369  int,
1370  int,
1371  << "This is a " << arg1 << 'x' << arg2 << " matrix, "
1372  << "but should be a " << arg3 << 'x' << arg4 << " matrix.");
1373 
1380  double,
1381  << "Least squares fit leaves a gap of " << arg1);
1382 
1389  int,
1390  int,
1391  << arg1 << " must be greater than " << arg2);
1392 } // namespace FETools
1393 
1394 
1395 #ifndef DOXYGEN
1396 
1397 namespace FETools
1398 {
1399  template <class FE>
1400  std::unique_ptr<FiniteElement<FE::dimension, FE::space_dimension>>
1401  FEFactory<FE>::get(const unsigned int degree) const
1402  {
1403  return std::make_unique<FE>(degree);
1404  }
1405 
1406  namespace Compositing
1407  {
1408  template <int dim, int spacedim>
1409  std::vector<bool>
1411  const std::initializer_list<
1412  std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>>
1413  &fe_systems)
1414  {
1415  std::vector<const FiniteElement<dim, spacedim> *> fes;
1416  std::vector<unsigned int> multiplicities;
1417 
1418  const auto extract =
1419  [&fes, &multiplicities](
1420  const std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>,
1421  unsigned int> &fe_system) {
1422  fes.push_back(fe_system.first.get());
1423  multiplicities.push_back(fe_system.second);
1424  };
1425 
1426  for (const auto &p : fe_systems)
1427  extract(p);
1428 
1429  return compute_restriction_is_additive_flags(fes, multiplicities);
1430  }
1431 
1432 
1433 
1434  template <int dim, int spacedim>
1437  const std::initializer_list<
1438  std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>>
1439  &fe_systems)
1440  {
1441  std::vector<const FiniteElement<dim, spacedim> *> fes;
1442  std::vector<unsigned int> multiplicities;
1443 
1444  const auto extract =
1445  [&fes, &multiplicities](
1446  const std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>,
1447  unsigned int> &fe_system) {
1448  fes.push_back(fe_system.first.get());
1449  multiplicities.push_back(fe_system.second);
1450  };
1451 
1452  for (const auto &p : fe_systems)
1453  extract(p);
1454 
1455  return multiply_dof_numbers(fes, multiplicities, true);
1456  }
1457 
1458 
1459 
1460  template <int dim, int spacedim>
1461  std::vector<ComponentMask>
1463  const std::initializer_list<
1464  std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>>
1465  &fe_systems)
1466  {
1467  std::vector<const FiniteElement<dim, spacedim> *> fes;
1468  std::vector<unsigned int> multiplicities;
1469 
1470  const auto extract =
1471  [&fes, &multiplicities](
1472  const std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>,
1473  unsigned int> &fe_system) {
1474  fes.push_back(fe_system.first.get());
1475  multiplicities.push_back(fe_system.second);
1476  };
1477 
1478  for (const auto &p : fe_systems)
1479  extract(p);
1480 
1481  return compute_nonzero_components(fes, multiplicities, true);
1482  }
1483  } // namespace Compositing
1484 } // namespace FETools
1485 
1486 #endif
1487 
1491 
1492 #endif /* dealii_fe_tools_H */
virtual ~FEFactoryBase() override=default
virtual std::unique_ptr< FiniteElement< dim, spacedim > > get(const unsigned int degree) const =0
virtual std::unique_ptr< FiniteElement< dim, spacedim > > get(const Quadrature< 1 > &quad) const =0
virtual std::unique_ptr< FiniteElement< FE::dimension, FE::space_dimension > > get(const unsigned int degree) const override
virtual std::unique_ptr< FiniteElement< FE::dimension, FE::space_dimension > > get(const Quadrature< 1 > &quad) const override
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcTriangulationMismatch()
static ::ExceptionBase & ExcNotGreaterThan(int arg1, int arg2)
static ::ExceptionBase & ExcGridNotRefinedAtLeastOnce()
#define DeclException0(Exception0)
Definition: exceptions.h:464
static ::ExceptionBase & ExcInvalidFEDimension(char arg1, int arg2)
static ::ExceptionBase & ExcInvalidFE()
static ::ExceptionBase & ExcMatrixDimensionMismatch(int arg1, int arg2, int arg3, int arg4)
static ::ExceptionBase & ExcInvalidFEName(std::string arg1)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:578
static ::ExceptionBase & ExcLeastSquaresError(double arg1)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
static ::ExceptionBase & ExcHangingNodesNotAllowed()
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
static ::ExceptionBase & ExcFENotPrimitive()
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
FiniteElementData< dim > multiply_dof_numbers(const std::vector< const FiniteElement< dim, spacedim > * > &fes, const std::vector< unsigned int > &multiplicities, const bool do_tensor_product=true)
std::vector< ComponentMask > compute_nonzero_components(const std::vector< const FiniteElement< dim, spacedim > * > &fes, const std::vector< unsigned int > &multiplicities, const bool do_tensor_product=true)
std::vector< bool > compute_restriction_is_additive_flags(const std::vector< const FiniteElement< dim, spacedim > * > &fes, const std::vector< unsigned int > &multiplicities)
void build_cell_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int >> &system_to_base_table, std::vector< std::pair< unsigned int, unsigned int >> &system_to_component_table, std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int >> &component_to_base_table, const FiniteElement< dim, spacedim > &finite_element, const bool do_tensor_product=true)
void build_face_tables(std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int >> &face_system_to_base_table, std::vector< std::pair< unsigned int, unsigned int >> &face_system_to_component_table, const FiniteElement< dim, spacedim > &finite_element, const bool do_tensor_product=true, const unsigned int face_no=0)
void interpolation_difference(const DoFHandler< dim, spacedim > &dof1, const InVector &z1, const FiniteElement< dim, spacedim > &fe2, OutVector &z1_difference)
std::vector< unsigned int > hierarchic_to_lexicographic_numbering(unsigned int degree)
void project_dg(const DoFHandler< dim, spacedim > &dof1, const InVector &u1, const DoFHandler< dim, spacedim > &dof2, OutVector &u2)
void get_back_interpolation_matrix(const FiniteElement< dim, spacedim > &fe1, const FiniteElement< dim, spacedim > &fe2, FullMatrix< number > &interpolation_matrix)
void compute_component_wise(const FiniteElement< dim, spacedim > &fe, std::vector< unsigned int > &renumbering, std::vector< std::vector< unsigned int >> &start_indices)
void get_projection_matrix(const FiniteElement< dim, spacedim > &fe1, const FiniteElement< dim, spacedim > &fe2, FullMatrix< number > &matrix)
void compute_embedding_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number >>> &matrices, const bool isotropic_only=false, const double threshold=1.e-12)
void compute_projection_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number >>> &matrices, const bool isotropic_only=false)
void add_fe_name(const std::string &name, const FEFactoryBase< dim, spacedim > *factory)
std::vector< unsigned int > lexicographic_to_hierarchic_numbering(unsigned int degree)
void compute_interpolation_to_quadrature_points_matrix(const FiniteElement< dim, spacedim > &fe, const Quadrature< dim > &quadrature, FullMatrix< double > &I_q)
FullMatrix< double > compute_node_matrix(const FiniteElement< dim, spacedim > &fe)
void convert_generalized_support_point_values_to_dof_values(const FiniteElement< dim, spacedim > &finite_element, const std::vector< Vector< number >> &support_point_values, std::vector< number > &dof_values)
void compute_projection_from_quadrature_points_matrix(const FiniteElement< dim, spacedim > &fe, const Quadrature< dim > &lhs_quadrature, const Quadrature< dim > &rhs_quadrature, FullMatrix< double > &X)
void get_interpolation_matrix(const FiniteElement< dim, spacedim > &fe1, const FiniteElement< dim, spacedim > &fe2, FullMatrix< number > &interpolation_matrix)
void extrapolate(const DoFHandler< dim, spacedim > &dof1, const InVector &z1, const DoFHandler< dim, spacedim > &dof2, OutVector &z2)
void compute_face_embedding_matrices(const FiniteElement< dim, spacedim > &fe, FullMatrix< number >(&matrices)[GeometryInfo< dim >::max_children_per_face], const unsigned int face_coarse, const unsigned int face_fine, const double threshold=1.e-12)
void interpolate(const DoFHandler< dim, spacedim > &dof1, const InVector &u1, const DoFHandler< dim, spacedim > &dof2, OutVector &u2)
std::unique_ptr< FiniteElement< dim, spacedim > > get_fe_by_name(const std::string &name)
void compute_block_renumbering(const FiniteElement< dim, spacedim > &fe, std::vector< types::global_dof_index > &renumbering, std::vector< types::global_dof_index > &block_data, bool return_start_indices=true)
void compute_projection_from_face_quadrature_points_matrix(const FiniteElement< dim, spacedim > &fe, const Quadrature< dim - 1 > &lhs_quadrature, const Quadrature< dim - 1 > &rhs_quadrature, const typename DoFHandler< dim, spacedim >::active_cell_iterator &cell, const unsigned int face, FullMatrix< double > &X)
void get_interpolation_difference_matrix(const FiniteElement< dim, spacedim > &fe1, const FiniteElement< dim, spacedim > &fe2, FullMatrix< number > &difference_matrix)
void compute_projection_from_quadrature_points(const FullMatrix< double > &projection_matrix, const std::vector< Tensor< 1, dim >> &vector_of_tensors_at_qp, std::vector< Tensor< 1, dim >> &vector_of_tensors_at_nodes)
void back_interpolate(const DoFHandler< dim, spacedim > &dof1, const InVector &u1, const FiniteElement< dim, spacedim > &fe2, OutVector &u1_interpolated)
@ matrix
Contents is actually a matrix.
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)