Reference documentation for deal.II version Git 409ee4b167 2020-08-14 09:46:12 -0400
\(\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 - 2020 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 
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->n_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 
86  // TODO: the implementation makes the assumption that all faces have the
87  // same number of dofs
88  AssertDimension(this->n_unique_faces(), 1);
89  const unsigned int face_no = 0;
90 
91  // Fill prolongation matrices with embedding operators
93  // TODO[TL]: for anisotropic refinement we will probably need a table of
94  // submatrices with an array for each refine case
96  for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
97  face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
98  this->n_dofs_per_face(face_no));
99  FETools::compute_face_embedding_matrices<dim, double>(*this,
100  face_embeddings,
101  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 (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
108  for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
109  {
110  for (unsigned int j = 0; j < face_embeddings[d].n(); ++j)
111  this->interface_constraints(target_row, j) = face_embeddings[d](i, j);
112  ++target_row;
113  }
114 }
115 
116 
117 
118 template <int dim>
119 std::string
121 {
122  // note that the
123  // FETools::get_fe_by_name
124  // function depends on the
125  // particular format of the string
126  // this function returns, so they
127  // have to be kept in synch
128 
129  // note that this->degree is the maximal
130  // polynomial degree and is thus one higher
131  // than the argument given to the
132  // constructor
133  std::ostringstream namebuf;
134  namebuf << "FE_RaviartThomasNodal<" << dim << ">(" << this->degree - 1 << ")";
135 
136  return namebuf.str();
137 }
138 
139 
140 template <int dim>
141 std::unique_ptr<FiniteElement<dim, dim>>
143 {
144  return std::make_unique<FE_RaviartThomasNodal<dim>>(*this);
145 }
146 
147 
148 //---------------------------------------------------------------------------
149 // Auxiliary and internal functions
150 //---------------------------------------------------------------------------
151 
152 
153 
154 template <int dim>
155 void
157 {
158  // TODO: the implementation makes the assumption that all faces have the
159  // same number of dofs
160  AssertDimension(this->n_unique_faces(), 1);
161  const unsigned int face_no = 0;
162 
163  this->generalized_support_points.resize(this->n_dofs_per_cell());
164  this->generalized_face_support_points[face_no].resize(
165  this->n_dofs_per_face(face_no));
166 
167  // Number of the point being entered
168  unsigned int current = 0;
169 
170  // On the faces, we choose as many
171  // Gauss points as necessary to
172  // determine the normal component
173  // uniquely. This is the deg of
174  // the Raviart-Thomas element plus
175  // one.
176  if (dim > 1)
177  {
178  QGauss<dim - 1> face_points(deg + 1);
179  Assert(face_points.size() == this->n_dofs_per_face(face_no),
180  ExcInternalError());
181  for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
182  this->generalized_face_support_points[face_no][k] =
183  face_points.point(k);
184  Quadrature<dim> faces =
186  face_points);
187  for (unsigned int k = 0; k < this->n_dofs_per_face(face_no) *
189  ++k)
190  this->generalized_support_points[k] =
192  this->reference_cell_type(),
193  0,
194  true,
195  false,
196  false,
197  this->n_dofs_per_face(face_no)));
198 
199  current =
200  this->n_dofs_per_face(face_no) * GeometryInfo<dim>::faces_per_cell;
201  }
202 
203  if (deg == 0)
204  return;
205  // In the interior, we need
206  // anisotropic Gauss quadratures,
207  // different for each direction.
208  QGauss<1> high(deg + 1);
209  QGauss<1> low(deg);
210 
211  for (unsigned int d = 0; d < dim; ++d)
212  {
213  std::unique_ptr<QAnisotropic<dim>> quadrature;
214  switch (dim)
215  {
216  case 1:
217  quadrature = std::make_unique<QAnisotropic<dim>>(high);
218  break;
219  case 2:
220  quadrature =
221  std::make_unique<QAnisotropic<dim>>(((d == 0) ? low : high),
222  ((d == 1) ? low : high));
223  break;
224  case 3:
225  quadrature =
226  std::make_unique<QAnisotropic<dim>>(((d == 0) ? low : high),
227  ((d == 1) ? low : high),
228  ((d == 2) ? low : high));
229  break;
230  default:
231  Assert(false, ExcNotImplemented());
232  }
233 
234  for (unsigned int k = 0; k < quadrature->size(); ++k)
235  this->generalized_support_points[current++] = quadrature->point(k);
236  }
237  Assert(current == this->n_dofs_per_cell(), ExcInternalError());
238 }
239 
240 
241 
242 template <int dim>
243 std::vector<unsigned int>
245 {
246  // the element is face-based and we have
247  // (deg+1)^(dim-1) DoFs per face
248  unsigned int dofs_per_face = 1;
249  for (unsigned int d = 1; d < dim; ++d)
250  dofs_per_face *= deg + 1;
251 
252  // and then there are interior dofs
253  const unsigned int interior_dofs = dim * deg * dofs_per_face;
254 
255  std::vector<unsigned int> dpo(dim + 1);
256  dpo[dim - 1] = dofs_per_face;
257  dpo[dim] = interior_dofs;
258 
259  return dpo;
260 }
261 
262 
263 
264 template <>
265 std::vector<bool>
267 {
268  Assert(false, ExcImpossibleInDim(1));
269  return std::vector<bool>();
270 }
271 
272 
273 
274 template <int dim>
275 std::vector<bool>
277 {
278  const unsigned int dofs_per_cell =
280  unsigned int dofs_per_face = deg + 1;
281  for (unsigned int d = 2; d < dim; ++d)
282  dofs_per_face *= deg + 1;
283  // all face dofs need to be
284  // non-additive, since they have
285  // continuity requirements.
286  // however, the interior dofs are
287  // made additive
288  std::vector<bool> ret_val(dofs_per_cell, false);
289  for (unsigned int i = GeometryInfo<dim>::faces_per_cell * dofs_per_face;
290  i < dofs_per_cell;
291  ++i)
292  ret_val[i] = true;
293 
294  return ret_val;
295 }
296 
297 
298 template <int dim>
299 bool
301  const unsigned int shape_index,
302  const unsigned int face_index) const
303 {
304  AssertIndexRange(shape_index, this->n_dofs_per_cell());
306 
307  // The first degrees of freedom are
308  // on the faces and each face has
309  // degree degrees.
310  const unsigned int support_face = shape_index / this->degree;
311 
312  // The only thing we know for sure
313  // is that shape functions with
314  // support on one face are zero on
315  // the opposite face.
316  if (support_face < GeometryInfo<dim>::faces_per_cell)
317  return (face_index != GeometryInfo<dim>::opposite_face[support_face]);
318 
319  // In all other cases, return true,
320  // which is safe
321  return true;
322 }
323 
324 
325 
326 template <int dim>
327 void
330  const std::vector<Vector<double>> &support_point_values,
331  std::vector<double> & nodal_values) const
332 {
333  Assert(support_point_values.size() == this->generalized_support_points.size(),
334  ExcDimensionMismatch(support_point_values.size(),
335  this->generalized_support_points.size()));
336  Assert(nodal_values.size() == this->n_dofs_per_cell(),
337  ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
338  Assert(support_point_values[0].size() == this->n_components(),
339  ExcDimensionMismatch(support_point_values[0].size(),
340  this->n_components()));
341 
342  // First do interpolation on
343  // faces. There, the component
344  // evaluated depends on the face
345  // direction and orientation.
346  unsigned int fbase = 0;
347  unsigned int f = 0;
348  for (; f < GeometryInfo<dim>::faces_per_cell;
349  ++f, fbase += this->n_dofs_per_face(f))
350  {
351  for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
352  {
353  nodal_values[fbase + i] = support_point_values[fbase + i](
355  }
356  }
357 
358  // The remaining points form dim
359  // chunks, one for each component.
360  const unsigned int istep = (this->n_dofs_per_cell() - fbase) / dim;
361  Assert((this->n_dofs_per_cell() - fbase) % dim == 0, ExcInternalError());
362 
363  f = 0;
364  while (fbase < this->n_dofs_per_cell())
365  {
366  for (unsigned int i = 0; i < istep; ++i)
367  {
368  nodal_values[fbase + i] = support_point_values[fbase + i](f);
369  }
370  fbase += istep;
371  ++f;
372  }
373  Assert(fbase == this->n_dofs_per_cell(), ExcInternalError());
374 }
375 
376 
377 
378 // TODO: There are tests that check that the following few functions don't
379 // produce assertion failures, but none that actually check whether they do the
380 // right thing. one example for such a test would be to project a function onto
381 // an hp space and make sure that the convergence order is correct with regard
382 // to the lowest used polynomial degree
383 
384 template <int dim>
385 bool
387 {
388  return true;
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 these
398  // identities if both FEs are
399  // FE_RaviartThomasNodals or the other is FE_Nothing.
400  // In either case, no dofs are assigned on the vertex,
401  // so we shouldn't be getting here at all.
402  if (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other) != nullptr)
403  return std::vector<std::pair<unsigned int, unsigned int>>();
404  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
405  return std::vector<std::pair<unsigned int, unsigned int>>();
406  else
407  {
408  Assert(false, ExcNotImplemented());
409  return std::vector<std::pair<unsigned int, unsigned int>>();
410  }
411 }
412 
413 
414 
415 template <int dim>
416 std::vector<std::pair<unsigned int, unsigned int>>
418  const FiniteElement<dim> &fe_other) const
419 {
420  // we can presently only compute
421  // these identities if both FEs are
422  // FE_RaviartThomasNodals or if the other
423  // one is FE_Nothing
424  if (const FE_RaviartThomasNodal<dim> *fe_q_other =
425  dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
426  {
427  // dofs are located on faces; these are
428  // only lines in 2d
429  if (dim != 2)
430  return std::vector<std::pair<unsigned int, unsigned int>>();
431 
432  // dofs are located along lines, so two
433  // dofs are identical only if in the
434  // following two cases (remember that
435  // the face support points are Gauss
436  // points):
437  // 1. this->degree = fe_q_other->degree,
438  // in the case, all the dofs on
439  // the line are identical
440  // 2. this->degree-1 and fe_q_other->degree-1
441  // are both even, i.e. this->dof_per_line
442  // and fe_q_other->dof_per_line are both odd,
443  // there exists only one point (the middle one)
444  // such that dofs are identical on this point
445  //
446  // to understand this, note that
447  // this->degree is the *maximal*
448  // polynomial degree, and is thus one
449  // higher than the argument given to
450  // the constructor
451  const unsigned int p = this->degree - 1;
452  const unsigned int q = fe_q_other->degree - 1;
453 
454  std::vector<std::pair<unsigned int, unsigned int>> identities;
455 
456  if (p == q)
457  for (unsigned int i = 0; i < p + 1; ++i)
458  identities.emplace_back(i, i);
459 
460  else if (p % 2 == 0 && q % 2 == 0)
461  identities.emplace_back(p / 2, q / 2);
462 
463  return identities;
464  }
465  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
466  {
467  // the FE_Nothing has no degrees of freedom, so there are no
468  // equivalencies to be recorded
469  return std::vector<std::pair<unsigned int, unsigned int>>();
470  }
471  else
472  {
473  Assert(false, ExcNotImplemented());
474  return std::vector<std::pair<unsigned int, unsigned int>>();
475  }
476 }
477 
478 
479 template <int dim>
480 std::vector<std::pair<unsigned int, unsigned int>>
482  const FiniteElement<dim> &fe_other,
483  const unsigned int face_no) const
484 {
485  // we can presently only compute
486  // these identities if both FEs are
487  // FE_RaviartThomasNodals or if the other
488  // one is FE_Nothing
489  if (const FE_RaviartThomasNodal<dim> *fe_q_other =
490  dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
491  {
492  // dofs are located on faces; these are
493  // only quads in 3d
494  if (dim != 3)
495  return std::vector<std::pair<unsigned int, unsigned int>>();
496 
497  // this works exactly like the line
498  // case above
499  const unsigned int p = this->n_dofs_per_quad(face_no);
500 
501  AssertDimension(fe_q_other->n_unique_faces(), 1);
502  const unsigned int q = fe_q_other->n_dofs_per_quad(0);
503 
504  std::vector<std::pair<unsigned int, unsigned int>> identities;
505 
506  if (p == q)
507  for (unsigned int i = 0; i < p; ++i)
508  identities.emplace_back(i, i);
509 
510  else if (p % 2 != 0 && q % 2 != 0)
511  identities.emplace_back(p / 2, q / 2);
512 
513  return identities;
514  }
515  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
516  {
517  // the FE_Nothing has no degrees of freedom, so there are no
518  // equivalencies to be recorded
519  return std::vector<std::pair<unsigned int, unsigned int>>();
520  }
521  else
522  {
523  Assert(false, ExcNotImplemented());
524  return std::vector<std::pair<unsigned int, unsigned int>>();
525  }
526 }
527 
528 
529 template <int dim>
532  const FiniteElement<dim> &fe_other,
533  const unsigned int codim) const
534 {
535  Assert(codim <= dim, ExcImpossibleInDim(dim));
536  (void)codim;
537 
538  // vertex/line/face/cell domination
539  // --------------------------------
540  if (const FE_RaviartThomasNodal<dim> *fe_rt_nodal_other =
541  dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
542  {
543  if (this->degree < fe_rt_nodal_other->degree)
545  else if (this->degree == fe_rt_nodal_other->degree)
547  else
549  }
550  else if (const FE_Nothing<dim> *fe_nothing =
551  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
552  {
553  if (fe_nothing->is_dominating())
555  else
556  // the FE_Nothing has no degrees of freedom and it is typically used
557  // in a context where we don't require any continuity along the
558  // interface
560  }
561 
562  Assert(false, ExcNotImplemented());
564 }
565 
566 
567 
568 template <>
569 void
571  const FiniteElement<1, 1> & /*x_source_fe*/,
572  FullMatrix<double> & /*interpolation_matrix*/,
573  const unsigned int) const
574 {
575  Assert(false, ExcImpossibleInDim(1));
576 }
577 
578 
579 template <>
580 void
582  const FiniteElement<1, 1> & /*x_source_fe*/,
583  const unsigned int /*subface*/,
584  FullMatrix<double> & /*interpolation_matrix*/,
585  const unsigned int) const
586 {
587  Assert(false, ExcImpossibleInDim(1));
588 }
589 
590 
591 
592 template <int dim>
593 void
595  const FiniteElement<dim> &x_source_fe,
596  FullMatrix<double> & interpolation_matrix,
597  const unsigned int face_no) const
598 {
599  // this is only implemented, if the
600  // source FE is also a
601  // RaviartThomasNodal element
602  AssertThrow((x_source_fe.get_name().find("FE_RaviartThomasNodal<") == 0) ||
603  (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(
604  &x_source_fe) != nullptr),
606 
607  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
608  ExcDimensionMismatch(interpolation_matrix.n(),
609  this->n_dofs_per_face(face_no)));
610  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
611  ExcDimensionMismatch(interpolation_matrix.m(),
612  x_source_fe.n_dofs_per_face(face_no)));
613 
614  // ok, source is a RaviartThomasNodal element, so
615  // we will be able to do the work
616  const FE_RaviartThomasNodal<dim> &source_fe =
617  dynamic_cast<const FE_RaviartThomasNodal<dim> &>(x_source_fe);
618 
619  // Make sure, that the element,
620  // for which the DoFs should be
621  // constrained is the one with
622  // the higher polynomial degree.
623  // Actually the procedure will work
624  // also if this assertion is not
625  // satisfied. But the matrices
626  // produced in that case might
627  // lead to problems in the
628  // hp procedures, which use this
629  // method.
630  Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
632 
633  // generate a quadrature
634  // with the generalized support points.
635  // This is later based as a
636  // basis for the QProjector,
637  // which returns the support
638  // points on the face.
639  Quadrature<dim - 1> quad_face_support(
640  source_fe.generalized_face_support_points[face_no]);
641 
642  // Rule of thumb for FP accuracy,
643  // that can be expected for a
644  // given polynomial degree.
645  // This value is used to cut
646  // off values close to zero.
647  double eps = 2e-13 * this->degree * (dim - 1);
648 
649  // compute the interpolation
650  // matrix by simply taking the
651  // value at the support points.
652  const Quadrature<dim> face_projection =
654  quad_face_support,
655  0);
656 
657  for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
658  {
659  const Point<dim> &p = face_projection.point(i);
660 
661  for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
662  {
663  double matrix_entry =
664  this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
665 
666  // Correct the interpolated
667  // value. I.e. if it is close
668  // to 1 or 0, make it exactly
669  // 1 or 0. Unfortunately, this
670  // is required to avoid problems
671  // with higher order elements.
672  if (std::fabs(matrix_entry - 1.0) < eps)
673  matrix_entry = 1.0;
674  if (std::fabs(matrix_entry) < eps)
675  matrix_entry = 0.0;
676 
677  interpolation_matrix(i, j) = matrix_entry;
678  }
679  }
680 
681  // make sure that the row sum of
682  // each of the matrices is 1 at
683  // this point. this must be so
684  // since the shape functions sum up
685  // to 1
686  for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
687  {
688  double sum = 0.;
689 
690  for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
691  sum += interpolation_matrix(j, i);
692 
693  Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
694  ExcInternalError());
695  }
696 }
697 
698 
699 template <int dim>
700 void
702  const FiniteElement<dim> &x_source_fe,
703  const unsigned int subface,
704  FullMatrix<double> & interpolation_matrix,
705  const unsigned int face_no) const
706 {
707  // this is only implemented, if the
708  // source FE is also a
709  // RaviartThomasNodal element
710  AssertThrow((x_source_fe.get_name().find("FE_RaviartThomasNodal<") == 0) ||
711  (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(
712  &x_source_fe) != nullptr),
714 
715  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
716  ExcDimensionMismatch(interpolation_matrix.n(),
717  this->n_dofs_per_face(face_no)));
718  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
719  ExcDimensionMismatch(interpolation_matrix.m(),
720  x_source_fe.n_dofs_per_face(face_no)));
721 
722  // ok, source is a RaviartThomasNodal element, so
723  // we will be able to do the work
724  const FE_RaviartThomasNodal<dim> &source_fe =
725  dynamic_cast<const FE_RaviartThomasNodal<dim> &>(x_source_fe);
726 
727  // Make sure, that the element,
728  // for which the DoFs should be
729  // constrained is the one with
730  // the higher polynomial degree.
731  // Actually the procedure will work
732  // also if this assertion is not
733  // satisfied. But the matrices
734  // produced in that case might
735  // lead to problems in the
736  // hp procedures, which use this
737  // method.
738  Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
740 
741  // generate a quadrature
742  // with the generalized support points.
743  // This is later based as a
744  // basis for the QProjector,
745  // which returns the support
746  // points on the face.
747  Quadrature<dim - 1> quad_face_support(
748  source_fe.generalized_face_support_points[face_no]);
749 
750  // Rule of thumb for FP accuracy,
751  // that can be expected for a
752  // given polynomial degree.
753  // This value is used to cut
754  // off values close to zero.
755  double eps = 2e-13 * this->degree * (dim - 1);
756 
757  // compute the interpolation
758  // matrix by simply taking the
759  // value at the support points.
760 
761  const Quadrature<dim> subface_projection =
763  quad_face_support,
764  0,
765  subface);
766 
767  for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
768  {
769  const Point<dim> &p = subface_projection.point(i);
770 
771  for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
772  {
773  double matrix_entry =
774  this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
775 
776  // Correct the interpolated
777  // value. I.e. if it is close
778  // to 1 or 0, make it exactly
779  // 1 or 0. Unfortunately, this
780  // is required to avoid problems
781  // with higher order elements.
782  if (std::fabs(matrix_entry - 1.0) < eps)
783  matrix_entry = 1.0;
784  if (std::fabs(matrix_entry) < eps)
785  matrix_entry = 0.0;
786 
787  interpolation_matrix(i, j) = matrix_entry;
788  }
789  }
790 
791  // make sure that the row sum of
792  // each of the matrices is 1 at
793  // this point. this must be so
794  // since the shape functions sum up
795  // to 1
796  for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
797  {
798  double sum = 0.;
799 
800  for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
801  sum += interpolation_matrix(j, i);
802 
803  Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
804  ExcInternalError());
805  }
806 }
807 
808 
809 
810 // explicit instantiations
811 #include "fe_raviart_thomas_nodal.inst"
812 
813 
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:1568
static unsigned int n_polynomials(const unsigned int degree)
std::vector< Point< dim > > generalized_support_points
Definition: fe.h:2452
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:2427
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)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1636
const unsigned int degree
Definition: fe_base.h:431
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:1521
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()
ReferenceCell::Type reference_cell_type() const
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:2415
T sum(const T &t, const MPI_Comm &mpi_communicator)
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)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
virtual std::string get_name() const =0
std::vector< std::vector< Point< dim - 1 > > > generalized_face_support_points
Definition: fe.h:2458
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:415
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
unsigned int n_dofs_per_cell() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
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:401
static ::ExceptionBase & ExcNotImplemented()
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
static ::ExceptionBase & ExcInternalError()
static std::vector< bool > get_ria_vector(const unsigned int degree)
std::vector< MappingKind > mapping_kind