Reference documentation for deal.II version 9.2.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.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 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 
19 
21 
22 #include <deal.II/fe/fe.h>
23 #include <deal.II/fe/fe_values.h>
24 #include <deal.II/fe/mapping.h>
25 
26 #include <deal.II/grid/tria.h>
28 
29 #include <algorithm>
30 #include <functional>
31 #include <numeric>
32 #include <typeinfo>
33 
35 
36 
37 /*------------------------------- FiniteElement ----------------------*/
38 
39 
40 template <int dim, int spacedim>
42  : update_each(update_default)
43 {}
44 
45 
46 
47 template <int dim, int spacedim>
48 std::size_t
50 {
51  return sizeof(*this);
52 }
53 
54 
55 
56 template <int dim, int spacedim>
58  const FiniteElementData<dim> & fe_data,
59  const std::vector<bool> & r_i_a_f,
60  const std::vector<ComponentMask> &nonzero_c)
61  : FiniteElementData<dim>(fe_data)
62  , adjust_quad_dof_index_for_face_orientation_table(dim == 3 ?
63  this->dofs_per_quad :
64  0,
65  dim == 3 ? 8 : 0)
66  , adjust_line_dof_index_for_line_orientation_table(
67  dim == 3 ? this->dofs_per_line : 0)
68  , system_to_base_table(this->dofs_per_cell)
69  , face_system_to_base_table(this->dofs_per_face)
70  , component_to_base_table(this->components,
71  std::make_pair(std::make_pair(0U, 0U), 0U))
72  ,
73 
74  // Special handling of vectors of length one: in this case, we
75  // assume that all entries were supposed to be equal
76  restriction_is_additive_flags(
77  r_i_a_f.size() == 1 ? std::vector<bool>(fe_data.dofs_per_cell, r_i_a_f[0]) :
78  r_i_a_f)
79  , nonzero_components(
80  nonzero_c.size() == 1 ?
81  std::vector<ComponentMask>(fe_data.dofs_per_cell, nonzero_c[0]) :
82  nonzero_c)
83  , n_nonzero_components_table(compute_n_nonzero_components(nonzero_components))
84  , cached_primitivity(std::find_if(n_nonzero_components_table.begin(),
85  n_nonzero_components_table.end(),
86  [](const unsigned int n_components) {
87  return n_components != 1U;
88  }) == n_nonzero_components_table.end())
89 {
92  this->dofs_per_cell));
94  for (unsigned int i = 0; i < nonzero_components.size(); ++i)
95  {
96  Assert(nonzero_components[i].size() == this->n_components(),
98  Assert(nonzero_components[i].n_selected_components() >= 1,
102  ExcInternalError());
103  }
104 
105  // initialize some tables in the default way, i.e. if there is only one
106  // (vector-)component; if the element is not primitive, leave these tables
107  // empty.
108  if (this->is_primitive())
109  {
112  for (unsigned int j = 0; j < this->dofs_per_cell; ++j)
113  system_to_component_table[j] = std::pair<unsigned, unsigned>(0, j);
114  for (unsigned int j = 0; j < this->dofs_per_face; ++j)
115  face_system_to_component_table[j] = std::pair<unsigned, unsigned>(0, j);
116  }
117 
118  for (unsigned int j = 0; j < this->dofs_per_cell; ++j)
119  system_to_base_table[j] = std::make_pair(std::make_pair(0U, 0U), j);
120  for (unsigned int j = 0; j < this->dofs_per_face; ++j)
121  face_system_to_base_table[j] = std::make_pair(std::make_pair(0U, 0U), j);
122 
123  // Fill with default value; may be changed by constructor of derived class.
125 
126  // initialize the restriction and prolongation matrices. the default
127  // constructor of FullMatrix<dim> initializes them with size zero
130  for (unsigned int ref = RefinementCase<dim>::cut_x;
131  ref < RefinementCase<dim>::isotropic_refinement + 1;
132  ++ref)
133  {
135  RefinementCase<dim>(ref)),
138  RefinementCase<dim>(ref)),
140  }
141 
143 }
144 
145 
146 
147 template <int dim, int spacedim>
148 std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>
149 FiniteElement<dim, spacedim>::operator^(const unsigned int multiplicity) const
150 {
151  return {this->clone(), multiplicity};
152 }
153 
154 
155 
156 template <int dim, int spacedim>
157 double
159  const Point<dim> &) const
160 {
162  return 0.;
163 }
164 
165 
166 
167 template <int dim, int spacedim>
168 double
170  const Point<dim> &,
171  const unsigned int) const
172 {
174  return 0.;
175 }
176 
177 
178 
179 template <int dim, int spacedim>
182  const Point<dim> &) const
183 {
185  return Tensor<1, dim>();
186 }
187 
188 
189 
190 template <int dim, int spacedim>
193  const Point<dim> &,
194  const unsigned int) const
195 {
197  return Tensor<1, dim>();
198 }
199 
200 
201 
202 template <int dim, int spacedim>
205  const Point<dim> &) const
206 {
208  return Tensor<2, dim>();
209 }
210 
211 
212 
213 template <int dim, int spacedim>
216  const unsigned int,
217  const Point<dim> &,
218  const unsigned int) const
219 {
221  return Tensor<2, dim>();
222 }
223 
224 
225 
226 template <int dim, int spacedim>
229  const Point<dim> &) const
230 {
232  return Tensor<3, dim>();
233 }
234 
235 
236 
237 template <int dim, int spacedim>
240  const unsigned int,
241  const Point<dim> &,
242  const unsigned int) const
243 {
245  return Tensor<3, dim>();
246 }
247 
248 
249 
250 template <int dim, int spacedim>
253  const Point<dim> &) const
254 {
256  return Tensor<4, dim>();
257 }
258 
259 
260 
261 template <int dim, int spacedim>
264  const unsigned int,
265  const Point<dim> &,
266  const unsigned int) const
267 {
269  return Tensor<4, dim>();
270 }
271 
272 
273 template <int dim, int spacedim>
274 void
276  const bool isotropic_restriction_only,
277  const bool isotropic_prolongation_only)
278 {
279  for (unsigned int ref_case = RefinementCase<dim>::cut_x;
280  ref_case <= RefinementCase<dim>::isotropic_refinement;
281  ++ref_case)
282  {
283  const unsigned int nc =
285 
286  for (unsigned int i = 0; i < nc; ++i)
287  {
288  if (this->restriction[ref_case - 1][i].m() != this->dofs_per_cell &&
289  (!isotropic_restriction_only ||
291  this->restriction[ref_case - 1][i].reinit(this->dofs_per_cell,
292  this->dofs_per_cell);
293  if (this->prolongation[ref_case - 1][i].m() != this->dofs_per_cell &&
294  (!isotropic_prolongation_only ||
296  this->prolongation[ref_case - 1][i].reinit(this->dofs_per_cell,
297  this->dofs_per_cell);
298  }
299  }
300 }
301 
302 
303 template <int dim, int spacedim>
304 const FullMatrix<double> &
306  const unsigned int child,
307  const RefinementCase<dim> &refinement_case) const
308 {
309  AssertIndexRange(refinement_case,
311  Assert(refinement_case != RefinementCase<dim>::no_refinement,
312  ExcMessage(
313  "Restriction matrices are only available for refined cells!"));
315  child, GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)));
316  // we use refinement_case-1 here. the -1 takes care of the origin of the
317  // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
318  // available and so the vector indices are shifted
319  Assert(restriction[refinement_case - 1][child].n() == this->dofs_per_cell,
321  return restriction[refinement_case - 1][child];
322 }
323 
324 
325 
326 template <int dim, int spacedim>
327 const FullMatrix<double> &
329  const unsigned int child,
330  const RefinementCase<dim> &refinement_case) const
331 {
332  AssertIndexRange(refinement_case,
334  Assert(refinement_case != RefinementCase<dim>::no_refinement,
335  ExcMessage(
336  "Prolongation matrices are only available for refined cells!"));
338  child, GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)));
339  // we use refinement_case-1 here. the -1 takes care
340  // of the origin of the vector, as for
341  // RefinementCase::no_refinement (=0) there is no
342  // data available and so the vector indices
343  // are shifted
344  Assert(prolongation[refinement_case - 1][child].n() == this->dofs_per_cell,
345  ExcEmbeddingVoid());
346  return prolongation[refinement_case - 1][child];
347 }
348 
349 
350 // TODO:[GK] This is probably not the most efficient way of doing this.
351 template <int dim, int spacedim>
352 unsigned int
354  const unsigned int index) const
355 {
356  AssertIndexRange(index, this->n_components());
357 
358  return first_block_of_base(component_to_base_table[index].first.first) +
359  component_to_base_table[index].second;
360 }
361 
362 
363 template <int dim, int spacedim>
366  const FEValuesExtractors::Scalar &scalar) const
367 {
368  AssertIndexRange(scalar.component, this->n_components());
369 
370  // TODO: it would be nice to verify that it is indeed possible
371  // to select this scalar component, i.e., that it is not part
372  // of a non-primitive element. unfortunately, there is no simple
373  // way to write such a condition...
374 
375  std::vector<bool> mask(this->n_components(), false);
376  mask[scalar.component] = true;
377  return mask;
378 }
379 
380 
381 template <int dim, int spacedim>
384  const FEValuesExtractors::Vector &vector) const
385 {
386  AssertIndexRange(vector.first_vector_component + dim - 1,
387  this->n_components());
388 
389  // TODO: it would be nice to verify that it is indeed possible
390  // to select these vector components, i.e., that they don't span
391  // beyond the beginning or end of anon-primitive element.
392  // unfortunately, there is no simple way to write such a condition...
393 
394  std::vector<bool> mask(this->n_components(), false);
395  for (unsigned int c = vector.first_vector_component;
396  c < vector.first_vector_component + dim;
397  ++c)
398  mask[c] = true;
399  return mask;
400 }
401 
402 
403 template <int dim, int spacedim>
406  const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const
407 {
410  this->n_components());
411 
412  // TODO: it would be nice to verify that it is indeed possible
413  // to select these vector components, i.e., that they don't span
414  // beyond the beginning or end of anon-primitive element.
415  // unfortunately, there is no simple way to write such a condition...
416 
417  std::vector<bool> mask(this->n_components(), false);
418  for (unsigned int c = sym_tensor.first_tensor_component;
419  c < sym_tensor.first_tensor_component +
421  ++c)
422  mask[c] = true;
423  return mask;
424 }
425 
426 
427 
428 template <int dim, int spacedim>
431 {
432  // if we get a block mask that represents all blocks, then
433  // do the same for the returned component mask
434  if (block_mask.represents_the_all_selected_mask())
435  return {};
436 
437  AssertDimension(block_mask.size(), this->n_blocks());
438 
439  std::vector<bool> component_mask(this->n_components(), false);
440  for (unsigned int c = 0; c < this->n_components(); ++c)
441  if (block_mask[component_to_block_index(c)] == true)
442  component_mask[c] = true;
443 
444  return component_mask;
445 }
446 
447 
448 
449 template <int dim, int spacedim>
450 BlockMask
452  const FEValuesExtractors::Scalar &scalar) const
453 {
454  // simply create the corresponding component mask (a simpler
455  // process) and then convert it to a block mask
456  return block_mask(component_mask(scalar));
457 }
458 
459 
460 template <int dim, int spacedim>
461 BlockMask
463  const FEValuesExtractors::Vector &vector) const
464 {
465  // simply create the corresponding component mask (a simpler
466  // process) and then convert it to a block mask
467  return block_mask(component_mask(vector));
468 }
469 
470 
471 template <int dim, int spacedim>
472 BlockMask
474  const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const
475 {
476  // simply create the corresponding component mask (a simpler
477  // process) and then convert it to a block mask
478  return block_mask(component_mask(sym_tensor));
479 }
480 
481 
482 
483 template <int dim, int spacedim>
484 BlockMask
486  const ComponentMask &component_mask) const
487 {
488  // if we get a component mask that represents all component, then
489  // do the same for the returned block mask
490  if (component_mask.represents_the_all_selected_mask())
491  return {};
492 
493  AssertDimension(component_mask.size(), this->n_components());
494 
495  // walk over all of the components
496  // of this finite element and see
497  // if we need to set the
498  // corresponding block. inside the
499  // block, walk over all the
500  // components that correspond to
501  // this block and make sure the
502  // component mask is set for all of
503  // them
504  std::vector<bool> block_mask(this->n_blocks(), false);
505  for (unsigned int c = 0; c < this->n_components();)
506  {
507  const unsigned int block = component_to_block_index(c);
508  if (component_mask[c] == true)
509  block_mask[block] = true;
510 
511  // now check all of the other
512  // components that correspond
513  // to this block
514  ++c;
515  while ((c < this->n_components()) &&
516  (component_to_block_index(c) == block))
517  {
518  Assert(component_mask[c] == block_mask[block],
519  ExcMessage(
520  "The component mask argument given to this function "
521  "is not a mask where the individual components belonging "
522  "to one block of the finite element are either all "
523  "selected or not selected. You can't call this function "
524  "with a component mask that splits blocks."));
525  ++c;
526  }
527  }
528 
529 
530  return block_mask;
531 }
532 
533 
534 
535 template <int dim, int spacedim>
536 unsigned int
538  const unsigned int face,
539  const bool face_orientation,
540  const bool face_flip,
541  const bool face_rotation) const
542 {
543  AssertIndexRange(face_index, this->dofs_per_face);
545 
546  // TODO: we could presumably solve the 3d case below using the
547  // adjust_quad_dof_index_for_face_orientation_table field. for the
548  // 2d case, we can't use adjust_line_dof_index_for_line_orientation_table
549  // since that array is empty (presumably because we thought that
550  // there are no flipped edges in 2d, but these can happen in
551  // DoFTools::make_periodicity_constraints, for example). so we
552  // would need to either fill this field, or rely on derived classes
553  // implementing this function, as we currently do
554 
555  // see the function's documentation for an explanation of this
556  // assertion -- in essence, derived classes have to implement
557  // an overloaded version of this function if we are to use any
558  // other than standard orientation
559  if ((face_orientation != true) || (face_flip != false) ||
560  (face_rotation != false))
561  Assert((this->dofs_per_line <= 1) && (this->dofs_per_quad <= 1),
562  ExcMessage(
563  "The function in this base class can not handle this case. "
564  "Rather, the derived class you are using must provide "
565  "an overloaded version but apparently hasn't done so. See "
566  "the documentation of this function for more information."));
567 
568  // we need to distinguish between DoFs on vertices, lines and in 3d quads.
569  // do so in a sequence of if-else statements
570  if (face_index < this->first_face_line_index)
571  // DoF is on a vertex
572  {
573  // get the number of the vertex on the face that corresponds to this DoF,
574  // along with the number of the DoF on this vertex
575  const unsigned int face_vertex = face_index / this->dofs_per_vertex;
576  const unsigned int dof_index_on_vertex =
577  face_index % this->dofs_per_vertex;
578 
579  // then get the number of this vertex on the cell and translate
580  // this to a DoF number on the cell
582  face, face_vertex, face_orientation, face_flip, face_rotation) *
583  this->dofs_per_vertex +
584  dof_index_on_vertex);
585  }
586  else if (face_index < this->first_face_quad_index)
587  // DoF is on a face
588  {
589  // do the same kind of translation as before. we need to only consider
590  // DoFs on the lines, i.e., ignoring those on the vertices
591  const unsigned int index = face_index - this->first_face_line_index;
592 
593  const unsigned int face_line = index / this->dofs_per_line;
594  const unsigned int dof_index_on_line = index % this->dofs_per_line;
595 
596  return (this->first_line_index +
598  face, face_line, face_orientation, face_flip, face_rotation) *
599  this->dofs_per_line +
600  dof_index_on_line);
601  }
602  else
603  // DoF is on a quad
604  {
605  Assert(dim >= 3, ExcInternalError());
606 
607  // ignore vertex and line dofs
608  const unsigned int index = face_index - this->first_face_quad_index;
609 
610  return (this->first_quad_index + face * this->dofs_per_quad + index);
611  }
612 }
613 
614 
615 
616 template <int dim, int spacedim>
617 unsigned int
619  const unsigned int index,
620  const bool face_orientation,
621  const bool face_flip,
622  const bool face_rotation) const
623 {
624  // general template for 1D and 2D: not
625  // implemented. in fact, the function
626  // shouldn't even be called unless we are
627  // in 3d, so throw an internal error
628  Assert(dim == 3, ExcInternalError());
629  if (dim < 3)
630  return index;
631 
632  // adjust dofs on 3d faces if the face is
633  // flipped. note that we query a table that
634  // derived elements need to have set up
635  // front. the exception are discontinuous
636  // elements for which there should be no
637  // face dofs anyway (i.e. dofs_per_quad==0
638  // in 3d), so we don't need the table, but
639  // the function should also not have been
640  // called
641  AssertIndexRange(index, this->dofs_per_quad);
643  8 * this->dofs_per_quad,
644  ExcInternalError());
646  index, 4 * face_orientation + 2 * face_flip + face_rotation);
647 }
648 
649 
650 
651 template <int dim, int spacedim>
652 unsigned int
654  const unsigned int index,
655  const bool line_orientation) const
656 {
657  // general template for 1D and 2D: do
658  // nothing. Do not throw an Assertion,
659  // however, in order to allow to call this
660  // function in 2D as well
661  if (dim < 3)
662  return index;
663 
664  AssertIndexRange(index, this->dofs_per_line);
666  this->dofs_per_line,
667  ExcInternalError());
668  if (line_orientation)
669  return index;
670  else
672 }
673 
674 
675 
676 template <int dim, int spacedim>
677 bool
679 {
680  for (unsigned int ref_case = RefinementCase<dim>::cut_x;
681  ref_case < RefinementCase<dim>::isotropic_refinement + 1;
682  ++ref_case)
683  for (unsigned int c = 0;
684  c < GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
685  ++c)
686  {
687  // make sure also the lazily initialized matrices are created
689  Assert((prolongation[ref_case - 1][c].m() == this->dofs_per_cell) ||
690  (prolongation[ref_case - 1][c].m() == 0),
691  ExcInternalError());
692  Assert((prolongation[ref_case - 1][c].n() == this->dofs_per_cell) ||
693  (prolongation[ref_case - 1][c].n() == 0),
694  ExcInternalError());
695  if ((prolongation[ref_case - 1][c].m() == 0) ||
696  (prolongation[ref_case - 1][c].n() == 0))
697  return false;
698  }
699  return true;
700 }
701 
702 
703 
704 template <int dim, int spacedim>
705 bool
707 {
708  for (unsigned int ref_case = RefinementCase<dim>::cut_x;
709  ref_case < RefinementCase<dim>::isotropic_refinement + 1;
710  ++ref_case)
711  for (unsigned int c = 0;
712  c < GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
713  ++c)
714  {
715  // make sure also the lazily initialized matrices are created
717  Assert((restriction[ref_case - 1][c].m() == this->dofs_per_cell) ||
718  (restriction[ref_case - 1][c].m() == 0),
719  ExcInternalError());
720  Assert((restriction[ref_case - 1][c].n() == this->dofs_per_cell) ||
721  (restriction[ref_case - 1][c].n() == 0),
722  ExcInternalError());
723  if ((restriction[ref_case - 1][c].m() == 0) ||
724  (restriction[ref_case - 1][c].n() == 0))
725  return false;
726  }
727  return true;
728 }
729 
730 
731 
732 template <int dim, int spacedim>
733 bool
735 {
736  const RefinementCase<dim> ref_case =
738 
739  for (unsigned int c = 0;
740  c < GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
741  ++c)
742  {
743  // make sure also the lazily initialized matrices are created
745  Assert((prolongation[ref_case - 1][c].m() == this->dofs_per_cell) ||
746  (prolongation[ref_case - 1][c].m() == 0),
747  ExcInternalError());
748  Assert((prolongation[ref_case - 1][c].n() == this->dofs_per_cell) ||
749  (prolongation[ref_case - 1][c].n() == 0),
750  ExcInternalError());
751  if ((prolongation[ref_case - 1][c].m() == 0) ||
752  (prolongation[ref_case - 1][c].n() == 0))
753  return false;
754  }
755  return true;
756 }
757 
758 
759 
760 template <int dim, int spacedim>
761 bool
763 {
764  const RefinementCase<dim> ref_case =
766 
767  for (unsigned int c = 0;
768  c < GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
769  ++c)
770  {
771  // make sure also the lazily initialized matrices are created
773  Assert((restriction[ref_case - 1][c].m() == this->dofs_per_cell) ||
774  (restriction[ref_case - 1][c].m() == 0),
775  ExcInternalError());
776  Assert((restriction[ref_case - 1][c].n() == this->dofs_per_cell) ||
777  (restriction[ref_case - 1][c].n() == 0),
778  ExcInternalError());
779  if ((restriction[ref_case - 1][c].m() == 0) ||
780  (restriction[ref_case - 1][c].n() == 0))
781  return false;
782  }
783  return true;
784 }
785 
786 
787 
788 template <int dim, int spacedim>
789 bool
791  const internal::SubfaceCase<dim> &subface_case) const
792 {
793  if (subface_case == internal::SubfaceCase<dim>::case_isotropic)
794  return (this->dofs_per_face == 0) || (interface_constraints.m() != 0);
795  else
796  return false;
797 }
798 
799 
800 
801 template <int dim, int spacedim>
802 bool
804 {
805  return false;
806 }
807 
808 
809 
810 template <int dim, int spacedim>
811 const FullMatrix<double> &
813  const internal::SubfaceCase<dim> &subface_case) const
814 {
815  (void)subface_case;
817  ExcMessage("Constraints for this element are only implemented "
818  "for the case that faces are refined isotropically "
819  "(which is always the case in 2d, and in 3d requires "
820  "that the neighboring cell of a coarse cell presents "
821  "exactly four children on the common face)."));
822  Assert((this->dofs_per_face == 0) || (interface_constraints.m() != 0),
823  ExcMessage("The finite element for which you try to obtain "
824  "hanging node constraints does not appear to "
825  "implement them."));
826 
827  if (dim == 1)
831 
832  return interface_constraints;
833 }
834 
835 
836 
837 template <int dim, int spacedim>
840 {
841  switch (dim)
842  {
843  case 1:
844  return {0U, 0U};
845  case 2:
846  return {this->dofs_per_vertex + 2 * this->dofs_per_line,
847  this->dofs_per_face};
848  case 3:
849  return {5 * this->dofs_per_vertex + 12 * this->dofs_per_line +
850  4 * this->dofs_per_quad,
851  this->dofs_per_face};
852  default:
853  Assert(false, ExcNotImplemented());
854  }
856 }
857 
858 
859 
860 template <int dim, int spacedim>
861 void
864  FullMatrix<double> &) const
865 {
866  // by default, no interpolation
867  // implemented. so throw exception,
868  // as documentation says
869  AssertThrow(
870  false,
872 }
873 
874 
875 
876 template <int dim, int spacedim>
877 void
880  FullMatrix<double> &) const
881 {
882  // by default, no interpolation
883  // implemented. so throw exception,
884  // as documentation says
885  AssertThrow(
886  false,
888 }
889 
890 
891 
892 template <int dim, int spacedim>
893 void
896  const unsigned int,
897  FullMatrix<double> &) const
898 {
899  // by default, no interpolation
900  // implemented. so throw exception,
901  // as documentation says
902  AssertThrow(
903  false,
905 }
906 
907 
908 
909 template <int dim, int spacedim>
910 std::vector<std::pair<unsigned int, unsigned int>>
912  const FiniteElement<dim, spacedim> &) const
913 {
914  Assert(false, ExcNotImplemented());
915  return std::vector<std::pair<unsigned int, unsigned int>>();
916 }
917 
918 
919 
920 template <int dim, int spacedim>
921 std::vector<std::pair<unsigned int, unsigned int>>
923  const FiniteElement<dim, spacedim> &) const
924 {
925  Assert(false, ExcNotImplemented());
926  return std::vector<std::pair<unsigned int, unsigned int>>();
927 }
928 
929 
930 
931 template <int dim, int spacedim>
932 std::vector<std::pair<unsigned int, unsigned int>>
934  const FiniteElement<dim, spacedim> &) const
935 {
936  Assert(false, ExcNotImplemented());
937  return std::vector<std::pair<unsigned int, unsigned int>>();
938 }
939 
940 
941 
942 template <int dim, int spacedim>
945  const FiniteElement<dim, spacedim> &fe_other) const
946 {
947  return this->compare_for_domination(fe_other, 1);
948 }
949 
950 
951 
952 template <int dim, int spacedim>
956  const unsigned int) const
957 {
958  Assert(false, ExcNotImplemented());
960 }
961 
962 
963 
964 template <int dim, int spacedim>
965 bool
968 {
969  // Compare fields in roughly increasing order of how expensive the
970  // comparison is
971  return ((typeid(*this) == typeid(f)) && (this->get_name() == f.get_name()) &&
972  (static_cast<const FiniteElementData<dim> &>(*this) ==
973  static_cast<const FiniteElementData<dim> &>(f)) &&
975 }
976 
977 
978 
979 template <int dim, int spacedim>
980 bool
983 {
984  return !(*this == f);
985 }
986 
987 
988 
989 template <int dim, int spacedim>
990 const std::vector<Point<dim>> &
992 {
993  // a finite element may define
994  // support points, but only if
995  // there are as many as there are
996  // degrees of freedom
997  Assert((unit_support_points.size() == 0) ||
998  (unit_support_points.size() == this->dofs_per_cell),
999  ExcInternalError());
1000  return unit_support_points;
1001 }
1002 
1003 
1004 
1005 template <int dim, int spacedim>
1006 bool
1008 {
1009  return (unit_support_points.size() != 0);
1010 }
1011 
1012 
1013 
1014 template <int dim, int spacedim>
1015 const std::vector<Point<dim>> &
1017 {
1018  // If the finite element implements generalized support points, return
1019  // those. Otherwise fall back to unit support points.
1020  return ((generalized_support_points.size() == 0) ?
1023 }
1024 
1025 
1026 
1027 template <int dim, int spacedim>
1028 bool
1030 {
1031  return (get_generalized_support_points().size() != 0);
1032 }
1033 
1034 
1035 
1036 template <int dim, int spacedim>
1037 Point<dim>
1039 {
1040  AssertIndexRange(index, this->dofs_per_cell);
1041  Assert(unit_support_points.size() == this->dofs_per_cell,
1043  return unit_support_points[index];
1044 }
1045 
1046 
1047 
1048 template <int dim, int spacedim>
1049 const std::vector<Point<dim - 1>> &
1051 {
1052  // a finite element may define
1053  // support points, but only if
1054  // there are as many as there are
1055  // degrees of freedom on a face
1056  Assert((unit_face_support_points.size() == 0) ||
1057  (unit_face_support_points.size() == this->dofs_per_face),
1058  ExcInternalError());
1059  return unit_face_support_points;
1060 }
1061 
1062 
1063 
1064 template <int dim, int spacedim>
1065 bool
1067 {
1068  return (unit_face_support_points.size() != 0);
1069 }
1070 
1071 
1072 
1073 template <int dim, int spacedim>
1074 Point<dim - 1>
1076  const unsigned int index) const
1077 {
1078  AssertIndexRange(index, this->dofs_per_face);
1081  return unit_face_support_points[index];
1082 }
1083 
1084 
1085 
1086 template <int dim, int spacedim>
1087 bool
1089  const unsigned int) const
1090 {
1091  return true;
1092 }
1093 
1094 
1095 
1096 template <int dim, int spacedim>
1099 {
1100  // Translate the ComponentMask into first_selected and n_components after
1101  // some error checking:
1102  const unsigned int n_total_components = this->n_components();
1103  Assert((n_total_components == mask.size()) || (mask.size() == 0),
1104  ExcMessage("The given ComponentMask has the wrong size."));
1105 
1106  const unsigned int n_selected =
1107  mask.n_selected_components(n_total_components);
1108  Assert(n_selected > 0,
1109  ExcMessage("You need at least one selected component."));
1110 
1111  const unsigned int first_selected =
1112  mask.first_selected_component(n_total_components);
1113 
1114 #ifdef DEBUG
1115  // check that it is contiguous:
1116  for (unsigned int c = 0; c < n_total_components; ++c)
1117  Assert((c < first_selected && (!mask[c])) ||
1118  (c >= first_selected && c < first_selected + n_selected &&
1119  mask[c]) ||
1120  (c >= first_selected + n_selected && !mask[c]),
1121  ExcMessage("Error: the given ComponentMask is not contiguous!"));
1122 #endif
1123 
1124  return get_sub_fe(first_selected, n_selected);
1125 }
1126 
1127 
1128 
1129 template <int dim, int spacedim>
1132  const unsigned int first_component,
1133  const unsigned int n_selected_components) const
1134 {
1135  // No complicated logic is needed here, because it is overridden in
1136  // FESystem<dim,spacedim>. Just make sure that what the user chose is valid:
1137  Assert(first_component == 0 && n_selected_components == this->n_components(),
1138  ExcMessage(
1139  "You can only select a whole FiniteElement, not a part of one."));
1140 
1141  (void)first_component;
1142  (void)n_selected_components;
1143  return *this;
1144 }
1145 
1146 
1147 
1148 template <int dim, int spacedim>
1149 std::pair<Table<2, bool>, std::vector<unsigned int>>
1151 {
1152  Assert(false, ExcNotImplemented());
1153  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
1154  Table<2, bool>(this->n_components(), this->dofs_per_cell),
1155  std::vector<unsigned int>(this->n_components()));
1156 }
1157 
1158 
1159 
1160 template <int dim, int spacedim>
1161 void
1164  const std::vector<Vector<double>> &,
1165  std::vector<double> &) const
1166 {
1168  ExcMessage("The element for which you are calling the current "
1169  "function does not have generalized support points (see "
1170  "the glossary for a definition of generalized support "
1171  "points). Consequently, the current function can not "
1172  "be defined and is not implemented by the element."));
1173  Assert(false, ExcNotImplemented());
1174 }
1175 
1176 
1177 
1178 template <int dim, int spacedim>
1179 std::size_t
1181 {
1182  return (
1183  sizeof(FiniteElementData<dim>) +
1195 }
1196 
1197 
1198 
1199 template <int dim, int spacedim>
1200 std::vector<unsigned int>
1202  const std::vector<ComponentMask> &nonzero_components)
1203 {
1204  std::vector<unsigned int> retval(nonzero_components.size());
1205  for (unsigned int i = 0; i < nonzero_components.size(); ++i)
1206  retval[i] = nonzero_components[i].n_selected_components();
1207  return retval;
1208 }
1209 
1210 
1211 
1212 /*------------------------------- FiniteElement ----------------------*/
1213 
1214 template <int dim, int spacedim>
1215 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1217  const UpdateFlags flags,
1218  const Mapping<dim, spacedim> &mapping,
1219  const Quadrature<dim - 1> & quadrature,
1221  spacedim>
1222  &output_data) const
1223 {
1224  return get_data(flags,
1225  mapping,
1227  output_data);
1228 }
1229 
1230 
1231 
1232 template <int dim, int spacedim>
1233 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1235  const UpdateFlags flags,
1236  const Mapping<dim, spacedim> &mapping,
1237  const Quadrature<dim - 1> & quadrature,
1239  spacedim>
1240  &output_data) const
1241 {
1242  return get_data(flags,
1243  mapping,
1245  output_data);
1246 }
1247 
1248 
1249 
1250 template <int dim, int spacedim>
1252 FiniteElement<dim, spacedim>::base_element(const unsigned int index) const
1253 {
1254  (void)index;
1255  AssertIndexRange(index, 1);
1256  // This function should not be
1257  // called for a system element
1259  return *this;
1260 }
1261 
1262 
1263 
1264 /*------------------------------- Explicit Instantiations -------------*/
1265 #include "fe.inst"
1266 
1267 
static ::ExceptionBase & ExcFEHasNoSupportPoints()
bool represents_the_all_selected_mask() const
size_type m() const
static const unsigned int invalid_unsigned_int
Definition: types.h:191
BlockMask block_mask(const FEValuesExtractors::Scalar &scalar) const
Definition: fe.cc:451
virtual Tensor< 2, dim > shape_grad_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
Definition: fe.cc:215
std::vector< std::vector< FullMatrix< double > > > restriction
Definition: fe.h:2414
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
virtual std::size_t memory_consumption() const
Definition: fe.cc:49
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const
Definition: fe.cc:878
virtual bool operator==(const FiniteElement< dim, spacedim > &fe) const
Definition: fe.cc:967
virtual Point< dim - 1 > unit_face_support_point(const unsigned int index) const
Definition: fe.cc:1075
std::vector< Point< dim > > generalized_support_points
Definition: fe.h:2465
virtual Tensor< 3, dim > shape_3rd_derivative(const unsigned int i, const Point< dim > &p) const
Definition: fe.cc:228
FullMatrix< double > interface_constraints
Definition: fe.h:2440
bool operator!=(const FiniteElement< dim, spacedim > &) const
Definition: fe.cc:982
FiniteElement(const FiniteElementData< dim > &fe_data, const std::vector< bool > &restriction_is_additive_flags, const std::vector< ComponentMask > &nonzero_components)
Definition: fe.cc:57
const FiniteElement< dim, spacedim > & get_sub_fe(const ComponentMask &mask) const
Definition: fe.cc:1098
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
Definition: fe.cc:169
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const
Definition: fe.cc:862
static ::ExceptionBase & ExcUnitShapeValuesDoNotExist()
virtual Point< dim > unit_support_point(const unsigned int index) const
Definition: fe.cc:1038
const std::vector< Point< dim - 1 > > & get_unit_face_support_points() const
Definition: fe.cc:1050
bool restriction_is_implemented() const
Definition: fe.cc:706
const std::vector< Point< dim > > & get_generalized_support_points() const
Definition: fe.cc:1016
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const
Definition: fe.cc:181
virtual Tensor< 3, dim > shape_3rd_derivative_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
Definition: fe.cc:239
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
const unsigned int dofs_per_quad
Definition: fe_base.h:237
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:933
bool isotropic_prolongation_is_implemented() const
Definition: fe.cc:734
std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > component_to_base_table
Definition: fe.h:2574
bool has_generalized_support_points() const
Definition: fe.cc:1029
bool represents_the_all_selected_mask() const
std::vector< Point< dim - 1 > > unit_face_support_points
Definition: fe.h:2459
STL namespace.
ComponentMask component_mask(const FEValuesExtractors::Scalar &scalar) const
Definition: fe.cc:365
static const char U
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const =0
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
static std::vector< unsigned int > compute_n_nonzero_components(const std::vector< ComponentMask > &nonzero_components)
Definition: fe.cc:1201
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
Definition: fe.cc:328
bool is_primitive() const
Definition: fe.h:3273
const std::vector< unsigned int > n_nonzero_components_table
Definition: fe.h:2599
const unsigned int dofs_per_line
Definition: fe_base.h:231
size_type n_elements() const
unsigned int n_blocks() const
bool has_face_support_points() const
Definition: fe.cc:1066
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2452
const unsigned int first_face_line_index
Definition: fe_base.h:263
bool has_support_points() const
Definition: fe.cc:1007
virtual std::unique_ptr< InternalDataBase > get_subface_data(const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const Quadrature< dim - 1 > &quadrature, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const
Definition: fe.cc:1234
unsigned int component_to_block_index(const unsigned int component) const
Definition: fe.cc:353
static ::ExceptionBase & ExcMessage(std::string arg1)
const FullMatrix< double > & constraints(const ::internal::SubfaceCase< dim > &subface_case=::internal::SubfaceCase< dim >::case_isotropic) const
Definition: fe.cc:812
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
Definition: fe.cc:192
const unsigned int first_quad_index
Definition: fe_base.h:253
virtual std::size_t memory_consumption() const
Definition: fe.cc:1180
size_type n() const
No update.
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2428
#define Assert(cond, exc)
Definition: exceptions.h:1419
UpdateFlags
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
Definition: fe.cc:305
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
virtual FiniteElementDomination::Domination compare_for_face_domination(const FiniteElement< dim, spacedim > &fe_other) const final
Definition: fe.cc:944
Abstract base class for mapping classes.
Definition: mapping.h:302
void reinit(const unsigned int n_blocks, const size_type n_elements_per_block)
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix) const
Definition: fe.cc:894
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
bool constraints_are_implemented(const ::internal::SubfaceCase< dim > &subface_case=::internal::SubfaceCase< dim >::case_isotropic) const
Definition: fe.cc:790
static ::ExceptionBase & ExcWrongInterfaceMatrixSize(int arg1, int arg2)
VectorType::value_type * end(VectorType &V)
unsigned int size() const
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
Definition: fe.cc:1252
virtual std::string get_name() const =0
static ::ExceptionBase & ExcProjectionVoid()
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:911
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
virtual void convert_generalized_support_point_values_to_dof_values(const std::vector< Vector< double >> &support_point_values, std::vector< double > &nodal_values) const
Definition: fe.cc:1163
const unsigned int dofs_per_cell
Definition: fe_base.h:282
Point< 2 > first
Definition: grid_out.cc:4352
types::global_dof_index first_block_of_base(const unsigned int b) const
Definition: fe.h:3193
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
unsigned int adjust_line_dof_index_for_line_orientation(const unsigned int index, const bool line_orientation) const
Definition: fe.cc:653
const std::vector< Point< dim > > & get_unit_support_points() const
Definition: fe.cc:991
virtual std::unique_ptr< InternalDataBase > get_face_data(const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const Quadrature< dim - 1 > &quadrature, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const
Definition: fe.cc:1216
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:922
Table< 2, int > adjust_quad_dof_index_for_face_orientation_table
Definition: fe.h:2488
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const
Definition: fe.cc:954
static ::ExceptionBase & ExcEmbeddingVoid()
virtual double shape_value(const unsigned int i, const Point< dim > &p) const
Definition: fe.cc:158
unsigned int n_components() const
const unsigned int first_face_quad_index
Definition: fe_base.h:268
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
VectorType::value_type * begin(VectorType &V)
std::vector< std::pair< unsigned int, unsigned int > > system_to_component_table
Definition: fe.h:2508
virtual bool hp_constraints_are_implemented() const
Definition: fe.cc:803
unsigned int first_selected_component(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const unsigned int dofs_per_face
Definition: fe_base.h:275
bool prolongation_is_implemented() const
Definition: fe.cc:678
const unsigned int first_line_index
Definition: fe_base.h:248
virtual Tensor< 4, dim > shape_4th_derivative_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
Definition: fe.cc:263
virtual Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const
Definition: fe.cc:204
const std::vector< ComponentMask > nonzero_components
Definition: fe.h:2590
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
Definition: fe.cc:275
std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > system_to_base_table
Definition: fe.h:2539
static ::ExceptionBase & ExcNotImplemented()
const unsigned int dofs_per_vertex
Definition: fe_base.h:225
Definition: table.h:699
std::vector< std::pair< unsigned int, unsigned int > > face_system_to_component_table
Definition: fe.h:2520
bool isotropic_restriction_is_implemented() const
Definition: fe.cc:762
BlockIndices base_to_block_indices
Definition: fe.h:2551
std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > face_system_to_base_table
Definition: fe.h:2545
virtual Tensor< 4, dim > shape_4th_derivative(const unsigned int i, const Point< dim > &p) const
Definition: fe.cc:252
unsigned int size() const
unsigned int size() const
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
unsigned int adjust_quad_dof_index_for_face_orientation(const unsigned int index, const bool face_orientation, const bool face_flip, const bool face_rotation) const
Definition: fe.cc:618
std::vector< int > adjust_line_dof_index_for_line_orientation_table
Definition: fe.h:2503
TableIndices< 2 > interface_constraints_size() const
Definition: fe.cc:839
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const
Definition: fe.cc:1150
void fill(InputIterator entries, const bool C_style_indexing=true)
const std::vector< bool > restriction_is_additive_flags
Definition: fe.h:2581
std::pair< std::unique_ptr< FiniteElement< dim, spacedim > >, unsigned int > operator^(const unsigned int multiplicity) const
Definition: fe.cc:149
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const
Definition: fe.cc:1088
virtual std::unique_ptr< InternalDataBase > get_data(const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const Quadrature< dim > &quadrature, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const =0
static ::ExceptionBase & ExcInternalError()