Reference documentation for deal.II version GIT relicensing-437-g81ec864850 2024-04-19 07:30:02+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1998 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
18
20
21#include <deal.II/fe/fe.h>
23#include <deal.II/fe/mapping.h>
24
25#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#ifndef DOXYGEN
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)
65 std::make_pair(std::make_pair(0U, 0U), 0U))
66 ,
67
68 // Special handling of vectors of length one: in this case, we
69 // assume that all entries were supposed to be equal
71 r_i_a_f.size() == 1 ?
72 std::vector<bool>(fe_data.n_dofs_per_cell(), r_i_a_f[0]) :
73 r_i_a_f)
75 nonzero_c.size() == 1 ?
76 std::vector<ComponentMask>(fe_data.n_dofs_per_cell(), nonzero_c[0]) :
77 nonzero_c)
81 [](const unsigned int n_components) {
82 return n_components != 1U;
83 }) == n_nonzero_components_table.end())
84{
85 Assert(restriction_is_additive_flags.size() == this->n_dofs_per_cell(),
86 ExcDimensionMismatch(restriction_is_additive_flags.size(),
87 this->n_dofs_per_cell()));
88 AssertDimension(nonzero_components.size(), this->n_dofs_per_cell());
89 for (unsigned int i = 0; i < nonzero_components.size(); ++i)
90 {
91 Assert(nonzero_components[i].size() == this->n_components(),
93 Assert(nonzero_components[i].n_selected_components() >= 1,
95 Assert(n_nonzero_components_table[i] >= 1, ExcInternalError());
96 Assert(n_nonzero_components_table[i] <= this->n_components(),
98 }
99
100 // initialize some tables in the default way, i.e. if there is only one
101 // (vector-)component; if the element is not primitive, leave these tables
102 // empty.
103 if (this->is_primitive())
104 {
105 system_to_component_table.resize(this->n_dofs_per_cell());
106 for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
107 system_to_component_table[j] = std::pair<unsigned, unsigned>(0, j);
108
109 face_system_to_component_table.resize(this->n_unique_faces());
110 for (unsigned int f = 0; f < this->n_unique_faces(); ++f)
111 {
112 face_system_to_component_table[f].resize(this->n_dofs_per_face(f));
113 for (unsigned int j = 0; j < this->n_dofs_per_face(f); ++j)
114 face_system_to_component_table[f][j] =
115 std::pair<unsigned, unsigned>(0, j);
116 }
117 }
118
119 for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
120 system_to_base_table[j] = std::make_pair(std::make_pair(0U, 0U), j);
121
122 face_system_to_base_table.resize(this->n_unique_faces());
123 for (unsigned int f = 0; f < this->n_unique_faces(); ++f)
124 {
125 face_system_to_base_table[f].resize(this->n_dofs_per_face(f));
126 for (unsigned int j = 0; j < this->n_dofs_per_face(f); ++j)
127 face_system_to_base_table[f][j] =
128 std::make_pair(std::make_pair(0U, 0U), j);
129 }
130
131 // Fill with default value; may be changed by constructor of derived class.
132 base_to_block_indices.reinit(1, 1);
133
134 // initialize the restriction and prolongation matrices. the default
135 // constructor of FullMatrix<dim> initializes them with size zero
138 for (const unsigned int ref_case :
139 RefinementCase<dim>::all_refinement_cases())
140 if (ref_case != RefinementCase<dim>::no_refinement)
141 {
142 prolongation[ref_case - 1].resize(GeometryInfo<dim>::n_children(
143 RefinementCase<dim>(ref_case)),
145 restriction[ref_case - 1].resize(GeometryInfo<dim>::n_children(
146 RefinementCase<dim>(ref_case)),
148 }
149
150
151 if (dim == 3)
152 {
153 adjust_quad_dof_index_for_face_orientation_table.resize(
154 this->n_unique_2d_subobjects());
155
156 for (unsigned int f = 0; f < this->n_unique_2d_subobjects(); ++f)
157 {
158 adjust_quad_dof_index_for_face_orientation_table[f] =
159 Table<2, int>(this->n_dofs_per_quad(f),
160 this->reference_cell().n_face_orientations(f));
161 adjust_quad_dof_index_for_face_orientation_table[f].fill(0);
162 }
163 }
164
165 unit_face_support_points.resize(this->n_unique_faces());
166 generalized_face_support_points.resize(this->n_unique_faces());
167}
168
169
170
171template <int dim, int spacedim>
172std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>
173FiniteElement<dim, spacedim>::operator^(const unsigned int multiplicity) const
174{
175 return {this->clone(), multiplicity};
176}
177
178
179
180template <int dim, int spacedim>
181double
183 const Point<dim> &) const
184{
185 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
186 return 0.;
187}
188
189
190
191template <int dim, int spacedim>
192double
194 const Point<dim> &,
195 const unsigned int) const
196{
197 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
198 return 0.;
199}
200
201
202
203template <int dim, int spacedim>
206 const Point<dim> &) const
207{
208 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
209 return Tensor<1, dim>();
210}
211
212
213
214template <int dim, int spacedim>
217 const Point<dim> &,
218 const unsigned int) const
219{
220 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
221 return Tensor<1, dim>();
222}
223
224
225
226template <int dim, int spacedim>
229 const Point<dim> &) const
230{
231 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
232 return Tensor<2, dim>();
233}
234
235
236
237template <int dim, int spacedim>
240 const unsigned int,
241 const Point<dim> &,
242 const unsigned int) const
243{
244 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
245 return Tensor<2, dim>();
246}
247
248
249
250template <int dim, int spacedim>
253 const Point<dim> &) const
254{
255 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
256 return Tensor<3, dim>();
257}
258
259
260
261template <int dim, int spacedim>
264 const unsigned int,
265 const Point<dim> &,
266 const unsigned int) const
267{
268 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
269 return Tensor<3, dim>();
270}
271
272
273
274template <int dim, int spacedim>
277 const Point<dim> &) const
278{
279 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
280 return Tensor<4, dim>();
281}
282
283
284
285template <int dim, int spacedim>
288 const unsigned int,
289 const Point<dim> &,
290 const unsigned int) const
291{
292 AssertThrow(false, ExcUnitShapeValuesDoNotExist());
293 return Tensor<4, dim>();
294}
295
296
297template <int dim, int spacedim>
298void
300 const bool isotropic_restriction_only,
301 const bool isotropic_prolongation_only)
302{
303 for (const unsigned int ref_case :
304 RefinementCase<dim>::all_refinement_cases())
305 if (ref_case != RefinementCase<dim>::no_refinement)
306 {
307 const unsigned int nc =
309
310 for (unsigned int i = 0; i < nc; ++i)
311 {
312 if (this->restriction[ref_case - 1][i].m() !=
313 this->n_dofs_per_cell() &&
314 (!isotropic_restriction_only ||
316 this->restriction[ref_case - 1][i].reinit(
317 this->n_dofs_per_cell(), this->n_dofs_per_cell());
318 if (this->prolongation[ref_case - 1][i].m() !=
319 this->n_dofs_per_cell() &&
320 (!isotropic_prolongation_only ||
322 this->prolongation[ref_case - 1][i].reinit(
323 this->n_dofs_per_cell(), this->n_dofs_per_cell());
324 }
325 }
326}
327
328
329template <int dim, int spacedim>
330const FullMatrix<double> &
332 const unsigned int child,
333 const RefinementCase<dim> &refinement_case) const
334{
335 AssertIndexRange(refinement_case,
339 "Restriction matrices are only available for refined cells!"));
341 child, GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)));
342 // we use refinement_case-1 here. the -1 takes care of the origin of the
343 // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
344 // available and so the vector indices are shifted
345 Assert(restriction[refinement_case - 1][child].n() == this->n_dofs_per_cell(),
346 ExcProjectionVoid());
347 return restriction[refinement_case - 1][child];
348}
349
350
351
352template <int dim, int spacedim>
353const FullMatrix<double> &
355 const unsigned int child,
356 const RefinementCase<dim> &refinement_case) const
357{
358 AssertIndexRange(refinement_case,
362 "Prolongation matrices are only available for refined cells!"));
364 child, GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)));
365 // we use refinement_case-1 here. the -1 takes care
366 // of the origin of the vector, as for
367 // RefinementCase::no_refinement (=0) there is no
368 // data available and so the vector indices
369 // are shifted
370 Assert(prolongation[refinement_case - 1][child].n() ==
371 this->n_dofs_per_cell(),
372 ExcEmbeddingVoid());
373 return prolongation[refinement_case - 1][child];
374}
375
376
377// TODO:[GK] This is probably not the most efficient way of doing this.
378template <int dim, int spacedim>
379unsigned int
381 const unsigned int index) const
382{
383 AssertIndexRange(index, this->n_components());
384
385 return first_block_of_base(component_to_base_table[index].first.first) +
386 component_to_base_table[index].second;
387}
388
389
390template <int dim, int spacedim>
393 const FEValuesExtractors::Scalar &scalar) const
394{
395 AssertIndexRange(scalar.component, this->n_components());
396
397 // TODO: it would be nice to verify that it is indeed possible
398 // to select this scalar component, i.e., that it is not part
399 // of a non-primitive element. unfortunately, there is no simple
400 // way to write such a condition...
401
402 std::vector<bool> mask(this->n_components(), false);
403 mask[scalar.component] = true;
404 return ComponentMask(mask);
405}
406
407
408template <int dim, int spacedim>
411 const FEValuesExtractors::Vector &vector) const
412{
414 this->n_components());
415
416 // TODO: it would be nice to verify that it is indeed possible
417 // to select these vector components, i.e., that they don't span
418 // beyond the beginning or end of anon-primitive element.
419 // unfortunately, there is no simple way to write such a condition...
420
421 std::vector<bool> mask(this->n_components(), false);
422 for (unsigned int c = vector.first_vector_component;
423 c < vector.first_vector_component + dim;
424 ++c)
425 mask[c] = true;
426 return ComponentMask(mask);
427}
428
429
430template <int dim, int spacedim>
433 const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const
434{
437 this->n_components());
438
439 // TODO: it would be nice to verify that it is indeed possible
440 // to select these vector components, i.e., that they don't span
441 // beyond the beginning or end of anon-primitive element.
442 // unfortunately, there is no simple way to write such a condition...
443
444 std::vector<bool> mask(this->n_components(), false);
445 for (unsigned int c = sym_tensor.first_tensor_component;
446 c < sym_tensor.first_tensor_component +
448 ++c)
449 mask[c] = true;
450 return ComponentMask(mask);
451}
452
453
454
455template <int dim, int spacedim>
458{
459 // if we get a block mask that represents all blocks, then
460 // do the same for the returned component mask
461 if (block_mask.represents_the_all_selected_mask())
462 return {};
463
464 AssertDimension(block_mask.size(), this->n_blocks());
465
466 std::vector<bool> component_mask(this->n_components(), false);
467 for (unsigned int c = 0; c < this->n_components(); ++c)
468 if (block_mask[component_to_block_index(c)] == true)
469 component_mask[c] = true;
470
471 return ComponentMask(component_mask);
472}
473
474
475
476template <int dim, int spacedim>
479 const FEValuesExtractors::Scalar &scalar) const
480{
481 // simply create the corresponding component mask (a simpler
482 // process) and then convert it to a block mask
483 return block_mask(component_mask(scalar));
484}
485
486
487template <int dim, int spacedim>
490 const FEValuesExtractors::Vector &vector) const
491{
492 // simply create the corresponding component mask (a simpler
493 // process) and then convert it to a block mask
494 return block_mask(component_mask(vector));
495}
496
497
498template <int dim, int spacedim>
501 const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const
502{
503 // simply create the corresponding component mask (a simpler
504 // process) and then convert it to a block mask
505 return block_mask(component_mask(sym_tensor));
506}
507
508
509
510template <int dim, int spacedim>
513 const ComponentMask &component_mask) const
514{
515 // if we get a component mask that represents all component, then
516 // do the same for the returned block mask
517 if (component_mask.represents_the_all_selected_mask())
518 return {};
519
520 AssertDimension(component_mask.size(), this->n_components());
521
522 // walk over all of the components
523 // of this finite element and see
524 // if we need to set the
525 // corresponding block. inside the
526 // block, walk over all the
527 // components that correspond to
528 // this block and make sure the
529 // component mask is set for all of
530 // them
531 std::vector<bool> block_mask(this->n_blocks(), false);
532 for (unsigned int c = 0; c < this->n_components();)
533 {
534 const unsigned int block = component_to_block_index(c);
535 if (component_mask[c] == true)
536 block_mask[block] = true;
537
538 // now check all of the other
539 // components that correspond
540 // to this block
541 ++c;
542 while ((c < this->n_components()) &&
543 (component_to_block_index(c) == block))
544 {
545 Assert(component_mask[c] == block_mask[block],
547 "The component mask argument given to this function "
548 "is not a mask where the individual components belonging "
549 "to one block of the finite element are either all "
550 "selected or not selected. You can't call this function "
551 "with a component mask that splits blocks."));
552 ++c;
553 }
554 }
555
556
557 return BlockMask(block_mask);
558}
559
560
561
562template <int dim, int spacedim>
563unsigned int
565 const unsigned int face_index,
566 const unsigned int face,
567 const unsigned char combined_orientation) const
568{
569 AssertIndexRange(face_index, this->n_dofs_per_face(face));
570 AssertIndexRange(face, this->reference_cell().n_faces());
571
572 // see the function's documentation for an explanation of this
573 // assertion -- in essence, derived classes have to implement
574 // an overloaded version of this function if we are to use any
575 // other than default (standard) orientation
576 if (combined_orientation !=
578 Assert((this->n_dofs_per_line() <= 1) && (this->n_dofs_per_quad(face) <= 1),
580 "The function in this base class can not handle this case. "
581 "Rather, the derived class you are using must provide "
582 "an overloaded version but apparently hasn't done so. See "
583 "the documentation of this function for more information."));
584
585 // we need to distinguish between DoFs on vertices, lines and (in 3d) quads.
586 // do so in a sequence of if-else statements
587 if (face_index < this->get_first_face_line_index(face))
588 // DoF is on a vertex
589 {
590 // get the number of the vertex on the face that corresponds to this DoF,
591 // along with the number of the DoF on this vertex
592 const unsigned int face_vertex = face_index / this->n_dofs_per_vertex();
593 const unsigned int dof_index_on_vertex =
594 face_index % this->n_dofs_per_vertex();
595
596 // then get the number of this vertex on the cell and translate
597 // this to a DoF number on the cell
598 return (this->reference_cell().face_to_cell_vertices(
599 face, face_vertex, combined_orientation) *
600 this->n_dofs_per_vertex() +
601 dof_index_on_vertex);
602 }
603 else if (face_index < this->get_first_face_quad_index(face))
604 // DoF is on a face
605 {
606 // do the same kind of translation as before. we need to only consider
607 // DoFs on the lines, i.e., ignoring those on the vertices
608 const unsigned int index =
609 face_index - this->get_first_face_line_index(face);
610
611 const unsigned int face_line = index / this->n_dofs_per_line();
612 const unsigned int dof_index_on_line = index % this->n_dofs_per_line();
613
614 return (this->get_first_line_index() +
615 this->reference_cell().face_to_cell_lines(face,
616 face_line,
617 combined_orientation) *
618 this->n_dofs_per_line() +
619 dof_index_on_line);
620 }
621 else
622 // DoF is on a quad
623 {
624 Assert(dim >= 3, ExcInternalError());
625
626 // ignore vertex and line dofs
627 const unsigned int index =
628 face_index - this->get_first_face_quad_index(face);
629
630 return (this->get_first_quad_index(face) + index);
631 }
632}
633
634
635
636template <int dim, int spacedim>
637unsigned int
639 const unsigned int index,
640 const unsigned int face,
641 const unsigned char combined_orientation) const
642{
643 // general template for 1d and 2d: not
644 // implemented. in fact, the function
645 // shouldn't even be called unless we are
646 // in 3d, so throw an internal error
647 Assert(dim == 3, ExcInternalError());
648 if (dim < 3)
649 return index;
650
651 // adjust dofs on 3d faces if the face is
652 // flipped. note that we query a table that
653 // derived elements need to have set up
654 // front. the exception are discontinuous
655 // elements for which there should be no
656 // face dofs anyway (i.e. dofs_per_quad==0
657 // in 3d), so we don't need the table, but
658 // the function should also not have been
659 // called
660 AssertIndexRange(index, this->n_dofs_per_quad(face));
661 const auto table_n = this->n_unique_2d_subobjects() == 1 ? 0 : face;
662 Assert(
663 adjust_quad_dof_index_for_face_orientation_table[table_n].n_elements() ==
664 (this->reference_cell().n_face_orientations(face)) *
665 this->n_dofs_per_quad(face),
667 return index + adjust_quad_dof_index_for_face_orientation_table[table_n](
668 index, combined_orientation);
669}
670
671
672
673template <int dim, int spacedim>
674unsigned int
676 const unsigned int index,
677 const unsigned char combined_orientation) const
678{
679 Assert(combined_orientation ==
681 combined_orientation ==
684
685 AssertIndexRange(index, this->n_dofs_per_line());
686 Assert(adjust_line_dof_index_for_line_orientation_table.size() ==
687 this->n_dofs_per_line(),
689 if (combined_orientation ==
691 return index;
692 else
693 return index + adjust_line_dof_index_for_line_orientation_table[index];
694}
695
696
697
698template <int dim, int spacedim>
699bool
701{
702 for (const unsigned int ref_case :
703 RefinementCase<dim>::all_refinement_cases())
704 if (ref_case != RefinementCase<dim>::no_refinement)
705 for (unsigned int c = 0;
706 c < GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
707 ++c)
708 {
709 // make sure also the lazily initialized matrices are created
710 get_prolongation_matrix(c, RefinementCase<dim>(ref_case));
711 Assert((prolongation[ref_case - 1][c].m() ==
712 this->n_dofs_per_cell()) ||
713 (prolongation[ref_case - 1][c].m() == 0),
715 Assert((prolongation[ref_case - 1][c].n() ==
716 this->n_dofs_per_cell()) ||
717 (prolongation[ref_case - 1][c].n() == 0),
719 if ((prolongation[ref_case - 1][c].m() == 0) ||
720 (prolongation[ref_case - 1][c].n() == 0))
721 return false;
722 }
723 return true;
724}
725
726
727
728template <int dim, int spacedim>
729bool
731{
732 for (const unsigned int ref_case :
733 RefinementCase<dim>::all_refinement_cases())
734 if (ref_case != RefinementCase<dim>::no_refinement)
735 for (unsigned int c = 0;
736 c < GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
737 ++c)
738 {
739 // make sure also the lazily initialized matrices are created
740 get_restriction_matrix(c, RefinementCase<dim>(ref_case));
741 Assert((restriction[ref_case - 1][c].m() ==
742 this->n_dofs_per_cell()) ||
743 (restriction[ref_case - 1][c].m() == 0),
745 Assert((restriction[ref_case - 1][c].n() ==
746 this->n_dofs_per_cell()) ||
747 (restriction[ref_case - 1][c].n() == 0),
749 if ((restriction[ref_case - 1][c].m() == 0) ||
750 (restriction[ref_case - 1][c].n() == 0))
751 return false;
752 }
753 return true;
754}
755
756
757
758template <int dim, int spacedim>
759bool
761{
762 const RefinementCase<dim> ref_case =
764
765 for (unsigned int c = 0;
766 c < GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
767 ++c)
768 {
769 // make sure also the lazily initialized matrices are created
770 get_prolongation_matrix(c, RefinementCase<dim>(ref_case));
771 Assert((prolongation[ref_case - 1][c].m() == this->n_dofs_per_cell()) ||
772 (prolongation[ref_case - 1][c].m() == 0),
774 Assert((prolongation[ref_case - 1][c].n() == this->n_dofs_per_cell()) ||
775 (prolongation[ref_case - 1][c].n() == 0),
777 if ((prolongation[ref_case - 1][c].m() == 0) ||
778 (prolongation[ref_case - 1][c].n() == 0))
779 return false;
780 }
781 return true;
782}
783
784
785
786template <int dim, int spacedim>
787bool
789{
790 const RefinementCase<dim> ref_case =
792
793 for (unsigned int c = 0;
794 c < GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
795 ++c)
796 {
797 // make sure also the lazily initialized matrices are created
798 get_restriction_matrix(c, RefinementCase<dim>(ref_case));
799 Assert((restriction[ref_case - 1][c].m() == this->n_dofs_per_cell()) ||
800 (restriction[ref_case - 1][c].m() == 0),
802 Assert((restriction[ref_case - 1][c].n() == this->n_dofs_per_cell()) ||
803 (restriction[ref_case - 1][c].n() == 0),
805 if ((restriction[ref_case - 1][c].m() == 0) ||
806 (restriction[ref_case - 1][c].n() == 0))
807 return false;
808 }
809 return true;
810}
811
812
813
814template <int dim, int spacedim>
815bool
817 const internal::SubfaceCase<dim> &subface_case) const
818{
820 {
821 unsigned int n_dofs_on_faces = 0;
822
823 for (const auto face_no : this->reference_cell().face_indices())
824 n_dofs_on_faces += this->n_dofs_per_face(face_no);
825
826 return (n_dofs_on_faces == 0) || (interface_constraints.m() != 0);
827 }
828 else
829 return false;
830}
831
832
833
834template <int dim, int spacedim>
835bool
837{
838 return false;
839}
840
841
842
843template <int dim, int spacedim>
844const FullMatrix<double> &
846 const internal::SubfaceCase<dim> &subface_case) const
847{
848 // TODO: the implementation makes the assumption that all faces have the
849 // same number of dofs
850 AssertDimension(this->n_unique_faces(), 1);
851 const unsigned int face_no = 0;
852 (void)face_no;
853
854 (void)subface_case;
856 ExcMessage("Constraints for this element are only implemented "
857 "for the case that faces are refined isotropically "
858 "(which is always the case in 2d, and in 3d requires "
859 "that the neighboring cell of a coarse cell presents "
860 "exactly four children on the common face)."));
861 Assert((this->n_dofs_per_face(face_no) == 0) ||
862 (interface_constraints.m() != 0),
863 ExcMessage("The finite element for which you try to obtain "
864 "hanging node constraints does not appear to "
865 "implement them."));
866
867 if (dim == 1)
868 Assert((interface_constraints.m() == 0) && (interface_constraints.n() == 0),
869 ExcWrongInterfaceMatrixSize(interface_constraints.m(),
870 interface_constraints.n()));
871
872 return interface_constraints;
873}
874
875
876
877template <int dim, int spacedim>
880{
881 // TODO: the implementation makes the assumption that all faces have the
882 // same number of dofs
883 AssertDimension(this->n_unique_faces(), 1);
884 const unsigned int face_no = 0;
885
886 switch (dim)
887 {
888 case 1:
889 return {0U, 0U};
890 case 2:
891 return {this->n_dofs_per_vertex() + 2 * this->n_dofs_per_line(),
892 this->n_dofs_per_face(face_no)};
893 case 3:
894 return {5 * this->n_dofs_per_vertex() + 12 * this->n_dofs_per_line() +
895 4 * this->n_dofs_per_quad(face_no),
896 this->n_dofs_per_face(face_no)};
897 default:
899 }
901}
902
903
904
905template <int dim, int spacedim>
906void
909 FullMatrix<double> &) const
910{
911 // by default, no interpolation
912 // implemented. so throw exception,
913 // as documentation says
915 false,
917}
918
919
920
921template <int dim, int spacedim>
922void
926 const unsigned int) 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
942 const unsigned int,
944 const unsigned int) const
945{
946 // by default, no interpolation
947 // implemented. so throw exception,
948 // as documentation says
950 false,
952}
953
954
955
956template <int dim, int spacedim>
957std::vector<std::pair<unsigned int, unsigned int>>
959 const FiniteElement<dim, spacedim> &) const
960{
962 return std::vector<std::pair<unsigned int, unsigned int>>();
963}
964
965
966
967template <int dim, int spacedim>
968std::vector<std::pair<unsigned int, unsigned int>>
970 const FiniteElement<dim, spacedim> &) const
971{
973 return std::vector<std::pair<unsigned int, unsigned int>>();
974}
975
976
977
978template <int dim, int spacedim>
979std::vector<std::pair<unsigned int, unsigned int>>
982 const unsigned int) const
983{
985 return std::vector<std::pair<unsigned int, unsigned int>>();
986}
987
988
989
990template <int dim, int spacedim>
994 const unsigned int) const
995{
998}
999
1000
1001
1002template <int dim, int spacedim>
1003bool
1005 const FiniteElement<dim, spacedim> &f) const
1006{
1007 // Compare fields in roughly increasing order of how expensive the
1008 // comparison is
1009 return ((typeid(*this) == typeid(f)) && (this->get_name() == f.get_name()) &&
1010 (static_cast<const FiniteElementData<dim> &>(*this) ==
1011 static_cast<const FiniteElementData<dim> &>(f)) &&
1012 (interface_constraints == f.interface_constraints));
1013}
1014
1015
1016
1017template <int dim, int spacedim>
1018bool
1020 const FiniteElement<dim, spacedim> &f) const
1021{
1022 return !(*this == f);
1023}
1024
1025
1026
1027template <int dim, int spacedim>
1028const std::vector<Point<dim>> &
1030{
1031 // a finite element may define
1032 // support points, but only if
1033 // there are as many as there are
1034 // degrees of freedom
1035 Assert((unit_support_points.empty()) ||
1036 (unit_support_points.size() == this->n_dofs_per_cell()),
1038 return unit_support_points;
1039}
1040
1041
1042
1043template <int dim, int spacedim>
1044bool
1046{
1047 return (unit_support_points.size() != 0);
1048}
1049
1050
1051
1052template <int dim, int spacedim>
1053const std::vector<Point<dim>> &
1055{
1056 // If the finite element implements generalized support points, return
1057 // those. Otherwise fall back to unit support points.
1058 return ((generalized_support_points.empty()) ? unit_support_points :
1059 generalized_support_points);
1060}
1061
1062
1063
1064template <int dim, int spacedim>
1065bool
1067{
1068 return (get_generalized_support_points().size() != 0);
1069}
1070
1071
1072
1073template <int dim, int spacedim>
1075FiniteElement<dim, spacedim>::unit_support_point(const unsigned int index) const
1076{
1077 AssertIndexRange(index, this->n_dofs_per_cell());
1078 Assert(unit_support_points.size() == this->n_dofs_per_cell(),
1079 ExcFEHasNoSupportPoints());
1080 return unit_support_points[index];
1081}
1082
1083
1084
1085template <int dim, int spacedim>
1086const std::vector<Point<dim - 1>> &
1088 const unsigned int face_no) const
1089{
1090 // a finite element may define
1091 // support points, but only if
1092 // there are as many as there are
1093 // degrees of freedom on a face
1094 Assert((unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
1095 .empty()) ||
1096 (unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
1097 .size() == this->n_dofs_per_face(face_no)),
1099 return unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no];
1100}
1101
1102
1103
1104template <int dim, int spacedim>
1105bool
1107 const unsigned int face_no) const
1108{
1109 return (unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
1110 .size() != 0);
1111}
1112
1113
1114
1115template <int dim, int spacedim>
1116Point<dim - 1>
1118 const unsigned int index,
1119 const unsigned int face_no) const
1120{
1121 AssertIndexRange(index, this->n_dofs_per_face(face_no));
1122 Assert(unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
1123 .size() == this->n_dofs_per_face(face_no),
1124 ExcFEHasNoSupportPoints());
1125 return unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
1126 [index];
1127}
1128
1129
1130
1131template <int dim, int spacedim>
1132bool
1134 const unsigned int) const
1135{
1136 return true;
1137}
1138
1139
1140
1141template <int dim, int spacedim>
1144{
1145 // Translate the ComponentMask into first_selected and n_components after
1146 // some error checking:
1147 const unsigned int n_total_components = this->n_components();
1148 Assert((n_total_components == mask.size()) || (mask.size() == 0),
1149 ExcMessage("The given ComponentMask has the wrong size."));
1150
1151 const unsigned int n_selected =
1152 mask.n_selected_components(n_total_components);
1153 Assert(n_selected > 0,
1154 ExcMessage("You need at least one selected component."));
1155
1156 const unsigned int first_selected =
1157 mask.first_selected_component(n_total_components);
1158
1159# ifdef DEBUG
1160 // check that it is contiguous:
1161 for (unsigned int c = 0; c < n_total_components; ++c)
1162 Assert((c < first_selected && (!mask[c])) ||
1163 (c >= first_selected && c < first_selected + n_selected &&
1164 mask[c]) ||
1165 (c >= first_selected + n_selected && !mask[c]),
1166 ExcMessage("Error: the given ComponentMask is not contiguous!"));
1167# endif
1168
1169 return get_sub_fe(first_selected, n_selected);
1170}
1171
1172
1173
1174template <int dim, int spacedim>
1177 const unsigned int first_component,
1178 const unsigned int n_selected_components) const
1179{
1180 // No complicated logic is needed here, because it is overridden in
1181 // FESystem<dim,spacedim>. Just make sure that what the user chose is valid:
1182 Assert(first_component == 0 && n_selected_components == this->n_components(),
1183 ExcMessage(
1184 "You can only select a whole FiniteElement, not a part of one."));
1185
1186 (void)first_component;
1187 (void)n_selected_components;
1188 return *this;
1189}
1190
1191
1192
1193template <int dim, int spacedim>
1194std::pair<Table<2, bool>, std::vector<unsigned int>>
1196{
1198 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
1199 Table<2, bool>(this->n_components(), this->n_dofs_per_cell()),
1200 std::vector<unsigned int>(this->n_components()));
1201}
1202
1203
1204
1205template <int dim, int spacedim>
1206void
1209 const std::vector<Vector<double>> &,
1210 std::vector<double> &) const
1211{
1212 Assert(has_generalized_support_points(),
1213 ExcMessage("The element for which you are calling the current "
1214 "function does not have generalized support points (see "
1215 "the glossary for a definition of generalized support "
1216 "points). Consequently, the current function can not "
1217 "be defined and is not implemented by the element."));
1219}
1220
1221
1222
1223template <int dim, int spacedim>
1224std::size_t
1226{
1227 return (
1228 sizeof(FiniteElementData<dim>) +
1231 MemoryConsumption::memory_consumption(interface_constraints) +
1232 MemoryConsumption::memory_consumption(system_to_component_table) +
1233 MemoryConsumption::memory_consumption(face_system_to_component_table) +
1234 MemoryConsumption::memory_consumption(system_to_base_table) +
1235 MemoryConsumption::memory_consumption(face_system_to_base_table) +
1236 MemoryConsumption::memory_consumption(component_to_base_table) +
1237 MemoryConsumption::memory_consumption(restriction_is_additive_flags) +
1238 MemoryConsumption::memory_consumption(nonzero_components) +
1239 MemoryConsumption::memory_consumption(n_nonzero_components_table));
1240}
1241
1242
1243
1244template <int dim, int spacedim>
1245std::vector<unsigned int>
1247 const std::vector<ComponentMask> &nonzero_components)
1248{
1249 std::vector<unsigned int> retval(nonzero_components.size());
1250 for (unsigned int i = 0; i < nonzero_components.size(); ++i)
1251 retval[i] = nonzero_components[i].n_selected_components();
1252 return retval;
1253}
1254
1255
1256
1257/*------------------------------- FiniteElement ----------------------*/
1258
1259# ifndef DOXYGEN
1260template <int dim, int spacedim>
1261std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1263 const UpdateFlags flags,
1264 const Mapping<dim, spacedim> &mapping,
1265 const hp::QCollection<dim - 1> &quadrature,
1267 spacedim>
1268 &output_data) const
1269{
1270 return get_data(flags,
1271 mapping,
1273 quadrature),
1274 output_data);
1275}
1276
1277
1278
1279template <int dim, int spacedim>
1280std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1282 const UpdateFlags flags,
1283 const Mapping<dim, spacedim> &mapping,
1284 const Quadrature<dim - 1> &quadrature,
1286 spacedim>
1287 &output_data) const
1288{
1289 return get_data(flags,
1290 mapping,
1292 quadrature),
1293 output_data);
1294}
1295
1296
1297
1298template <int dim, int spacedim>
1299inline void
1302 const unsigned int face_no,
1303 const hp::QCollection<dim - 1> &quadrature,
1304 const Mapping<dim, spacedim> &mapping,
1305 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
1307 &mapping_data,
1308 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1310 spacedim>
1311 &output_data) const
1312{
1313 // base class version, implement overridden function in derived classes
1314 AssertDimension(quadrature.size(), 1);
1315 fill_fe_face_values(cell,
1316 face_no,
1317 quadrature[0],
1318 mapping,
1319 mapping_internal,
1320 mapping_data,
1321 fe_internal,
1322 output_data);
1323}
1324
1325
1326
1327template <int dim, int spacedim>
1328inline void
1331 const unsigned int face_no,
1332 const Quadrature<dim - 1> &quadrature,
1333 const Mapping<dim, spacedim> &mapping,
1334 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
1336 &mapping_data,
1337 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1339 spacedim>
1340 &output_data) const
1341{
1342 Assert(false,
1343 ExcMessage("Use of a deprecated interface, please implement "
1344 "fill_fe_face_values taking a hp::QCollection argument"));
1345 (void)cell;
1346 (void)face_no;
1347 (void)quadrature;
1348 (void)mapping;
1349 (void)mapping_internal;
1350 (void)mapping_data;
1351 (void)fe_internal;
1352 (void)output_data;
1353}
1354# endif
1355
1356
1357
1358template <int dim, int spacedim>
1359std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
1361 const UpdateFlags flags,
1362 const Mapping<dim, spacedim> &mapping,
1363 const Quadrature<dim - 1> &quadrature,
1365 spacedim>
1366 &output_data) const
1367{
1368 return get_data(flags,
1369 mapping,
1371 this->reference_cell(), quadrature),
1372 output_data);
1373}
1374
1375
1376
1377template <int dim, int spacedim>
1379FiniteElement<dim, spacedim>::base_element(const unsigned int index) const
1380{
1381 (void)index;
1382 AssertIndexRange(index, 1);
1383 // This function should not be
1384 // called for a system element
1385 Assert(base_to_block_indices.size() == 1, ExcInternalError());
1386 return *this;
1387}
1388
1389
1390#endif
1391/*------------------------------- Explicit Instantiations -------------*/
1392#include "fe.inst"
1393
1394
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:446
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:2582
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:2589
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
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const unsigned char combined_orientation=ReferenceCell::default_combined_face_orientation()) 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
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:2485
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) 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:2521
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)
unsigned int adjust_quad_dof_index_for_face_orientation(const unsigned int index, const unsigned int face_no, const unsigned char combined_orientation) const
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:2564
unsigned int adjust_line_dof_index_for_line_orientation(const unsigned int index, const unsigned char combined_orientation) const
std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > component_to_base_table
Definition fe.h:2557
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:2424
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
const std::vector< ComponentMask > nonzero_components
Definition fe.h:2573
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:318
Definition point.h:111
static constexpr unsigned char default_combined_face_orientation()
static constexpr unsigned char reversed_combined_line_orientation()
unsigned int size() const
Definition collection.h:264
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 2 > first
Definition grid_out.cc:4623
#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.
#define DEAL_II_NOT_IMPLEMENTED()
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
if(marked_vertices.size() !=0) for(auto it
for(unsigned int j=best_vertex+1;j< vertices.size();++j) if(vertices_to_use[j])
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition operators.h:49
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
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:220
const Iterator const std_cxx20::type_identity_t< Iterator > & end
Definition parallel.h:610
const Iterator & begin
Definition parallel.h:609
STL namespace.
static unsigned int n_children(const RefinementCase< dim > &refinement_case)