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