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