Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
fe.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1998 - 2023 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)
63 dim == 3 ? this->n_dofs_per_line() : 0)
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
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)
76 nonzero_c.size() == 1 ?
77 std::vector<ComponentMask>(fe_data.n_dofs_per_cell(), nonzero_c[0]) :
78 nonzero_c)
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_2d_subobjects());
156
157 for (unsigned int f = 0; f < this->n_unique_2d_subobjects(); ++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().n_face_orientations(f));
162 adjust_quad_dof_index_for_face_orientation_table[f].fill(0);
163 }
164 }
165
166 unit_face_support_points.resize(this->n_unique_faces());
167 generalized_face_support_points.resize(this->n_unique_faces());
168}
169
170
171
172template <int dim, int spacedim>
173std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>
174FiniteElement<dim, spacedim>::operator^(const unsigned int multiplicity) const
175{
176 return {this->clone(), multiplicity};
177}
178
179
180
181template <int dim, int spacedim>
182double
184 const Point<dim> &) const
185{
186 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
187 return 0.;
188}
189
190
191
192template <int dim, int spacedim>
193double
195 const Point<dim> &,
196 const unsigned int) const
197{
198 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
199 return 0.;
200}
201
202
203
204template <int dim, int spacedim>
207 const Point<dim> &) const
208{
209 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
210 return Tensor<1, dim>();
211}
212
213
214
215template <int dim, int spacedim>
218 const Point<dim> &,
219 const unsigned int) const
220{
221 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
222 return Tensor<1, dim>();
223}
224
225
226
227template <int dim, int spacedim>
230 const Point<dim> &) const
231{
232 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
233 return Tensor<2, dim>();
234}
235
236
237
238template <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<2, dim>();
247}
248
249
250
251template <int dim, int spacedim>
254 const Point<dim> &) const
255{
256 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
257 return Tensor<3, dim>();
258}
259
260
261
262template <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<3, dim>();
271}
272
273
274
275template <int dim, int spacedim>
278 const Point<dim> &) const
279{
280 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
281 return Tensor<4, dim>();
282}
283
284
285
286template <int dim, int spacedim>
289 const unsigned int,
290 const Point<dim> &,
291 const unsigned int) const
292{
293 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
294 return Tensor<4, dim>();
295}
296
297
298template <int dim, int spacedim>
299void
301 const bool isotropic_restriction_only,
302 const bool isotropic_prolongation_only)
303{
304 for (unsigned int ref_case = RefinementCase<dim>::cut_x;
305 ref_case <= RefinementCase<dim>::isotropic_refinement;
306 ++ref_case)
307 {
308 const unsigned int nc =
310
311 for (unsigned int i = 0; i < nc; ++i)
312 {
313 if (this->restriction[ref_case - 1][i].m() !=
314 this->n_dofs_per_cell() &&
315 (!isotropic_restriction_only ||
317 this->restriction[ref_case - 1][i].reinit(this->n_dofs_per_cell(),
318 this->n_dofs_per_cell());
319 if (this->prolongation[ref_case - 1][i].m() !=
320 this->n_dofs_per_cell() &&
321 (!isotropic_prolongation_only ||
323 this->prolongation[ref_case - 1][i].reinit(this->n_dofs_per_cell(),
324 this->n_dofs_per_cell());
325 }
326 }
327}
328
329
330template <int dim, int spacedim>
331const FullMatrix<double> &
333 const unsigned int child,
334 const RefinementCase<dim> &refinement_case) const
335{
336 AssertIndexRange(refinement_case,
340 "Restriction matrices are only available for refined cells!"));
342 child, GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)));
343 // we use refinement_case-1 here. the -1 takes care of the origin of the
344 // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
345 // available and so the vector indices are shifted
346 Assert(restriction[refinement_case - 1][child].n() == this->n_dofs_per_cell(),
347 ExcProjectionVoid());
348 return restriction[refinement_case - 1][child];
349}
350
351
352
353template <int dim, int spacedim>
354const FullMatrix<double> &
356 const unsigned int child,
357 const RefinementCase<dim> &refinement_case) const
358{
359 AssertIndexRange(refinement_case,
363 "Prolongation matrices are only available for refined cells!"));
365 child, GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)));
366 // we use refinement_case-1 here. the -1 takes care
367 // of the origin of the vector, as for
368 // RefinementCase::no_refinement (=0) there is no
369 // data available and so the vector indices
370 // are shifted
371 Assert(prolongation[refinement_case - 1][child].n() ==
372 this->n_dofs_per_cell(),
373 ExcEmbeddingVoid());
374 return prolongation[refinement_case - 1][child];
375}
376
377
378// TODO:[GK] This is probably not the most efficient way of doing this.
379template <int dim, int spacedim>
380unsigned int
382 const unsigned int index) const
383{
384 AssertIndexRange(index, this->n_components());
385
386 return first_block_of_base(component_to_base_table[index].first.first) +
387 component_to_base_table[index].second;
388}
389
390
391template <int dim, int spacedim>
394 const FEValuesExtractors::Scalar &scalar) const
395{
396 AssertIndexRange(scalar.component, this->n_components());
397
398 // TODO: it would be nice to verify that it is indeed possible
399 // to select this scalar component, i.e., that it is not part
400 // of a non-primitive element. unfortunately, there is no simple
401 // way to write such a condition...
402
403 std::vector<bool> mask(this->n_components(), false);
404 mask[scalar.component] = true;
405 return mask;
406}
407
408
409template <int dim, int spacedim>
412 const FEValuesExtractors::Vector &vector) const
413{
415 this->n_components());
416
417 // TODO: it would be nice to verify that it is indeed possible
418 // to select these vector components, i.e., that they don't span
419 // beyond the beginning or end of anon-primitive element.
420 // unfortunately, there is no simple way to write such a condition...
421
422 std::vector<bool> mask(this->n_components(), false);
423 for (unsigned int c = vector.first_vector_component;
424 c < vector.first_vector_component + dim;
425 ++c)
426 mask[c] = true;
427 return mask;
428}
429
430
431template <int dim, int spacedim>
434 const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const
435{
438 this->n_components());
439
440 // TODO: it would be nice to verify that it is indeed possible
441 // to select these vector components, i.e., that they don't span
442 // beyond the beginning or end of anon-primitive element.
443 // unfortunately, there is no simple way to write such a condition...
444
445 std::vector<bool> mask(this->n_components(), false);
446 for (unsigned int c = sym_tensor.first_tensor_component;
447 c < sym_tensor.first_tensor_component +
449 ++c)
450 mask[c] = true;
451 return mask;
452}
453
454
455
456template <int dim, int spacedim>
459{
460 // if we get a block mask that represents all blocks, then
461 // do the same for the returned component mask
462 if (block_mask.represents_the_all_selected_mask())
463 return {};
464
465 AssertDimension(block_mask.size(), this->n_blocks());
466
467 std::vector<bool> component_mask(this->n_components(), false);
468 for (unsigned int c = 0; c < this->n_components(); ++c)
469 if (block_mask[component_to_block_index(c)] == true)
470 component_mask[c] = true;
471
472 return component_mask;
473}
474
475
476
477template <int dim, int spacedim>
480 const FEValuesExtractors::Scalar &scalar) const
481{
482 // simply create the corresponding component mask (a simpler
483 // process) and then convert it to a block mask
484 return block_mask(component_mask(scalar));
485}
486
487
488template <int dim, int spacedim>
491 const FEValuesExtractors::Vector &vector) const
492{
493 // simply create the corresponding component mask (a simpler
494 // process) and then convert it to a block mask
495 return block_mask(component_mask(vector));
496}
497
498
499template <int dim, int spacedim>
502 const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const
503{
504 // simply create the corresponding component mask (a simpler
505 // process) and then convert it to a block mask
506 return block_mask(component_mask(sym_tensor));
507}
508
509
510
511template <int dim, int spacedim>
514 const ComponentMask &component_mask) const
515{
516 // if we get a component mask that represents all component, then
517 // do the same for the returned block mask
518 if (component_mask.represents_the_all_selected_mask())
519 return {};
520
521 AssertDimension(component_mask.size(), this->n_components());
522
523 // walk over all of the components
524 // of this finite element and see
525 // if we need to set the
526 // corresponding block. inside the
527 // block, walk over all the
528 // components that correspond to
529 // this block and make sure the
530 // component mask is set for all of
531 // them
532 std::vector<bool> block_mask(this->n_blocks(), false);
533 for (unsigned int c = 0; c < this->n_components();)
534 {
535 const unsigned int block = component_to_block_index(c);
536 if (component_mask[c] == true)
537 block_mask[block] = true;
538
539 // now check all of the other
540 // components that correspond
541 // to this block
542 ++c;
543 while ((c < this->n_components()) &&
544 (component_to_block_index(c) == block))
545 {
546 Assert(component_mask[c] == block_mask[block],
548 "The component mask argument given to this function "
549 "is not a mask where the individual components belonging "
550 "to one block of the finite element are either all "
551 "selected or not selected. You can't call this function "
552 "with a component mask that splits blocks."));
553 ++c;
554 }
555 }
556
557
558 return block_mask;
559}
560
561
562
563template <int dim, int spacedim>
564unsigned int
565FiniteElement<dim, spacedim>::face_to_cell_index(const unsigned int face_index,
566 const unsigned int face,
567 const bool face_orientation,
568 const bool face_flip,
569 const bool face_rotation) const
570{
571 AssertIndexRange(face_index, this->n_dofs_per_face(face));
572 AssertIndexRange(face, this->reference_cell().n_faces());
573
574 // TODO: we could presumably solve the 3d case below using the
575 // adjust_quad_dof_index_for_face_orientation_table field. for the
576 // 2d case, we can't use adjust_line_dof_index_for_line_orientation_table
577 // since that array is empty (presumably because we thought that
578 // there are no flipped edges in 2d, but these can happen in
579 // DoFTools::make_periodicity_constraints, for example). so we
580 // would need to either fill this field, or rely on derived classes
581 // implementing this function, as we currently do
582
583 // see the function's documentation for an explanation of this
584 // assertion -- in essence, derived classes have to implement
585 // an overloaded version of this function if we are to use any
586 // other than standard orientation
587 if ((face_orientation != true) || (face_flip != false) ||
588 (face_rotation != false))
589 Assert((this->n_dofs_per_line() <= 1) && (this->n_dofs_per_quad(face) <= 1),
591 "The function in this base class can not handle this case. "
592 "Rather, the derived class you are using must provide "
593 "an overloaded version but apparently hasn't done so. See "
594 "the documentation of this function for more information."));
595
596 // we need to distinguish between DoFs on vertices, lines and in 3d quads.
597 // do so in a sequence of if-else statements
598 if (face_index < this->get_first_face_line_index(face))
599 // DoF is on a vertex
600 {
601 // get the number of the vertex on the face that corresponds to this DoF,
602 // along with the number of the DoF on this vertex
603 const unsigned int face_vertex = face_index / this->n_dofs_per_vertex();
604 const unsigned int dof_index_on_vertex =
605 face_index % this->n_dofs_per_vertex();
606
607 // then get the number of this vertex on the cell and translate
608 // this to a DoF number on the cell
609 return (this->reference_cell().face_to_cell_vertices(
610 face,
611 face_vertex,
612 (face_orientation ? 1 : 0) + (face_rotation ? 2 : 0) +
613 (face_flip ? 4 : 0)) *
614 this->n_dofs_per_vertex() +
615 dof_index_on_vertex);
616 }
617 else if (face_index < this->get_first_face_quad_index(face))
618 // DoF is on a face
619 {
620 // do the same kind of translation as before. we need to only consider
621 // DoFs on the lines, i.e., ignoring those on the vertices
622 const unsigned int index =
623 face_index - this->get_first_face_line_index(face);
624
625 const unsigned int face_line = index / this->n_dofs_per_line();
626 const unsigned int dof_index_on_line = index % this->n_dofs_per_line();
627
628 return (
629 this->get_first_line_index() +
630 this->reference_cell().face_to_cell_lines(face,
631 face_line,
632 (face_orientation ? 1 : 0) +
633 (face_rotation ? 2 : 0) +
634 (face_flip ? 4 : 0)) *
635 this->n_dofs_per_line() +
636 dof_index_on_line);
637 }
638 else
639 // DoF is on a quad
640 {
641 Assert(dim >= 3, ExcInternalError());
642
643 // ignore vertex and line dofs
644 const unsigned int index =
645 face_index - this->get_first_face_quad_index(face);
646
647 return (this->get_first_quad_index(face) + index);
648 }
649}
650
651
652
653template <int dim, int spacedim>
654unsigned int
656 const unsigned int index,
657 const unsigned int face,
658 const bool face_orientation,
659 const bool face_flip,
660 const bool face_rotation) const
661{
662 // general template for 1d and 2d: not
663 // implemented. in fact, the function
664 // shouldn't even be called unless we are
665 // in 3d, so throw an internal error
666 Assert(dim == 3, ExcInternalError());
667 if (dim < 3)
668 return index;
669
670 // adjust dofs on 3d faces if the face is
671 // flipped. note that we query a table that
672 // derived elements need to have set up
673 // front. the exception are discontinuous
674 // elements for which there should be no
675 // face dofs anyway (i.e. dofs_per_quad==0
676 // in 3d), so we don't need the table, but
677 // the function should also not have been
678 // called
679 AssertIndexRange(index, this->n_dofs_per_quad(face));
680 const auto table_n = this->n_unique_2d_subobjects() == 1 ? 0 : face;
681 Assert(
682 adjust_quad_dof_index_for_face_orientation_table[table_n].n_elements() ==
683 (this->reference_cell().n_face_orientations(face)) *
684 this->n_dofs_per_quad(face),
686 return index + adjust_quad_dof_index_for_face_orientation_table[table_n](
687 index,
688 (face_orientation ? 4 : 0) + (face_flip ? 2 : 0) +
689 (face_rotation ? 1 : 0));
690}
691
692
693
694template <int dim, int spacedim>
695unsigned int
697 const unsigned int index,
698 const bool line_orientation) const
699{
700 // general template for 1d and 2d: do
701 // nothing. Do not throw an Assertion,
702 // however, in order to allow to call this
703 // function in 2d as well
704 if (dim < 3)
705 return index;
706
707 AssertIndexRange(index, this->n_dofs_per_line());
708 Assert(adjust_line_dof_index_for_line_orientation_table.size() ==
709 this->n_dofs_per_line(),
711 if (line_orientation)
712 return index;
713 else
714 return index + adjust_line_dof_index_for_line_orientation_table[index];
715}
716
717
718
719template <int dim, int spacedim>
720bool
722{
723 for (unsigned int ref_case = RefinementCase<dim>::cut_x;
724 ref_case < RefinementCase<dim>::isotropic_refinement + 1;
725 ++ref_case)
726 for (unsigned int c = 0;
727 c < GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
728 ++c)
729 {
730 // make sure also the lazily initialized matrices are created
731 get_prolongation_matrix(c, RefinementCase<dim>(ref_case));
732 Assert((prolongation[ref_case - 1][c].m() == this->n_dofs_per_cell()) ||
733 (prolongation[ref_case - 1][c].m() == 0),
735 Assert((prolongation[ref_case - 1][c].n() == this->n_dofs_per_cell()) ||
736 (prolongation[ref_case - 1][c].n() == 0),
738 if ((prolongation[ref_case - 1][c].m() == 0) ||
739 (prolongation[ref_case - 1][c].n() == 0))
740 return false;
741 }
742 return true;
743}
744
745
746
747template <int dim, int spacedim>
748bool
750{
751 for (unsigned int ref_case = RefinementCase<dim>::cut_x;
752 ref_case < RefinementCase<dim>::isotropic_refinement + 1;
753 ++ref_case)
754 for (unsigned int c = 0;
755 c < GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
756 ++c)
757 {
758 // make sure also the lazily initialized matrices are created
759 get_restriction_matrix(c, RefinementCase<dim>(ref_case));
760 Assert((restriction[ref_case - 1][c].m() == this->n_dofs_per_cell()) ||
761 (restriction[ref_case - 1][c].m() == 0),
763 Assert((restriction[ref_case - 1][c].n() == this->n_dofs_per_cell()) ||
764 (restriction[ref_case - 1][c].n() == 0),
766 if ((restriction[ref_case - 1][c].m() == 0) ||
767 (restriction[ref_case - 1][c].n() == 0))
768 return false;
769 }
770 return true;
771}
772
773
774
775template <int dim, int spacedim>
776bool
778{
779 const RefinementCase<dim> ref_case =
781
782 for (unsigned int c = 0;
783 c < GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
784 ++c)
785 {
786 // make sure also the lazily initialized matrices are created
787 get_prolongation_matrix(c, RefinementCase<dim>(ref_case));
788 Assert((prolongation[ref_case - 1][c].m() == this->n_dofs_per_cell()) ||
789 (prolongation[ref_case - 1][c].m() == 0),
791 Assert((prolongation[ref_case - 1][c].n() == this->n_dofs_per_cell()) ||
792 (prolongation[ref_case - 1][c].n() == 0),
794 if ((prolongation[ref_case - 1][c].m() == 0) ||
795 (prolongation[ref_case - 1][c].n() == 0))
796 return false;
797 }
798 return true;
799}
800
801
802
803template <int dim, int spacedim>
804bool
806{
807 const RefinementCase<dim> ref_case =
809
810 for (unsigned int c = 0;
811 c < GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
812 ++c)
813 {
814 // make sure also the lazily initialized matrices are created
815 get_restriction_matrix(c, RefinementCase<dim>(ref_case));
816 Assert((restriction[ref_case - 1][c].m() == this->n_dofs_per_cell()) ||
817 (restriction[ref_case - 1][c].m() == 0),
819 Assert((restriction[ref_case - 1][c].n() == this->n_dofs_per_cell()) ||
820 (restriction[ref_case - 1][c].n() == 0),
822 if ((restriction[ref_case - 1][c].m() == 0) ||
823 (restriction[ref_case - 1][c].n() == 0))
824 return false;
825 }
826 return true;
827}
828
829
830
831template <int dim, int spacedim>
832bool
834 const internal::SubfaceCase<dim> &subface_case) const
835{
837 {
838 unsigned int n_dofs_on_faces = 0;
839
840 for (const auto face_no : this->reference_cell().face_indices())
841 n_dofs_on_faces += this->n_dofs_per_face(face_no);
842
843 return (n_dofs_on_faces == 0) || (interface_constraints.m() != 0);
844 }
845 else
846 return false;
847}
848
849
850
851template <int dim, int spacedim>
852bool
854{
855 return false;
856}
857
858
859
860template <int dim, int spacedim>
861const FullMatrix<double> &
863 const internal::SubfaceCase<dim> &subface_case) const
864{
865 // TODO: the implementation makes the assumption that all faces have the
866 // same number of dofs
867 AssertDimension(this->n_unique_faces(), 1);
868 const unsigned int face_no = 0;
869 (void)face_no;
870
871 (void)subface_case;
873 ExcMessage("Constraints for this element are only implemented "
874 "for the case that faces are refined isotropically "
875 "(which is always the case in 2d, and in 3d requires "
876 "that the neighboring cell of a coarse cell presents "
877 "exactly four children on the common face)."));
878 Assert((this->n_dofs_per_face(face_no) == 0) ||
879 (interface_constraints.m() != 0),
880 ExcMessage("The finite element for which you try to obtain "
881 "hanging node constraints does not appear to "
882 "implement them."));
883
884 if (dim == 1)
885 Assert((interface_constraints.m() == 0) && (interface_constraints.n() == 0),
886 ExcWrongInterfaceMatrixSize(interface_constraints.m(),
887 interface_constraints.n()));
888
889 return interface_constraints;
890}
891
892
893
894template <int dim, int spacedim>
897{
898 // TODO: the implementation makes the assumption that all faces have the
899 // same number of dofs
900 AssertDimension(this->n_unique_faces(), 1);
901 const unsigned int face_no = 0;
902
903 switch (dim)
904 {
905 case 1:
906 return {0U, 0U};
907 case 2:
908 return {this->n_dofs_per_vertex() + 2 * this->n_dofs_per_line(),
909 this->n_dofs_per_face(face_no)};
910 case 3:
911 return {5 * this->n_dofs_per_vertex() + 12 * this->n_dofs_per_line() +
912 4 * this->n_dofs_per_quad(face_no),
913 this->n_dofs_per_face(face_no)};
914 default:
915 Assert(false, ExcNotImplemented());
916 }
918}
919
920
921
922template <int dim, int spacedim>
923void
926 FullMatrix<double> &) const
927{
928 // by default, no interpolation
929 // implemented. so throw exception,
930 // as documentation says
932 false,
934}
935
936
937
938template <int dim, int spacedim>
939void
943 const unsigned int) const
944{
945 // by default, no interpolation
946 // implemented. so throw exception,
947 // as documentation says
949 false,
951}
952
953
954
955template <int dim, int spacedim>
956void
959 const unsigned int,
961 const unsigned int) const
962{
963 // by default, no interpolation
964 // implemented. so throw exception,
965 // as documentation says
967 false,
969}
970
971
972
973template <int dim, int spacedim>
974std::vector<std::pair<unsigned int, unsigned int>>
976 const FiniteElement<dim, spacedim> &) const
977{
978 Assert(false, ExcNotImplemented());
979 return std::vector<std::pair<unsigned int, unsigned int>>();
980}
981
982
983
984template <int dim, int spacedim>
985std::vector<std::pair<unsigned int, unsigned int>>
987 const FiniteElement<dim, spacedim> &) const
988{
989 Assert(false, ExcNotImplemented());
990 return std::vector<std::pair<unsigned int, unsigned int>>();
991}
992
993
994
995template <int dim, int spacedim>
996std::vector<std::pair<unsigned int, unsigned int>>
999 const unsigned int) const
1000{
1001 Assert(false, ExcNotImplemented());
1002 return std::vector<std::pair<unsigned int, unsigned int>>();
1003}
1004
1005
1006
1007template <int dim, int spacedim>
1011 const unsigned int) const
1012{
1013 Assert(false, ExcNotImplemented());
1015}
1016
1017
1018
1019template <int dim, int spacedim>
1020bool
1022 const FiniteElement<dim, spacedim> &f) const
1023{
1024 // Compare fields in roughly increasing order of how expensive the
1025 // comparison is
1026 return ((typeid(*this) == typeid(f)) && (this->get_name() == f.get_name()) &&
1027 (static_cast<const FiniteElementData<dim> &>(*this) ==
1028 static_cast<const FiniteElementData<dim> &>(f)) &&
1029 (interface_constraints == f.interface_constraints));
1030}
1031
1032
1033
1034template <int dim, int spacedim>
1035bool
1037 const FiniteElement<dim, spacedim> &f) const
1038{
1039 return !(*this == f);
1040}
1041
1042
1043
1044template <int dim, int spacedim>
1045const std::vector<Point<dim>> &
1047{
1048 // a finite element may define
1049 // support points, but only if
1050 // there are as many as there are
1051 // degrees of freedom
1052 Assert((unit_support_points.size() == 0) ||
1053 (unit_support_points.size() == this->n_dofs_per_cell()),
1055 return unit_support_points;
1056}
1057
1058
1059
1060template <int dim, int spacedim>
1061bool
1063{
1064 return (unit_support_points.size() != 0);
1065}
1066
1067
1068
1069template <int dim, int spacedim>
1070const std::vector<Point<dim>> &
1072{
1073 // If the finite element implements generalized support points, return
1074 // those. Otherwise fall back to unit support points.
1075 return ((generalized_support_points.size() == 0) ?
1076 unit_support_points :
1077 generalized_support_points);
1078}
1079
1080
1081
1082template <int dim, int spacedim>
1083bool
1085{
1086 return (get_generalized_support_points().size() != 0);
1087}
1088
1089
1090
1091template <int dim, int spacedim>
1093FiniteElement<dim, spacedim>::unit_support_point(const unsigned int index) const
1094{
1095 AssertIndexRange(index, this->n_dofs_per_cell());
1096 Assert(unit_support_points.size() == this->n_dofs_per_cell(),
1097 ExcFEHasNoSupportPoints());
1098 return unit_support_points[index];
1099}
1100
1101
1102
1103template <int dim, int spacedim>
1104const std::vector<Point<dim - 1>> &
1106 const unsigned int face_no) const
1107{
1108 // a finite element may define
1109 // support points, but only if
1110 // there are as many as there are
1111 // degrees of freedom on a face
1112 Assert((unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
1113 .size() == 0) ||
1114 (unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
1115 .size() == this->n_dofs_per_face(face_no)),
1117 return unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no];
1118}
1119
1120
1121
1122template <int dim, int spacedim>
1123bool
1125 const unsigned int face_no) const
1126{
1127 return (unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
1128 .size() != 0);
1129}
1130
1131
1132
1133template <int dim, int spacedim>
1134Point<dim - 1>
1136 const unsigned int index,
1137 const unsigned int face_no) const
1138{
1139 AssertIndexRange(index, this->n_dofs_per_face(face_no));
1140 Assert(unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
1141 .size() == this->n_dofs_per_face(face_no),
1142 ExcFEHasNoSupportPoints());
1143 return unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
1144 [index];
1145}
1146
1147
1148
1149template <int dim, int spacedim>
1150bool
1152 const unsigned int) const
1153{
1154 return true;
1155}
1156
1157
1158
1159template <int dim, int spacedim>
1162{
1163 // Translate the ComponentMask into first_selected and n_components after
1164 // some error checking:
1165 const unsigned int n_total_components = this->n_components();
1166 Assert((n_total_components == mask.size()) || (mask.size() == 0),
1167 ExcMessage("The given ComponentMask has the wrong size."));
1168
1169 const unsigned int n_selected =
1170 mask.n_selected_components(n_total_components);
1171 Assert(n_selected > 0,
1172 ExcMessage("You need at least one selected component."));
1173
1174 const unsigned int first_selected =
1175 mask.first_selected_component(n_total_components);
1176
1177#ifdef DEBUG
1178 // check that it is contiguous:
1179 for (unsigned int c = 0; c < n_total_components; ++c)
1180 Assert((c < first_selected && (!mask[c])) ||
1181 (c >= first_selected && c < first_selected + n_selected &&
1182 mask[c]) ||
1183 (c >= first_selected + n_selected && !mask[c]),
1184 ExcMessage("Error: the given ComponentMask is not contiguous!"));
1185#endif
1186
1187 return get_sub_fe(first_selected, n_selected);
1188}
1189
1190
1191
1192template <int dim, int spacedim>
1195 const unsigned int first_component,
1196 const unsigned int n_selected_components) const
1197{
1198 // No complicated logic is needed here, because it is overridden in
1199 // FESystem<dim,spacedim>. Just make sure that what the user chose is valid:
1200 Assert(first_component == 0 && n_selected_components == this->n_components(),
1201 ExcMessage(
1202 "You can only select a whole FiniteElement, not a part of one."));
1203
1204 (void)first_component;
1205 (void)n_selected_components;
1206 return *this;
1207}
1208
1209
1210
1211template <int dim, int spacedim>
1212std::pair<Table<2, bool>, std::vector<unsigned int>>
1214{
1215 Assert(false, ExcNotImplemented());
1216 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
1217 Table<2, bool>(this->n_components(), this->n_dofs_per_cell()),
1218 std::vector<unsigned int>(this->n_components()));
1219}
1220
1221
1222
1223template <int dim, int spacedim>
1224void
1227 const std::vector<Vector<double>> &,
1228 std::vector<double> &) const
1229{
1230 Assert(has_generalized_support_points(),
1231 ExcMessage("The element for which you are calling the current "
1232 "function does not have generalized support points (see "
1233 "the glossary for a definition of generalized support "
1234 "points). Consequently, the current function can not "
1235 "be defined and is not implemented by the element."));
1236 Assert(false, ExcNotImplemented());
1237}
1238
1239
1240
1241template <int dim, int spacedim>
1242std::size_t
1244{
1245 return (
1246 sizeof(FiniteElementData<dim>) +
1249 MemoryConsumption::memory_consumption(interface_constraints) +
1250 MemoryConsumption::memory_consumption(system_to_component_table) +
1251 MemoryConsumption::memory_consumption(face_system_to_component_table) +
1252 MemoryConsumption::memory_consumption(system_to_base_table) +
1253 MemoryConsumption::memory_consumption(face_system_to_base_table) +
1254 MemoryConsumption::memory_consumption(component_to_base_table) +
1255 MemoryConsumption::memory_consumption(restriction_is_additive_flags) +
1256 MemoryConsumption::memory_consumption(nonzero_components) +
1257 MemoryConsumption::memory_consumption(n_nonzero_components_table));
1258}
1259
1260
1261
1262template <int dim, int spacedim>
1263std::vector<unsigned int>
1265 const std::vector<ComponentMask> &nonzero_components)
1266{
1267 std::vector<unsigned int> retval(nonzero_components.size());
1268 for (unsigned int i = 0; i < nonzero_components.size(); ++i)
1269 retval[i] = nonzero_components[i].n_selected_components();
1270 return retval;
1271}
1272
1273
1274
1275/*------------------------------- FiniteElement ----------------------*/
1276
1277#ifndef DOXYGEN
1278template <int dim, int spacedim>
1279std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1281 const UpdateFlags flags,
1282 const Mapping<dim, spacedim> & mapping,
1283 const hp::QCollection<dim - 1> &quadrature,
1285 spacedim>
1286 &output_data) const
1287{
1288 return get_data(flags,
1289 mapping,
1291 quadrature),
1292 output_data);
1293}
1294
1295
1296
1297template <int dim, int spacedim>
1298std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1300 const UpdateFlags flags,
1301 const Mapping<dim, spacedim> &mapping,
1302 const Quadrature<dim - 1> & quadrature,
1304 spacedim>
1305 &output_data) const
1306{
1307 return get_data(flags,
1308 mapping,
1310 quadrature),
1311 output_data);
1312}
1313
1314
1315
1316template <int dim, int spacedim>
1317inline void
1320 const unsigned int face_no,
1321 const hp::QCollection<dim - 1> & quadrature,
1322 const Mapping<dim, spacedim> & mapping,
1323 const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
1325 & mapping_data,
1326 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1328 spacedim>
1329 &output_data) const
1330{
1331 // base class version, implement overridden function in derived classes
1332 AssertDimension(quadrature.size(), 1);
1333 fill_fe_face_values(cell,
1334 face_no,
1335 quadrature[0],
1336 mapping,
1337 mapping_internal,
1338 mapping_data,
1339 fe_internal,
1340 output_data);
1341}
1342
1343
1344
1345template <int dim, int spacedim>
1346inline void
1349 const unsigned int face_no,
1350 const Quadrature<dim - 1> & quadrature,
1351 const Mapping<dim, spacedim> & mapping,
1352 const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
1354 & mapping_data,
1355 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1357 spacedim>
1358 &output_data) const
1359{
1360 Assert(false,
1361 ExcMessage("Use of a deprecated interface, please implement "
1362 "fill_fe_face_values taking a hp::QCollection argument"));
1363 (void)cell;
1364 (void)face_no;
1365 (void)quadrature;
1366 (void)mapping;
1367 (void)mapping_internal;
1368 (void)mapping_data;
1369 (void)fe_internal;
1370 (void)output_data;
1371}
1372#endif
1373
1374
1375
1376template <int dim, int spacedim>
1377std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1379 const UpdateFlags flags,
1380 const Mapping<dim, spacedim> &mapping,
1381 const Quadrature<dim - 1> & quadrature,
1383 spacedim>
1384 &output_data) const
1385{
1386 return get_data(flags,
1387 mapping,
1389 this->reference_cell(), quadrature),
1390 output_data);
1391}
1392
1393
1394
1395template <int dim, int spacedim>
1397FiniteElement<dim, spacedim>::base_element(const unsigned int index) const
1398{
1399 (void)index;
1400 AssertIndexRange(index, 1);
1401 // This function should not be
1402 // called for a system element
1403 Assert(base_to_block_indices.size() == 1, ExcInternalError());
1404 return *this;
1405}
1406
1407
1408
1409/*------------------------------- Explicit Instantiations -------------*/
1410#include "fe.inst"
1411
1412
bool represents_the_all_selected_mask() const
unsigned int size() const
bool represents_the_all_selected_mask() const
unsigned int size() const
const unsigned int components
Definition fe_data.h:447
unsigned int n_dofs_per_cell() const
unsigned int n_dofs_per_line() const
unsigned int n_components() 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:2588
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
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
const bool cached_primitivity
Definition fe.h:2595
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
bool has_generalized_support_points() const
virtual Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const
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
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
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:2491
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
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:2527
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)
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:2570
std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > component_to_base_table
Definition fe.h:2563
TableIndices< 2 > interface_constraints_size() const
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
FullMatrix< double > interface_constraints
Definition fe.h:2428
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
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:2579
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
Abstract base class for mapping classes.
Definition mapping.h:317
Definition point.h:112
unsigned int size() const
Definition collection.h:265
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Point< 2 > first
Definition grid_out.cc:4615
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
UpdateFlags
@ update_default
No update.
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition operators.h:50
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 > > &line_support_points, const std::vector< unsigned int > &renumbering)
static const unsigned int invalid_unsigned_int
Definition types.h:213
STL namespace.
static unsigned int n_children(const RefinementCase< dim > &refinement_case)