Reference documentation for deal.II version GIT relicensing-439-g5fda5c893d 2024-04-20 06:50: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_dgp_nonparametric.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2002 - 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
15
17
19
20#include <deal.II/fe/fe.h>
24#include <deal.II/fe/mapping.h>
25
26#include <deal.II/grid/tria.h>
28
29#include <memory>
30#include <sstream>
31
32
34
35template <int dim, int spacedim>
37 const unsigned int degree)
38 : FiniteElement<dim, spacedim>(
39 FiniteElementData<dim>(get_dpo_vector(degree),
40 1,
41 degree,
42 FiniteElementData<dim>::L2),
43 std::vector<bool>(
44 FiniteElementData<dim>(get_dpo_vector(degree), 1, degree)
45 .n_dofs_per_cell(),
46 true),
47 std::vector<ComponentMask>(
48 FiniteElementData<dim>(get_dpo_vector(degree), 1, degree)
49 .n_dofs_per_cell(),
50 ComponentMask(std::vector<bool>(1, true))))
51 , polynomial_space(Polynomials::Legendre::generate_complete_basis(degree))
52{
53 const unsigned int n_dofs = this->n_dofs_per_cell();
54 for (const unsigned int ref_case :
57 {
58 if (dim != 2 && ref_case != RefinementCase<dim>::isotropic_refinement)
59 // do nothing, as anisotropic
60 // refinement is not
61 // implemented so far
62 continue;
63
64 const unsigned int nc =
66 for (unsigned int i = 0; i < nc; ++i)
67 {
68 this->prolongation[ref_case - 1][i].reinit(n_dofs, n_dofs);
69 // Fill prolongation matrices with
70 // embedding operators
71 for (unsigned int j = 0; j < n_dofs; ++j)
72 this->prolongation[ref_case - 1][i](j, j) = 1.;
73 }
74 }
75
76 // restriction can be defined
77 // through projection for
78 // discontinuous elements, but is
79 // presently not implemented for DGPNonparametric
80 // elements.
81 //
82 // if it were, then the following
83 // snippet would be the right code
84 // if ((degree < Matrices::n_projection_matrices) &&
85 // (Matrices::projection_matrices[degree] != 0))
86 // {
87 // restriction[0].fill (Matrices::projection_matrices[degree]);
88 // }
89 // else
90 // // matrix undefined, set size to zero
91 // for (unsigned int i=0;i<GeometryInfo<dim>::max_children_per_cell;++i)
92 // restriction[i].reinit(0, 0);
93 // since not implemented, set to
94 // "empty". however, that is done in the
95 // default constructor already, so do nothing
96 // for (unsigned int i=0;i<GeometryInfo<dim>::max_children_per_cell;++i)
97 // this->restriction[i].reinit(0, 0);
98
99 // note further, that these
100 // elements have neither support
101 // nor face-support points, so
102 // leave these fields empty
103}
104
105
106
107template <int dim, int spacedim>
108std::string
110{
111 // note that the
112 // FETools::get_fe_by_name
113 // function depends on the
114 // particular format of the string
115 // this function returns, so they
116 // have to be kept in synch
117
118 std::ostringstream namebuf;
119 namebuf << "FE_DGPNonparametric<" << Utilities::dim_string(dim, spacedim)
120 << ">(" << this->degree << ")";
121
122 return namebuf.str();
123}
124
125
126
127template <int dim, int spacedim>
128std::unique_ptr<FiniteElement<dim, spacedim>>
130{
131 return std::make_unique<FE_DGPNonparametric<dim, spacedim>>(*this);
132}
133
134
135
136template <int dim, int spacedim>
137double
139 const Point<dim> &p) const
140{
141 (void)i;
142 (void)p;
143 AssertIndexRange(i, this->n_dofs_per_cell());
144 AssertThrow(false,
146 return 0;
147}
148
149
150
151template <int dim, int spacedim>
152double
154 const unsigned int i,
155 const Point<dim> &p,
156 const unsigned int component) const
157{
158 (void)i;
159 (void)p;
160 (void)component;
161 AssertIndexRange(i, this->n_dofs_per_cell());
162 AssertIndexRange(component, 1);
163 AssertThrow(false,
165 return 0;
166}
167
168
169
170template <int dim, int spacedim>
173 const Point<dim> &p) const
174{
175 (void)i;
176 (void)p;
177 AssertIndexRange(i, this->n_dofs_per_cell());
178 AssertThrow(false,
180 return Tensor<1, dim>();
181}
182
183
184template <int dim, int spacedim>
187 const unsigned int i,
188 const Point<dim> &p,
189 const unsigned int component) const
190{
191 (void)i;
192 (void)p;
193 (void)component;
194 AssertIndexRange(i, this->n_dofs_per_cell());
195 AssertIndexRange(component, 1);
196 AssertThrow(false,
198 return Tensor<1, dim>();
199}
200
201
202
203template <int dim, int spacedim>
206 const Point<dim> &p) const
207{
208 (void)i;
209 (void)p;
210 AssertIndexRange(i, this->n_dofs_per_cell());
211 AssertThrow(false,
213 return Tensor<2, dim>();
214}
215
216
217
218template <int dim, int spacedim>
221 const unsigned int i,
222 const Point<dim> &p,
223 const unsigned int component) const
224{
225 (void)i;
226 (void)p;
227 (void)component;
228 AssertIndexRange(i, this->n_dofs_per_cell());
229 AssertIndexRange(component, 1);
230 AssertThrow(false,
232 return Tensor<2, dim>();
233}
234
235
236//---------------------------------------------------------------------------
237// Auxiliary functions
238//---------------------------------------------------------------------------
239
240
241template <int dim, int spacedim>
242std::vector<unsigned int>
244{
245 std::vector<unsigned int> dpo(dim + 1, static_cast<unsigned int>(0));
246 dpo[dim] = deg + 1;
247 for (unsigned int i = 1; i < dim; ++i)
248 {
249 dpo[dim] *= deg + 1 + i;
250 dpo[dim] /= i + 1;
251 }
252 return dpo;
253}
254
255
256
257template <int dim, int spacedim>
260 const UpdateFlags flags) const
261{
262 UpdateFlags out = flags;
263
266
267 return out;
268}
269
270
271
272//---------------------------------------------------------------------------
273// Data field initialization
274//---------------------------------------------------------------------------
275
276template <int dim, int spacedim>
277std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
279 const UpdateFlags update_flags,
281 const Quadrature<dim> &,
283 spacedim>
284 & /*output_data*/) const
285{
286 // generate a new data object
287 auto data_ptr =
288 std::make_unique<typename FiniteElement<dim, spacedim>::InternalDataBase>();
289 data_ptr->update_each = requires_update_flags(update_flags);
290
291 // other than that, there is nothing we can add here as discussed
292 // in the general documentation of this class
293
294 return data_ptr;
295}
296
297
298
299//---------------------------------------------------------------------------
300// Fill data of FEValues
301//---------------------------------------------------------------------------
302
303template <int dim, int spacedim>
304void
308 const Quadrature<dim> &,
312 &mapping_data,
313 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
315 spacedim>
316 &output_data) const
317{
320
321 const unsigned int n_q_points = mapping_data.quadrature_points.size();
322
323 std::vector<double> values(
324 (fe_internal.update_each & update_values) ? this->n_dofs_per_cell() : 0);
325 std::vector<Tensor<1, dim>> grads(
326 (fe_internal.update_each & update_gradients) ? this->n_dofs_per_cell() : 0);
327 std::vector<Tensor<2, dim>> grad_grads(
328 (fe_internal.update_each & update_hessians) ? this->n_dofs_per_cell() : 0);
329 std::vector<Tensor<3, dim>> empty_vector_of_3rd_order_tensors;
330 std::vector<Tensor<4, dim>> empty_vector_of_4th_order_tensors;
331
332 if (fe_internal.update_each & (update_values | update_gradients))
333 for (unsigned int i = 0; i < n_q_points; ++i)
334 {
335 polynomial_space.evaluate(mapping_data.quadrature_points[i],
336 values,
337 grads,
338 grad_grads,
339 empty_vector_of_3rd_order_tensors,
340 empty_vector_of_4th_order_tensors);
341
342 if (fe_internal.update_each & update_values)
343 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
344 output_data.shape_values[k][i] = values[k];
345
346 if (fe_internal.update_each & update_gradients)
347 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
348 output_data.shape_gradients[k][i] = grads[k];
349
350 if (fe_internal.update_each & update_hessians)
351 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
352 output_data.shape_hessians[k][i] = grad_grads[k];
353 }
354}
355
356
357
358template <int dim, int spacedim>
359void
362 const unsigned int,
367 &mapping_data,
368 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
370 spacedim>
371 &output_data) const
372{
375
376 const unsigned int n_q_points = mapping_data.quadrature_points.size();
377
378 std::vector<double> values(
379 (fe_internal.update_each & update_values) ? this->n_dofs_per_cell() : 0);
380 std::vector<Tensor<1, dim>> grads(
381 (fe_internal.update_each & update_gradients) ? this->n_dofs_per_cell() : 0);
382 std::vector<Tensor<2, dim>> grad_grads(
383 (fe_internal.update_each & update_hessians) ? this->n_dofs_per_cell() : 0);
384 std::vector<Tensor<3, dim>> empty_vector_of_3rd_order_tensors;
385 std::vector<Tensor<4, dim>> empty_vector_of_4th_order_tensors;
386
387 if (fe_internal.update_each & (update_values | update_gradients))
388 for (unsigned int i = 0; i < n_q_points; ++i)
389 {
390 polynomial_space.evaluate(mapping_data.quadrature_points[i],
391 values,
392 grads,
393 grad_grads,
394 empty_vector_of_3rd_order_tensors,
395 empty_vector_of_4th_order_tensors);
396
397 if (fe_internal.update_each & update_values)
398 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
399 output_data.shape_values[k][i] = values[k];
400
401 if (fe_internal.update_each & update_gradients)
402 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
403 output_data.shape_gradients[k][i] = grads[k];
404
405 if (fe_internal.update_each & update_hessians)
406 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
407 output_data.shape_hessians[k][i] = grad_grads[k];
408 }
409}
410
411
412
413template <int dim, int spacedim>
414void
417 const unsigned int,
418 const unsigned int,
419 const Quadrature<dim - 1> &,
423 &mapping_data,
424 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
426 spacedim>
427 &output_data) const
428{
431
432 const unsigned int n_q_points = mapping_data.quadrature_points.size();
433
434 std::vector<double> values(
435 (fe_internal.update_each & update_values) ? this->n_dofs_per_cell() : 0);
436 std::vector<Tensor<1, dim>> grads(
437 (fe_internal.update_each & update_gradients) ? this->n_dofs_per_cell() : 0);
438 std::vector<Tensor<2, dim>> grad_grads(
439 (fe_internal.update_each & update_hessians) ? this->n_dofs_per_cell() : 0);
440 std::vector<Tensor<3, dim>> empty_vector_of_3rd_order_tensors;
441 std::vector<Tensor<4, dim>> empty_vector_of_4th_order_tensors;
442
443 if (fe_internal.update_each & (update_values | update_gradients))
444 for (unsigned int i = 0; i < n_q_points; ++i)
445 {
446 polynomial_space.evaluate(mapping_data.quadrature_points[i],
447 values,
448 grads,
449 grad_grads,
450 empty_vector_of_3rd_order_tensors,
451 empty_vector_of_4th_order_tensors);
452
453 if (fe_internal.update_each & update_values)
454 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
455 output_data.shape_values[k][i] = values[k];
456
457 if (fe_internal.update_each & update_gradients)
458 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
459 output_data.shape_gradients[k][i] = grads[k];
460
461 if (fe_internal.update_each & update_hessians)
462 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
463 output_data.shape_hessians[k][i] = grad_grads[k];
464 }
465}
466
467
468
469template <int dim, int spacedim>
470void
472 const FiniteElement<dim, spacedim> &x_source_fe,
473 FullMatrix<double> &interpolation_matrix,
474 const unsigned int) const
475{
476 // this is only implemented, if the source
477 // FE is also a DGPNonparametric element. in that case,
478 // both elements have no dofs on their
479 // faces and the face interpolation matrix
480 // is necessarily empty -- i.e. there isn't
481 // much we need to do here.
482 (void)interpolation_matrix;
484 (x_source_fe.get_name().find("FE_DGPNonparametric<") == 0) ||
485 (dynamic_cast<const FE_DGPNonparametric<dim, spacedim> *>(&x_source_fe) !=
486 nullptr),
488
489 Assert(interpolation_matrix.m() == 0,
490 ExcDimensionMismatch(interpolation_matrix.m(), 0));
491 Assert(interpolation_matrix.n() == 0,
492 ExcDimensionMismatch(interpolation_matrix.n(), 0));
493}
494
495
496
497template <int dim, int spacedim>
498void
500 const FiniteElement<dim, spacedim> &x_source_fe,
501 const unsigned int,
502 FullMatrix<double> &interpolation_matrix,
503 const unsigned int) const
504{
505 // this is only implemented, if the source
506 // FE is also a DGPNonparametric element. in that case,
507 // both elements have no dofs on their
508 // faces and the face interpolation matrix
509 // is necessarily empty -- i.e. there isn't
510 // much we need to do here.
511 (void)interpolation_matrix;
513 (x_source_fe.get_name().find("FE_DGPNonparametric<") == 0) ||
514 (dynamic_cast<const FE_DGPNonparametric<dim, spacedim> *>(&x_source_fe) !=
515 nullptr),
517
518 Assert(interpolation_matrix.m() == 0,
519 ExcDimensionMismatch(interpolation_matrix.m(), 0));
520 Assert(interpolation_matrix.n() == 0,
521 ExcDimensionMismatch(interpolation_matrix.n(), 0));
522}
523
524
525
526template <int dim, int spacedim>
527bool
532
533
534
535template <int dim, int spacedim>
536std::vector<std::pair<unsigned int, unsigned int>>
538 const FiniteElement<dim, spacedim> &fe_other) const
539{
540 // there are no such constraints for DGPNonparametric
541 // elements at all
542 if (dynamic_cast<const FE_DGPNonparametric<dim, spacedim> *>(&fe_other) !=
543 nullptr)
544 return std::vector<std::pair<unsigned int, unsigned int>>();
545 else
546 {
548 return std::vector<std::pair<unsigned int, unsigned int>>();
549 }
550}
551
552
553
554template <int dim, int spacedim>
555std::vector<std::pair<unsigned int, unsigned int>>
557 const FiniteElement<dim, spacedim> &fe_other) const
558{
559 // there are no such constraints for DGPNonparametric
560 // elements at all
561 if (dynamic_cast<const FE_DGPNonparametric<dim, spacedim> *>(&fe_other) !=
562 nullptr)
563 return std::vector<std::pair<unsigned int, unsigned int>>();
564 else
565 {
567 return std::vector<std::pair<unsigned int, unsigned int>>();
568 }
569}
570
571
572
573template <int dim, int spacedim>
574std::vector<std::pair<unsigned int, unsigned int>>
576 const FiniteElement<dim, spacedim> &fe_other,
577 const unsigned int) const
578{
579 // there are no such constraints for DGPNonparametric
580 // elements at all
581 if (dynamic_cast<const FE_DGPNonparametric<dim, spacedim> *>(&fe_other) !=
582 nullptr)
583 return std::vector<std::pair<unsigned int, unsigned int>>();
584 else
585 {
587 return std::vector<std::pair<unsigned int, unsigned int>>();
588 }
589}
590
591
592
593template <int dim, int spacedim>
596 const FiniteElement<dim, spacedim> &fe_other,
597 const unsigned int codim) const
598{
599 Assert(codim <= dim, ExcImpossibleInDim(dim));
600
601 // vertex/line/face domination
602 // ---------------------------
603 if (codim > 0)
604 // this is a discontinuous element, so by definition there will
605 // be no constraints wherever this element comes together with
606 // any other kind of element
608
609 // cell domination
610 // ---------------
611 if (const FE_DGPNonparametric<dim, spacedim> *fe_nonparametric_other =
612 dynamic_cast<const FE_DGPNonparametric<dim, spacedim> *>(&fe_other))
613 {
614 if (this->degree < fe_nonparametric_other->degree)
616 else if (this->degree == fe_nonparametric_other->degree)
618 else
620 }
621 else if (const FE_Nothing<dim> *fe_nothing =
622 dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
623 {
624 if (fe_nothing->is_dominating())
626 else
627 // the FE_Nothing has no degrees of freedom and it is typically used
628 // in a context where we don't require any continuity along the
629 // interface
631 }
632
635}
636
637
638
639template <int dim, int spacedim>
640bool
642 const unsigned int,
643 const unsigned int) const
644{
645 return true;
646}
647
648
649
650template <int dim, int spacedim>
651std::size_t
657
658
659
660template <int dim, int spacedim>
661unsigned int
663{
664 return this->degree;
665}
666
667
668
669// explicit instantiations
670#include "fe_dgp_nonparametric.inst"
671
672
virtual std::unique_ptr< typename FiniteElement< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const Quadrature< dim > &quadrature, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const override
unsigned int get_degree() const
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
virtual void fill_fe_subface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature< 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 typename FiniteElement< dim, spacedim >::InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const override
virtual Tensor< 2, dim > shape_grad_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
virtual Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const override
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 typename FiniteElement< dim, spacedim >::InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const override
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const override
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
virtual double shape_value(const unsigned int i, const Point< dim > &p) const override
virtual bool hp_constraints_are_implemented() const override
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 override
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
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 override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
virtual std::string get_name() const override
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
virtual void fill_fe_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const Quadrature< dim > &quadrature, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &mapping_data, const typename FiniteElement< dim, spacedim >::InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const override
virtual std::size_t memory_consumption() const override
unsigned int n_dofs_per_cell() const
virtual std::string get_name() const =0
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition fe.h:2427
size_type n() const
size_type m() const
Abstract base class for mapping classes.
Definition mapping.h:318
Definition point.h:111
static std::array< RefinementCase< dim >, n_refinement_cases > all_refinement_cases()
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define AssertThrow(cond, exc)
UpdateFlags
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
#define DEAL_II_NOT_IMPLEMENTED()
std::string dim_string(const int dim, const int spacedim)
Definition utilities.cc:555
STL namespace.
static unsigned int n_children(const RefinementCase< dim > &refinement_case)