Reference documentation for deal.II version 9.3.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_raviart_thomas_nodal.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2003 - 2021 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 
19 
21 
22 #include <deal.II/fe/fe.h>
23 #include <deal.II/fe/fe_nothing.h>
25 #include <deal.II/fe/fe_tools.h>
26 #include <deal.II/fe/fe_values.h>
27 #include <deal.II/fe/mapping.h>
28 
29 #include <deal.II/grid/tria.h>
31 
32 #include <memory>
33 #include <sstream>
34 
35 
37 
38 // TODO: implement the adjust_quad_dof_index_for_face_orientation_table and
39 // adjust_line_dof_index_for_line_orientation_table fields, and write tests
40 // similar to bits/face_orientation_and_fe_q_*
41 
42 template <int dim>
44  : FE_PolyTensor<dim>(PolynomialsRaviartThomas<dim>(deg),
46  dim,
47  deg + 1,
48  FiniteElementData<dim>::Hdiv),
49  get_ria_vector(deg),
50  std::vector<ComponentMask>(
51  PolynomialsRaviartThomas<dim>::n_polynomials(deg),
52  std::vector<bool>(dim, true)))
53 {
54  Assert(dim >= 2, ExcImpossibleInDim(dim));
55  const unsigned int n_dofs = this->n_dofs_per_cell();
56 
58  // First, initialize the
59  // generalized support points and
60  // quadrature weights, since they
61  // are required for interpolation.
63 
64  // Now compute the inverse node matrix, generating the correct
65  // basis functions from the raw ones. For a discussion of what
66  // exactly happens here, see FETools::compute_node_matrix.
68  this->inverse_node_matrix.reinit(n_dofs, n_dofs);
69  this->inverse_node_matrix.invert(M);
70  // From now on, the shape functions provided by FiniteElement::shape_value
71  // and similar functions will be the correct ones, not
72  // the raw shape functions from the polynomial space anymore.
73 
74  // Reinit the vectors of
75  // prolongation matrices to the
76  // right sizes. There are no
77  // restriction matrices implemented
78  for (unsigned int ref_case = RefinementCase<dim>::cut_x;
79  ref_case < RefinementCase<dim>::isotropic_refinement + 1;
80  ++ref_case)
81  {
82  const unsigned int nc =
84 
85  for (unsigned int i = 0; i < nc; ++i)
86  this->prolongation[ref_case - 1][i].reinit(n_dofs, n_dofs);
87  }
88 
89  // TODO: the implementation makes the assumption that all faces have the
90  // same number of dofs
91  AssertDimension(this->n_unique_faces(), 1);
92  const unsigned int face_no = 0;
93 
94  // Fill prolongation matrices with embedding operators
96  // TODO[TL]: for anisotropic refinement we will probably need a table of
97  // submatrices with an array for each refine case
99  for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
100  face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
101  this->n_dofs_per_face(face_no));
102  FETools::compute_face_embedding_matrices<dim, double>(*this,
103  face_embeddings,
104  0,
105  0);
106  this->interface_constraints.reinit((1 << (dim - 1)) *
107  this->n_dofs_per_face(face_no),
108  this->n_dofs_per_face(face_no));
109  unsigned int target_row = 0;
110  for (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
111  for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
112  {
113  for (unsigned int j = 0; j < face_embeddings[d].n(); ++j)
114  this->interface_constraints(target_row, j) = face_embeddings[d](i, j);
115  ++target_row;
116  }
117 
118  // We need to initialize the dof permutation table and the one for the sign
119  // change.
121 }
122 
123 
124 
125 template <int dim>
126 std::string
128 {
129  // note that the
130  // FETools::get_fe_by_name
131  // function depends on the
132  // particular format of the string
133  // this function returns, so they
134  // have to be kept in synch
135 
136  // note that this->degree is the maximal
137  // polynomial degree and is thus one higher
138  // than the argument given to the
139  // constructor
140  std::ostringstream namebuf;
141  namebuf << "FE_RaviartThomasNodal<" << dim << ">(" << this->degree - 1 << ")";
142 
143  return namebuf.str();
144 }
145 
146 
147 template <int dim>
148 std::unique_ptr<FiniteElement<dim, dim>>
150 {
151  return std::make_unique<FE_RaviartThomasNodal<dim>>(*this);
152 }
153 
154 
155 //---------------------------------------------------------------------------
156 // Auxiliary and internal functions
157 //---------------------------------------------------------------------------
158 
159 
160 
161 template <int dim>
162 void
164 {
165  // TODO: the implementation makes the assumption that all faces have the
166  // same number of dofs
167  AssertDimension(this->n_unique_faces(), 1);
168  const unsigned int face_no = 0;
169 
170  this->generalized_support_points.resize(this->n_dofs_per_cell());
171  this->generalized_face_support_points[face_no].resize(
172  this->n_dofs_per_face(face_no));
173 
174  // Number of the point being entered
175  unsigned int current = 0;
176 
177  // On the faces, we choose as many
178  // Gauss points as necessary to
179  // determine the normal component
180  // uniquely. This is the deg of
181  // the Raviart-Thomas element plus
182  // one.
183  if (dim > 1)
184  {
185  QGauss<dim - 1> face_points(deg + 1);
186  Assert(face_points.size() == this->n_dofs_per_face(face_no),
187  ExcInternalError());
188  for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
189  this->generalized_face_support_points[face_no][k] =
190  face_points.point(k);
191  Quadrature<dim> faces =
193  face_points);
194  for (unsigned int k = 0; k < this->n_dofs_per_face(face_no) *
196  ++k)
197  this->generalized_support_points[k] = faces.point(
198  k + QProjector<dim>::DataSetDescriptor::face(this->reference_cell(),
199  0,
200  true,
201  false,
202  false,
203  this->n_dofs_per_face(
204  face_no)));
205 
206  current =
207  this->n_dofs_per_face(face_no) * GeometryInfo<dim>::faces_per_cell;
208  }
209 
210  if (deg == 0)
211  return;
212  // In the interior, we need
213  // anisotropic Gauss quadratures,
214  // different for each direction.
215  QGauss<1> high(deg + 1);
216  QGauss<1> low(deg);
217 
218  for (unsigned int d = 0; d < dim; ++d)
219  {
220  std::unique_ptr<QAnisotropic<dim>> quadrature;
221  switch (dim)
222  {
223  case 1:
224  quadrature = std::make_unique<QAnisotropic<dim>>(high);
225  break;
226  case 2:
227  quadrature =
228  std::make_unique<QAnisotropic<dim>>(((d == 0) ? low : high),
229  ((d == 1) ? low : high));
230  break;
231  case 3:
232  quadrature =
233  std::make_unique<QAnisotropic<dim>>(((d == 0) ? low : high),
234  ((d == 1) ? low : high),
235  ((d == 2) ? low : high));
236  break;
237  default:
238  Assert(false, ExcNotImplemented());
239  }
240 
241  for (unsigned int k = 0; k < quadrature->size(); ++k)
242  this->generalized_support_points[current++] = quadrature->point(k);
243  }
244  Assert(current == this->n_dofs_per_cell(), ExcInternalError());
245 }
246 
247 
248 
249 template <int dim>
250 void
253 {
254  // for 1D and 2D, do nothing
255  if (dim < 3)
256  return;
257 
258  // TODO: Implement this for this class
259  return;
260 }
261 
262 
263 
264 template <int dim>
265 std::vector<unsigned int>
267 {
268  // the element is face-based and we have
269  // (deg+1)^(dim-1) DoFs per face
270  unsigned int dofs_per_face = 1;
271  for (unsigned int d = 1; d < dim; ++d)
272  dofs_per_face *= deg + 1;
273 
274  // and then there are interior dofs
275  const unsigned int interior_dofs = dim * deg * dofs_per_face;
276 
277  std::vector<unsigned int> dpo(dim + 1);
278  dpo[dim - 1] = dofs_per_face;
279  dpo[dim] = interior_dofs;
280 
281  return dpo;
282 }
283 
284 
285 
286 template <>
287 std::vector<bool>
289 {
290  Assert(false, ExcImpossibleInDim(1));
291  return std::vector<bool>();
292 }
293 
294 
295 
296 template <int dim>
297 std::vector<bool>
299 {
300  const unsigned int dofs_per_cell =
302  unsigned int dofs_per_face = deg + 1;
303  for (unsigned int d = 2; d < dim; ++d)
304  dofs_per_face *= deg + 1;
305  // all face dofs need to be
306  // non-additive, since they have
307  // continuity requirements.
308  // however, the interior dofs are
309  // made additive
310  std::vector<bool> ret_val(dofs_per_cell, false);
311  for (unsigned int i = GeometryInfo<dim>::faces_per_cell * dofs_per_face;
312  i < dofs_per_cell;
313  ++i)
314  ret_val[i] = true;
315 
316  return ret_val;
317 }
318 
319 
320 template <int dim>
321 bool
323  const unsigned int shape_index,
324  const unsigned int face_index) const
325 {
326  AssertIndexRange(shape_index, this->n_dofs_per_cell());
328 
329  // The first degrees of freedom are
330  // on the faces and each face has
331  // degree degrees.
332  const unsigned int support_face = shape_index / this->degree;
333 
334  // The only thing we know for sure
335  // is that shape functions with
336  // support on one face are zero on
337  // the opposite face.
338  if (support_face < GeometryInfo<dim>::faces_per_cell)
339  return (face_index != GeometryInfo<dim>::opposite_face[support_face]);
340 
341  // In all other cases, return true,
342  // which is safe
343  return true;
344 }
345 
346 
347 
348 template <int dim>
349 void
352  const std::vector<Vector<double>> &support_point_values,
353  std::vector<double> & nodal_values) const
354 {
355  Assert(support_point_values.size() == this->generalized_support_points.size(),
356  ExcDimensionMismatch(support_point_values.size(),
357  this->generalized_support_points.size()));
358  Assert(nodal_values.size() == this->n_dofs_per_cell(),
359  ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
360  Assert(support_point_values[0].size() == this->n_components(),
361  ExcDimensionMismatch(support_point_values[0].size(),
362  this->n_components()));
363 
364  // First do interpolation on
365  // faces. There, the component
366  // evaluated depends on the face
367  // direction and orientation.
368  unsigned int fbase = 0;
369  unsigned int f = 0;
370  for (; f < GeometryInfo<dim>::faces_per_cell;
371  ++f, fbase += this->n_dofs_per_face(f))
372  {
373  for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
374  {
375  nodal_values[fbase + i] = support_point_values[fbase + i](
377  }
378  }
379 
380  // The remaining points form dim
381  // chunks, one for each component.
382  const unsigned int istep = (this->n_dofs_per_cell() - fbase) / dim;
383  Assert((this->n_dofs_per_cell() - fbase) % dim == 0, ExcInternalError());
384 
385  f = 0;
386  while (fbase < this->n_dofs_per_cell())
387  {
388  for (unsigned int i = 0; i < istep; ++i)
389  {
390  nodal_values[fbase + i] = support_point_values[fbase + i](f);
391  }
392  fbase += istep;
393  ++f;
394  }
395  Assert(fbase == this->n_dofs_per_cell(), ExcInternalError());
396 }
397 
398 
399 
400 // TODO: There are tests that check that the following few functions don't
401 // produce assertion failures, but none that actually check whether they do the
402 // right thing. one example for such a test would be to project a function onto
403 // an hp-space and make sure that the convergence order is correct with regard
404 // to the lowest used polynomial degree
405 
406 template <int dim>
407 bool
409 {
410  return true;
411 }
412 
413 
414 template <int dim>
415 std::vector<std::pair<unsigned int, unsigned int>>
417  const FiniteElement<dim> &fe_other) const
418 {
419  // we can presently only compute these
420  // identities if both FEs are
421  // FE_RaviartThomasNodals or the other is FE_Nothing.
422  // In either case, no dofs are assigned on the vertex,
423  // so we shouldn't be getting here at all.
424  if (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other) != nullptr)
425  return std::vector<std::pair<unsigned int, unsigned int>>();
426  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
427  return std::vector<std::pair<unsigned int, unsigned int>>();
428  else
429  {
430  Assert(false, ExcNotImplemented());
431  return std::vector<std::pair<unsigned int, unsigned int>>();
432  }
433 }
434 
435 
436 
437 template <int dim>
438 std::vector<std::pair<unsigned int, unsigned int>>
440  const FiniteElement<dim> &fe_other) const
441 {
442  // we can presently only compute
443  // these identities if both FEs are
444  // FE_RaviartThomasNodals or if the other
445  // one is FE_Nothing
446  if (const FE_RaviartThomasNodal<dim> *fe_q_other =
447  dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
448  {
449  // dofs are located on faces; these are
450  // only lines in 2d
451  if (dim != 2)
452  return std::vector<std::pair<unsigned int, unsigned int>>();
453 
454  // dofs are located along lines, so two
455  // dofs are identical only if in the
456  // following two cases (remember that
457  // the face support points are Gauss
458  // points):
459  // 1. this->degree = fe_q_other->degree,
460  // in the case, all the dofs on
461  // the line are identical
462  // 2. this->degree-1 and fe_q_other->degree-1
463  // are both even, i.e. this->dof_per_line
464  // and fe_q_other->dof_per_line are both odd,
465  // there exists only one point (the middle one)
466  // such that dofs are identical on this point
467  //
468  // to understand this, note that
469  // this->degree is the *maximal*
470  // polynomial degree, and is thus one
471  // higher than the argument given to
472  // the constructor
473  const unsigned int p = this->degree - 1;
474  const unsigned int q = fe_q_other->degree - 1;
475 
476  std::vector<std::pair<unsigned int, unsigned int>> identities;
477 
478  if (p == q)
479  for (unsigned int i = 0; i < p + 1; ++i)
480  identities.emplace_back(i, i);
481 
482  else if (p % 2 == 0 && q % 2 == 0)
483  identities.emplace_back(p / 2, q / 2);
484 
485  return identities;
486  }
487  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
488  {
489  // the FE_Nothing has no degrees of freedom, so there are no
490  // equivalencies to be recorded
491  return std::vector<std::pair<unsigned int, unsigned int>>();
492  }
493  else
494  {
495  Assert(false, ExcNotImplemented());
496  return std::vector<std::pair<unsigned int, unsigned int>>();
497  }
498 }
499 
500 
501 template <int dim>
502 std::vector<std::pair<unsigned int, unsigned int>>
504  const FiniteElement<dim> &fe_other,
505  const unsigned int face_no) const
506 {
507  // we can presently only compute
508  // these identities if both FEs are
509  // FE_RaviartThomasNodals or if the other
510  // one is FE_Nothing
511  if (const FE_RaviartThomasNodal<dim> *fe_q_other =
512  dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
513  {
514  // dofs are located on faces; these are
515  // only quads in 3d
516  if (dim != 3)
517  return std::vector<std::pair<unsigned int, unsigned int>>();
518 
519  // this works exactly like the line
520  // case above
521  const unsigned int p = this->n_dofs_per_quad(face_no);
522 
523  AssertDimension(fe_q_other->n_unique_faces(), 1);
524  const unsigned int q = fe_q_other->n_dofs_per_quad(0);
525 
526  std::vector<std::pair<unsigned int, unsigned int>> identities;
527 
528  if (p == q)
529  for (unsigned int i = 0; i < p; ++i)
530  identities.emplace_back(i, i);
531 
532  else if (p % 2 != 0 && q % 2 != 0)
533  identities.emplace_back(p / 2, q / 2);
534 
535  return identities;
536  }
537  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
538  {
539  // the FE_Nothing has no degrees of freedom, so there are no
540  // equivalencies to be recorded
541  return std::vector<std::pair<unsigned int, unsigned int>>();
542  }
543  else
544  {
545  Assert(false, ExcNotImplemented());
546  return std::vector<std::pair<unsigned int, unsigned int>>();
547  }
548 }
549 
550 
551 template <int dim>
554  const FiniteElement<dim> &fe_other,
555  const unsigned int codim) const
556 {
557  Assert(codim <= dim, ExcImpossibleInDim(dim));
558  (void)codim;
559 
560  // vertex/line/face/cell domination
561  // --------------------------------
562  if (const FE_RaviartThomasNodal<dim> *fe_rt_nodal_other =
563  dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
564  {
565  if (this->degree < fe_rt_nodal_other->degree)
567  else if (this->degree == fe_rt_nodal_other->degree)
569  else
571  }
572  else if (const FE_Nothing<dim> *fe_nothing =
573  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
574  {
575  if (fe_nothing->is_dominating())
577  else
578  // the FE_Nothing has no degrees of freedom and it is typically used
579  // in a context where we don't require any continuity along the
580  // interface
582  }
583 
584  Assert(false, ExcNotImplemented());
586 }
587 
588 
589 
590 template <>
591 void
593  const FiniteElement<1, 1> & /*x_source_fe*/,
594  FullMatrix<double> & /*interpolation_matrix*/,
595  const unsigned int) const
596 {
597  Assert(false, ExcImpossibleInDim(1));
598 }
599 
600 
601 template <>
602 void
604  const FiniteElement<1, 1> & /*x_source_fe*/,
605  const unsigned int /*subface*/,
606  FullMatrix<double> & /*interpolation_matrix*/,
607  const unsigned int) const
608 {
609  Assert(false, ExcImpossibleInDim(1));
610 }
611 
612 
613 
614 template <int dim>
615 void
617  const FiniteElement<dim> &x_source_fe,
618  FullMatrix<double> & interpolation_matrix,
619  const unsigned int face_no) const
620 {
621  // this is only implemented, if the
622  // source FE is also a
623  // RaviartThomasNodal element
624  AssertThrow((x_source_fe.get_name().find("FE_RaviartThomasNodal<") == 0) ||
625  (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(
626  &x_source_fe) != nullptr),
628 
629  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
630  ExcDimensionMismatch(interpolation_matrix.n(),
631  this->n_dofs_per_face(face_no)));
632  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
633  ExcDimensionMismatch(interpolation_matrix.m(),
634  x_source_fe.n_dofs_per_face(face_no)));
635 
636  // ok, source is a RaviartThomasNodal element, so
637  // we will be able to do the work
638  const FE_RaviartThomasNodal<dim> &source_fe =
639  dynamic_cast<const FE_RaviartThomasNodal<dim> &>(x_source_fe);
640 
641  // Make sure, that the element,
642  // for which the DoFs should be
643  // constrained is the one with
644  // the higher polynomial degree.
645  // Actually the procedure will work
646  // also if this assertion is not
647  // satisfied. But the matrices
648  // produced in that case might
649  // lead to problems in the
650  // hp-procedures, which use this
651  // method.
652  Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
654 
655  // generate a quadrature
656  // with the generalized support points.
657  // This is later based as a
658  // basis for the QProjector,
659  // which returns the support
660  // points on the face.
661  Quadrature<dim - 1> quad_face_support(
662  source_fe.generalized_face_support_points[face_no]);
663 
664  // Rule of thumb for FP accuracy,
665  // that can be expected for a
666  // given polynomial degree.
667  // This value is used to cut
668  // off values close to zero.
669  double eps = 2e-13 * this->degree * (dim - 1);
670 
671  // compute the interpolation
672  // matrix by simply taking the
673  // value at the support points.
674  const Quadrature<dim> face_projection =
676  quad_face_support,
677  0);
678 
679  for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
680  {
681  const Point<dim> &p = face_projection.point(i);
682 
683  for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
684  {
685  double matrix_entry =
686  this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
687 
688  // Correct the interpolated
689  // value. I.e. if it is close
690  // to 1 or 0, make it exactly
691  // 1 or 0. Unfortunately, this
692  // is required to avoid problems
693  // with higher order elements.
694  if (std::fabs(matrix_entry - 1.0) < eps)
695  matrix_entry = 1.0;
696  if (std::fabs(matrix_entry) < eps)
697  matrix_entry = 0.0;
698 
699  interpolation_matrix(i, j) = matrix_entry;
700  }
701  }
702 
703  // make sure that the row sum of
704  // each of the matrices is 1 at
705  // this point. this must be so
706  // since the shape functions sum up
707  // to 1
708  for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
709  {
710  double sum = 0.;
711 
712  for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
713  sum += interpolation_matrix(j, i);
714 
715  Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
716  ExcInternalError());
717  }
718 }
719 
720 
721 template <int dim>
722 void
724  const FiniteElement<dim> &x_source_fe,
725  const unsigned int subface,
726  FullMatrix<double> & interpolation_matrix,
727  const unsigned int face_no) const
728 {
729  // this is only implemented, if the
730  // source FE is also a
731  // RaviartThomasNodal element
732  AssertThrow((x_source_fe.get_name().find("FE_RaviartThomasNodal<") == 0) ||
733  (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(
734  &x_source_fe) != nullptr),
736 
737  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
738  ExcDimensionMismatch(interpolation_matrix.n(),
739  this->n_dofs_per_face(face_no)));
740  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
741  ExcDimensionMismatch(interpolation_matrix.m(),
742  x_source_fe.n_dofs_per_face(face_no)));
743 
744  // ok, source is a RaviartThomasNodal element, so
745  // we will be able to do the work
746  const FE_RaviartThomasNodal<dim> &source_fe =
747  dynamic_cast<const FE_RaviartThomasNodal<dim> &>(x_source_fe);
748 
749  // Make sure, that the element,
750  // for which the DoFs should be
751  // constrained is the one with
752  // the higher polynomial degree.
753  // Actually the procedure will work
754  // also if this assertion is not
755  // satisfied. But the matrices
756  // produced in that case might
757  // lead to problems in the
758  // hp-procedures, which use this
759  // method.
760  Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
762 
763  // generate a quadrature
764  // with the generalized support points.
765  // This is later based as a
766  // basis for the QProjector,
767  // which returns the support
768  // points on the face.
769  Quadrature<dim - 1> quad_face_support(
770  source_fe.generalized_face_support_points[face_no]);
771 
772  // Rule of thumb for FP accuracy,
773  // that can be expected for a
774  // given polynomial degree.
775  // This value is used to cut
776  // off values close to zero.
777  double eps = 2e-13 * this->degree * (dim - 1);
778 
779  // compute the interpolation
780  // matrix by simply taking the
781  // value at the support points.
782 
783  const Quadrature<dim> subface_projection =
785  quad_face_support,
786  0,
787  subface);
788 
789  for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
790  {
791  const Point<dim> &p = subface_projection.point(i);
792 
793  for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
794  {
795  double matrix_entry =
796  this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
797 
798  // Correct the interpolated
799  // value. I.e. if it is close
800  // to 1 or 0, make it exactly
801  // 1 or 0. Unfortunately, this
802  // is required to avoid problems
803  // with higher order elements.
804  if (std::fabs(matrix_entry - 1.0) < eps)
805  matrix_entry = 1.0;
806  if (std::fabs(matrix_entry) < eps)
807  matrix_entry = 0.0;
808 
809  interpolation_matrix(i, j) = matrix_entry;
810  }
811  }
812 
813  // make sure that the row sum of
814  // each of the matrices is 1 at
815  // this point. this must be so
816  // since the shape functions sum up
817  // to 1
818  for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
819  {
820  double sum = 0.;
821 
822  for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
823  sum += interpolation_matrix(j, i);
824 
825  Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
826  ExcInternalError());
827  }
828 }
829 
830 
831 
832 // explicit instantiations
833 #include "fe_raviart_thomas_nodal.inst"
834 
835 
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim > &fe_other, const unsigned int face_no=0) const override
size_type m() const
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
static unsigned int n_polynomials(const unsigned int degree)
std::vector< Point< dim > > generalized_support_points
Definition: fe.h:2455
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 > interface_constraints
Definition: fe.h:2430
static void project_to_subface(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)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
static Quadrature< dim > project_to_all_faces(const Quadrature< dim - 1 > &quadrature)
Definition: qprojector.h:579
FE_RaviartThomasNodal(const unsigned int p)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
const unsigned int degree
Definition: fe_base.h:435
const Point< dim > & point(const unsigned int i) const
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
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
void invert(const FullMatrix< number2 > &M)
STL namespace.
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
unsigned int n_dofs_per_quad(unsigned int face_no=0) const
FullMatrix< double > inverse_node_matrix
FullMatrix< double > compute_node_matrix(const FiniteElement< dim, spacedim > &fe)
virtual void get_face_interpolation_matrix(const FiniteElement< dim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
static ::ExceptionBase & ExcInterpolationNotImplemented()
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
virtual bool hp_constraints_are_implemented() const override
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
size_type n() const
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2418
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void initialize_support_points(const unsigned int rt_degree)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
virtual std::string get_name() const =0
std::vector< std::vector< Point< dim - 1 > > > generalized_face_support_points
Definition: fe.h:2461
Expression fabs(const Expression &x)
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false) const
Definition: fe.cc:568
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_unique_faces() const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim > &fe_other) const override
const unsigned int dofs_per_cell
Definition: fe_base.h:419
virtual std::string get_name() const override
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
static void project_to_face(const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim >> &q_points)
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim > &fe_other, const unsigned int codim=0) const override final
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
unsigned int n_components() const
unsigned int n_dofs_per_cell() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim > &fe_other) const override
const unsigned int dofs_per_face
Definition: fe_base.h:405
static ::ExceptionBase & ExcNotImplemented()
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
virtual void get_subface_interpolation_matrix(const FiniteElement< dim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
void initialize_quad_dof_index_permutation_and_sign_change()
ReferenceCell reference_cell() const
static ::ExceptionBase & ExcInternalError()
static std::vector< bool > get_ria_vector(const unsigned int degree)
std::vector< MappingKind > mapping_kind