Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
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/thread_management.h>
26 
27 #include <deal.II/fe/fe.h>
28 
29 #include <deal.II/lac/full_matrix.h>
30 
31 
32 DEAL_II_NAMESPACE_OPEN
33 
143 template <class PolynomialType, int dim, int spacedim = dim>
144 class FE_PolyTensor : public FiniteElement<dim, spacedim>
145 {
146 public:
153  FE_PolyTensor(const unsigned int degree,
154  const FiniteElementData<dim> & fe_data,
155  const std::vector<bool> & restriction_is_additive_flags,
156  const std::vector<ComponentMask> &nonzero_components);
157 
158  // for documentation, see the FiniteElement base class
159  virtual UpdateFlags
160  requires_update_flags(const UpdateFlags update_flags) const override;
161 
168  virtual double
169  shape_value(const unsigned int i, const Point<dim> &p) const override;
170 
171  // documentation inherited from the base class
172  virtual double
173  shape_value_component(const unsigned int i,
174  const Point<dim> & p,
175  const unsigned int component) const override;
176 
183  virtual Tensor<1, dim>
184  shape_grad(const unsigned int i, const Point<dim> &p) const override;
185 
186  // documentation inherited from the base class
187  virtual Tensor<1, dim>
188  shape_grad_component(const unsigned int i,
189  const Point<dim> & p,
190  const unsigned int component) const override;
191 
198  virtual Tensor<2, dim>
199  shape_grad_grad(const unsigned int i, const Point<dim> &p) const override;
200 
201  // documentation inherited from the base class
202  virtual Tensor<2, dim>
203  shape_grad_grad_component(const unsigned int i,
204  const Point<dim> & p,
205  const unsigned int component) const override;
206 
207 protected:
213 
214 
215  /* NOTE: The following function has its definition inlined into the class
216  declaration because we otherwise run into a compiler error with MS Visual
217  Studio. */
218  virtual std::unique_ptr<
221  const UpdateFlags update_flags,
222  const Mapping<dim, spacedim> & /*mapping*/,
223  const Quadrature<dim> &quadrature,
225  spacedim>
226  & /*output_data*/) const override
227  {
228  // generate a new data object and
229  // initialize some fields
230  std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
231  data_ptr = std_cxx14::make_unique<InternalData>();
232  auto &data = dynamic_cast<InternalData &>(*data_ptr);
233  data.update_each = requires_update_flags(update_flags);
234 
235  const unsigned int n_q_points = quadrature.size();
236 
237  // some scratch arrays
238  std::vector<Tensor<1, dim>> values(0);
239  std::vector<Tensor<2, dim>> grads(0);
240  std::vector<Tensor<3, dim>> grad_grads(0);
241  std::vector<Tensor<4, dim>> third_derivatives(0);
242  std::vector<Tensor<5, dim>> fourth_derivatives(0);
243 
244  if (update_flags & (update_values | update_gradients | update_hessians))
245  data.sign_change.resize(this->dofs_per_cell);
246 
247  // initialize fields only if really
248  // necessary. otherwise, don't
249  // allocate memory
250  if (update_flags & update_values)
251  {
252  values.resize(this->dofs_per_cell);
253  data.shape_values.reinit(this->dofs_per_cell, n_q_points);
254  if (mapping_type != mapping_none)
255  data.transformed_shape_values.resize(n_q_points);
256  }
257 
258  if (update_flags & update_gradients)
259  {
260  grads.resize(this->dofs_per_cell);
261  data.shape_grads.reinit(this->dofs_per_cell, n_q_points);
262  data.transformed_shape_grads.resize(n_q_points);
263 
268  data.untransformed_shape_grads.resize(n_q_points);
269  }
270 
271  if (update_flags & update_hessians)
272  {
273  grad_grads.resize(this->dofs_per_cell);
274  data.shape_grad_grads.reinit(this->dofs_per_cell, n_q_points);
275  data.transformed_shape_hessians.resize(n_q_points);
276  if (mapping_type != mapping_none)
277  data.untransformed_shape_hessian_tensors.resize(n_q_points);
278  }
279 
280  // Compute shape function values
281  // and derivatives and hessians on
282  // the reference cell.
283  // Make sure, that for the
284  // node values N_i holds
285  // N_i(v_j)=\delta_ij for all basis
286  // functions v_j
287  if (update_flags & (update_values | update_gradients))
288  for (unsigned int k = 0; k < n_q_points; ++k)
289  {
290  poly_space.compute(quadrature.point(k),
291  values,
292  grads,
293  grad_grads,
294  third_derivatives,
295  fourth_derivatives);
296 
297  if (update_flags & update_values)
298  {
299  if (inverse_node_matrix.n_cols() == 0)
300  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
301  data.shape_values[i][k] = values[i];
302  else
303  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
304  {
305  Tensor<1, dim> add_values;
306  for (unsigned int j = 0; j < this->dofs_per_cell; ++j)
307  add_values += inverse_node_matrix(j, i) * values[j];
308  data.shape_values[i][k] = add_values;
309  }
310  }
311 
312  if (update_flags & update_gradients)
313  {
314  if (inverse_node_matrix.n_cols() == 0)
315  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
316  data.shape_grads[i][k] = grads[i];
317  else
318  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
319  {
320  Tensor<2, dim> add_grads;
321  for (unsigned int j = 0; j < this->dofs_per_cell; ++j)
322  add_grads += inverse_node_matrix(j, i) * grads[j];
323  data.shape_grads[i][k] = add_grads;
324  }
325  }
326 
327  if (update_flags & update_hessians)
328  {
329  if (inverse_node_matrix.n_cols() == 0)
330  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
331  data.shape_grad_grads[i][k] = grad_grads[i];
332  else
333  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
334  {
335  Tensor<3, dim> add_grad_grads;
336  for (unsigned int j = 0; j < this->dofs_per_cell; ++j)
337  add_grad_grads +=
338  inverse_node_matrix(j, i) * grad_grads[j];
339  data.shape_grad_grads[i][k] = add_grad_grads;
340  }
341  }
342  }
343  return data_ptr;
344  }
345 
346  virtual void
347  fill_fe_values(
348  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
349  const CellSimilarity::Similarity cell_similarity,
350  const Quadrature<dim> & quadrature,
351  const Mapping<dim, spacedim> & mapping,
352  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
353  const ::internal::FEValuesImplementation::MappingRelatedData<dim,
354  spacedim>
355  & mapping_data,
356  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
358  spacedim>
359  &output_data) const override;
360 
361  virtual void
362  fill_fe_face_values(
363  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
364  const unsigned int face_no,
365  const Quadrature<dim - 1> & quadrature,
366  const Mapping<dim, spacedim> & mapping,
367  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
368  const ::internal::FEValuesImplementation::MappingRelatedData<dim,
369  spacedim>
370  & mapping_data,
371  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
373  spacedim>
374  &output_data) const override;
375 
376  virtual void
377  fill_fe_subface_values(
378  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
379  const unsigned int face_no,
380  const unsigned int sub_no,
381  const Quadrature<dim - 1> & quadrature,
382  const Mapping<dim, spacedim> & mapping,
383  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
384  const ::internal::FEValuesImplementation::MappingRelatedData<dim,
385  spacedim>
386  & mapping_data,
387  const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
389  spacedim>
390  &output_data) const override;
391 
401  class InternalData : public FiniteElement<dim, spacedim>::InternalDataBase
402  {
403  public:
409 
416 
423 
427  mutable std::vector<double> sign_change;
428  mutable std::vector<Tensor<1, spacedim>> transformed_shape_values;
429  // for shape_gradient computations
430  mutable std::vector<Tensor<2, spacedim>> transformed_shape_grads;
431  mutable std::vector<Tensor<2, dim>> untransformed_shape_grads;
432  // for shape_hessian computations
433  mutable std::vector<Tensor<3, spacedim>> transformed_shape_hessians;
434  mutable std::vector<Tensor<3, dim>> untransformed_shape_hessian_tensors;
435  };
436 
437 
438 
443  PolynomialType poly_space;
444 
457 
462 
469 
473  mutable std::vector<Tensor<1, dim>> cached_values;
474 
478  mutable std::vector<Tensor<2, dim>> cached_grads;
479 
484  mutable std::vector<Tensor<3, dim>> cached_grad_grads;
485 };
486 
487 DEAL_II_NAMESPACE_CLOSE
488 
489 #endif
Shape function values.
Threads::Mutex cache_mutex
Point< dim > cached_point
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
MappingType
Definition: mapping.h:61
MappingType mapping_type
const unsigned int degree
Definition: fe_base.h:298
FE_PolyTensor(const unsigned int degree, const FiniteElementData< dim > &fe_data, const std::vector< bool > &restriction_is_additive_flags, const std::vector< ComponentMask > &nonzero_components)
const Point< dim > & point(const unsigned int i) const
std::vector< Tensor< 1, dim > > cached_values
std::vector< Tensor< 2, dim > > cached_grads
PolynomialType poly_space
virtual Tensor< 2, dim > shape_grad_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const override
UpdateFlags
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
Abstract base class for mapping classes.
Definition: dof_tools.h:57
virtual double shape_value(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
std::vector< double > sign_change
std::vector< Tensor< 3, dim > > cached_grad_grads
Second derivatives of shape functions.
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
unsigned int size() const
const unsigned int dofs_per_cell
Definition: fe_base.h:282
Table< 2, Tensor< 1, dim > > shape_values
Shape function gradients.
const std::vector< ComponentMask > nonzero_components
Definition: fe.h:2618
Definition: table.h:37
Table< 2, DerivativeForm< 1, dim, spacedim > > shape_grads
FullMatrix< double > inverse_node_matrix
const std::vector< bool > restriction_is_additive_flags
Definition: fe.h:2609
Table< 2, DerivativeForm< 2, dim, spacedim > > shape_grad_grads