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_dgp_nonparametric.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2002 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
18
20
21#include <deal.II/fe/fe.h>
25#include <deal.II/fe/mapping.h>
26
27#include <deal.II/grid/tria.h>
29
30#include <memory>
31#include <sstream>
32
33
35
36template <int dim, int spacedim>
38 const unsigned int degree)
39 : FiniteElement<dim, spacedim>(
40 FiniteElementData<dim>(get_dpo_vector(degree),
41 1,
42 degree,
43 FiniteElementData<dim>::L2),
44 std::vector<bool>(
45 FiniteElementData<dim>(get_dpo_vector(degree), 1, degree)
46 .n_dofs_per_cell(),
47 true),
48 std::vector<ComponentMask>(
49 FiniteElementData<dim>(get_dpo_vector(degree), 1, degree)
50 .n_dofs_per_cell(),
51 std::vector<bool>(1, true)))
52 , polynomial_space(Polynomials::Legendre::generate_complete_basis(degree))
53{
54 const unsigned int n_dofs = this->n_dofs_per_cell();
55 for (unsigned int ref_case = RefinementCase<dim>::cut_x;
56 ref_case < RefinementCase<dim>::isotropic_refinement + 1;
57 ++ref_case)
58 {
59 if (dim != 2 && ref_case != RefinementCase<dim>::isotropic_refinement)
60 // do nothing, as anisotropic
61 // refinement is not
62 // implemented so far
63 continue;
64
65 const unsigned int nc =
67 for (unsigned int i = 0; i < nc; ++i)
68 {
69 this->prolongation[ref_case - 1][i].reinit(n_dofs, n_dofs);
70 // Fill prolongation matrices with
71 // embedding operators
72 for (unsigned int j = 0; j < n_dofs; ++j)
73 this->prolongation[ref_case - 1][i](j, j) = 1.;
74 }
75 }
76
77 // restriction can be defined
78 // through projection for
79 // discontinuous elements, but is
80 // presently not implemented for DGPNonparametric
81 // elements.
82 //
83 // if it were, then the following
84 // snippet would be the right code
85 // if ((degree < Matrices::n_projection_matrices) &&
86 // (Matrices::projection_matrices[degree] != 0))
87 // {
88 // restriction[0].fill (Matrices::projection_matrices[degree]);
89 // }
90 // else
91 // // matrix undefined, set size to zero
92 // for (unsigned int i=0;i<GeometryInfo<dim>::max_children_per_cell;++i)
93 // restriction[i].reinit(0, 0);
94 // since not implemented, set to
95 // "empty". however, that is done in the
96 // default constructor already, so do nothing
97 // for (unsigned int i=0;i<GeometryInfo<dim>::max_children_per_cell;++i)
98 // this->restriction[i].reinit(0, 0);
99
100 // note further, that these
101 // elements have neither support
102 // nor face-support points, so
103 // leave these fields empty
104}
105
106
107
108template <int dim, int spacedim>
109std::string
111{
112 // note that the
113 // FETools::get_fe_by_name
114 // function depends on the
115 // particular format of the string
116 // this function returns, so they
117 // have to be kept in synch
118
119 std::ostringstream namebuf;
120 namebuf << "FE_DGPNonparametric<" << Utilities::dim_string(dim, spacedim)
121 << ">(" << this->degree << ")";
122
123 return namebuf.str();
124}
125
126
127
128template <int dim, int spacedim>
129std::unique_ptr<FiniteElement<dim, spacedim>>
131{
132 return std::make_unique<FE_DGPNonparametric<dim, spacedim>>(*this);
133}
134
135
136
137template <int dim, int spacedim>
138double
140 const Point<dim> & p) const
141{
142 (void)i;
143 (void)p;
144 AssertIndexRange(i, this->n_dofs_per_cell());
145 AssertThrow(false,
147 return 0;
148}
149
150
151
152template <int dim, int spacedim>
153double
155 const unsigned int i,
156 const Point<dim> & p,
157 const unsigned int component) const
158{
159 (void)i;
160 (void)p;
161 (void)component;
162 AssertIndexRange(i, this->n_dofs_per_cell());
163 AssertIndexRange(component, 1);
164 AssertThrow(false,
166 return 0;
167}
168
169
170
171template <int dim, int spacedim>
174 const Point<dim> & p) const
175{
176 (void)i;
177 (void)p;
178 AssertIndexRange(i, this->n_dofs_per_cell());
179 AssertThrow(false,
181 return Tensor<1, dim>();
182}
183
184
185template <int dim, int spacedim>
188 const unsigned int i,
189 const Point<dim> & p,
190 const unsigned int component) const
191{
192 (void)i;
193 (void)p;
194 (void)component;
195 AssertIndexRange(i, this->n_dofs_per_cell());
196 AssertIndexRange(component, 1);
197 AssertThrow(false,
199 return Tensor<1, dim>();
200}
201
202
203
204template <int dim, int spacedim>
207 const Point<dim> & p) const
208{
209 (void)i;
210 (void)p;
211 AssertIndexRange(i, this->n_dofs_per_cell());
212 AssertThrow(false,
214 return Tensor<2, dim>();
215}
216
217
218
219template <int dim, int spacedim>
222 const unsigned int i,
223 const Point<dim> & p,
224 const unsigned int component) const
225{
226 (void)i;
227 (void)p;
228 (void)component;
229 AssertIndexRange(i, this->n_dofs_per_cell());
230 AssertIndexRange(component, 1);
231 AssertThrow(false,
233 return Tensor<2, dim>();
234}
235
236
237//---------------------------------------------------------------------------
238// Auxiliary functions
239//---------------------------------------------------------------------------
240
241
242template <int dim, int spacedim>
243std::vector<unsigned int>
245{
246 std::vector<unsigned int> dpo(dim + 1, static_cast<unsigned int>(0));
247 dpo[dim] = deg + 1;
248 for (unsigned int i = 1; i < dim; ++i)
249 {
250 dpo[dim] *= deg + 1 + i;
251 dpo[dim] /= i + 1;
252 }
253 return dpo;
254}
255
256
257
258template <int dim, int spacedim>
261 const UpdateFlags flags) const
262{
263 UpdateFlags out = flags;
264
267
268 return out;
269}
270
271
272
273//---------------------------------------------------------------------------
274// Data field initialization
275//---------------------------------------------------------------------------
276
277template <int dim, int spacedim>
278std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
280 const UpdateFlags update_flags,
282 const Quadrature<dim> &,
284 spacedim>
285 & /*output_data*/) const
286{
287 // generate a new data object
288 auto data_ptr =
289 std::make_unique<typename FiniteElement<dim, spacedim>::InternalDataBase>();
290 data_ptr->update_each = requires_update_flags(update_flags);
291
292 // other than that, there is nothing we can add here as discussed
293 // in the general documentation of this class
294
295 return data_ptr;
296}
297
298
299
300//---------------------------------------------------------------------------
301// Fill data of FEValues
302//---------------------------------------------------------------------------
303
304template <int dim, int spacedim>
305void
309 const Quadrature<dim> &,
313 & mapping_data,
314 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
316 spacedim>
317 &output_data) const
318{
321
322 const unsigned int n_q_points = mapping_data.quadrature_points.size();
323
324 std::vector<double> values(
325 (fe_internal.update_each & update_values) ? this->n_dofs_per_cell() : 0);
326 std::vector<Tensor<1, dim>> grads(
327 (fe_internal.update_each & update_gradients) ? this->n_dofs_per_cell() : 0);
328 std::vector<Tensor<2, dim>> grad_grads(
329 (fe_internal.update_each & update_hessians) ? this->n_dofs_per_cell() : 0);
330 std::vector<Tensor<3, dim>> empty_vector_of_3rd_order_tensors;
331 std::vector<Tensor<4, dim>> empty_vector_of_4th_order_tensors;
332
333 if (fe_internal.update_each & (update_values | update_gradients))
334 for (unsigned int i = 0; i < n_q_points; ++i)
335 {
336 polynomial_space.evaluate(mapping_data.quadrature_points[i],
337 values,
338 grads,
339 grad_grads,
340 empty_vector_of_3rd_order_tensors,
341 empty_vector_of_4th_order_tensors);
342
343 if (fe_internal.update_each & update_values)
344 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
345 output_data.shape_values[k][i] = values[k];
346
347 if (fe_internal.update_each & update_gradients)
348 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
349 output_data.shape_gradients[k][i] = grads[k];
350
351 if (fe_internal.update_each & update_hessians)
352 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
353 output_data.shape_hessians[k][i] = grad_grads[k];
354 }
355}
356
357
358
359template <int dim, int spacedim>
360void
363 const unsigned int,
368 & mapping_data,
369 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
371 spacedim>
372 &output_data) const
373{
376
377 const unsigned int n_q_points = mapping_data.quadrature_points.size();
378
379 std::vector<double> values(
380 (fe_internal.update_each & update_values) ? this->n_dofs_per_cell() : 0);
381 std::vector<Tensor<1, dim>> grads(
382 (fe_internal.update_each & update_gradients) ? this->n_dofs_per_cell() : 0);
383 std::vector<Tensor<2, dim>> grad_grads(
384 (fe_internal.update_each & update_hessians) ? this->n_dofs_per_cell() : 0);
385 std::vector<Tensor<3, dim>> empty_vector_of_3rd_order_tensors;
386 std::vector<Tensor<4, dim>> empty_vector_of_4th_order_tensors;
387
388 if (fe_internal.update_each & (update_values | update_gradients))
389 for (unsigned int i = 0; i < n_q_points; ++i)
390 {
391 polynomial_space.evaluate(mapping_data.quadrature_points[i],
392 values,
393 grads,
394 grad_grads,
395 empty_vector_of_3rd_order_tensors,
396 empty_vector_of_4th_order_tensors);
397
398 if (fe_internal.update_each & update_values)
399 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
400 output_data.shape_values[k][i] = values[k];
401
402 if (fe_internal.update_each & update_gradients)
403 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
404 output_data.shape_gradients[k][i] = grads[k];
405
406 if (fe_internal.update_each & update_hessians)
407 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
408 output_data.shape_hessians[k][i] = grad_grads[k];
409 }
410}
411
412
413
414template <int dim, int spacedim>
415void
418 const unsigned int,
419 const unsigned int,
420 const Quadrature<dim - 1> &,
424 & mapping_data,
425 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
427 spacedim>
428 &output_data) const
429{
432
433 const unsigned int n_q_points = mapping_data.quadrature_points.size();
434
435 std::vector<double> values(
436 (fe_internal.update_each & update_values) ? this->n_dofs_per_cell() : 0);
437 std::vector<Tensor<1, dim>> grads(
438 (fe_internal.update_each & update_gradients) ? this->n_dofs_per_cell() : 0);
439 std::vector<Tensor<2, dim>> grad_grads(
440 (fe_internal.update_each & update_hessians) ? this->n_dofs_per_cell() : 0);
441 std::vector<Tensor<3, dim>> empty_vector_of_3rd_order_tensors;
442 std::vector<Tensor<4, dim>> empty_vector_of_4th_order_tensors;
443
444 if (fe_internal.update_each & (update_values | update_gradients))
445 for (unsigned int i = 0; i < n_q_points; ++i)
446 {
447 polynomial_space.evaluate(mapping_data.quadrature_points[i],
448 values,
449 grads,
450 grad_grads,
451 empty_vector_of_3rd_order_tensors,
452 empty_vector_of_4th_order_tensors);
453
454 if (fe_internal.update_each & update_values)
455 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
456 output_data.shape_values[k][i] = values[k];
457
458 if (fe_internal.update_each & update_gradients)
459 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
460 output_data.shape_gradients[k][i] = grads[k];
461
462 if (fe_internal.update_each & update_hessians)
463 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
464 output_data.shape_hessians[k][i] = grad_grads[k];
465 }
466}
467
468
469
470template <int dim, int spacedim>
471void
473 const FiniteElement<dim, spacedim> &x_source_fe,
474 FullMatrix<double> & interpolation_matrix,
475 const unsigned int) const
476{
477 // this is only implemented, if the source
478 // FE is also a DGPNonparametric element. in that case,
479 // both elements have no dofs on their
480 // faces and the face interpolation matrix
481 // is necessarily empty -- i.e. there isn't
482 // much we need to do here.
483 (void)interpolation_matrix;
485 (x_source_fe.get_name().find("FE_DGPNonparametric<") == 0) ||
486 (dynamic_cast<const FE_DGPNonparametric<dim, spacedim> *>(&x_source_fe) !=
487 nullptr),
489
490 Assert(interpolation_matrix.m() == 0,
491 ExcDimensionMismatch(interpolation_matrix.m(), 0));
492 Assert(interpolation_matrix.n() == 0,
493 ExcDimensionMismatch(interpolation_matrix.n(), 0));
494}
495
496
497
498template <int dim, int spacedim>
499void
501 const FiniteElement<dim, spacedim> &x_source_fe,
502 const unsigned int,
503 FullMatrix<double> &interpolation_matrix,
504 const unsigned int) const
505{
506 // this is only implemented, if the source
507 // FE is also a DGPNonparametric element. in that case,
508 // both elements have no dofs on their
509 // faces and the face interpolation matrix
510 // is necessarily empty -- i.e. there isn't
511 // much we need to do here.
512 (void)interpolation_matrix;
514 (x_source_fe.get_name().find("FE_DGPNonparametric<") == 0) ||
515 (dynamic_cast<const FE_DGPNonparametric<dim, spacedim> *>(&x_source_fe) !=
516 nullptr),
518
519 Assert(interpolation_matrix.m() == 0,
520 ExcDimensionMismatch(interpolation_matrix.m(), 0));
521 Assert(interpolation_matrix.n() == 0,
522 ExcDimensionMismatch(interpolation_matrix.n(), 0));
523}
524
525
526
527template <int dim, int spacedim>
528bool
530{
531 return true;
532}
533
534
535
536template <int dim, int spacedim>
537std::vector<std::pair<unsigned int, unsigned int>>
539 const FiniteElement<dim, spacedim> &fe_other) const
540{
541 // there are no such constraints for DGPNonparametric
542 // elements at all
543 if (dynamic_cast<const FE_DGPNonparametric<dim, spacedim> *>(&fe_other) !=
544 nullptr)
545 return std::vector<std::pair<unsigned int, unsigned int>>();
546 else
547 {
548 Assert(false, ExcNotImplemented());
549 return std::vector<std::pair<unsigned int, unsigned int>>();
550 }
551}
552
553
554
555template <int dim, int spacedim>
556std::vector<std::pair<unsigned int, unsigned int>>
558 const FiniteElement<dim, spacedim> &fe_other) const
559{
560 // there are no such constraints for DGPNonparametric
561 // elements at all
562 if (dynamic_cast<const FE_DGPNonparametric<dim, spacedim> *>(&fe_other) !=
563 nullptr)
564 return std::vector<std::pair<unsigned int, unsigned int>>();
565 else
566 {
567 Assert(false, ExcNotImplemented());
568 return std::vector<std::pair<unsigned int, unsigned int>>();
569 }
570}
571
572
573
574template <int dim, int spacedim>
575std::vector<std::pair<unsigned int, unsigned int>>
577 const FiniteElement<dim, spacedim> &fe_other,
578 const unsigned int) const
579{
580 // there are no such constraints for DGPNonparametric
581 // elements at all
582 if (dynamic_cast<const FE_DGPNonparametric<dim, spacedim> *>(&fe_other) !=
583 nullptr)
584 return std::vector<std::pair<unsigned int, unsigned int>>();
585 else
586 {
587 Assert(false, ExcNotImplemented());
588 return std::vector<std::pair<unsigned int, unsigned int>>();
589 }
590}
591
592
593
594template <int dim, int spacedim>
597 const FiniteElement<dim, spacedim> &fe_other,
598 const unsigned int codim) const
599{
600 Assert(codim <= dim, ExcImpossibleInDim(dim));
601
602 // vertex/line/face domination
603 // ---------------------------
604 if (codim > 0)
605 // this is a discontinuous element, so by definition there will
606 // be no constraints wherever this element comes together with
607 // any other kind of element
609
610 // cell domination
611 // ---------------
612 if (const FE_DGPNonparametric<dim, spacedim> *fe_nonparametric_other =
613 dynamic_cast<const FE_DGPNonparametric<dim, spacedim> *>(&fe_other))
614 {
615 if (this->degree < fe_nonparametric_other->degree)
617 else if (this->degree == fe_nonparametric_other->degree)
619 else
621 }
622 else if (const FE_Nothing<dim> *fe_nothing =
623 dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
624 {
625 if (fe_nothing->is_dominating())
627 else
628 // the FE_Nothing has no degrees of freedom and it is typically used
629 // in a context where we don't require any continuity along the
630 // interface
632 }
633
634 Assert(false, ExcNotImplemented());
636}
637
638
639
640template <int dim, int spacedim>
641bool
643 const unsigned int,
644 const unsigned int) const
645{
646 return true;
647}
648
649
650
651template <int dim, int spacedim>
652std::size_t
654{
655 Assert(false, ExcNotImplemented());
656 return 0;
657}
658
659
660
661template <int dim, int spacedim>
662unsigned int
664{
665 return this->degree;
666}
667
668
669
670// explicit instantiations
671#include "fe_dgp_nonparametric.inst"
672
673
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:2416
size_type n() const
size_type m() const
Abstract base class for mapping classes.
Definition mapping.h:317
Definition point.h:112
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#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.
std::string dim_string(const int dim, const int spacedim)
Definition utilities.cc:556
STL namespace.
static unsigned int n_children(const RefinementCase< dim > &refinement_case)