Reference documentation for deal.II version GIT 9042b9283b 2023-12-02 14:50: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_abf.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2003 - 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 
16 
20 #include <deal.II/base/table.h>
21 #include <deal.II/base/utilities.h>
22 
24 
25 #include <deal.II/fe/fe.h>
26 #include <deal.II/fe/fe_abf.h>
27 #include <deal.II/fe/fe_tools.h>
28 #include <deal.II/fe/fe_values.h>
29 #include <deal.II/fe/mapping.h>
30 
31 #include <deal.II/grid/tria.h>
33 
34 #include <iostream>
35 #include <memory>
36 #include <sstream>
37 
39 
40 // TODO: implement the adjust_quad_dof_index_for_face_orientation_table and
41 // adjust_line_dof_index_for_line_orientation_table fields, and write tests
42 // similar to bits/face_orientation_and_fe_q_*
43 
44 template <int dim>
45 FE_ABF<dim>::FE_ABF(const unsigned int deg)
46  : FE_PolyTensor<dim>(
47  PolynomialsABF<dim>(deg),
49  dim,
50  deg + 2,
51  FiniteElementData<dim>::Hdiv),
52  std::vector<bool>(PolynomialsABF<dim>::n_polynomials(deg), true),
53  std::vector<ComponentMask>(PolynomialsABF<dim>::n_polynomials(deg),
54  ComponentMask(std::vector<bool>(dim, true))))
55  , rt_order(deg)
56 {
57  Assert(dim >= 2, ExcImpossibleInDim(dim));
58  const unsigned int n_dofs = this->n_dofs_per_cell();
59 
61  // First, initialize the
62  // generalized support points and
63  // quadrature weights, since they
64  // are required for interpolation.
66 
67  // Now compute the inverse node matrix, generating the correct
68  // basis functions from the raw ones. For a discussion of what
69  // exactly happens here, see FETools::compute_node_matrix.
71  this->inverse_node_matrix.reinit(n_dofs, n_dofs);
72  this->inverse_node_matrix.invert(M);
73  // From now on, the shape functions provided by FiniteElement::shape_value
74  // and similar functions will be the correct ones, not
75  // the raw shape functions from the polynomial space anymore.
76 
77  // Reinit the vectors of
78  // restriction and prolongation
79  // matrices to the right sizes.
80  // Restriction only for isotropic
81  // refinement
83  // Fill prolongation matrices with embedding operators
84  FETools::compute_embedding_matrices(*this, this->prolongation, false, 1.e-10);
85 
87 
88  // TODO: the implementation makes the assumption that all faces have the
89  // same number of dofs
90  AssertDimension(this->n_unique_faces(), 1);
91  const unsigned int face_no = 0;
92 
93  // TODO[TL]: for anisotropic refinement we will probably need a table of
94  // submatrices with an array for each refine case
95  std::vector<FullMatrix<double>> face_embeddings(
96  1 << (dim - 1),
97  FullMatrix<double>(this->n_dofs_per_face(face_no),
98  this->n_dofs_per_face(face_no)));
99  // TODO: Something goes wrong there. The error of the least squares fit
100  // is to large ...
101  // FETools::compute_face_embedding_matrices(*this, face_embeddings.data(), 0,
102  // 0);
103  this->interface_constraints.reinit((1 << (dim - 1)) *
104  this->n_dofs_per_face(face_no),
105  this->n_dofs_per_face(face_no));
106  unsigned int target_row = 0;
107  for (const auto &face_embedding : face_embeddings)
108  for (unsigned int i = 0; i < face_embedding.m(); ++i)
109  {
110  for (unsigned int j = 0; j < face_embedding.n(); ++j)
111  this->interface_constraints(target_row, j) = face_embedding(i, j);
112  ++target_row;
113  }
114 
115  // We need to initialize the dof permutation table and the one for the sign
116  // change.
118 }
119 
120 
121 template <int dim>
122 void
124 {
125  // for 1d and 2d, do nothing
126  if (dim < 3)
127  return;
128 
129  // TODO: Implement this for this class
130  return;
131 }
132 
133 
134 template <int dim>
135 std::string
137 {
138  // note that the
139  // FETools::get_fe_by_name
140  // function depends on the
141  // particular format of the string
142  // this function returns, so they
143  // have to be kept in synch
144 
145  std::ostringstream namebuf;
146 
147  namebuf << "FE_ABF<" << dim << ">(" << rt_order << ")";
148 
149  return namebuf.str();
150 }
151 
152 
153 
154 template <int dim>
155 std::unique_ptr<FiniteElement<dim, dim>>
157 {
158  return std::make_unique<FE_ABF<dim>>(rt_order);
159 }
160 
161 
162 //---------------------------------------------------------------------------
163 // Auxiliary and internal functions
164 //---------------------------------------------------------------------------
165 
166 
167 
168 // Version for 2d and higher. See above for 1d version
169 template <int dim>
170 void
172 {
173  QGauss<dim> cell_quadrature(deg + 2);
174  const unsigned int n_interior_points = cell_quadrature.size();
175 
176  // TODO: the implementation makes the assumption that all faces have the
177  // same number of dofs
178  AssertDimension(this->n_unique_faces(), 1);
179  const unsigned int face_no = 0;
180 
181  unsigned int n_face_points = (dim > 1) ? 1 : 0;
182  // compute (deg+1)^(dim-1)
183  for (unsigned int d = 1; d < dim; ++d)
184  n_face_points *= deg + 1;
185 
186  this->generalized_support_points.resize(
187  GeometryInfo<dim>::faces_per_cell * n_face_points + n_interior_points);
188  this->generalized_face_support_points[face_no].resize(n_face_points);
189 
190 
191  // These might be required when the faces contribution is computed
192  // Therefore they will be initialized at this point.
193  std::array<std::unique_ptr<AnisotropicPolynomials<dim>>, dim> polynomials_abf;
194 
195  // Generate x_1^{i} x_2^{r+1} ...
196  for (unsigned int dd = 0; dd < dim; ++dd)
197  {
198  std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
199  for (unsigned int d = 0; d < dim; ++d)
200  poly[d].push_back(Polynomials::Monomial<double>(deg + 1));
202 
203  polynomials_abf[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
204  }
205 
206  // Number of the point being entered
207  unsigned int current = 0;
208 
209  if (dim > 1)
210  {
211  QGauss<dim - 1> face_points(deg + 1);
214 
215  boundary_weights.reinit(n_face_points, legendre.n());
216 
217  // Assert (face_points.size() == this->n_dofs_per_face(),
218  // ExcInternalError());
219 
220  for (unsigned int k = 0; k < n_face_points; ++k)
221  {
222  this->generalized_face_support_points[face_no][k] =
223  face_points.point(k);
224  // Compute its quadrature
225  // contribution for each
226  // moment.
227  for (unsigned int i = 0; i < legendre.n(); ++i)
228  {
229  boundary_weights(k, i) =
230  face_points.weight(k) *
231  legendre.compute_value(i, face_points.point(k));
232  }
233  }
234 
235  Quadrature<dim> faces =
237  face_points);
238  for (; current < GeometryInfo<dim>::faces_per_cell * n_face_points;
239  ++current)
240  {
241  // Enter the support point
242  // into the vector
243  this->generalized_support_points[current] = faces.point(current);
244  }
245 
246 
247  // Now initialize edge interior weights for the ABF elements.
248  // These are completely independent from the usual edge moments. They
249  // stem from applying the Gauss theorem to the nodal values, which
250  // was necessary to cast the ABF elements into the deal.II framework
251  // for vector valued elements.
252  boundary_weights_abf.reinit(faces.size(), polynomials_abf[0]->n() * dim);
253  for (unsigned int k = 0; k < faces.size(); ++k)
254  {
255  for (unsigned int i = 0; i < polynomials_abf[0]->n() * dim; ++i)
256  {
257  boundary_weights_abf(k, i) =
258  polynomials_abf[i % dim]->compute_value(i / dim,
259  faces.point(k)) *
260  faces.weight(k);
261  }
262  }
263  }
264 
265  // Create Legendre basis for the
266  // space D_xi Q_k
267  if (deg > 0)
268  {
269  std::array<std::unique_ptr<AnisotropicPolynomials<dim>>, dim> polynomials;
270 
271  for (unsigned int dd = 0; dd < dim; ++dd)
272  {
273  std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
274  for (unsigned int d = 0; d < dim; ++d)
277 
278  polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
279  }
280 
281  interior_weights.reinit(
282  TableIndices<3>(n_interior_points, polynomials[0]->n(), dim));
283 
284  for (unsigned int k = 0; k < cell_quadrature.size(); ++k)
285  {
286  for (unsigned int i = 0; i < polynomials[0]->n(); ++i)
287  for (unsigned int d = 0; d < dim; ++d)
288  interior_weights(k, i, d) =
289  cell_quadrature.weight(k) *
290  polynomials[d]->compute_value(i, cell_quadrature.point(k));
291  }
292  }
293 
294 
295  // Decouple the creation of the generalized support points
296  // from computation of interior weights.
297  for (unsigned int k = 0; k < cell_quadrature.size(); ++k)
298  this->generalized_support_points[current++] = cell_quadrature.point(k);
299 
300  // Additional functionality for the ABF elements
301  // TODO: Here the canonical extension of the principle
302  // behind the ABF elements is implemented. It is unclear,
303  // if this really leads to the ABF spaces in 3d!
304  interior_weights_abf.reinit(TableIndices<3>(cell_quadrature.size(),
305  polynomials_abf[0]->n() * dim,
306  dim));
307  Tensor<1, dim> poly_grad;
308 
309  for (unsigned int k = 0; k < cell_quadrature.size(); ++k)
310  {
311  for (unsigned int i = 0; i < polynomials_abf[0]->n() * dim; ++i)
312  {
313  poly_grad =
314  polynomials_abf[i % dim]->compute_grad(i / dim,
315  cell_quadrature.point(k)) *
316  cell_quadrature.weight(k);
317  // The minus sign comes from the use of the Gauss theorem to replace
318  // the divergence.
319  for (unsigned int d = 0; d < dim; ++d)
320  interior_weights_abf(k, i, d) = -poly_grad[d];
321  }
322  }
323 
324  Assert(current == this->generalized_support_points.size(),
325  ExcInternalError());
326 }
327 
328 
329 
330 // This function is the same Raviart-Thomas interpolation performed by
331 // interpolate. Still, we cannot use interpolate, since it was written
332 // for smooth functions. The functions interpolated here are not
333 // smooth, maybe even not continuous. Therefore, we must double the
334 // number of quadrature points in each direction in order to integrate
335 // only smooth functions.
336 
337 // Then again, the interpolated function is chosen such that the
338 // moments coincide with the function to be interpolated.
339 
340 template <int dim>
341 void
343 {
344  if (dim == 1)
345  {
346  unsigned int iso = RefinementCase<dim>::isotropic_refinement - 1;
347  for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
348  ++i)
349  this->restriction[iso][i].reinit(0, 0);
350  return;
351  }
352  unsigned int iso = RefinementCase<dim>::isotropic_refinement - 1;
353  QGauss<dim - 1> q_base(rt_order + 1);
354  const unsigned int n_face_points = q_base.size();
355  // First, compute interpolation on
356  // subfaces
357  for (const unsigned int face : GeometryInfo<dim>::face_indices())
358  {
359  // The shape functions of the
360  // child cell are evaluated
361  // in the quadrature points
362  // of a full face.
363  Quadrature<dim> q_face =
364  QProjector<dim>::project_to_face(this->reference_cell(), q_base, face);
365  // Store shape values, since the
366  // evaluation suffers if not
367  // ordered by point
368  Table<2, double> cached_values_face(this->n_dofs_per_cell(),
369  q_face.size());
370  for (unsigned int k = 0; k < q_face.size(); ++k)
371  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
372  cached_values_face(i, k) = this->shape_value_component(
373  i, q_face.point(k), GeometryInfo<dim>::unit_normal_direction[face]);
374 
375  for (unsigned int sub = 0; sub < GeometryInfo<dim>::max_children_per_face;
376  ++sub)
377  {
378  // The weight functions for
379  // the coarse face are
380  // evaluated on the subface
381  // only.
383  this->reference_cell(), q_base, face, sub);
384  const unsigned int child = GeometryInfo<dim>::child_cell_on_face(
386 
387  // On a certain face, we must
388  // compute the moments of ALL
389  // fine level functions with
390  // the coarse level weight
391  // functions belonging to
392  // that face. Due to the
393  // orthogonalization process
394  // when building the shape
395  // functions, these weights
396  // are equal to the
397  // corresponding shape
398  // functions.
399  for (unsigned int k = 0; k < n_face_points; ++k)
400  for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
401  ++i_child)
402  for (unsigned int i_face = 0;
403  i_face < this->n_dofs_per_face(face);
404  ++i_face)
405  {
406  // The quadrature
407  // weights on the
408  // subcell are NOT
409  // transformed, so we
410  // have to do it here.
411  this->restriction[iso][child](
412  face * this->n_dofs_per_face(face) + i_face, i_child) +=
413  Utilities::fixed_power<dim - 1>(.5) * q_sub.weight(k) *
414  cached_values_face(i_child, k) *
415  this->shape_value_component(
416  face * this->n_dofs_per_face(face) + i_face,
417  q_sub.point(k),
419  }
420  }
421  }
422 
423  if (rt_order == 0)
424  return;
425 
426  // Create Legendre basis for the
427  // space D_xi Q_k. Here, we cannot
428  // use the shape functions
429  std::array<std::unique_ptr<AnisotropicPolynomials<dim>>, dim> polynomials;
430  for (unsigned int dd = 0; dd < dim; ++dd)
431  {
432  std::vector<std::vector<Polynomials::Polynomial<double>>> poly(dim);
433  for (unsigned int d = 0; d < dim; ++d)
435  poly[dd] = Polynomials::Legendre::generate_complete_basis(rt_order - 1);
436 
437  polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
438  }
439 
440  // TODO: the implementation makes the assumption that all faces have the
441  // same number of dofs
442  AssertDimension(this->n_unique_faces(), 1);
443  const unsigned int face_no = 0;
444 
445  QGauss<dim> q_cell(rt_order + 1);
446  const unsigned int start_cell_dofs =
447  GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face(face_no);
448 
449  // Store shape values, since the
450  // evaluation suffers if not
451  // ordered by point
452  Table<3, double> cached_values_cell(this->n_dofs_per_cell(),
453  q_cell.size(),
454  dim);
455  for (unsigned int k = 0; k < q_cell.size(); ++k)
456  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
457  for (unsigned int d = 0; d < dim; ++d)
458  cached_values_cell(i, k, d) =
459  this->shape_value_component(i, q_cell.point(k), d);
460 
461  for (unsigned int child = 0; child < GeometryInfo<dim>::max_children_per_cell;
462  ++child)
463  {
464  Quadrature<dim> q_sub =
466  q_cell,
467  child);
468 
469  for (unsigned int k = 0; k < q_sub.size(); ++k)
470  for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
471  ++i_child)
472  for (unsigned int d = 0; d < dim; ++d)
473  for (unsigned int i_weight = 0; i_weight < polynomials[d]->n();
474  ++i_weight)
475  {
476  this->restriction[iso][child](start_cell_dofs + i_weight * dim +
477  d,
478  i_child) +=
479  q_sub.weight(k) * cached_values_cell(i_child, k, d) *
480  polynomials[d]->compute_value(i_weight, q_sub.point(k));
481  }
482  }
483 }
484 
485 
486 
487 template <int dim>
488 std::vector<unsigned int>
489 FE_ABF<dim>::get_dpo_vector(const unsigned int rt_order)
490 {
491  if (dim == 1)
492  {
493  Assert(false, ExcImpossibleInDim(1));
494  return std::vector<unsigned int>();
495  }
496 
497  // the element is face-based (not
498  // to be confused with George
499  // W. Bush's Faith Based
500  // Initiative...), and we have
501  // (rt_order+1)^(dim-1) DoFs per face
502  unsigned int dofs_per_face = 1;
503  for (unsigned int d = 0; d < dim - 1; ++d)
504  dofs_per_face *= rt_order + 1;
505 
506  // and then there are interior dofs
507  const unsigned int interior_dofs = dim * (rt_order + 1) * dofs_per_face;
508 
509  std::vector<unsigned int> dpo(dim + 1);
510  dpo[dim - 1] = dofs_per_face;
511  dpo[dim] = interior_dofs;
512 
513  return dpo;
514 }
515 
516 //---------------------------------------------------------------------------
517 // Data field initialization
518 //---------------------------------------------------------------------------
519 
520 template <int dim>
521 bool
522 FE_ABF<dim>::has_support_on_face(const unsigned int shape_index,
523  const unsigned int face_index) const
524 {
525  AssertIndexRange(shape_index, this->n_dofs_per_cell());
527 
528  // Return computed values if we
529  // know them easily. Otherwise, it
530  // is always safe to return true.
531  switch (rt_order)
532  {
533  case 0:
534  {
535  switch (dim)
536  {
537  case 2:
538  {
539  // only on the one
540  // non-adjacent face
541  // are the values
542  // actually zero. list
543  // these in a table
544  return (face_index !=
545  GeometryInfo<dim>::opposite_face[shape_index]);
546  }
547 
548  default:
549  return true;
550  }
551  }
552 
553  default: // other rt_order
554  return true;
555  }
556 
557  return true;
558 }
559 
560 
561 
562 template <int dim>
563 void
565  const std::vector<Vector<double>> &support_point_values,
566  std::vector<double> &nodal_values) const
567 {
568  Assert(support_point_values.size() == this->generalized_support_points.size(),
569  ExcDimensionMismatch(support_point_values.size(),
570  this->generalized_support_points.size()));
571  Assert(support_point_values[0].size() == this->n_components(),
572  ExcDimensionMismatch(support_point_values[0].size(),
573  this->n_components()));
574  Assert(nodal_values.size() == this->n_dofs_per_cell(),
575  ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
576 
577  std::fill(nodal_values.begin(), nodal_values.end(), 0.);
578 
579  const unsigned int n_face_points = boundary_weights.size(0);
580  for (const unsigned int face : GeometryInfo<dim>::face_indices())
581  for (unsigned int k = 0; k < n_face_points; ++k)
582  for (unsigned int i = 0; i < boundary_weights.size(1); ++i)
583  {
584  nodal_values[i + face * this->n_dofs_per_face(face)] +=
585  boundary_weights(k, i) *
586  support_point_values[face * n_face_points + k][GeometryInfo<
587  dim>::unit_normal_direction[face]];
588  }
589 
590  // TODO: the implementation makes the assumption that all faces have the
591  // same number of dofs
592  AssertDimension(this->n_unique_faces(), 1);
593  const unsigned int face_no = 0;
594 
595  const unsigned int start_cell_dofs =
596  GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face(face_no);
597  const unsigned int start_cell_points =
598  GeometryInfo<dim>::faces_per_cell * n_face_points;
599 
600  for (unsigned int k = 0; k < interior_weights.size(0); ++k)
601  for (unsigned int i = 0; i < interior_weights.size(1); ++i)
602  for (unsigned int d = 0; d < dim; ++d)
603  nodal_values[start_cell_dofs + i * dim + d] +=
604  interior_weights(k, i, d) *
605  support_point_values[k + start_cell_points][d];
606 
607  const unsigned int start_abf_dofs =
608  start_cell_dofs + interior_weights.size(1) * dim;
609 
610  // Cell integral of ABF terms
611  for (unsigned int k = 0; k < interior_weights_abf.size(0); ++k)
612  for (unsigned int i = 0; i < interior_weights_abf.size(1); ++i)
613  for (unsigned int d = 0; d < dim; ++d)
614  nodal_values[start_abf_dofs + i] +=
615  interior_weights_abf(k, i, d) *
616  support_point_values[k + start_cell_points][d];
617 
618  // Face integral of ABF terms
619  for (const unsigned int face : GeometryInfo<dim>::face_indices())
620  {
621  const double n_orient = GeometryInfo<dim>::unit_normal_orientation[face];
622  for (unsigned int fp = 0; fp < n_face_points; ++fp)
623  {
624  // TODO: Check what the face_orientation, face_flip and face_rotation
625  // have to be in 3d
627  this->reference_cell(), face, false, false, false, n_face_points);
628  for (unsigned int i = 0; i < boundary_weights_abf.size(1); ++i)
629  nodal_values[start_abf_dofs + i] +=
630  n_orient * boundary_weights_abf(k + fp, i) *
631  support_point_values[face * n_face_points + fp][GeometryInfo<
632  dim>::unit_normal_direction[face]];
633  }
634  }
635 
636  // TODO: Check if this "correction" can be removed.
637  for (unsigned int i = 0; i < boundary_weights_abf.size(1); ++i)
638  if (std::fabs(nodal_values[start_abf_dofs + i]) < 1.0e-16)
639  nodal_values[start_abf_dofs + i] = 0.0;
640 }
641 
642 
643 
644 template <int dim>
645 std::size_t
647 {
648  Assert(false, ExcNotImplemented());
649  return 0;
650 }
651 
652 
653 
654 /*-------------- Explicit Instantiations -------------------------------*/
655 #include "fe_abf.inst"
656 
void initialize_quad_dof_index_permutation_and_sign_change()
Definition: fe_abf.cc:123
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
Definition: fe_abf.cc:156
virtual std::size_t memory_consumption() const override
Definition: fe_abf.cc:646
void initialize_restriction()
Definition: fe_abf.cc:342
virtual void convert_generalized_support_point_values_to_dof_values(const std::vector< Vector< double >> &support_point_values, std::vector< double > &nodal_values) const override
Definition: fe_abf.cc:564
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition: fe_abf.cc:522
virtual std::string get_name() const override
Definition: fe_abf.cc:136
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
Definition: fe_abf.cc:489
void initialize_support_points(const unsigned int rt_degree)
Definition: fe_abf.cc:171
friend class FE_ABF
Definition: fe_abf.h:256
FullMatrix< double > inverse_node_matrix
std::vector< MappingKind > mapping_kind
unsigned int n_dofs_per_cell() const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_unique_faces() const
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
FullMatrix< double > interface_constraints
Definition: fe.h:2428
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2416
void invert(const FullMatrix< number2 > &M)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:730
static std::vector< Polynomial< number > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:551
static DataSetDescriptor face(const ReferenceCell &reference_cell, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
Definition: qprojector.cc:1319
static void project_to_subface(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, const unsigned int subface_no, std::vector< Point< dim >> &q_points, const RefinementCase< dim - 1 > &ref_case=RefinementCase< dim - 1 >::isotropic_refinement)
static Quadrature< dim > project_to_child(const ReferenceCell &reference_cell, const Quadrature< dim > &quadrature, const unsigned int child_no)
Definition: qprojector.cc:1227
static Quadrature< dim > project_to_all_faces(const ReferenceCell &reference_cell, const hp::QCollection< dim - 1 > &quadrature)
static void project_to_face(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim >> &q_points)
double weight(const unsigned int i) const
const Point< dim > & point(const unsigned int i) const
unsigned int size() const
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
@ mapping_raviart_thomas
Definition: mapping.h:133
Expression fabs(const Expression &x)
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)
FullMatrix< double > compute_node_matrix(const FiniteElement< dim, spacedim > &fe)
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
double legendre(unsigned int l, double x)
Definition: cmath.h:75
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement)