Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
fe_dgp_nonparametric.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2002 - 2019 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 
17 #include <deal.II/base/quadrature.h>
18 #include <deal.II/base/std_cxx14/memory.h>
19 
20 #include <deal.II/dofs/dof_accessor.h>
21 
22 #include <deal.II/fe/fe.h>
23 #include <deal.II/fe/fe_dgp_nonparametric.h>
24 #include <deal.II/fe/fe_nothing.h>
25 #include <deal.II/fe/fe_values.h>
26 #include <deal.II/fe/mapping.h>
27 
28 #include <deal.II/grid/tria.h>
29 #include <deal.II/grid/tria_iterator.h>
30 
31 #include <sstream>
32 
33 
34 DEAL_II_NAMESPACE_OPEN
35 
36 template <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).dofs_per_cell,
46  true),
47  std::vector<ComponentMask>(
48  FiniteElementData<dim>(get_dpo_vector(degree), 1, degree).dofs_per_cell,
49  std::vector<bool>(1, true)))
50  , polynomial_space(Polynomials::Legendre::generate_complete_basis(degree))
51 {
52  const unsigned int n_dofs = this->dofs_per_cell;
53  for (unsigned int ref_case = RefinementCase<dim>::cut_x;
54  ref_case < RefinementCase<dim>::isotropic_refinement + 1;
55  ++ref_case)
56  {
57  if (dim != 2 && ref_case != RefinementCase<dim>::isotropic_refinement)
58  // do nothing, as anisotropic
59  // refinement is not
60  // implemented so far
61  continue;
62 
63  const unsigned int nc =
65  for (unsigned int i = 0; i < nc; ++i)
66  {
67  this->prolongation[ref_case - 1][i].reinit(n_dofs, n_dofs);
68  // Fill prolongation matrices with
69  // embedding operators
70  for (unsigned int j = 0; j < n_dofs; ++j)
71  this->prolongation[ref_case - 1][i](j, j) = 1.;
72  }
73  }
74 
75  // restriction can be defined
76  // through projection for
77  // discontinuous elements, but is
78  // presently not implemented for DGPNonparametric
79  // elements.
80  //
81  // if it were, then the following
82  // snippet would be the right code
83  // if ((degree < Matrices::n_projection_matrices) &&
84  // (Matrices::projection_matrices[degree] != 0))
85  // {
86  // restriction[0].fill (Matrices::projection_matrices[degree]);
87  // }
88  // else
89  // // matrix undefined, set size to zero
90  // for (unsigned int i=0;i<GeometryInfo<dim>::max_children_per_cell;++i)
91  // restriction[i].reinit(0, 0);
92  // since not implemented, set to
93  // "empty". however, that is done in the
94  // default constructor already, so do nothing
95  // for (unsigned int i=0;i<GeometryInfo<dim>::max_children_per_cell;++i)
96  // this->restriction[i].reinit(0, 0);
97 
98  // note further, that these
99  // elements have neither support
100  // nor face-support points, so
101  // leave these fields empty
102 }
103 
104 
105 
106 template <int dim, int spacedim>
107 std::string
109 {
110  // note that the
111  // FETools::get_fe_by_name
112  // function depends on the
113  // particular format of the string
114  // this function returns, so they
115  // have to be kept in synch
116 
117  std::ostringstream namebuf;
118  namebuf << "FE_DGPNonparametric<" << Utilities::dim_string(dim, spacedim)
119  << ">(" << this->degree << ")";
120 
121  return namebuf.str();
122 }
123 
124 
125 
126 template <int dim, int spacedim>
127 std::unique_ptr<FiniteElement<dim, spacedim>>
129 {
130  return std_cxx14::make_unique<FE_DGPNonparametric<dim, spacedim>>(*this);
131 }
132 
133 
134 
135 template <int dim, int spacedim>
136 double
138  const Point<dim> & p) const
139 {
140  (void)i;
141  (void)p;
142  Assert(i < this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
143  AssertThrow(false,
145  return 0;
146 }
147 
148 
149 
150 template <int dim, int spacedim>
151 double
153  const unsigned int i,
154  const Point<dim> & p,
155  const unsigned int component) const
156 {
157  (void)i;
158  (void)p;
159  (void)component;
160  Assert(i < this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
161  Assert(component == 0, ExcIndexRange(component, 0, 1));
162  AssertThrow(false,
164  return 0;
165 }
166 
167 
168 
169 template <int dim, int spacedim>
172  const Point<dim> & p) const
173 {
174  (void)i;
175  (void)p;
176  Assert(i < this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
177  AssertThrow(false,
179  return Tensor<1, dim>();
180 }
181 
182 
183 template <int dim, int spacedim>
186  const unsigned int i,
187  const Point<dim> & p,
188  const unsigned int component) const
189 {
190  (void)i;
191  (void)p;
192  (void)component;
193  Assert(i < this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
194  Assert(component == 0, ExcIndexRange(component, 0, 1));
195  AssertThrow(false,
197  return Tensor<1, dim>();
198 }
199 
200 
201 
202 template <int dim, int spacedim>
205  const Point<dim> & p) const
206 {
207  (void)i;
208  (void)p;
209  Assert(i < this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
210  AssertThrow(false,
212  return Tensor<2, dim>();
213 }
214 
215 
216 
217 template <int dim, int spacedim>
220  const unsigned int i,
221  const Point<dim> & p,
222  const unsigned int component) const
223 {
224  (void)i;
225  (void)p;
226  (void)component;
227  Assert(i < this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
228  Assert(component == 0, ExcIndexRange(component, 0, 1));
229  AssertThrow(false,
231  return Tensor<2, dim>();
232 }
233 
234 
235 //---------------------------------------------------------------------------
236 // Auxiliary functions
237 //---------------------------------------------------------------------------
238 
239 
240 template <int dim, int spacedim>
241 std::vector<unsigned int>
243 {
244  std::vector<unsigned int> dpo(dim + 1, static_cast<unsigned int>(0));
245  dpo[dim] = deg + 1;
246  for (unsigned int i = 1; i < dim; ++i)
247  {
248  dpo[dim] *= deg + 1 + i;
249  dpo[dim] /= i + 1;
250  }
251  return dpo;
252 }
253 
254 
255 
256 template <int dim, int spacedim>
259  const UpdateFlags flags) const
260 {
261  UpdateFlags out = flags;
262 
265 
266  return out;
267 }
268 
269 
270 
271 //---------------------------------------------------------------------------
272 // Data field initialization
273 //---------------------------------------------------------------------------
274 
275 template <int dim, int spacedim>
276 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
278  const UpdateFlags update_flags,
279  const Mapping<dim, spacedim> &,
280  const Quadrature<dim> &,
282  spacedim>
283  & /*output_data*/) const
284 {
285  // generate a new data object
286  auto data_ptr = std_cxx14::make_unique<
288  data_ptr->update_each = requires_update_flags(update_flags);
289 
290  // other than that, there is nothing we can add here as discussed
291  // in the general documentation of this class
292 
293  return data_ptr;
294 }
295 
296 
297 
298 //---------------------------------------------------------------------------
299 // Fill data of FEValues
300 //---------------------------------------------------------------------------
301 
302 template <int dim, int spacedim>
303 void
307  const Quadrature<dim> &,
308  const Mapping<dim, spacedim> &,
310  const ::internal::FEValuesImplementation::MappingRelatedData<dim,
311  spacedim>
312  & mapping_data,
313  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
315  spacedim>
316  &output_data) const
317 {
319  ExcInternalError());
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->dofs_per_cell : 0);
325  std::vector<Tensor<1, dim>> grads(
326  (fe_internal.update_each & update_gradients) ? this->dofs_per_cell : 0);
327  std::vector<Tensor<2, dim>> grad_grads(
328  (fe_internal.update_each & update_hessians) ? this->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.compute(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->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->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->dofs_per_cell; ++k)
352  output_data.shape_hessians[k][i] = grad_grads[k];
353  }
354 }
355 
356 
357 
358 template <int dim, int spacedim>
359 void
362  const unsigned int,
363  const Quadrature<dim - 1> &,
364  const Mapping<dim, spacedim> &,
366  const ::internal::FEValuesImplementation::MappingRelatedData<dim,
367  spacedim>
368  & mapping_data,
369  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
371  spacedim>
372  &output_data) const
373 {
375  ExcInternalError());
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->dofs_per_cell : 0);
381  std::vector<Tensor<1, dim>> grads(
382  (fe_internal.update_each & update_gradients) ? this->dofs_per_cell : 0);
383  std::vector<Tensor<2, dim>> grad_grads(
384  (fe_internal.update_each & update_hessians) ? this->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.compute(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->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->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->dofs_per_cell; ++k)
408  output_data.shape_hessians[k][i] = grad_grads[k];
409  }
410 }
411 
412 
413 
414 template <int dim, int spacedim>
415 void
418  const unsigned int,
419  const unsigned int,
420  const Quadrature<dim - 1> &,
421  const Mapping<dim, spacedim> &,
423  const ::internal::FEValuesImplementation::MappingRelatedData<dim,
424  spacedim>
425  & mapping_data,
426  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
428  spacedim>
429  &output_data) const
430 {
432  ExcInternalError());
433 
434  const unsigned int n_q_points = mapping_data.quadrature_points.size();
435 
436  std::vector<double> values(
437  (fe_internal.update_each & update_values) ? this->dofs_per_cell : 0);
438  std::vector<Tensor<1, dim>> grads(
439  (fe_internal.update_each & update_gradients) ? this->dofs_per_cell : 0);
440  std::vector<Tensor<2, dim>> grad_grads(
441  (fe_internal.update_each & update_hessians) ? this->dofs_per_cell : 0);
442  std::vector<Tensor<3, dim>> empty_vector_of_3rd_order_tensors;
443  std::vector<Tensor<4, dim>> empty_vector_of_4th_order_tensors;
444 
445  if (fe_internal.update_each & (update_values | update_gradients))
446  for (unsigned int i = 0; i < n_q_points; ++i)
447  {
448  polynomial_space.compute(mapping_data.quadrature_points[i],
449  values,
450  grads,
451  grad_grads,
452  empty_vector_of_3rd_order_tensors,
453  empty_vector_of_4th_order_tensors);
454 
455  if (fe_internal.update_each & update_values)
456  for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
457  output_data.shape_values[k][i] = values[k];
458 
459  if (fe_internal.update_each & update_gradients)
460  for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
461  output_data.shape_gradients[k][i] = grads[k];
462 
463  if (fe_internal.update_each & update_hessians)
464  for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
465  output_data.shape_hessians[k][i] = grad_grads[k];
466  }
467 }
468 
469 
470 
471 template <int dim, int spacedim>
472 void
474  const FiniteElement<dim, spacedim> &x_source_fe,
475  FullMatrix<double> & interpolation_matrix) 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;
484  AssertThrow(
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 
498 template <int dim, int spacedim>
499 void
501  const FiniteElement<dim, spacedim> &x_source_fe,
502  const unsigned int,
503  FullMatrix<double> &interpolation_matrix) 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;
512  AssertThrow(
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 
526 template <int dim, int spacedim>
527 bool
529 {
530  return true;
531 }
532 
533 
534 
535 template <int dim, int spacedim>
536 std::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  {
547  Assert(false, ExcNotImplemented());
548  return std::vector<std::pair<unsigned int, unsigned int>>();
549  }
550 }
551 
552 
553 
554 template <int dim, int spacedim>
555 std::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  {
566  Assert(false, ExcNotImplemented());
567  return std::vector<std::pair<unsigned int, unsigned int>>();
568  }
569 }
570 
571 
572 
573 template <int dim, int spacedim>
574 std::vector<std::pair<unsigned int, unsigned int>>
576  const FiniteElement<dim, spacedim> &fe_other) const
577 {
578  // there are no such constraints for DGPNonparametric
579  // elements at all
580  if (dynamic_cast<const FE_DGPNonparametric<dim, spacedim> *>(&fe_other) !=
581  nullptr)
582  return std::vector<std::pair<unsigned int, unsigned int>>();
583  else
584  {
585  Assert(false, ExcNotImplemented());
586  return std::vector<std::pair<unsigned int, unsigned int>>();
587  }
588 }
589 
590 
591 
592 template <int dim, int spacedim>
595  const FiniteElement<dim, spacedim> &fe_other,
596  const unsigned int codim) const
597 {
598  Assert(codim <= dim, ExcImpossibleInDim(dim));
599 
600  // vertex/line/face domination
601  // ---------------------------
602  if (codim > 0)
603  // this is a discontinuous element, so by definition there will
604  // be no constraints wherever this element comes together with
605  // any other kind of element
607 
608  // cell domination
609  // ---------------
610  if (const FE_DGPNonparametric<dim, spacedim> *fe_nonparametric_other =
611  dynamic_cast<const FE_DGPNonparametric<dim, spacedim> *>(&fe_other))
612  {
613  if (this->degree < fe_nonparametric_other->degree)
615  else if (this->degree == fe_nonparametric_other->degree)
617  else
619  }
620  else if (const FE_Nothing<dim> *fe_nothing =
621  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
622  {
623  if (fe_nothing->is_dominating())
625  else
626  // the FE_Nothing has no degrees of freedom and it is typically used
627  // in a context where we don't require any continuity along the
628  // interface
630  }
631 
632  Assert(false, ExcNotImplemented());
634 }
635 
636 
637 
638 template <int dim, int spacedim>
639 bool
641  const unsigned int,
642  const unsigned int) const
643 {
644  return true;
645 }
646 
647 
648 
649 template <int dim, int spacedim>
650 std::size_t
652 {
653  Assert(false, ExcNotImplemented());
654  return 0;
655 }
656 
657 
658 
659 template <int dim, int spacedim>
660 unsigned int
662 {
663  return this->degree;
664 }
665 
666 
667 
668 // explicit instantiations
669 #include "fe_dgp_nonparametric.inst"
670 
671 
672 DEAL_II_NAMESPACE_CLOSE
Shape function values.
size_type m() const
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const override
friend class FE_DGPNonparametric
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
virtual std::size_t memory_consumption() const override
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix) const override
virtual double shape_value(const unsigned int i, const Point< dim > &p) const override
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
STL namespace.
Transformed quadrature points.
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
virtual bool hp_constraints_are_implemented() const override
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
size_type n() const
unsigned int get_degree() const
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2456
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
#define Assert(cond, exc)
Definition: exceptions.h:1407
UpdateFlags
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Abstract base class for mapping classes.
Definition: dof_tools.h:57
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
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
virtual std::string get_name() const =0
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
Second derivatives of shape functions.
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:457
const unsigned int dofs_per_cell
Definition: fe_base.h:282
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
virtual Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const override
Shape function gradients.
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
virtual std::string get_name() const override
virtual Tensor< 2, dim > shape_grad_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
static ::ExceptionBase & ExcNotImplemented()
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
UpdateFlags update_each
Definition: fe.h:715
static ::ExceptionBase & ExcInternalError()