Reference documentation for deal.II version GIT 0582167846 2023-09-24 21:20: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\}}\)
fe_poly_tensor.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_fe_poly_tensor_h
17 #define dealii_fe_poly_tensor_h
18 
19 
20 #include <deal.II/base/config.h>
21 
23 #include <deal.II/base/mutex.h>
26 
27 #include <deal.II/fe/fe.h>
28 
30 
31 #include <memory>
32 
34 
140 template <int dim, int spacedim = dim>
141 class FE_PolyTensor : public FiniteElement<dim, spacedim>
142 {
143 public:
151  const FiniteElementData<dim> &fe_data,
152  const std::vector<bool> &restriction_is_additive_flags,
153  const std::vector<ComponentMask> &nonzero_components);
154 
155 
160 
161  // for documentation, see the FiniteElement base class
162  virtual UpdateFlags
163  requires_update_flags(const UpdateFlags update_flags) const override;
164 
171  virtual double
172  shape_value(const unsigned int i, const Point<dim> &p) const override;
173 
174  // documentation inherited from the base class
175  virtual double
176  shape_value_component(const unsigned int i,
177  const Point<dim> &p,
178  const unsigned int component) const override;
179 
186  virtual Tensor<1, dim>
187  shape_grad(const unsigned int i, const Point<dim> &p) const override;
188 
189  // documentation inherited from the base class
190  virtual Tensor<1, dim>
191  shape_grad_component(const unsigned int i,
192  const Point<dim> &p,
193  const unsigned int component) const override;
194 
201  virtual Tensor<2, dim>
202  shape_grad_grad(const unsigned int i, const Point<dim> &p) const override;
203 
204  // documentation inherited from the base class
205  virtual Tensor<2, dim>
206  shape_grad_grad_component(const unsigned int i,
207  const Point<dim> &p,
208  const unsigned int component) const override;
209 
210 protected:
218  std::vector<MappingKind> mapping_kind;
219 
224  bool
226 
240  bool
242  const unsigned int face_no,
243  const bool face_orientation,
244  const bool face_flip,
245  const bool face_rotation) const;
246 
266 
271  get_mapping_kind(const unsigned int i) const;
272 
273  /* NOTE: The following function has its definition inlined into the class
274  declaration because we otherwise run into a compiler error with MS Visual
275  Studio. */
276  virtual std::unique_ptr<
279  const UpdateFlags update_flags,
280  const Mapping<dim, spacedim> &mapping,
281  const Quadrature<dim> &quadrature,
283  spacedim>
284  &output_data) const override
285  {
286  (void)mapping;
287  (void)output_data;
288  // generate a new data object and
289  // initialize some fields
290  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
291  data_ptr = std::make_unique<InternalData>();
292  auto &data = dynamic_cast<InternalData &>(*data_ptr);
293  data.update_each = requires_update_flags(update_flags);
294 
295  const unsigned int n_q_points = quadrature.size();
296 
297  // some scratch arrays
298  std::vector<Tensor<1, dim>> values(0);
299  std::vector<Tensor<2, dim>> grads(0);
300  std::vector<Tensor<3, dim>> grad_grads(0);
301  std::vector<Tensor<4, dim>> third_derivatives(0);
302  std::vector<Tensor<5, dim>> fourth_derivatives(0);
303 
304  if (update_flags & (update_values | update_gradients | update_hessians))
305  data.dof_sign_change.resize(this->dofs_per_cell);
306 
307  // initialize fields only if really
308  // necessary. otherwise, don't
309  // allocate memory
310 
311  const bool update_transformed_shape_values =
312  std::any_of(this->mapping_kind.begin(),
313  this->mapping_kind.end(),
314  [](const MappingKind t) { return t != mapping_none; });
315 
316  const bool update_transformed_shape_grads =
317  std::any_of(this->mapping_kind.begin(),
318  this->mapping_kind.end(),
319  [](const MappingKind t) {
320  return (t == mapping_raviart_thomas || t == mapping_piola ||
321  t == mapping_nedelec || t == mapping_contravariant);
322  });
323 
324  const bool update_transformed_shape_hessian_tensors =
325  update_transformed_shape_values;
326 
327  if (update_flags & update_values)
328  {
329  values.resize(this->n_dofs_per_cell());
330  data.shape_values.reinit(this->n_dofs_per_cell(), n_q_points);
331  if (update_transformed_shape_values)
332  data.transformed_shape_values.resize(n_q_points);
333  }
334 
335  if (update_flags & update_gradients)
336  {
337  grads.resize(this->n_dofs_per_cell());
338  data.shape_grads.reinit(this->n_dofs_per_cell(), n_q_points);
339  data.transformed_shape_grads.resize(n_q_points);
340 
341  if (update_transformed_shape_grads)
342  data.untransformed_shape_grads.resize(n_q_points);
343  }
344 
345  if (update_flags & update_hessians)
346  {
347  grad_grads.resize(this->n_dofs_per_cell());
348  data.shape_grad_grads.reinit(this->n_dofs_per_cell(), n_q_points);
349  data.transformed_shape_hessians.resize(n_q_points);
350  if (update_transformed_shape_hessian_tensors)
351  data.untransformed_shape_hessian_tensors.resize(n_q_points);
352  }
353 
354  // Compute shape function values
355  // and derivatives and hessians on
356  // the reference cell.
357  // Make sure, that for the
358  // node values N_i holds
359  // N_i(v_j)=\delta_ij for all basis
360  // functions v_j
361  if (update_flags & (update_values | update_gradients))
362  for (unsigned int k = 0; k < n_q_points; ++k)
363  {
364  poly_space->evaluate(quadrature.point(k),
365  values,
366  grads,
367  grad_grads,
368  third_derivatives,
369  fourth_derivatives);
370 
371  if (update_flags & update_values)
372  {
373  if (inverse_node_matrix.n_cols() == 0)
374  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
375  data.shape_values[i][k] = values[i];
376  else
377  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
378  {
379  Tensor<1, dim> add_values;
380  for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
381  add_values += inverse_node_matrix(j, i) * values[j];
382  data.shape_values[i][k] = add_values;
383  }
384  }
385 
386  if (update_flags & update_gradients)
387  {
388  if (inverse_node_matrix.n_cols() == 0)
389  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
390  data.shape_grads[i][k] = grads[i];
391  else
392  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
393  {
394  Tensor<2, dim> add_grads;
395  for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
396  add_grads += inverse_node_matrix(j, i) * grads[j];
397  data.shape_grads[i][k] = add_grads;
398  }
399  }
400 
401  if (update_flags & update_hessians)
402  {
403  if (inverse_node_matrix.n_cols() == 0)
404  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
405  data.shape_grad_grads[i][k] = grad_grads[i];
406  else
407  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
408  {
409  Tensor<3, dim> add_grad_grads;
410  for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
411  add_grad_grads +=
412  inverse_node_matrix(j, i) * grad_grads[j];
413  data.shape_grad_grads[i][k] = add_grad_grads;
414  }
415  }
416  }
417  return data_ptr;
418  }
419 
420  virtual void
422  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
423  const CellSimilarity::Similarity cell_similarity,
424  const Quadrature<dim> &quadrature,
425  const Mapping<dim, spacedim> &mapping,
426  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
428  &mapping_data,
429  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
431  spacedim>
432  &output_data) const override;
433 
435 
436  virtual void
438  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
439  const unsigned int face_no,
440  const hp::QCollection<dim - 1> &quadrature,
441  const Mapping<dim, spacedim> &mapping,
442  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
444  &mapping_data,
445  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
447  spacedim>
448  &output_data) const override;
449 
450  virtual void
452  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
453  const unsigned int face_no,
454  const unsigned int sub_no,
455  const Quadrature<dim - 1> &quadrature,
456  const Mapping<dim, spacedim> &mapping,
457  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
459  &mapping_data,
460  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
462  spacedim>
463  &output_data) const override;
464 
474  class InternalData : public FiniteElement<dim, spacedim>::InternalDataBase
475  {
476  public:
482 
489 
496 
500  mutable std::vector<double> dof_sign_change;
501  mutable std::vector<Tensor<1, spacedim>> transformed_shape_values;
502  // for shape_gradient computations
503  mutable std::vector<Tensor<2, spacedim>> transformed_shape_grads;
504  mutable std::vector<Tensor<2, dim>> untransformed_shape_grads;
505  // for shape_hessian computations
506  mutable std::vector<Tensor<3, spacedim>> transformed_shape_hessians;
507  mutable std::vector<Tensor<3, dim>> untransformed_shape_hessian_tensors;
508  };
509 
510 
511 
516  const std::unique_ptr<const TensorPolynomialsBase<dim>> poly_space;
517 
532 
537 
544 
548  mutable std::vector<Tensor<1, dim>> cached_values;
549 
553  mutable std::vector<Tensor<2, dim>> cached_grads;
554 
559  mutable std::vector<Tensor<3, dim>> cached_grad_grads;
560 };
561 
563 
564 #endif
std::vector< Tensor< 3, spacedim > > transformed_shape_hessians
Table< 2, DerivativeForm< 2, dim, spacedim > > shape_grad_grads
Table< 2, DerivativeForm< 1, dim, spacedim > > shape_grads
std::vector< Tensor< 2, spacedim > > transformed_shape_grads
std::vector< Tensor< 3, dim > > untransformed_shape_hessian_tensors
std::vector< Tensor< 2, dim > > untransformed_shape_grads
Table< 2, Tensor< 1, dim > > shape_values
std::vector< Tensor< 1, spacedim > > transformed_shape_values
std::vector< double > dof_sign_change
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
virtual Tensor< 2, dim > shape_grad_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
std::vector< Table< 2, bool > > adjust_quad_dof_sign_for_face_orientation_table
std::vector< Tensor< 2, dim > > cached_grads
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
FullMatrix< double > inverse_node_matrix
FE_PolyTensor(const FE_PolyTensor &fe)
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
Point< dim > cached_point
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 Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const override
Threads::Mutex cache_mutex
virtual double shape_value(const unsigned int i, const Point< dim > &p) 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
bool adjust_quad_dof_sign_for_face_orientation(const unsigned int index, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation) const
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
std::vector< MappingKind > mapping_kind
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
MappingKind get_mapping_kind(const unsigned int i) const
FE_PolyTensor(const TensorPolynomialsBase< dim > &polynomials, const FiniteElementData< dim > &fe_data, const std::vector< bool > &restriction_is_additive_flags, const std::vector< ComponentMask > &nonzero_components)
std::vector< Tensor< 1, dim > > cached_values
bool single_mapping_kind() const
const std::unique_ptr< const TensorPolynomialsBase< dim > > poly_space
std::vector< Tensor< 3, dim > > cached_grad_grads
unsigned int n_dofs_per_cell() const
const unsigned int dofs_per_cell
Definition: fe_data.h:437
const std::vector< bool > restriction_is_additive_flags
Definition: fe.h:2570
friend class InternalDataBase
Definition: fe.h:3061
const std::vector< ComponentMask > nonzero_components
Definition: fe.h:2579
Definition: point.h:112
const Point< dim > & point(const unsigned int i) const
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
UpdateFlags
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_gradients
Shape function gradients.
MappingKind
Definition: mapping.h:78