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