Reference documentation for deal.II version Git a39983f358 2020-01-20 20:47:10 -0500
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
fe_poly_tensor.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 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 #ifndef dealii_fe_poly_tensor_h
17 #define dealii_fe_poly_tensor_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/derivative_form.h>
23 #include <deal.II/base/quadrature.h>
24 #include <deal.II/base/std_cxx14/memory.h>
25 #include <deal.II/base/tensor_polynomials_base.h>
26 #include <deal.II/base/thread_management.h>
27 
28 #include <deal.II/fe/fe.h>
29 
30 #include <deal.II/lac/full_matrix.h>
31 
32 
33 DEAL_II_NAMESPACE_OPEN
34 
142 template <int dim, int spacedim = dim>
143 class FE_PolyTensor : public FiniteElement<dim, spacedim>
144 {
145 public:
152  FE_PolyTensor(const TensorPolynomialsBase<dim> &polynomials,
153  const FiniteElementData<dim> & fe_data,
154  const std::vector<bool> & restriction_is_additive_flags,
155  const std::vector<ComponentMask> &nonzero_components);
156 
157 
161  FE_PolyTensor(const FE_PolyTensor &fe);
162 
163  // for documentation, see the FiniteElement base class
164  virtual UpdateFlags
165  requires_update_flags(const UpdateFlags update_flags) const override;
166 
173  virtual double
174  shape_value(const unsigned int i, const Point<dim> &p) const override;
175 
176  // documentation inherited from the base class
177  virtual double
178  shape_value_component(const unsigned int i,
179  const Point<dim> & p,
180  const unsigned int component) const override;
181 
188  virtual Tensor<1, dim>
189  shape_grad(const unsigned int i, const Point<dim> &p) const override;
190 
191  // documentation inherited from the base class
192  virtual Tensor<1, dim>
193  shape_grad_component(const unsigned int i,
194  const Point<dim> & p,
195  const unsigned int component) const override;
196 
203  virtual Tensor<2, dim>
204  shape_grad_grad(const unsigned int i, const Point<dim> &p) const override;
205 
206  // documentation inherited from the base class
207  virtual Tensor<2, dim>
208  shape_grad_grad_component(const unsigned int i,
209  const Point<dim> & p,
210  const unsigned int component) const override;
211 
212 protected:
220  std::vector<MappingKind> mapping_kind;
221 
226  bool
227  single_mapping_kind() const;
228 
233  get_mapping_kind(const unsigned int i) const;
234 
235  /* NOTE: The following function has its definition inlined into the class
236  declaration because we otherwise run into a compiler error with MS Visual
237  Studio. */
238  virtual std::unique_ptr<
241  const UpdateFlags update_flags,
242  const Mapping<dim, spacedim> & /*mapping*/,
243  const Quadrature<dim> &quadrature,
245  spacedim>
246  & /*output_data*/) const override
247  {
248  // generate a new data object and
249  // initialize some fields
250  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
251  data_ptr = std_cxx14::make_unique<InternalData>();
252  auto &data = dynamic_cast<InternalData &>(*data_ptr);
253  data.update_each = requires_update_flags(update_flags);
254 
255  const unsigned int n_q_points = quadrature.size();
256 
257  // some scratch arrays
258  std::vector<Tensor<1, dim>> values(0);
259  std::vector<Tensor<2, dim>> grads(0);
260  std::vector<Tensor<3, dim>> grad_grads(0);
261  std::vector<Tensor<4, dim>> third_derivatives(0);
262  std::vector<Tensor<5, dim>> fourth_derivatives(0);
263 
264  if (update_flags & (update_values | update_gradients | update_hessians))
265  data.sign_change.resize(this->dofs_per_cell);
266 
267  // initialize fields only if really
268  // necessary. otherwise, don't
269  // allocate memory
270 
271  const bool update_transformed_shape_values =
272  std::any_of(this->mapping_kind.begin(),
273  this->mapping_kind.end(),
274  [](const MappingKind t) { return t != mapping_none; });
275 
276  const bool update_transformed_shape_grads =
277  std::any_of(this->mapping_kind.begin(),
278  this->mapping_kind.end(),
279  [](const MappingKind t) {
280  return (t == mapping_raviart_thomas || t == mapping_piola ||
282  });
283 
284  const bool update_transformed_shape_hessian_tensors =
285  update_transformed_shape_values;
286 
287  if (update_flags & update_values)
288  {
289  values.resize(this->dofs_per_cell);
290  data.shape_values.reinit(this->dofs_per_cell, n_q_points);
291  if (update_transformed_shape_values)
292  data.transformed_shape_values.resize(n_q_points);
293  }
294 
295  if (update_flags & update_gradients)
296  {
297  grads.resize(this->dofs_per_cell);
298  data.shape_grads.reinit(this->dofs_per_cell, n_q_points);
299  data.transformed_shape_grads.resize(n_q_points);
300 
301  if (update_transformed_shape_grads)
302  data.untransformed_shape_grads.resize(n_q_points);
303  }
304 
305  if (update_flags & update_hessians)
306  {
307  grad_grads.resize(this->dofs_per_cell);
308  data.shape_grad_grads.reinit(this->dofs_per_cell, n_q_points);
309  data.transformed_shape_hessians.resize(n_q_points);
310  if (update_transformed_shape_hessian_tensors)
311  data.untransformed_shape_hessian_tensors.resize(n_q_points);
312  }
313 
314  // Compute shape function values
315  // and derivatives and hessians on
316  // the reference cell.
317  // Make sure, that for the
318  // node values N_i holds
319  // N_i(v_j)=\delta_ij for all basis
320  // functions v_j
321  if (update_flags & (update_values | update_gradients))
322  for (unsigned int k = 0; k < n_q_points; ++k)
323  {
324  poly_space->evaluate(quadrature.point(k),
325  values,
326  grads,
327  grad_grads,
328  third_derivatives,
329  fourth_derivatives);
330 
331  if (update_flags & update_values)
332  {
333  if (inverse_node_matrix.n_cols() == 0)
334  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
335  data.shape_values[i][k] = values[i];
336  else
337  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
338  {
339  Tensor<1, dim> add_values;
340  for (unsigned int j = 0; j < this->dofs_per_cell; ++j)
341  add_values += inverse_node_matrix(j, i) * values[j];
342  data.shape_values[i][k] = add_values;
343  }
344  }
345 
346  if (update_flags & update_gradients)
347  {
348  if (inverse_node_matrix.n_cols() == 0)
349  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
350  data.shape_grads[i][k] = grads[i];
351  else
352  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
353  {
354  Tensor<2, dim> add_grads;
355  for (unsigned int j = 0; j < this->dofs_per_cell; ++j)
356  add_grads += inverse_node_matrix(j, i) * grads[j];
357  data.shape_grads[i][k] = add_grads;
358  }
359  }
360 
361  if (update_flags & update_hessians)
362  {
363  if (inverse_node_matrix.n_cols() == 0)
364  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
365  data.shape_grad_grads[i][k] = grad_grads[i];
366  else
367  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
368  {
369  Tensor<3, dim> add_grad_grads;
370  for (unsigned int j = 0; j < this->dofs_per_cell; ++j)
371  add_grad_grads +=
372  inverse_node_matrix(j, i) * grad_grads[j];
373  data.shape_grad_grads[i][k] = add_grad_grads;
374  }
375  }
376  }
377  return data_ptr;
378  }
379 
380  virtual void
381  fill_fe_values(
382  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
383  const CellSimilarity::Similarity cell_similarity,
384  const Quadrature<dim> & quadrature,
385  const Mapping<dim, spacedim> & mapping,
386  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
387  const ::internal::FEValuesImplementation::MappingRelatedData<dim,
388  spacedim>
389  & mapping_data,
390  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
392  spacedim>
393  &output_data) const override;
394 
395  virtual void
396  fill_fe_face_values(
397  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
398  const unsigned int face_no,
399  const Quadrature<dim - 1> & quadrature,
400  const Mapping<dim, spacedim> & mapping,
401  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
402  const ::internal::FEValuesImplementation::MappingRelatedData<dim,
403  spacedim>
404  & mapping_data,
405  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
407  spacedim>
408  &output_data) const override;
409 
410  virtual void
411  fill_fe_subface_values(
412  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
413  const unsigned int face_no,
414  const unsigned int sub_no,
415  const Quadrature<dim - 1> & quadrature,
416  const Mapping<dim, spacedim> & mapping,
417  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
418  const ::internal::FEValuesImplementation::MappingRelatedData<dim,
419  spacedim>
420  & mapping_data,
421  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
423  spacedim>
424  &output_data) const override;
425 
435  class InternalData : public FiniteElement<dim, spacedim>::InternalDataBase
436  {
437  public:
443 
450 
457 
461  mutable std::vector<double> sign_change;
462  mutable std::vector<Tensor<1, spacedim>> transformed_shape_values;
463  // for shape_gradient computations
464  mutable std::vector<Tensor<2, spacedim>> transformed_shape_grads;
465  mutable std::vector<Tensor<2, dim>> untransformed_shape_grads;
466  // for shape_hessian computations
467  mutable std::vector<Tensor<3, spacedim>> transformed_shape_hessians;
468  mutable std::vector<Tensor<3, dim>> untransformed_shape_hessian_tensors;
469  };
470 
471 
472 
477  const std::unique_ptr<const TensorPolynomialsBase<dim>> poly_space;
478 
491 
495  mutable std::mutex cache_mutex;
496 
503 
507  mutable std::vector<Tensor<1, dim>> cached_values;
508 
512  mutable std::vector<Tensor<2, dim>> cached_grads;
513 
518  mutable std::vector<Tensor<3, dim>> cached_grad_grads;
519 };
520 
521 DEAL_II_NAMESPACE_CLOSE
522 
523 #endif
Shape function values.
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
std::vector< Tensor< 2, dim > > cached_grads
MappingKind get_mapping_kind(const unsigned int i) const
bool single_mapping_kind() const
std::vector< Tensor< 1, dim > > cached_values
std::vector< double > sign_change
const Point< dim > & point(const unsigned int i) const
FullMatrix< double > inverse_node_matrix
Table< 2, DerivativeForm< 1, dim, spacedim > > shape_grads
MappingKind
Definition: mapping.h:62
Table< 2, DerivativeForm< 2, dim, spacedim > > shape_grad_grads
virtual Tensor< 2, dim > shape_grad_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
UpdateFlags
virtual std::unique_ptr< typename FiniteElement< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags update_flags, const Mapping< dim, spacedim > &, const Quadrature< dim > &quadrature, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &) const override
Abstract base class for mapping classes.
Definition: mapping.h:302
std::mutex cache_mutex
Point< dim > cached_point
std::vector< Tensor< 3, dim > > cached_grad_grads
virtual Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const override
virtual double shape_value(const unsigned int i, const Point< dim > &p) const override
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const override
Second derivatives of shape functions.
unsigned int size() const
const unsigned int dofs_per_cell
Definition: fe_base.h:282
const std::unique_ptr< const TensorPolynomialsBase< dim > > poly_space
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
Shape function gradients.
Table< 2, Tensor< 1, dim > > shape_values
const std::vector< ComponentMask > nonzero_components
Definition: fe.h:2618
Definition: table.h:699
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)
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
const std::vector< bool > restriction_is_additive_flags
Definition: fe.h:2609
std::vector< MappingKind > mapping_kind