deal.II version GIT relicensing-3683-g2da9d4bdba 2025-07-07 18:10:00+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\}}\)
Loading...
Searching...
No Matches
fe_poly_tensor.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2005 - 2025 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_fe_poly_tensor_h
16#define dealii_fe_poly_tensor_h
17
18
19#include <deal.II/base/config.h>
20
22#include <deal.II/base/mutex.h>
26
27#include <deal.II/fe/fe.h>
28
30
31#include <memory>
32
34
140template <int dim, int spacedim = dim>
141class FE_PolyTensor : public FiniteElement<dim, spacedim>
142{
143public:
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
210protected:
211#ifndef DOXYGEN
212 class InternalData;
213#endif
214
219 void
222 const typename QProjector<dim>::DataSetDescriptor &offset,
223 const unsigned int n_q_points,
224 const Mapping<dim, spacedim> &mapping,
225 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
227 &mapping_data,
228 const typename FE_PolyTensor<dim, spacedim>::InternalData &fe_internal,
230 &output_data) const;
231
239 std::vector<MappingKind> mapping_kind;
240
245 bool
247
258 bool
260 const unsigned int index,
261 const unsigned int face_no,
262 const types::geometric_orientation combined_orientation) const;
263
283
288 get_mapping_kind(const unsigned int i) const;
289
290 /* NOTE: The following function has its definition inlined into the class
291 declaration because we otherwise run into a compiler error with MS Visual
292 Studio. */
293 virtual std::unique_ptr<
296 const UpdateFlags update_flags,
297 const Mapping<dim, spacedim> &mapping,
298 const Quadrature<dim> &quadrature,
300 spacedim>
301 &output_data) const override
302 {
303 (void)mapping;
304 (void)output_data;
305 // generate a new data object and
306 // initialize some fields
307 std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
308 data_ptr = std::make_unique<InternalData>();
309 auto &data = dynamic_cast<InternalData &>(*data_ptr);
310 data.update_each = requires_update_flags(update_flags);
311
312 const unsigned int n_q_points = quadrature.size();
313
314 // some scratch arrays
315 std::vector<Tensor<1, dim>> values(0);
316 std::vector<Tensor<2, dim>> grads(0);
317 std::vector<Tensor<3, dim>> grad_grads(0);
318 std::vector<Tensor<4, dim>> third_derivatives(0);
319 std::vector<Tensor<5, dim>> fourth_derivatives(0);
320
321 if (update_flags & (update_values | update_gradients | update_hessians))
322 data.dof_sign_change.resize(this->dofs_per_cell);
323
324 // initialize fields only if really
325 // necessary. otherwise, don't
326 // allocate memory
327
328 const bool update_transformed_shape_values =
329 std::any_of(this->mapping_kind.begin(),
330 this->mapping_kind.end(),
331 [](const MappingKind t) { return t != mapping_none; });
332
333 const bool update_transformed_shape_grads =
334 std::any_of(this->mapping_kind.begin(),
335 this->mapping_kind.end(),
336 [](const MappingKind t) {
337 return (t == mapping_raviart_thomas || t == mapping_piola ||
338 t == mapping_nedelec || t == mapping_contravariant);
339 });
340
341 const bool update_transformed_shape_hessian_tensors =
342 update_transformed_shape_values;
343
344 if (update_flags & update_values)
345 {
346 values.resize(this->n_dofs_per_cell());
347 data.shape_values.reinit(this->n_dofs_per_cell(), n_q_points);
348 if (update_transformed_shape_values)
349 data.transformed_shape_values.resize(n_q_points);
350 }
351
352 if (update_flags & update_gradients)
353 {
354 grads.resize(this->n_dofs_per_cell());
355 data.shape_grads.reinit(this->n_dofs_per_cell(), n_q_points);
356 data.transformed_shape_grads.resize(n_q_points);
357
358 if (update_transformed_shape_grads)
359 data.untransformed_shape_grads.resize(n_q_points);
360 }
361
362 if (update_flags & update_hessians)
363 {
364 grad_grads.resize(this->n_dofs_per_cell());
365 data.shape_grad_grads.reinit(this->n_dofs_per_cell(), n_q_points);
366 data.transformed_shape_hessians.resize(n_q_points);
367 if (update_transformed_shape_hessian_tensors)
368 data.untransformed_shape_hessian_tensors.resize(n_q_points);
369 }
370
371 // Compute shape function values
372 // and derivatives and hessians on
373 // the reference cell.
374 // Make sure, that for the
375 // node values N_i holds
376 // N_i(v_j)=\delta_ij for all basis
377 // functions v_j
378 if (update_flags & (update_values | update_gradients))
379 for (unsigned int k = 0; k < n_q_points; ++k)
380 {
381 poly_space->evaluate(quadrature.point(k),
382 values,
383 grads,
384 grad_grads,
385 third_derivatives,
386 fourth_derivatives);
387
388 if (update_flags & update_values)
389 {
390 if (inverse_node_matrix.n_cols() == 0)
391 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
392 data.shape_values[i][k] = values[i];
393 else
394 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
395 {
396 Tensor<1, dim> add_values;
397 for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
398 add_values += inverse_node_matrix(j, i) * values[j];
399 data.shape_values[i][k] = add_values;
400 }
401 }
402
403 if (update_flags & update_gradients)
404 {
405 if (inverse_node_matrix.n_cols() == 0)
406 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
407 data.shape_grads[i][k] = grads[i];
408 else
409 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
410 {
411 Tensor<2, dim> add_grads;
412 for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
413 add_grads += inverse_node_matrix(j, i) * grads[j];
414 data.shape_grads[i][k] = add_grads;
415 }
416 }
417
418 if (update_flags & update_hessians)
419 {
420 if (inverse_node_matrix.n_cols() == 0)
421 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
422 data.shape_grad_grads[i][k] = grad_grads[i];
423 else
424 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
425 {
426 Tensor<3, dim> add_grad_grads;
427 for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
428 add_grad_grads +=
429 inverse_node_matrix(j, i) * grad_grads[j];
430 data.shape_grad_grads[i][k] = add_grad_grads;
431 }
432 }
433 }
434 return data_ptr;
435 }
436
437 virtual void
440 const CellSimilarity::Similarity cell_similarity,
441 const Quadrature<dim> &quadrature,
442 const Mapping<dim, spacedim> &mapping,
443 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
445 &mapping_data,
446 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
448 spacedim>
449 &output_data) const override;
450
451 using FiniteElement<dim, spacedim>::fill_fe_face_values;
452
453 virtual void
456 const unsigned int face_no,
457 const hp::QCollection<dim - 1> &quadrature,
458 const Mapping<dim, spacedim> &mapping,
459 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
461 &mapping_data,
462 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
464 spacedim>
465 &output_data) const override;
466
467 virtual void
470 const unsigned int face_no,
471 const unsigned int sub_no,
472 const Quadrature<dim - 1> &quadrature,
473 const Mapping<dim, spacedim> &mapping,
474 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
476 &mapping_data,
477 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
479 spacedim>
480 &output_data) const override;
481
491 class InternalData : public FiniteElement<dim, spacedim>::InternalDataBase
492 {
493 public:
499
506
513
517 mutable std::vector<double> dof_sign_change;
518 mutable std::vector<Tensor<1, spacedim>> transformed_shape_values;
519 // for shape_gradient computations
520 mutable std::vector<Tensor<2, spacedim>> transformed_shape_grads;
521 mutable std::vector<Tensor<2, dim>> untransformed_shape_grads;
522 // for shape_hessian computations
523 mutable std::vector<Tensor<3, spacedim>> transformed_shape_hessians;
524 mutable std::vector<Tensor<3, dim>> untransformed_shape_hessian_tensors;
525 };
526
527
528
533 const std::unique_ptr<const TensorPolynomialsBase<dim>> poly_space;
534
549
554
561
565 mutable std::vector<Tensor<1, dim>> cached_values;
566
570 mutable std::vector<Tensor<2, dim>> cached_grads;
571
576 mutable std::vector<Tensor<3, dim>> cached_grad_grads;
577};
578
580
581#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
std::vector< Table< 2, bool > > adjust_quad_dof_sign_for_face_orientation_table
std::vector< Tensor< 2, dim > > cached_grads
bool adjust_quad_dof_sign_for_face_orientation(const unsigned int index, const unsigned int face_no, const types::geometric_orientation combined_orientation) const
FullMatrix< double > inverse_node_matrix
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const override
virtual Tensor< 2, dim > shape_grad_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
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 Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const override
void compute_fill(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const typename QProjector< dim >::DataSetDescriptor &offset, const unsigned int n_q_points, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &mapping_data, const typename FE_PolyTensor< dim, spacedim >::InternalData &fe_internal, internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const
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
Threads::Mutex cache_mutex
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 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
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
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
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
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:436
const std::vector< bool > restriction_is_additive_flags
Definition fe.h:2690
const std::vector< ComponentMask > nonzero_components
Definition fe.h:2699
Abstract base class for mapping classes.
Definition mapping.h:320
Definition point.h:113
Class storing the offset index into a Quadrature rule created by project_to_all_faces() or project_to...
Definition qprojector.h:334
const Point< dim > & point(const unsigned int i) const
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
UpdateFlags
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_gradients
Shape function gradients.
MappingKind
Definition mapping.h:81
std::vector< index_type > data
Definition mpi.cc:750
unsigned char geometric_orientation
Definition types.h:40