Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.3.3
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
fe.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1998 - 2021 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>
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
40template <int dim, int spacedim>
42 : update_each(update_default)
43{}
44
45
46
47template <int dim, int spacedim>
48std::size_t
50{
51 return sizeof(*this);
52}
53
54
55
56template <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_line_dof_index_for_line_orientation_table(
63 dim == 3 ? this->n_dofs_per_line() : 0)
64 , system_to_base_table(this->n_dofs_per_cell())
65 , component_to_base_table(this->components,
66 std::make_pair(std::make_pair(0U, 0U), 0U))
67 ,
68
69 // Special handling of vectors of length one: in this case, we
70 // assume that all entries were supposed to be equal
71 restriction_is_additive_flags(
72 r_i_a_f.size() == 1 ?
73 std::vector<bool>(fe_data.n_dofs_per_cell(), r_i_a_f[0]) :
74 r_i_a_f)
75 , nonzero_components(
76 nonzero_c.size() == 1 ?
77 std::vector<ComponentMask>(fe_data.n_dofs_per_cell(), nonzero_c[0]) :
78 nonzero_c)
79 , n_nonzero_components_table(compute_n_nonzero_components(nonzero_components))
80 , cached_primitivity(std::find_if(n_nonzero_components_table.begin(),
81 n_nonzero_components_table.end(),
82 [](const unsigned int n_components) {
83 return n_components != 1U;
84 }) == n_nonzero_components_table.end())
85{
86 Assert(restriction_is_additive_flags.size() == this->n_dofs_per_cell(),
87 ExcDimensionMismatch(restriction_is_additive_flags.size(),
88 this->n_dofs_per_cell()));
89 AssertDimension(nonzero_components.size(), this->n_dofs_per_cell());
90 for (unsigned int i = 0; i < nonzero_components.size(); ++i)
91 {
92 Assert(nonzero_components[i].size() == this->n_components(),
94 Assert(nonzero_components[i].n_selected_components() >= 1,
96 Assert(n_nonzero_components_table[i] >= 1, ExcInternalError());
97 Assert(n_nonzero_components_table[i] <= this->n_components(),
99 }
100
101 // initialize some tables in the default way, i.e. if there is only one
102 // (vector-)component; if the element is not primitive, leave these tables
103 // empty.
104 if (this->is_primitive())
105 {
106 system_to_component_table.resize(this->n_dofs_per_cell());
107 for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
108 system_to_component_table[j] = std::pair<unsigned, unsigned>(0, j);
109
110 face_system_to_component_table.resize(this->n_unique_faces());
111 for (unsigned int f = 0; f < this->n_unique_faces(); ++f)
112 {
113 face_system_to_component_table[f].resize(this->n_dofs_per_face(f));
114 for (unsigned int j = 0; j < this->n_dofs_per_face(f); ++j)
115 face_system_to_component_table[f][j] =
116 std::pair<unsigned, unsigned>(0, j);
117 }
118 }
119
120 for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
121 system_to_base_table[j] = std::make_pair(std::make_pair(0U, 0U), j);
122
123 face_system_to_base_table.resize(this->n_unique_faces());
124 for (unsigned int f = 0; f < this->n_unique_faces(); ++f)
125 {
126 face_system_to_base_table[f].resize(this->n_dofs_per_face(f));
127 for (unsigned int j = 0; j < this->n_dofs_per_face(f); ++j)
128 face_system_to_base_table[f][j] =
129 std::make_pair(std::make_pair(0U, 0U), j);
130 }
131
132 // Fill with default value; may be changed by constructor of derived class.
133 base_to_block_indices.reinit(1, 1);
134
135 // initialize the restriction and prolongation matrices. the default
136 // constructor of FullMatrix<dim> initializes them with size zero
139 for (unsigned int ref = RefinementCase<dim>::cut_x;
140 ref < RefinementCase<dim>::isotropic_refinement + 1;
141 ++ref)
142 {
143 prolongation[ref - 1].resize(GeometryInfo<dim>::n_children(
146 restriction[ref - 1].resize(GeometryInfo<dim>::n_children(
149 }
150
151
152 if (dim == 3)
153 {
154 adjust_quad_dof_index_for_face_orientation_table.resize(
155 this->n_unique_quads());
156
157 for (unsigned int f = 0; f < this->n_unique_quads(); ++f)
158 {
159 adjust_quad_dof_index_for_face_orientation_table[f] =
160 Table<2, int>(this->n_dofs_per_quad(f),
161 this->reference_cell().face_reference_cell(f) ==
163 8 :
164 6);
165 adjust_quad_dof_index_for_face_orientation_table[f].fill(0);
166 }
167 }
168
169 unit_face_support_points.resize(this->n_unique_faces());
170 generalized_face_support_points.resize(this->n_unique_faces());
171}
172
173
174
175template <int dim, int spacedim>
176std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>
177FiniteElement<dim, spacedim>::operator^(const unsigned int multiplicity) const
178{
179 return {this->clone(), multiplicity};
180}
181
182
183
184template <int dim, int spacedim>
185double
187 const Point<dim> &) const
188{
190 return 0.;
191}
192
193
194
195template <int dim, int spacedim>
196double
198 const Point<dim> &,
199 const unsigned int) const
200{
202 return 0.;
203}
204
205
206
207template <int dim, int spacedim>
210 const Point<dim> &) const
211{
213 return Tensor<1, dim>();
214}
215
216
217
218template <int dim, int spacedim>
221 const Point<dim> &,
222 const unsigned int) const
223{
225 return Tensor<1, dim>();
226}
227
228
229
230template <int dim, int spacedim>
233 const Point<dim> &) const
234{
236 return Tensor<2, dim>();
237}
238
239
240
241template <int dim, int spacedim>
244 const unsigned int,
245 const Point<dim> &,
246 const unsigned int) const
247{
249 return Tensor<2, dim>();
250}
251
252
253
254template <int dim, int spacedim>
257 const Point<dim> &) const
258{
260 return Tensor<3, dim>();
261}
262
263
264
265template <int dim, int spacedim>
268 const unsigned int,
269 const Point<dim> &,
270 const unsigned int) const
271{
273 return Tensor<3, dim>();
274}
275
276
277
278template <int dim, int spacedim>
281 const Point<dim> &) const
282{
284 return Tensor<4, dim>();
285}
286
287
288
289template <int dim, int spacedim>
292 const unsigned int,
293 const Point<dim> &,
294 const unsigned int) const
295{
297 return Tensor<4, dim>();
298}
299
300
301template <int dim, int spacedim>
302void
304 const bool isotropic_restriction_only,
305 const bool isotropic_prolongation_only)
306{
307 for (unsigned int ref_case = RefinementCase<dim>::cut_x;
308 ref_case <= RefinementCase<dim>::isotropic_refinement;
309 ++ref_case)
310 {
311 const unsigned int nc =
313
314 for (unsigned int i = 0; i < nc; ++i)
315 {
316 if (this->restriction[ref_case - 1][i].m() !=
317 this->n_dofs_per_cell() &&
318 (!isotropic_restriction_only ||
320 this->restriction[ref_case - 1][i].reinit(this->n_dofs_per_cell(),
321 this->n_dofs_per_cell());
322 if (this->prolongation[ref_case - 1][i].m() !=
323 this->n_dofs_per_cell() &&
324 (!isotropic_prolongation_only ||
326 this->prolongation[ref_case - 1][i].reinit(this->n_dofs_per_cell(),
327 this->n_dofs_per_cell());
328 }
329 }
330}
331
332
333template <int dim, int spacedim>
334const FullMatrix<double> &
336 const unsigned int child,
337 const RefinementCase<dim> &refinement_case) const
338{
339 AssertIndexRange(refinement_case,
343 "Restriction matrices are only available for refined cells!"));
345 child, GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)));
346 // we use refinement_case-1 here. the -1 takes care of the origin of the
347 // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
348 // available and so the vector indices are shifted
349 Assert(restriction[refinement_case - 1][child].n() == this->n_dofs_per_cell(),
351 return restriction[refinement_case - 1][child];
352}
353
354
355
356template <int dim, int spacedim>
357const FullMatrix<double> &
359 const unsigned int child,
360 const RefinementCase<dim> &refinement_case) const
361{
362 AssertIndexRange(refinement_case,
366 "Prolongation matrices are only available for refined cells!"));
368 child, GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)));
369 // we use refinement_case-1 here. the -1 takes care
370 // of the origin of the vector, as for
371 // RefinementCase::no_refinement (=0) there is no
372 // data available and so the vector indices
373 // are shifted
374 Assert(prolongation[refinement_case - 1][child].n() ==
375 this->n_dofs_per_cell(),
377 return prolongation[refinement_case - 1][child];
378}
379
380
381// TODO:[GK] This is probably not the most efficient way of doing this.
382template <int dim, int spacedim>
383unsigned int
385 const unsigned int index) const
386{
387 AssertIndexRange(index, this->n_components());
388
390 component_to_base_table[index].second;
391}
392
393
394template <int dim, int spacedim>
397 const FEValuesExtractors::Scalar &scalar) const
398{
399 AssertIndexRange(scalar.component, this->n_components());
400
401 // TODO: it would be nice to verify that it is indeed possible
402 // to select this scalar component, i.e., that it is not part
403 // of a non-primitive element. unfortunately, there is no simple
404 // way to write such a condition...
405
406 std::vector<bool> mask(this->n_components(), false);
407 mask[scalar.component] = true;
408 return mask;
409}
410
411
412template <int dim, int spacedim>
415 const FEValuesExtractors::Vector &vector) const
416{
418 this->n_components());
419
420 // TODO: it would be nice to verify that it is indeed possible
421 // to select these vector components, i.e., that they don't span
422 // beyond the beginning or end of anon-primitive element.
423 // unfortunately, there is no simple way to write such a condition...
424
425 std::vector<bool> mask(this->n_components(), false);
426 for (unsigned int c = vector.first_vector_component;
427 c < vector.first_vector_component + dim;
428 ++c)
429 mask[c] = true;
430 return mask;
431}
432
433
434template <int dim, int spacedim>
437 const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const
438{
441 this->n_components());
442
443 // TODO: it would be nice to verify that it is indeed possible
444 // to select these vector components, i.e., that they don't span
445 // beyond the beginning or end of anon-primitive element.
446 // unfortunately, there is no simple way to write such a condition...
447
448 std::vector<bool> mask(this->n_components(), false);
449 for (unsigned int c = sym_tensor.first_tensor_component;
450 c < sym_tensor.first_tensor_component +
452 ++c)
453 mask[c] = true;
454 return mask;
455}
456
457
458
459template <int dim, int spacedim>
462{
463 // if we get a block mask that represents all blocks, then
464 // do the same for the returned component mask
466 return {};
467
468 AssertDimension(block_mask.size(), this->n_blocks());
469
470 std::vector<bool> component_mask(this->n_components(), false);
471 for (unsigned int c = 0; c < this->n_components(); ++c)
472 if (block_mask[component_to_block_index(c)] == true)
473 component_mask[c] = true;
474
475 return component_mask;
476}
477
478
479
480template <int dim, int spacedim>
483 const FEValuesExtractors::Scalar &scalar) const
484{
485 // simply create the corresponding component mask (a simpler
486 // process) and then convert it to a block mask
487 return block_mask(component_mask(scalar));
488}
489
490
491template <int dim, int spacedim>
494 const FEValuesExtractors::Vector &vector) const
495{
496 // simply create the corresponding component mask (a simpler
497 // process) and then convert it to a block mask
498 return block_mask(component_mask(vector));
499}
500
501
502template <int dim, int spacedim>
505 const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const
506{
507 // simply create the corresponding component mask (a simpler
508 // process) and then convert it to a block mask
509 return block_mask(component_mask(sym_tensor));
510}
511
512
513
514template <int dim, int spacedim>
517 const ComponentMask &component_mask) const
518{
519 // if we get a component mask that represents all component, then
520 // do the same for the returned block mask
522 return {};
523
524 AssertDimension(component_mask.size(), this->n_components());
525
526 // walk over all of the components
527 // of this finite element and see
528 // if we need to set the
529 // corresponding block. inside the
530 // block, walk over all the
531 // components that correspond to
532 // this block and make sure the
533 // component mask is set for all of
534 // them
535 std::vector<bool> block_mask(this->n_blocks(), false);
536 for (unsigned int c = 0; c < this->n_components();)
537 {
538 const unsigned int block = component_to_block_index(c);
539 if (component_mask[c] == true)
540 block_mask[block] = true;
541
542 // now check all of the other
543 // components that correspond
544 // to this block
545 ++c;
546 while ((c < this->n_components()) &&
547 (component_to_block_index(c) == block))
548 {
549 Assert(component_mask[c] == block_mask[block],
551 "The component mask argument given to this function "
552 "is not a mask where the individual components belonging "
553 "to one block of the finite element are either all "
554 "selected or not selected. You can't call this function "
555 "with a component mask that splits blocks."));
556 ++c;
557 }
558 }
559
560
561 return block_mask;
562}
563
564
565
566template <int dim, int spacedim>
567unsigned int
568FiniteElement<dim, spacedim>::face_to_cell_index(const unsigned int face_index,
569 const unsigned int face,
570 const bool face_orientation,
571 const bool face_flip,
572 const bool face_rotation) const
573{
574 AssertIndexRange(face_index, this->n_dofs_per_face(face));
575 AssertIndexRange(face, this->reference_cell().n_faces());
576
577 // TODO: we could presumably solve the 3d case below using the
578 // adjust_quad_dof_index_for_face_orientation_table field. for the
579 // 2d case, we can't use adjust_line_dof_index_for_line_orientation_table
580 // since that array is empty (presumably because we thought that
581 // there are no flipped edges in 2d, but these can happen in
582 // DoFTools::make_periodicity_constraints, for example). so we
583 // would need to either fill this field, or rely on derived classes
584 // implementing this function, as we currently do
585
586 // see the function's documentation for an explanation of this
587 // assertion -- in essence, derived classes have to implement
588 // an overloaded version of this function if we are to use any
589 // other than standard orientation
590 if ((face_orientation != true) || (face_flip != false) ||
591 (face_rotation != false))
592 Assert((this->n_dofs_per_line() <= 1) && (this->n_dofs_per_quad(face) <= 1),
594 "The function in this base class can not handle this case. "
595 "Rather, the derived class you are using must provide "
596 "an overloaded version but apparently hasn't done so. See "
597 "the documentation of this function for more information."));
598
599 // we need to distinguish between DoFs on vertices, lines and in 3d quads.
600 // do so in a sequence of if-else statements
601 if (face_index < this->get_first_face_line_index(face))
602 // DoF is on a vertex
603 {
604 // get the number of the vertex on the face that corresponds to this DoF,
605 // along with the number of the DoF on this vertex
606 const unsigned int face_vertex = face_index / this->n_dofs_per_vertex();
607 const unsigned int dof_index_on_vertex =
608 face_index % this->n_dofs_per_vertex();
609
610 // then get the number of this vertex on the cell and translate
611 // this to a DoF number on the cell
612 return (this->reference_cell().face_to_cell_vertices(face,
613 face_vertex,
614 face_orientation +
615 2 * face_rotation +
616 4 * face_flip) *
617 this->n_dofs_per_vertex() +
618 dof_index_on_vertex);
619 }
620 else if (face_index < this->get_first_face_quad_index(face))
621 // DoF is on a face
622 {
623 // do the same kind of translation as before. we need to only consider
624 // DoFs on the lines, i.e., ignoring those on the vertices
625 const unsigned int index =
626 face_index - this->get_first_face_line_index(face);
627
628 const unsigned int face_line = index / this->n_dofs_per_line();
629 const unsigned int dof_index_on_line = index % this->n_dofs_per_line();
630
631 return (this->get_first_line_index() +
632 this->reference_cell().face_to_cell_lines(face,
633 face_line,
634 face_orientation +
635 2 * face_rotation +
636 4 * face_flip) *
637 this->n_dofs_per_line() +
638 dof_index_on_line);
639 }
640 else
641 // DoF is on a quad
642 {
643 Assert(dim >= 3, ExcInternalError());
644
645 // ignore vertex and line dofs
646 const unsigned int index =
647 face_index - this->get_first_face_quad_index(face);
648
649 return (this->get_first_quad_index(face) + index);
650 }
651}
652
653
654
655template <int dim, int spacedim>
656unsigned int
658 const unsigned int index,
659 const unsigned int face,
660 const bool face_orientation,
661 const bool face_flip,
662 const bool face_rotation) const
663{
664 // general template for 1D and 2D: not
665 // implemented. in fact, the function
666 // shouldn't even be called unless we are
667 // in 3d, so throw an internal error
668 Assert(dim == 3, ExcInternalError());
669 if (dim < 3)
670 return index;
671
672 // adjust dofs on 3d faces if the face is
673 // flipped. note that we query a table that
674 // derived elements need to have set up
675 // front. the exception are discontinuous
676 // elements for which there should be no
677 // face dofs anyway (i.e. dofs_per_quad==0
678 // in 3d), so we don't need the table, but
679 // the function should also not have been
680 // called
681 AssertIndexRange(index, this->n_dofs_per_quad(face));
683 [this->n_unique_quads() == 1 ? 0 : face]
684 .n_elements() == (this->reference_cell().face_reference_cell(
686 8 :
687 6) *
688 this->n_dofs_per_quad(face),
690 return index +
692 [this->n_unique_quads() == 1 ? 0 : face](
693 index, 4 * face_orientation + 2 * face_flip + face_rotation);
694}
695
696
697
698template <int dim, int spacedim>
699unsigned int
701 const unsigned int index,
702 const bool line_orientation) const
703{
704 // general template for 1D and 2D: do
705 // nothing. Do not throw an Assertion,
706 // however, in order to allow to call this
707 // function in 2D as well
708 if (dim < 3)
709 return index;
710
711 AssertIndexRange(index, this->n_dofs_per_line());
713 this->n_dofs_per_line(),
715 if (line_orientation)
716 return index;
717 else
719}
720
721
722
723template <int dim, int spacedim>
724bool
726{
727 for (unsigned int ref_case = RefinementCase<dim>::cut_x;
728 ref_case < RefinementCase<dim>::isotropic_refinement + 1;
729 ++ref_case)
730 for (unsigned int c = 0;
731 c < GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
732 ++c)
733 {
734 // make sure also the lazily initialized matrices are created
736 Assert((prolongation[ref_case - 1][c].m() == this->n_dofs_per_cell()) ||
737 (prolongation[ref_case - 1][c].m() == 0),
739 Assert((prolongation[ref_case - 1][c].n() == this->n_dofs_per_cell()) ||
740 (prolongation[ref_case - 1][c].n() == 0),
742 if ((prolongation[ref_case - 1][c].m() == 0) ||
743 (prolongation[ref_case - 1][c].n() == 0))
744 return false;
745 }
746 return true;
747}
748
749
750
751template <int dim, int spacedim>
752bool
754{
755 for (unsigned int ref_case = RefinementCase<dim>::cut_x;
756 ref_case < RefinementCase<dim>::isotropic_refinement + 1;
757 ++ref_case)
758 for (unsigned int c = 0;
759 c < GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
760 ++c)
761 {
762 // make sure also the lazily initialized matrices are created
764 Assert((restriction[ref_case - 1][c].m() == this->n_dofs_per_cell()) ||
765 (restriction[ref_case - 1][c].m() == 0),
767 Assert((restriction[ref_case - 1][c].n() == this->n_dofs_per_cell()) ||
768 (restriction[ref_case - 1][c].n() == 0),
770 if ((restriction[ref_case - 1][c].m() == 0) ||
771 (restriction[ref_case - 1][c].n() == 0))
772 return false;
773 }
774 return true;
775}
776
777
778
779template <int dim, int spacedim>
780bool
782{
783 const RefinementCase<dim> ref_case =
785
786 for (unsigned int c = 0;
787 c < GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
788 ++c)
789 {
790 // make sure also the lazily initialized matrices are created
792 Assert((prolongation[ref_case - 1][c].m() == this->n_dofs_per_cell()) ||
793 (prolongation[ref_case - 1][c].m() == 0),
795 Assert((prolongation[ref_case - 1][c].n() == this->n_dofs_per_cell()) ||
796 (prolongation[ref_case - 1][c].n() == 0),
798 if ((prolongation[ref_case - 1][c].m() == 0) ||
799 (prolongation[ref_case - 1][c].n() == 0))
800 return false;
801 }
802 return true;
803}
804
805
806
807template <int dim, int spacedim>
808bool
810{
811 const RefinementCase<dim> ref_case =
813
814 for (unsigned int c = 0;
815 c < GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
816 ++c)
817 {
818 // make sure also the lazily initialized matrices are created
820 Assert((restriction[ref_case - 1][c].m() == this->n_dofs_per_cell()) ||
821 (restriction[ref_case - 1][c].m() == 0),
823 Assert((restriction[ref_case - 1][c].n() == this->n_dofs_per_cell()) ||
824 (restriction[ref_case - 1][c].n() == 0),
826 if ((restriction[ref_case - 1][c].m() == 0) ||
827 (restriction[ref_case - 1][c].n() == 0))
828 return false;
829 }
830 return true;
831}
832
833
834
835template <int dim, int spacedim>
836bool
838 const internal::SubfaceCase<dim> &subface_case) const
839{
840 // TODO: the implementation makes the assumption that all faces have the
841 // same number of dofs
843 const unsigned int face_no = 0;
844
846 return (this->n_dofs_per_face(face_no) == 0) ||
847 (interface_constraints.m() != 0);
848 else
849 return false;
850}
851
852
853
854template <int dim, int spacedim>
855bool
857{
858 return false;
859}
860
861
862
863template <int dim, int spacedim>
864const FullMatrix<double> &
866 const internal::SubfaceCase<dim> &subface_case) const
867{
868 // TODO: the implementation makes the assumption that all faces have the
869 // same number of dofs
871 const unsigned int face_no = 0;
872 (void)face_no;
873
874 (void)subface_case;
876 ExcMessage("Constraints for this element are only implemented "
877 "for the case that faces are refined isotropically "
878 "(which is always the case in 2d, and in 3d requires "
879 "that the neighboring cell of a coarse cell presents "
880 "exactly four children on the common face)."));
881 Assert((this->n_dofs_per_face(face_no) == 0) ||
882 (interface_constraints.m() != 0),
883 ExcMessage("The finite element for which you try to obtain "
884 "hanging node constraints does not appear to "
885 "implement them."));
886
887 if (dim == 1)
891
893}
894
895
896
897template <int dim, int spacedim>
900{
901 // TODO: the implementation makes the assumption that all faces have the
902 // same number of dofs
904 const unsigned int face_no = 0;
905
906 switch (dim)
907 {
908 case 1:
909 return {0U, 0U};
910 case 2:
911 return {this->n_dofs_per_vertex() + 2 * this->n_dofs_per_line(),
912 this->n_dofs_per_face(face_no)};
913 case 3:
914 return {5 * this->n_dofs_per_vertex() + 12 * this->n_dofs_per_line() +
915 4 * this->n_dofs_per_quad(face_no),
916 this->n_dofs_per_face(face_no)};
917 default:
918 Assert(false, ExcNotImplemented());
919 }
921}
922
923
924
925template <int dim, int spacedim>
926void
929 FullMatrix<double> &) const
930{
931 // by default, no interpolation
932 // implemented. so throw exception,
933 // as documentation says
935 false,
937}
938
939
940
941template <int dim, int spacedim>
942void
946 const unsigned int) const
947{
948 // by default, no interpolation
949 // implemented. so throw exception,
950 // as documentation says
952 false,
954}
955
956
957
958template <int dim, int spacedim>
959void
962 const unsigned int,
964 const unsigned int) const
965{
966 // by default, no interpolation
967 // implemented. so throw exception,
968 // as documentation says
970 false,
972}
973
974
975
976template <int dim, int spacedim>
977std::vector<std::pair<unsigned int, unsigned int>>
979 const FiniteElement<dim, spacedim> &) const
980{
981 Assert(false, ExcNotImplemented());
982 return std::vector<std::pair<unsigned int, unsigned int>>();
983}
984
985
986
987template <int dim, int spacedim>
988std::vector<std::pair<unsigned int, unsigned int>>
990 const FiniteElement<dim, spacedim> &) const
991{
992 Assert(false, ExcNotImplemented());
993 return std::vector<std::pair<unsigned int, unsigned int>>();
994}
995
996
997
998template <int dim, int spacedim>
999std::vector<std::pair<unsigned int, unsigned int>>
1002 const unsigned int) const
1003{
1004 Assert(false, ExcNotImplemented());
1005 return std::vector<std::pair<unsigned int, unsigned int>>();
1006}
1007
1008
1009
1010template <int dim, int spacedim>
1014 const unsigned int) const
1015{
1016 Assert(false, ExcNotImplemented());
1018}
1019
1020
1021
1022template <int dim, int spacedim>
1023bool
1026{
1027 // Compare fields in roughly increasing order of how expensive the
1028 // comparison is
1029 return ((typeid(*this) == typeid(f)) && (this->get_name() == f.get_name()) &&
1030 (static_cast<const FiniteElementData<dim> &>(*this) ==
1031 static_cast<const FiniteElementData<dim> &>(f)) &&
1033}
1034
1035
1036
1037template <int dim, int spacedim>
1038bool
1041{
1042 return !(*this == f);
1043}
1044
1045
1046
1047template <int dim, int spacedim>
1048const std::vector<Point<dim>> &
1050{
1051 // a finite element may define
1052 // support points, but only if
1053 // there are as many as there are
1054 // degrees of freedom
1055 Assert((unit_support_points.size() == 0) ||
1056 (unit_support_points.size() == this->n_dofs_per_cell()),
1058 return unit_support_points;
1059}
1060
1061
1062
1063template <int dim, int spacedim>
1064bool
1066{
1067 return (unit_support_points.size() != 0);
1068}
1069
1070
1071
1072template <int dim, int spacedim>
1073const std::vector<Point<dim>> &
1075{
1076 // If the finite element implements generalized support points, return
1077 // those. Otherwise fall back to unit support points.
1078 return ((generalized_support_points.size() == 0) ?
1081}
1082
1083
1084
1085template <int dim, int spacedim>
1086bool
1088{
1089 return (get_generalized_support_points().size() != 0);
1090}
1091
1092
1093
1094template <int dim, int spacedim>
1096FiniteElement<dim, spacedim>::unit_support_point(const unsigned int index) const
1097{
1098 AssertIndexRange(index, this->n_dofs_per_cell());
1099 Assert(unit_support_points.size() == this->n_dofs_per_cell(),
1101 return unit_support_points[index];
1102}
1103
1104
1105
1106template <int dim, int spacedim>
1107const std::vector<Point<dim - 1>> &
1109 const unsigned int face_no) const
1110{
1111 // a finite element may define
1112 // support points, but only if
1113 // there are as many as there are
1114 // degrees of freedom on a face
1115 Assert((unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
1116 .size() == 0) ||
1117 (unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
1118 .size() == this->n_dofs_per_face(face_no)),
1120 return unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no];
1121}
1122
1123
1124
1125template <int dim, int spacedim>
1126bool
1128 const unsigned int face_no) const
1129{
1130 return (unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
1131 .size() != 0);
1132}
1133
1134
1135
1136template <int dim, int spacedim>
1137Point<dim - 1>
1139 const unsigned int index,
1140 const unsigned int face_no) const
1141{
1142 AssertIndexRange(index, this->n_dofs_per_face(face_no));
1143 Assert(unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
1144 .size() == this->n_dofs_per_face(face_no),
1146 return unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
1147 [index];
1148}
1149
1150
1151
1152template <int dim, int spacedim>
1153bool
1155 const unsigned int) const
1156{
1157 return true;
1158}
1159
1160
1161
1162template <int dim, int spacedim>
1165{
1166 // Translate the ComponentMask into first_selected and n_components after
1167 // some error checking:
1168 const unsigned int n_total_components = this->n_components();
1169 Assert((n_total_components == mask.size()) || (mask.size() == 0),
1170 ExcMessage("The given ComponentMask has the wrong size."));
1171
1172 const unsigned int n_selected =
1173 mask.n_selected_components(n_total_components);
1174 Assert(n_selected > 0,
1175 ExcMessage("You need at least one selected component."));
1176
1177 const unsigned int first_selected =
1178 mask.first_selected_component(n_total_components);
1179
1180#ifdef DEBUG
1181 // check that it is contiguous:
1182 for (unsigned int c = 0; c < n_total_components; ++c)
1183 Assert((c < first_selected && (!mask[c])) ||
1184 (c >= first_selected && c < first_selected + n_selected &&
1185 mask[c]) ||
1186 (c >= first_selected + n_selected && !mask[c]),
1187 ExcMessage("Error: the given ComponentMask is not contiguous!"));
1188#endif
1189
1190 return get_sub_fe(first_selected, n_selected);
1191}
1192
1193
1194
1195template <int dim, int spacedim>
1198 const unsigned int first_component,
1199 const unsigned int n_selected_components) const
1200{
1201 // No complicated logic is needed here, because it is overridden in
1202 // FESystem<dim,spacedim>. Just make sure that what the user chose is valid:
1203 Assert(first_component == 0 && n_selected_components == this->n_components(),
1204 ExcMessage(
1205 "You can only select a whole FiniteElement, not a part of one."));
1206
1207 (void)first_component;
1208 (void)n_selected_components;
1209 return *this;
1210}
1211
1212
1213
1214template <int dim, int spacedim>
1215std::pair<Table<2, bool>, std::vector<unsigned int>>
1217{
1218 Assert(false, ExcNotImplemented());
1219 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
1221 std::vector<unsigned int>(this->n_components()));
1222}
1223
1224
1225
1226template <int dim, int spacedim>
1227void
1230 const std::vector<Vector<double>> &,
1231 std::vector<double> &) const
1232{
1234 ExcMessage("The element for which you are calling the current "
1235 "function does not have generalized support points (see "
1236 "the glossary for a definition of generalized support "
1237 "points). Consequently, the current function can not "
1238 "be defined and is not implemented by the element."));
1239 Assert(false, ExcNotImplemented());
1240}
1241
1242
1243
1244template <int dim, int spacedim>
1245std::size_t
1247{
1248 return (
1249 sizeof(FiniteElementData<dim>) +
1261}
1262
1263
1264
1265template <int dim, int spacedim>
1266std::vector<unsigned int>
1268 const std::vector<ComponentMask> &nonzero_components)
1269{
1270 std::vector<unsigned int> retval(nonzero_components.size());
1271 for (unsigned int i = 0; i < nonzero_components.size(); ++i)
1272 retval[i] = nonzero_components[i].n_selected_components();
1273 return retval;
1274}
1275
1276
1277
1278/*------------------------------- FiniteElement ----------------------*/
1279
1280#ifndef DOXYGEN
1281template <int dim, int spacedim>
1282std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1284 const UpdateFlags flags,
1285 const Mapping<dim, spacedim> & mapping,
1286 const hp::QCollection<dim - 1> &quadrature,
1288 spacedim>
1289 &output_data) const
1290{
1291 return get_data(flags,
1292 mapping,
1294 quadrature),
1295 output_data);
1296}
1297
1298
1299
1300template <int dim, int spacedim>
1301std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1303 const UpdateFlags flags,
1304 const Mapping<dim, spacedim> &mapping,
1305 const Quadrature<dim - 1> & quadrature,
1307 spacedim>
1308 &output_data) const
1309{
1310 return get_data(flags,
1311 mapping,
1313 quadrature),
1314 output_data);
1315}
1316
1317
1318
1319template <int dim, int spacedim>
1320inline void
1323 const unsigned int face_no,
1324 const hp::QCollection<dim - 1> & quadrature,
1325 const Mapping<dim, spacedim> & mapping,
1326 const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
1327 const ::internal::FEValuesImplementation::MappingRelatedData<dim,
1328 spacedim>
1329 & mapping_data,
1330 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1332 spacedim>
1333 &output_data) const
1334{
1335 // base class version, implement overridden function in derived classes
1336 AssertDimension(quadrature.size(), 1);
1338 face_no,
1339 quadrature[0],
1340 mapping,
1341 mapping_internal,
1342 mapping_data,
1343 fe_internal,
1344 output_data);
1345}
1346
1347
1348
1349template <int dim, int spacedim>
1350inline void
1353 const unsigned int face_no,
1354 const Quadrature<dim - 1> & quadrature,
1355 const Mapping<dim, spacedim> & mapping,
1356 const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
1357 const ::internal::FEValuesImplementation::MappingRelatedData<dim,
1358 spacedim>
1359 & mapping_data,
1360 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1362 spacedim>
1363 &output_data) const
1364{
1365 Assert(false,
1366 ExcMessage("Use of a deprecated interface, please implement "
1367 "fill_fe_face_values taking a hp::QCollection argument"));
1368 (void)cell;
1369 (void)face_no;
1370 (void)quadrature;
1371 (void)mapping;
1372 (void)mapping_internal;
1373 (void)mapping_data;
1374 (void)fe_internal;
1375 (void)output_data;
1376}
1377#endif
1378
1379
1380
1381template <int dim, int spacedim>
1382std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1384 const UpdateFlags flags,
1385 const Mapping<dim, spacedim> &mapping,
1386 const Quadrature<dim - 1> & quadrature,
1388 spacedim>
1389 &output_data) const
1390{
1391 return get_data(flags,
1392 mapping,
1394 this->reference_cell(), quadrature),
1395 output_data);
1396}
1397
1398
1399
1400template <int dim, int spacedim>
1402FiniteElement<dim, spacedim>::base_element(const unsigned int index) const
1403{
1404 (void)index;
1405 AssertIndexRange(index, 1);
1406 // This function should not be
1407 // called for a system element
1409 return *this;
1410}
1411
1412
1413
1414/*------------------------------- Explicit Instantiations -------------*/
1415#include "fe.inst"
1416
1417
unsigned int size() const
bool represents_the_all_selected_mask() const
unsigned int size() const
bool represents_the_all_selected_mask() const
unsigned int size() const
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
unsigned int first_selected_component(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
unsigned int get_first_line_index() const
unsigned int n_dofs_per_vertex() const
unsigned int n_dofs_per_cell() const
unsigned int n_dofs_per_line() const
unsigned int get_first_quad_index(const unsigned int quad_no=0) const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_blocks() const
unsigned int n_components() const
unsigned int n_unique_quads() const
ReferenceCell reference_cell() const
unsigned int get_first_face_line_index(const unsigned int face_no=0) const
unsigned int get_first_face_quad_index(const unsigned int face_no=0) const
unsigned int n_unique_faces() const
unsigned int n_dofs_per_quad(unsigned int face_no=0) const
virtual std::size_t memory_consumption() const
bool constraints_are_implemented(const ::internal::SubfaceCase< dim > &subface_case=::internal::SubfaceCase< dim >::case_isotropic) const
bool isotropic_prolongation_is_implemented() const
virtual std::string get_name() const =0
virtual Point< dim > unit_support_point(const unsigned int index) const
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const
const std::vector< unsigned int > n_nonzero_components_table
Definition: fe.h:2590
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
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
std::vector< std::vector< Point< dim - 1 > > > unit_face_support_points
Definition: fe.h:2449
bool prolongation_is_implemented() const
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const
virtual std::unique_ptr< InternalDataBase > get_face_data(const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const hp::QCollection< dim - 1 > &quadrature, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const
virtual Point< dim - 1 > unit_face_support_point(const unsigned int index, const unsigned int face_no=0) const
static std::vector< unsigned int > compute_n_nonzero_components(const std::vector< ComponentMask > &nonzero_components)
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
bool has_face_support_points(const unsigned int face_no=0) const
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other, const unsigned int face_no=0) const
ComponentMask component_mask(const FEValuesExtractors::Scalar &scalar) const
bool operator!=(const FiniteElement< dim, spacedim > &) const
bool has_support_points() const
std::vector< std::vector< FullMatrix< double > > > restriction
Definition: fe.h:2404
bool has_generalized_support_points() const
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const =0
virtual Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const
std::vector< Table< 2, int > > adjust_quad_dof_index_for_face_orientation_table
Definition: fe.h:2478
virtual bool operator==(const FiniteElement< dim, spacedim > &fe) const
unsigned int adjust_line_dof_index_for_line_orientation(const unsigned int index, const bool line_orientation) const
const std::vector< Point< dim > > & get_unit_support_points() const
BlockIndices base_to_block_indices
Definition: fe.h:2542
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
const std::vector< Point< dim - 1 > > & get_unit_face_support_points(const unsigned int face_no=0) const
bool isotropic_restriction_is_implemented() const
std::vector< int > adjust_line_dof_index_for_line_orientation_table
Definition: fe.h:2493
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const
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
const FiniteElement< dim, spacedim > & get_sub_fe(const ComponentMask &mask) const
virtual Tensor< 4, dim > shape_4th_derivative_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const
unsigned int component_to_block_index(const unsigned int component) const
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const hp::QCollection< dim - 1 > &quadrature, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &mapping_data, const InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const
BlockMask block_mask(const FEValuesExtractors::Scalar &scalar) const
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const
virtual Tensor< 3, dim > shape_3rd_derivative(const unsigned int i, const Point< dim > &p) const
std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > system_to_base_table
Definition: fe.h:2529
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
virtual void convert_generalized_support_point_values_to_dof_values(const std::vector< Vector< double > > &support_point_values, std::vector< double > &nodal_values) const
FiniteElement(const FiniteElementData< dim > &fe_data, const std::vector< bool > &restriction_is_additive_flags, const std::vector< ComponentMask > &nonzero_components)
std::vector< std::pair< unsigned int, unsigned int > > system_to_component_table
Definition: fe.h:2498
virtual Tensor< 3, dim > shape_3rd_derivative_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
const FullMatrix< double > & constraints(const ::internal::SubfaceCase< dim > &subface_case=::internal::SubfaceCase< dim >::case_isotropic) const
std::pair< std::unique_ptr< FiniteElement< dim, spacedim > >, unsigned int > operator^(const unsigned int multiplicity) const
const std::vector< bool > restriction_is_additive_flags
Definition: fe.h:2572
std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > component_to_base_table
Definition: fe.h:2565
TableIndices< 2 > interface_constraints_size() const
std::vector< std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > > face_system_to_base_table
Definition: fe.h:2536
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2442
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
bool restriction_is_implemented() const
std::vector< std::vector< std::pair< unsigned int, unsigned int > > > face_system_to_component_table
Definition: fe.h:2510
FullMatrix< double > interface_constraints
Definition: fe.h:2430
const std::vector< Point< dim > > & get_generalized_support_points() const
virtual Tensor< 2, dim > shape_grad_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
virtual double shape_value(const unsigned int i, const Point< dim > &p) const
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const
virtual std::size_t memory_consumption() const
virtual Tensor< 4, dim > shape_4th_derivative(const unsigned int i, const Point< dim > &p) const
std::vector< Point< dim > > generalized_support_points
Definition: fe.h:2455
unsigned int adjust_quad_dof_index_for_face_orientation(const unsigned int index, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation) const
const std::vector< ComponentMask > nonzero_components
Definition: fe.h:2581
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
virtual bool hp_constraints_are_implemented() const
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2418
types::global_dof_index first_block_of_base(const unsigned int b) const
size_type n() const
size_type m() const
Definition: point.h:111
Definition: vector.h:110
unsigned int size() const
Definition: collection.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
UpdateFlags
@ update_default
No update.
Point< 2 > first
Definition: grid_out.cc:4587
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcWrongInterfaceMatrixSize(int arg1, int arg2)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
static ::ExceptionBase & ExcUnitShapeValuesDoNotExist()
static ::ExceptionBase & ExcProjectionVoid()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcEmbeddingVoid()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcFEHasNoSupportPoints()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
static const char U
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
constexpr const ReferenceCell Quadrilateral
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
static const unsigned int invalid_unsigned_int
Definition: types.h:196
STL namespace.
static unsigned int n_children(const RefinementCase< dim > &refinement_case)