Reference documentation for deal.II version 9.5.0
\(\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.cc
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
17#include <deal.II/base/array_view.h>
18#include <deal.II/base/derivative_form.h>
19#include <deal.II/base/polynomials_abf.h>
20#include <deal.II/base/polynomials_bdm.h>
21#include <deal.II/base/polynomials_bernardi_raugel.h>
22#include <deal.II/base/polynomials_nedelec.h>
23#include <deal.II/base/polynomials_raviart_thomas.h>
24#include <deal.II/base/polynomials_rt_bubbles.h>
25#include <deal.II/base/qprojector.h>
26
27#include <deal.II/fe/fe_poly_tensor.h>
28#include <deal.II/fe/fe_values.h>
29#include <deal.II/fe/mapping_cartesian.h>
30
31#include <deal.II/grid/cell_id.h>
32#include <deal.II/grid/tria.h>
33
34DEAL_II_NAMESPACE_OPEN
35
36namespace internal
37{
38 namespace FE_PolyTensor
39 {
40 namespace
41 {
42 template <int spacedim>
43 void
44 get_dof_sign_change_h_div(
45 const typename ::Triangulation<1, spacedim>::cell_iterator &,
46 const FiniteElement<1, spacedim> &,
47 const std::vector<MappingKind> &,
48 std::vector<double> &)
49 {
50 // Nothing to do in 1d.
51 }
52
53
54
55 // TODO: This function is not a consistent fix of the orientation issue
56 // like in 3d. It is rather kept not to break legacy behavior in 2d but
57 // should be replaced. See also the implementation of
58 // FE_RaviartThomas<dim>::initialize_quad_dof_index_permutation_and_sign_change()
59 // or other H(div) conforming elements such as FE_ABF<dim> and
60 // FE_BDM<dim>.
61 template <int spacedim>
62 void
63 get_dof_sign_change_h_div(
64 const typename ::Triangulation<2, spacedim>::cell_iterator &cell,
65 const FiniteElement<2, spacedim> & fe,
66 const std::vector<MappingKind> &mapping_kind,
67 std::vector<double> & face_sign)
68 {
69 const unsigned int dim = 2;
70 // const unsigned int spacedim = 2;
71
72 const CellId this_cell_id = cell->id();
73
74 for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
75 {
76 typename ::Triangulation<dim, spacedim>::face_iterator face =
77 cell->face(f);
78 if (!face->at_boundary())
79 {
80 const unsigned int nn = cell->neighbor_face_no(f);
81 const typename ::Triangulation<dim,
82 spacedim>::cell_iterator
83 neighbor_cell_at_face = cell->neighbor(f);
84 const CellId neigbor_cell_id = neighbor_cell_at_face->id();
85
86 // Only fix sign if the orientation is opposite and only do so
87 // on the face dofs on the cell with smaller cell_id.
88 if (((nn + f) % 2 == 0) && this_cell_id < neigbor_cell_id)
89 for (unsigned int j = 0; j < fe.n_dofs_per_face(f); ++j)
90 {
91 const unsigned int cell_j = fe.face_to_cell_index(j, f);
92
93 Assert(f * fe.n_dofs_per_face(f) + j < face_sign.size(),
94 ExcInternalError());
95 Assert(mapping_kind.size() == 1 ||
96 cell_j < mapping_kind.size(),
97 ExcInternalError());
98
99 // TODO: This is probably only going to work for those
100 // elements for which all dofs are face dofs
101 if ((mapping_kind.size() > 1 ?
102 mapping_kind[cell_j] :
103 mapping_kind[0]) == mapping_raviart_thomas)
104 face_sign[f * fe.n_dofs_per_face(f) + j] = -1.0;
105 }
106 }
107 }
108 }
109
110
111
112 template <int spacedim>
113 void
114 get_dof_sign_change_h_div(
115 const typename ::Triangulation<3, spacedim>::cell_iterator
116 & /*cell*/,
117 const FiniteElement<3, spacedim> & /*fe*/,
118 const std::vector<MappingKind> & /*mapping_kind*/,
119 std::vector<double> & /*face_sign*/)
120 {
121 // Nothing to do. In 3d we take care of it through the
122 // adjust_quad_dof_sign_for_face_orientation_table
123 }
124
125 template <int spacedim>
126 void
127 get_dof_sign_change_nedelec(
128 const typename ::Triangulation<1, spacedim>::cell_iterator
129 & /*cell*/,
130 const FiniteElement<1, spacedim> & /*fe*/,
131 const std::vector<MappingKind> & /*mapping_kind*/,
132 std::vector<double> & /*line_dof_sign*/)
133 {
134 // nothing to do in 1d
135 }
136
137
138
139 template <int spacedim>
140 void
141 get_dof_sign_change_nedelec(
142 const typename ::Triangulation<2, spacedim>::cell_iterator &cell,
143 const FiniteElement<2, spacedim> & /*fe*/,
144 const std::vector<MappingKind> &mapping_kind,
145 std::vector<double> & line_dof_sign)
146 {
147 const unsigned int dim = 2;
148 // TODO: This fixes only lowest order
149 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
150 if (!(cell->line_orientation(l)) &&
151 mapping_kind[0] == mapping_nedelec)
152 line_dof_sign[l] = -1.0;
153 }
154
155
156 template <int spacedim>
157 void
158 get_dof_sign_change_nedelec(
159 const typename ::Triangulation<3, spacedim>::cell_iterator &cell,
160 const FiniteElement<3, spacedim> & /*fe*/,
161 const std::vector<MappingKind> &mapping_kind,
162 std::vector<double> & line_dof_sign)
163 {
164 const unsigned int dim = 3;
165 // TODO: This is probably only going to work for those elements for
166 // which all dofs are face dofs
167 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
168 if (!(cell->line_orientation(l)) &&
169 mapping_kind[0] == mapping_nedelec)
170 line_dof_sign[l] = -1.0;
171 }
172 } // namespace
173 } // namespace FE_PolyTensor
174} // namespace internal
175
176
177template <int dim, int spacedim>
178FE_PolyTensor<dim, spacedim>::FE_PolyTensor(
179 const TensorPolynomialsBase<dim> &polynomials,
180 const FiniteElementData<dim> & fe_data,
181 const std::vector<bool> & restriction_is_additive_flags,
182 const std::vector<ComponentMask> &nonzero_components)
183 : FiniteElement<dim, spacedim>(fe_data,
184 restriction_is_additive_flags,
185 nonzero_components)
186 , mapping_kind({MappingKind::mapping_none})
187 , poly_space(polynomials.clone())
188{
189 cached_point(0) = -1;
190 // Set up the table converting
191 // components to base
192 // components. Since we have only
193 // one base element, everything
194 // remains zero except the
195 // component in the base, which is
196 // the component itself
197 for (unsigned int comp = 0; comp < this->n_components(); ++comp)
198 this->component_to_base_table[comp].first.second = comp;
199
200 if (dim == 3)
201 {
202 adjust_quad_dof_sign_for_face_orientation_table.resize(
203 this->n_unique_2d_subobjects());
204
205 for (unsigned int f = 0; f < this->n_unique_2d_subobjects(); ++f)
206 {
207 adjust_quad_dof_sign_for_face_orientation_table[f] =
208 Table<2, bool>(this->n_dofs_per_quad(f),
209 this->reference_cell().n_face_orientations(f));
210 adjust_quad_dof_sign_for_face_orientation_table[f].fill(false);
211 }
212 }
213}
214
215
216
217template <int dim, int spacedim>
218FE_PolyTensor<dim, spacedim>::FE_PolyTensor(const FE_PolyTensor &fe)
219 : FiniteElement<dim, spacedim>(fe)
220 , mapping_kind(fe.mapping_kind)
221 , adjust_quad_dof_sign_for_face_orientation_table(
222 fe.adjust_quad_dof_sign_for_face_orientation_table)
223 , poly_space(fe.poly_space->clone())
224 , inverse_node_matrix(fe.inverse_node_matrix)
225{}
226
227
228
229template <int dim, int spacedim>
230bool
231FE_PolyTensor<dim, spacedim>::single_mapping_kind() const
232{
233 return mapping_kind.size() == 1;
234}
235
236
237template <int dim, int spacedim>
238bool
239FE_PolyTensor<dim, spacedim>::adjust_quad_dof_sign_for_face_orientation(
240 const unsigned int index,
241 const unsigned int face,
242 const bool face_orientation,
243 const bool face_flip,
244 const bool face_rotation) const
245{
246 // do nothing in 1d and 2d
247 if (dim < 3)
248 return false;
249
250 // The exception are discontinuous
251 // elements for which there should be no
252 // face dofs anyway (i.e. dofs_per_quad==0
253 // in 3d), so we don't need the table, but
254 // the function should also not have been
255 // called
256 AssertIndexRange(index, this->n_dofs_per_quad(face));
257 Assert(adjust_quad_dof_sign_for_face_orientation_table
258 [this->n_unique_2d_subobjects() == 1 ? 0 : face]
259 .n_elements() ==
260 this->reference_cell().n_face_orientations(face) *
261 this->n_dofs_per_quad(face),
262 ExcInternalError());
263
264 return adjust_quad_dof_sign_for_face_orientation_table
265 [this->n_unique_2d_subobjects() == 1 ? 0 : face](
266 index,
267 4 * static_cast<int>(face_orientation) + 2 * static_cast<int>(face_flip) +
268 static_cast<int>(face_rotation));
269}
270
271
272template <int dim, int spacedim>
273MappingKind
274FE_PolyTensor<dim, spacedim>::get_mapping_kind(const unsigned int i) const
275{
276 if (single_mapping_kind())
277 return mapping_kind[0];
278
279 AssertIndexRange(i, mapping_kind.size());
280 return mapping_kind[i];
281}
282
283
284
285template <int dim, int spacedim>
286double
287FE_PolyTensor<dim, spacedim>::shape_value(const unsigned int,
288 const Point<dim> &) const
289
290{
291 Assert(false, (typename FiniteElement<dim, spacedim>::ExcFENotPrimitive()));
292 return 0.;
293}
294
295
296
297template <int dim, int spacedim>
298double
299FE_PolyTensor<dim, spacedim>::shape_value_component(
300 const unsigned int i,
301 const Point<dim> & p,
302 const unsigned int component) const
303{
304 AssertIndexRange(i, this->n_dofs_per_cell());
305 AssertIndexRange(component, dim);
306
307 std::lock_guard<std::mutex> lock(cache_mutex);
308
309 if (cached_point != p || cached_values.size() == 0)
310 {
311 cached_point = p;
312 cached_values.resize(poly_space->n());
313
314 std::vector<Tensor<4, dim>> dummy1;
315 std::vector<Tensor<5, dim>> dummy2;
316 poly_space->evaluate(
317 p, cached_values, cached_grads, cached_grad_grads, dummy1, dummy2);
318 }
319
320 double s = 0;
321 if (inverse_node_matrix.n_cols() == 0)
322 return cached_values[i][component];
323 else
324 for (unsigned int j = 0; j < inverse_node_matrix.n_cols(); ++j)
325 s += inverse_node_matrix(j, i) * cached_values[j][component];
326 return s;
327}
328
329
330
331template <int dim, int spacedim>
332Tensor<1, dim>
333FE_PolyTensor<dim, spacedim>::shape_grad(const unsigned int,
334 const Point<dim> &) const
335{
336 Assert(false, (typename FiniteElement<dim, spacedim>::ExcFENotPrimitive()));
337 return Tensor<1, dim>();
338}
339
340
341
342template <int dim, int spacedim>
343Tensor<1, dim>
344FE_PolyTensor<dim, spacedim>::shape_grad_component(
345 const unsigned int i,
346 const Point<dim> & p,
347 const unsigned int component) const
348{
349 AssertIndexRange(i, this->n_dofs_per_cell());
350 AssertIndexRange(component, dim);
351
352 std::lock_guard<std::mutex> lock(cache_mutex);
353
354 if (cached_point != p || cached_grads.size() == 0)
355 {
356 cached_point = p;
357 cached_grads.resize(poly_space->n());
358
359 std::vector<Tensor<4, dim>> dummy1;
360 std::vector<Tensor<5, dim>> dummy2;
361 poly_space->evaluate(
362 p, cached_values, cached_grads, cached_grad_grads, dummy1, dummy2);
363 }
364
365 Tensor<1, dim> s;
366 if (inverse_node_matrix.n_cols() == 0)
367 return cached_grads[i][component];
368 else
369 for (unsigned int j = 0; j < inverse_node_matrix.n_cols(); ++j)
370 s += inverse_node_matrix(j, i) * cached_grads[j][component];
371
372 return s;
373}
374
375
376
377template <int dim, int spacedim>
378Tensor<2, dim>
379FE_PolyTensor<dim, spacedim>::shape_grad_grad(const unsigned int,
380 const Point<dim> &) const
381{
382 Assert(false, (typename FiniteElement<dim, spacedim>::ExcFENotPrimitive()));
383 return Tensor<2, dim>();
384}
385
386
387
388template <int dim, int spacedim>
389Tensor<2, dim>
390FE_PolyTensor<dim, spacedim>::shape_grad_grad_component(
391 const unsigned int i,
392 const Point<dim> & p,
393 const unsigned int component) const
394{
395 AssertIndexRange(i, this->n_dofs_per_cell());
396 AssertIndexRange(component, dim);
397
398 std::lock_guard<std::mutex> lock(cache_mutex);
399
400 if (cached_point != p || cached_grad_grads.size() == 0)
401 {
402 cached_point = p;
403 cached_grad_grads.resize(poly_space->n());
404
405 std::vector<Tensor<4, dim>> dummy1;
406 std::vector<Tensor<5, dim>> dummy2;
407 poly_space->evaluate(
408 p, cached_values, cached_grads, cached_grad_grads, dummy1, dummy2);
409 }
410
411 Tensor<2, dim> s;
412 if (inverse_node_matrix.n_cols() == 0)
413 return cached_grad_grads[i][component];
414 else
415 for (unsigned int j = 0; j < inverse_node_matrix.n_cols(); ++j)
416 s += inverse_node_matrix(i, j) * cached_grad_grads[j][component];
417
418 return s;
419}
420
421
422//---------------------------------------------------------------------------
423// Fill data of FEValues
424//---------------------------------------------------------------------------
425
426template <int dim, int spacedim>
427void
428FE_PolyTensor<dim, spacedim>::fill_fe_values(
429 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
430 const CellSimilarity::Similarity cell_similarity,
431 const Quadrature<dim> & quadrature,
432 const Mapping<dim, spacedim> & mapping,
433 const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
434 const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
435 & mapping_data,
436 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
437 ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
438 spacedim>
439 &output_data) const
440{
441 // convert data object to internal
442 // data for this class. fails with
443 // an exception if that is not
444 // possible
445 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
446 ExcInternalError());
447 const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
448
449 const unsigned int n_q_points = quadrature.size();
450
451 Assert(!(fe_data.update_each & update_values) ||
452 fe_data.shape_values.size()[0] == this->n_dofs_per_cell(),
453 ExcDimensionMismatch(fe_data.shape_values.size()[0],
454 this->n_dofs_per_cell()));
455 Assert(!(fe_data.update_each & update_values) ||
456 fe_data.shape_values.size()[1] == n_q_points,
457 ExcDimensionMismatch(fe_data.shape_values.size()[1], n_q_points));
458
459 // TODO: The dof_sign_change only affects Nedelec elements and is not the
460 // correct thing on complicated meshes for higher order Nedelec elements.
461 // Something similar to FE_Q should be done to permute dofs and to change the
462 // dof signs. A static way using tables (as done in the RaviartThomas<dim>
463 // class) is preferable.
464 std::fill(fe_data.dof_sign_change.begin(),
465 fe_data.dof_sign_change.end(),
466 1.0);
467 internal::FE_PolyTensor::get_dof_sign_change_nedelec(cell,
468 *this,
469 this->mapping_kind,
470 fe_data.dof_sign_change);
471
472 // TODO: This, similarly to the Nedelec case, is just a legacy function in 2d
473 // and affects only face_dofs of H(div) conformal FEs. It does nothing in 1d.
474 // Also nothing in 3d since we take care of it by using the
475 // adjust_quad_dof_sign_for_face_orientation_table.
476 internal::FE_PolyTensor::get_dof_sign_change_h_div(cell,
477 *this,
478 this->mapping_kind,
479 fe_data.dof_sign_change);
480
481 // What is the first dof_index on a quad?
482 const unsigned int first_quad_index = this->get_first_quad_index();
483 // How many dofs per quad and how many quad dofs do we have at all?
484 const unsigned int n_dofs_per_quad = this->n_dofs_per_quad();
485 const unsigned int n_quad_dofs =
486 n_dofs_per_quad * GeometryInfo<dim>::faces_per_cell;
487
488 for (unsigned int dof_index = 0; dof_index < this->n_dofs_per_cell();
489 ++dof_index)
490 {
491 /*
492 * This assumes that the dofs are ordered by first vertices, lines, quads
493 * and volume dofs. Note that in 2d this always gives false.
494 */
495 const bool is_quad_dof =
496 (dim == 2 ? false :
497 (first_quad_index <= dof_index) &&
498 (dof_index < first_quad_index + n_quad_dofs));
499
500 // TODO: This hack is not pretty and it is only here to handle the 2d
501 // case and the Nedelec legacy case. In 2d dof_sign of a face_dof is never
502 // handled by the
503 // >>if(is_quad_dof){...}<< but still a possible dof sign change must be
504 // handled, also for line_dofs in 3d such as in Nedelec. In these cases
505 // this is encoded in the array fe_data.dof_sign_change[dof_index]. In 3d
506 // it is handles with a table. This array is allocated in
507 // fe_poly_tensor.h.
508 double dof_sign = 1.0;
509 // under some circumstances fe_data.dof_sign_change is not allocated
510 if (fe_data.update_each & update_values)
511 dof_sign = fe_data.dof_sign_change[dof_index];
512
513 if (is_quad_dof)
514 {
515 /*
516 * Find the face belonging to this dof_index. This is integer
517 * division.
518 */
519 const unsigned int face_index_from_dof_index =
520 (dof_index - first_quad_index) / (n_dofs_per_quad);
521
522 const unsigned int local_quad_dof_index = dof_index % n_dofs_per_quad;
523
524 // Correct the dof_sign if necessary
525 if (adjust_quad_dof_sign_for_face_orientation(
526 local_quad_dof_index,
527 face_index_from_dof_index,
528 cell->face_orientation(face_index_from_dof_index),
529 cell->face_flip(face_index_from_dof_index),
530 cell->face_rotation(face_index_from_dof_index)))
531 dof_sign = -1.0;
532 }
533
534 const MappingKind mapping_kind = get_mapping_kind(dof_index);
535
536 const unsigned int first =
537 output_data.shape_function_to_row_table
538 [dof_index * this->n_components() +
539 this->get_nonzero_components(dof_index).first_selected_component()];
540
541 // update the shape function values as necessary
542 //
543 // we only need to do this if the current cell is not a translation of
544 // the previous one; or, even if it is a translation, if we use
545 // mappings other than the standard mappings that require us to
546 // recompute values and derivatives because of possible sign changes
547 if (fe_data.update_each & update_values &&
548 ((cell_similarity != CellSimilarity::translation) ||
549 ((mapping_kind == mapping_piola) ||
550 (mapping_kind == mapping_raviart_thomas) ||
551 (mapping_kind == mapping_nedelec))))
552 {
553 switch (mapping_kind)
554 {
555 case mapping_none:
556 {
557 for (unsigned int k = 0; k < n_q_points; ++k)
558 for (unsigned int d = 0; d < dim; ++d)
559 output_data.shape_values(first + d, k) =
560 fe_data.shape_values[dof_index][k][d];
561 break;
562 }
563
564 case mapping_covariant:
565 case mapping_contravariant:
566 {
567 mapping.transform(
568 make_array_view(fe_data.shape_values, dof_index),
569 mapping_kind,
570 mapping_internal,
571 make_array_view(fe_data.transformed_shape_values));
572
573 for (unsigned int k = 0; k < n_q_points; ++k)
574 for (unsigned int d = 0; d < dim; ++d)
575 output_data.shape_values(first + d, k) =
576 fe_data.transformed_shape_values[k][d];
577
578 break;
579 }
580
581 case mapping_raviart_thomas:
582 case mapping_piola:
583 {
584 mapping.transform(
585 make_array_view(fe_data.shape_values, dof_index),
586 mapping_piola,
587 mapping_internal,
588 make_array_view(fe_data.transformed_shape_values));
589 for (unsigned int k = 0; k < n_q_points; ++k)
590 for (unsigned int d = 0; d < dim; ++d)
591 output_data.shape_values(first + d, k) =
592 dof_sign * fe_data.transformed_shape_values[k][d];
593 break;
594 }
595
596 case mapping_nedelec:
597 {
598 mapping.transform(
599 make_array_view(fe_data.shape_values, dof_index),
600 mapping_covariant,
601 mapping_internal,
602 make_array_view(fe_data.transformed_shape_values));
603
604 for (unsigned int k = 0; k < n_q_points; ++k)
605 for (unsigned int d = 0; d < dim; ++d)
606 output_data.shape_values(first + d, k) =
607 dof_sign * fe_data.transformed_shape_values[k][d];
608
609 break;
610 }
611
612 default:
613 Assert(false, ExcNotImplemented());
614 }
615 }
616
617 // update gradients. apply the same logic as above
618 if (fe_data.update_each & update_gradients &&
619 ((cell_similarity != CellSimilarity::translation) ||
620 ((mapping_kind == mapping_piola) ||
621 (mapping_kind == mapping_raviart_thomas) ||
622 (mapping_kind == mapping_nedelec))))
623
624 {
625 switch (mapping_kind)
626 {
627 case mapping_none:
628 {
629 mapping.transform(
630 make_array_view(fe_data.shape_grads, dof_index),
631 mapping_covariant,
632 mapping_internal,
633 make_array_view(fe_data.transformed_shape_grads));
634 for (unsigned int k = 0; k < n_q_points; ++k)
635 for (unsigned int d = 0; d < dim; ++d)
636 output_data.shape_gradients[first + d][k] =
637 fe_data.transformed_shape_grads[k][d];
638 break;
639 }
640 case mapping_covariant:
641 {
642 mapping.transform(
643 make_array_view(fe_data.shape_grads, dof_index),
644 mapping_covariant_gradient,
645 mapping_internal,
646 make_array_view(fe_data.transformed_shape_grads));
647
648 for (unsigned int k = 0; k < n_q_points; ++k)
649 for (unsigned int d = 0; d < spacedim; ++d)
650 for (unsigned int n = 0; n < spacedim; ++n)
651 fe_data.transformed_shape_grads[k][d] -=
652 output_data.shape_values(first + n, k) *
653 mapping_data.jacobian_pushed_forward_grads[k][n][d];
654
655 for (unsigned int k = 0; k < n_q_points; ++k)
656 for (unsigned int d = 0; d < dim; ++d)
657 output_data.shape_gradients[first + d][k] =
658 fe_data.transformed_shape_grads[k][d];
659
660 break;
661 }
662 case mapping_contravariant:
663 {
664 for (unsigned int k = 0; k < n_q_points; ++k)
665 fe_data.untransformed_shape_grads[k] =
666 fe_data.shape_grads[dof_index][k];
667 mapping.transform(
668 make_array_view(fe_data.untransformed_shape_grads),
669 mapping_contravariant_gradient,
670 mapping_internal,
671 make_array_view(fe_data.transformed_shape_grads));
672
673 for (unsigned int k = 0; k < n_q_points; ++k)
674 for (unsigned int d = 0; d < spacedim; ++d)
675 for (unsigned int n = 0; n < spacedim; ++n)
676 fe_data.transformed_shape_grads[k][d] +=
677 output_data.shape_values(first + n, k) *
678 mapping_data.jacobian_pushed_forward_grads[k][d][n];
679
680
681 for (unsigned int k = 0; k < n_q_points; ++k)
682 for (unsigned int d = 0; d < dim; ++d)
683 output_data.shape_gradients[first + d][k] =
684 fe_data.transformed_shape_grads[k][d];
685
686 break;
687 }
688 case mapping_raviart_thomas:
689 case mapping_piola:
690 {
691 for (unsigned int k = 0; k < n_q_points; ++k)
692 fe_data.untransformed_shape_grads[k] =
693 fe_data.shape_grads[dof_index][k];
694 mapping.transform(
695 make_array_view(fe_data.untransformed_shape_grads),
696 mapping_piola_gradient,
697 mapping_internal,
698 make_array_view(fe_data.transformed_shape_grads));
699
700 for (unsigned int k = 0; k < n_q_points; ++k)
701 for (unsigned int d = 0; d < spacedim; ++d)
702 for (unsigned int n = 0; n < spacedim; ++n)
703 fe_data.transformed_shape_grads[k][d] +=
704 (output_data.shape_values(first + n, k) *
705 mapping_data
706 .jacobian_pushed_forward_grads[k][d][n]) -
707 (output_data.shape_values(first + d, k) *
708 mapping_data.jacobian_pushed_forward_grads[k][n][n]);
709
710 for (unsigned int k = 0; k < n_q_points; ++k)
711 for (unsigned int d = 0; d < dim; ++d)
712 output_data.shape_gradients[first + d][k] =
713 dof_sign * fe_data.transformed_shape_grads[k][d];
714
715 break;
716 }
717
718 case mapping_nedelec:
719 {
720 // treat the gradients of
721 // this particular shape
722 // function at all
723 // q-points. if Dv is the
724 // gradient of the shape
725 // function on the unit
726 // cell, then
727 // (J^-T)Dv(J^-1) is the
728 // value we want to have on
729 // the real cell.
730 for (unsigned int k = 0; k < n_q_points; ++k)
731 fe_data.untransformed_shape_grads[k] =
732 fe_data.shape_grads[dof_index][k];
733
734 mapping.transform(
735 make_array_view(fe_data.untransformed_shape_grads),
736 mapping_covariant_gradient,
737 mapping_internal,
738 make_array_view(fe_data.transformed_shape_grads));
739
740 for (unsigned int k = 0; k < n_q_points; ++k)
741 for (unsigned int d = 0; d < spacedim; ++d)
742 for (unsigned int n = 0; n < spacedim; ++n)
743 fe_data.transformed_shape_grads[k][d] -=
744 output_data.shape_values(first + n, k) *
745 mapping_data.jacobian_pushed_forward_grads[k][n][d];
746
747 for (unsigned int k = 0; k < n_q_points; ++k)
748 for (unsigned int d = 0; d < dim; ++d)
749 output_data.shape_gradients[first + d][k] =
750 dof_sign * fe_data.transformed_shape_grads[k][d];
751
752 break;
753 }
754
755 default:
756 Assert(false, ExcNotImplemented());
757 }
758 }
759
760 // update hessians. apply the same logic as above
761 if (fe_data.update_each & update_hessians &&
762 ((cell_similarity != CellSimilarity::translation) ||
763 ((mapping_kind == mapping_piola) ||
764 (mapping_kind == mapping_raviart_thomas) ||
765 (mapping_kind == mapping_nedelec))))
766
767 {
768 switch (mapping_kind)
769 {
770 case mapping_none:
771 {
772 mapping.transform(
773 make_array_view(fe_data.shape_grad_grads, dof_index),
774 mapping_covariant_gradient,
775 mapping_internal,
776 make_array_view(fe_data.transformed_shape_hessians));
777
778 for (unsigned int k = 0; k < n_q_points; ++k)
779 for (unsigned int d = 0; d < spacedim; ++d)
780 for (unsigned int n = 0; n < spacedim; ++n)
781 fe_data.transformed_shape_hessians[k][d] -=
782 output_data.shape_gradients[first + d][k][n] *
783 mapping_data.jacobian_pushed_forward_grads[k][n];
784
785 for (unsigned int k = 0; k < n_q_points; ++k)
786 for (unsigned int d = 0; d < dim; ++d)
787 output_data.shape_hessians[first + d][k] =
788 fe_data.transformed_shape_hessians[k][d];
789
790 break;
791 }
792 case mapping_covariant:
793 {
794 for (unsigned int k = 0; k < n_q_points; ++k)
795 fe_data.untransformed_shape_hessian_tensors[k] =
796 fe_data.shape_grad_grads[dof_index][k];
797
798 mapping.transform(
799 make_array_view(
800 fe_data.untransformed_shape_hessian_tensors),
801 mapping_covariant_hessian,
802 mapping_internal,
803 make_array_view(fe_data.transformed_shape_hessians));
804
805 for (unsigned int k = 0; k < n_q_points; ++k)
806 for (unsigned int d = 0; d < spacedim; ++d)
807 for (unsigned int n = 0; n < spacedim; ++n)
808 for (unsigned int i = 0; i < spacedim; ++i)
809 for (unsigned int j = 0; j < spacedim; ++j)
810 {
811 fe_data.transformed_shape_hessians[k][d][i][j] -=
812 (output_data.shape_values(first + n, k) *
813 mapping_data
814 .jacobian_pushed_forward_2nd_derivatives
815 [k][n][d][i][j]) +
816 (output_data.shape_gradients[first + d][k][n] *
817 mapping_data
818 .jacobian_pushed_forward_grads[k][n][i][j]) +
819 (output_data.shape_gradients[first + n][k][i] *
820 mapping_data
821 .jacobian_pushed_forward_grads[k][n][d][j]) +
822 (output_data.shape_gradients[first + n][k][j] *
823 mapping_data
824 .jacobian_pushed_forward_grads[k][n][i][d]);
825 }
826
827 for (unsigned int k = 0; k < n_q_points; ++k)
828 for (unsigned int d = 0; d < dim; ++d)
829 output_data.shape_hessians[first + d][k] =
830 fe_data.transformed_shape_hessians[k][d];
831
832 break;
833 }
834 case mapping_contravariant:
835 {
836 for (unsigned int k = 0; k < n_q_points; ++k)
837 fe_data.untransformed_shape_hessian_tensors[k] =
838 fe_data.shape_grad_grads[dof_index][k];
839
840 mapping.transform(
841 make_array_view(
842 fe_data.untransformed_shape_hessian_tensors),
843 mapping_contravariant_hessian,
844 mapping_internal,
845 make_array_view(fe_data.transformed_shape_hessians));
846
847 for (unsigned int k = 0; k < n_q_points; ++k)
848 for (unsigned int d = 0; d < spacedim; ++d)
849 for (unsigned int n = 0; n < spacedim; ++n)
850 for (unsigned int i = 0; i < spacedim; ++i)
851 for (unsigned int j = 0; j < spacedim; ++j)
852 {
853 fe_data.transformed_shape_hessians[k][d][i][j] +=
854 (output_data.shape_values(first + n, k) *
855 mapping_data
856 .jacobian_pushed_forward_2nd_derivatives
857 [k][d][n][i][j]) +
858 (output_data.shape_gradients[first + n][k][i] *
859 mapping_data
860 .jacobian_pushed_forward_grads[k][d][n][j]) +
861 (output_data.shape_gradients[first + n][k][j] *
862 mapping_data
863 .jacobian_pushed_forward_grads[k][d][i][n]) -
864 (output_data.shape_gradients[first + d][k][n] *
865 mapping_data
866 .jacobian_pushed_forward_grads[k][n][i][j]);
867 for (unsigned int m = 0; m < spacedim; ++m)
868 fe_data
869 .transformed_shape_hessians[k][d][i][j] -=
870 (mapping_data
871 .jacobian_pushed_forward_grads[k][d][i]
872 [m] *
873 mapping_data
874 .jacobian_pushed_forward_grads[k][m][n]
875 [j] *
876 output_data.shape_values(first + n, k)) +
877 (mapping_data
878 .jacobian_pushed_forward_grads[k][d][m]
879 [j] *
880 mapping_data
881 .jacobian_pushed_forward_grads[k][m][i]
882 [n] *
883 output_data.shape_values(first + n, k));
884 }
885
886 for (unsigned int k = 0; k < n_q_points; ++k)
887 for (unsigned int d = 0; d < dim; ++d)
888 output_data.shape_hessians[first + d][k] =
889 fe_data.transformed_shape_hessians[k][d];
890
891 break;
892 }
893 case mapping_raviart_thomas:
894 case mapping_piola:
895 {
896 for (unsigned int k = 0; k < n_q_points; ++k)
897 fe_data.untransformed_shape_hessian_tensors[k] =
898 fe_data.shape_grad_grads[dof_index][k];
899
900 mapping.transform(
901 make_array_view(
902 fe_data.untransformed_shape_hessian_tensors),
903 mapping_piola_hessian,
904 mapping_internal,
905 make_array_view(fe_data.transformed_shape_hessians));
906
907 for (unsigned int k = 0; k < n_q_points; ++k)
908 for (unsigned int d = 0; d < spacedim; ++d)
909 for (unsigned int n = 0; n < spacedim; ++n)
910 for (unsigned int i = 0; i < spacedim; ++i)
911 for (unsigned int j = 0; j < spacedim; ++j)
912 {
913 fe_data.transformed_shape_hessians[k][d][i][j] +=
914 (output_data.shape_values(first + n, k) *
915 mapping_data
916 .jacobian_pushed_forward_2nd_derivatives
917 [k][d][n][i][j]) +
918 (output_data.shape_gradients[first + n][k][i] *
919 mapping_data
920 .jacobian_pushed_forward_grads[k][d][n][j]) +
921 (output_data.shape_gradients[first + n][k][j] *
922 mapping_data
923 .jacobian_pushed_forward_grads[k][d][i][n]) -
924 (output_data.shape_gradients[first + d][k][n] *
925 mapping_data
926 .jacobian_pushed_forward_grads[k][n][i][j]);
927
928 fe_data.transformed_shape_hessians[k][d][i][j] -=
929 (output_data.shape_values(first + d, k) *
930 mapping_data
931 .jacobian_pushed_forward_2nd_derivatives
932 [k][n][n][i][j]) +
933 (output_data.shape_gradients[first + d][k][i] *
934 mapping_data
935 .jacobian_pushed_forward_grads[k][n][n][j]) +
936 (output_data.shape_gradients[first + d][k][j] *
937 mapping_data
938 .jacobian_pushed_forward_grads[k][n][n][i]);
939
940 for (unsigned int m = 0; m < spacedim; ++m)
941 {
942 fe_data
943 .transformed_shape_hessians[k][d][i][j] -=
944 (mapping_data
945 .jacobian_pushed_forward_grads[k][d][i]
946 [m] *
947 mapping_data
948 .jacobian_pushed_forward_grads[k][m][n]
949 [j] *
950 output_data.shape_values(first + n, k)) +
951 (mapping_data
952 .jacobian_pushed_forward_grads[k][d][m]
953 [j] *
954 mapping_data
955 .jacobian_pushed_forward_grads[k][m][i]
956 [n] *
957 output_data.shape_values(first + n, k));
958
959 fe_data
960 .transformed_shape_hessians[k][d][i][j] +=
961 (mapping_data
962 .jacobian_pushed_forward_grads[k][n][i]
963 [m] *
964 mapping_data
965 .jacobian_pushed_forward_grads[k][m][n]
966 [j] *
967 output_data.shape_values(first + d, k)) +
968 (mapping_data
969 .jacobian_pushed_forward_grads[k][n][m]
970 [j] *
971 mapping_data
972 .jacobian_pushed_forward_grads[k][m][i]
973 [n] *
974 output_data.shape_values(first + d, k));
975 }
976 }
977
978 for (unsigned int k = 0; k < n_q_points; ++k)
979 for (unsigned int d = 0; d < dim; ++d)
980 output_data.shape_hessians[first + d][k] =
981 dof_sign * fe_data.transformed_shape_hessians[k][d];
982
983 break;
984 }
985
986 case mapping_nedelec:
987 {
988 for (unsigned int k = 0; k < n_q_points; ++k)
989 fe_data.untransformed_shape_hessian_tensors[k] =
990 fe_data.shape_grad_grads[dof_index][k];
991
992 mapping.transform(
993 make_array_view(
994 fe_data.untransformed_shape_hessian_tensors),
995 mapping_covariant_hessian,
996 mapping_internal,
997 make_array_view(fe_data.transformed_shape_hessians));
998
999 for (unsigned int k = 0; k < n_q_points; ++k)
1000 for (unsigned int d = 0; d < spacedim; ++d)
1001 for (unsigned int n = 0; n < spacedim; ++n)
1002 for (unsigned int i = 0; i < spacedim; ++i)
1003 for (unsigned int j = 0; j < spacedim; ++j)
1004 {
1005 fe_data.transformed_shape_hessians[k][d][i][j] -=
1006 (output_data.shape_values(first + n, k) *
1007 mapping_data
1008 .jacobian_pushed_forward_2nd_derivatives
1009 [k][n][d][i][j]) +
1010 (output_data.shape_gradients[first + d][k][n] *
1011 mapping_data
1012 .jacobian_pushed_forward_grads[k][n][i][j]) +
1013 (output_data.shape_gradients[first + n][k][i] *
1014 mapping_data
1015 .jacobian_pushed_forward_grads[k][n][d][j]) +
1016 (output_data.shape_gradients[first + n][k][j] *
1017 mapping_data
1018 .jacobian_pushed_forward_grads[k][n][i][d]);
1019 }
1020
1021 for (unsigned int k = 0; k < n_q_points; ++k)
1022 for (unsigned int d = 0; d < dim; ++d)
1023 output_data.shape_hessians[first + d][k] =
1024 dof_sign * fe_data.transformed_shape_hessians[k][d];
1025
1026 break;
1027 }
1028
1029 default:
1030 Assert(false, ExcNotImplemented());
1031 }
1032 }
1033
1034 // third derivatives are not implemented
1035 if (fe_data.update_each & update_3rd_derivatives &&
1036 ((cell_similarity != CellSimilarity::translation) ||
1037 ((mapping_kind == mapping_piola) ||
1038 (mapping_kind == mapping_raviart_thomas) ||
1039 (mapping_kind == mapping_nedelec))))
1040 {
1041 Assert(false, ExcNotImplemented())
1042 }
1043 }
1044}
1045
1046
1047
1048template <int dim, int spacedim>
1049void
1050FE_PolyTensor<dim, spacedim>::fill_fe_face_values(
1051 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1052 const unsigned int face_no,
1053 const hp::QCollection<dim - 1> & quadrature,
1054 const Mapping<dim, spacedim> & mapping,
1055 const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
1056 const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1057 & mapping_data,
1058 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1059 ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
1060 spacedim>
1061 &output_data) const
1062{
1063 AssertDimension(quadrature.size(), 1);
1064
1065 // convert data object to internal
1066 // data for this class. fails with
1067 // an exception if that is not
1068 // possible
1069 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
1070 ExcInternalError());
1071 const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
1072
1073 const unsigned int n_q_points = quadrature[0].size();
1074 // offset determines which data set
1075 // to take (all data sets for all
1076 // faces are stored contiguously)
1077
1078 const auto offset =
1079 QProjector<dim>::DataSetDescriptor::face(this->reference_cell(),
1080 face_no,
1081 cell->face_orientation(face_no),
1082 cell->face_flip(face_no),
1083 cell->face_rotation(face_no),
1084 n_q_points);
1085
1086 // TODO: Size assertions
1087
1088 // TODO: The dof_sign_change only affects Nedelec elements and is not the
1089 // correct thing on complicated meshes for higher order Nedelec elements.
1090 // Something similar to FE_Q should be done to permute dofs and to change the
1091 // dof signs. A static way using tables (as done in the RaviartThomas<dim>
1092 // class) is preferable.
1093 std::fill(fe_data.dof_sign_change.begin(),
1094 fe_data.dof_sign_change.end(),
1095 1.0);
1096 internal::FE_PolyTensor::get_dof_sign_change_nedelec(cell,
1097 *this,
1098 this->mapping_kind,
1099 fe_data.dof_sign_change);
1100
1101 // TODO: This, similarly to the Nedelec case, is just a legacy function in 2d
1102 // and affects only face_dofs of H(div) conformal FEs. It does nothing in 1d.
1103 // Also nothing in 3d since we take care of it by using the
1104 // adjust_quad_dof_sign_for_face_orientation_table.
1105 internal::FE_PolyTensor::get_dof_sign_change_h_div(cell,
1106 *this,
1107 this->mapping_kind,
1108 fe_data.dof_sign_change);
1109
1110 // What is the first dof_index on a quad?
1111 const unsigned int first_quad_index = this->get_first_quad_index();
1112 // How many dofs per quad and how many quad dofs do we have at all?
1113 const unsigned int n_dofs_per_quad = this->n_dofs_per_quad();
1114 const unsigned int n_quad_dofs =
1115 n_dofs_per_quad * GeometryInfo<dim>::faces_per_cell;
1116
1117 for (unsigned int dof_index = 0; dof_index < this->n_dofs_per_cell();
1118 ++dof_index)
1119 {
1120 /*
1121 * This assumes that the dofs are ordered by first vertices, lines, quads
1122 * and volume dofs. Note that in 2d this always gives false.
1123 */
1124 const bool is_quad_dof =
1125 (dim == 2 ? false :
1126 (first_quad_index <= dof_index) &&
1127 (dof_index < first_quad_index + n_quad_dofs));
1128
1129 // TODO: This hack is not pretty and it is only here to handle the 2d
1130 // case and the Nedelec legacy case. In 2d dof_sign of a face_dof is never
1131 // handled by the
1132 // >>if(is_quad_dof){...}<< but still a possible dof sign change must be
1133 // handled, also for line_dofs in 3d such as in Nedelec. In these cases
1134 // this is encoded in the array fe_data.dof_sign_change[dof_index]. In 3d
1135 // it is handles with a table. This array is allocated in
1136 // fe_poly_tensor.h.
1137 double dof_sign = 1.0;
1138 // under some circumstances fe_data.dof_sign_change is not allocated
1139 if (fe_data.update_each & update_values)
1140 dof_sign = fe_data.dof_sign_change[dof_index];
1141
1142 if (is_quad_dof)
1143 {
1144 /*
1145 * Find the face belonging to this dof_index. This is integer
1146 * division.
1147 */
1148 unsigned int face_index_from_dof_index =
1149 (dof_index - first_quad_index) / (n_dofs_per_quad);
1150
1151 unsigned int local_quad_dof_index = dof_index % n_dofs_per_quad;
1152
1153 // Correct the dof_sign if necessary
1154 if (adjust_quad_dof_sign_for_face_orientation(
1155 local_quad_dof_index,
1156 face_index_from_dof_index,
1157 cell->face_orientation(face_index_from_dof_index),
1158 cell->face_flip(face_index_from_dof_index),
1159 cell->face_rotation(face_index_from_dof_index)))
1160 dof_sign = -1.0;
1161 }
1162
1163 const MappingKind mapping_kind = get_mapping_kind(dof_index);
1164
1165 const unsigned int first =
1166 output_data.shape_function_to_row_table
1167 [dof_index * this->n_components() +
1168 this->get_nonzero_components(dof_index).first_selected_component()];
1169
1170 if (fe_data.update_each & update_values)
1171 {
1172 switch (mapping_kind)
1173 {
1174 case mapping_none:
1175 {
1176 for (unsigned int k = 0; k < n_q_points; ++k)
1177 for (unsigned int d = 0; d < dim; ++d)
1178 output_data.shape_values(first + d, k) =
1179 fe_data.shape_values[dof_index][k + offset][d];
1180 break;
1181 }
1182
1183 case mapping_covariant:
1184 case mapping_contravariant:
1185 {
1186 const ArrayView<Tensor<1, spacedim>>
1187 transformed_shape_values =
1188 make_array_view(fe_data.transformed_shape_values,
1189 offset,
1190 n_q_points);
1191 mapping.transform(make_array_view(fe_data.shape_values,
1192 dof_index,
1193 offset,
1194 n_q_points),
1195 mapping_kind,
1196 mapping_internal,
1197 transformed_shape_values);
1198
1199 for (unsigned int k = 0; k < n_q_points; ++k)
1200 for (unsigned int d = 0; d < dim; ++d)
1201 output_data.shape_values(first + d, k) =
1202 transformed_shape_values[k][d];
1203
1204 break;
1205 }
1206 case mapping_raviart_thomas:
1207 case mapping_piola:
1208 {
1209 const ArrayView<Tensor<1, spacedim>>
1210 transformed_shape_values =
1211 make_array_view(fe_data.transformed_shape_values,
1212 offset,
1213 n_q_points);
1214 mapping.transform(make_array_view(fe_data.shape_values,
1215 dof_index,
1216 offset,
1217 n_q_points),
1218 mapping_piola,
1219 mapping_internal,
1220 transformed_shape_values);
1221 for (unsigned int k = 0; k < n_q_points; ++k)
1222 for (unsigned int d = 0; d < dim; ++d)
1223 output_data.shape_values(first + d, k) =
1224 dof_sign * transformed_shape_values[k][d];
1225 break;
1226 }
1227
1228 case mapping_nedelec:
1229 {
1230 const ArrayView<Tensor<1, spacedim>>
1231 transformed_shape_values =
1232 make_array_view(fe_data.transformed_shape_values,
1233 offset,
1234 n_q_points);
1235 mapping.transform(make_array_view(fe_data.shape_values,
1236 dof_index,
1237 offset,
1238 n_q_points),
1239 mapping_covariant,
1240 mapping_internal,
1241 transformed_shape_values);
1242
1243 for (unsigned int k = 0; k < n_q_points; ++k)
1244 for (unsigned int d = 0; d < dim; ++d)
1245 output_data.shape_values(first + d, k) =
1246 dof_sign * transformed_shape_values[k][d];
1247
1248 break;
1249 }
1250
1251 default:
1252 Assert(false, ExcNotImplemented());
1253 }
1254 }
1255
1256 if (fe_data.update_each & update_gradients)
1257 {
1258 switch (mapping_kind)
1259 {
1260 case mapping_none:
1261 {
1262 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1263 make_array_view(fe_data.transformed_shape_grads,
1264 offset,
1265 n_q_points);
1266 mapping.transform(make_array_view(fe_data.shape_grads,
1267 dof_index,
1268 offset,
1269 n_q_points),
1270 mapping_covariant,
1271 mapping_internal,
1272 transformed_shape_grads);
1273 for (unsigned int k = 0; k < n_q_points; ++k)
1274 for (unsigned int d = 0; d < dim; ++d)
1275 output_data.shape_gradients[first + d][k] =
1276 transformed_shape_grads[k][d];
1277 break;
1278 }
1279
1280 case mapping_covariant:
1281 {
1282 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1283 make_array_view(fe_data.transformed_shape_grads,
1284 offset,
1285 n_q_points);
1286 mapping.transform(make_array_view(fe_data.shape_grads,
1287 dof_index,
1288 offset,
1289 n_q_points),
1290 mapping_covariant_gradient,
1291 mapping_internal,
1292 transformed_shape_grads);
1293
1294 for (unsigned int k = 0; k < n_q_points; ++k)
1295 for (unsigned int d = 0; d < spacedim; ++d)
1296 for (unsigned int n = 0; n < spacedim; ++n)
1297 transformed_shape_grads[k][d] -=
1298 output_data.shape_values(first + n, k) *
1299 mapping_data.jacobian_pushed_forward_grads[k][n][d];
1300
1301 for (unsigned int k = 0; k < n_q_points; ++k)
1302 for (unsigned int d = 0; d < dim; ++d)
1303 output_data.shape_gradients[first + d][k] =
1304 transformed_shape_grads[k][d];
1305 break;
1306 }
1307
1308 case mapping_contravariant:
1309 {
1310 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1311 make_array_view(fe_data.transformed_shape_grads,
1312 offset,
1313 n_q_points);
1314 for (unsigned int k = 0; k < n_q_points; ++k)
1315 fe_data.untransformed_shape_grads[k + offset] =
1316 fe_data.shape_grads[dof_index][k + offset];
1317 mapping.transform(
1318 make_array_view(fe_data.untransformed_shape_grads,
1319 offset,
1320 n_q_points),
1321 mapping_contravariant_gradient,
1322 mapping_internal,
1323 transformed_shape_grads);
1324
1325 for (unsigned int k = 0; k < n_q_points; ++k)
1326 for (unsigned int d = 0; d < spacedim; ++d)
1327 for (unsigned int n = 0; n < spacedim; ++n)
1328 transformed_shape_grads[k][d] +=
1329 output_data.shape_values(first + n, k) *
1330 mapping_data.jacobian_pushed_forward_grads[k][d][n];
1331
1332 for (unsigned int k = 0; k < n_q_points; ++k)
1333 for (unsigned int d = 0; d < dim; ++d)
1334 output_data.shape_gradients[first + d][k] =
1335 transformed_shape_grads[k][d];
1336
1337 break;
1338 }
1339
1340 case mapping_raviart_thomas:
1341 case mapping_piola:
1342 {
1343 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1344 make_array_view(fe_data.transformed_shape_grads,
1345 offset,
1346 n_q_points);
1347 for (unsigned int k = 0; k < n_q_points; ++k)
1348 fe_data.untransformed_shape_grads[k + offset] =
1349 fe_data.shape_grads[dof_index][k + offset];
1350 mapping.transform(
1351 make_array_view(fe_data.untransformed_shape_grads,
1352 offset,
1353 n_q_points),
1354 mapping_piola_gradient,
1355 mapping_internal,
1356 transformed_shape_grads);
1357
1358 for (unsigned int k = 0; k < n_q_points; ++k)
1359 for (unsigned int d = 0; d < spacedim; ++d)
1360 for (unsigned int n = 0; n < spacedim; ++n)
1361 transformed_shape_grads[k][d] +=
1362 (output_data.shape_values(first + n, k) *
1363 mapping_data
1364 .jacobian_pushed_forward_grads[k][d][n]) -
1365 (output_data.shape_values(first + d, k) *
1366 mapping_data.jacobian_pushed_forward_grads[k][n][n]);
1367
1368 for (unsigned int k = 0; k < n_q_points; ++k)
1369 for (unsigned int d = 0; d < dim; ++d)
1370 output_data.shape_gradients[first + d][k] =
1371 dof_sign * transformed_shape_grads[k][d];
1372
1373 break;
1374 }
1375
1376 case mapping_nedelec:
1377 {
1378 // treat the gradients of
1379 // this particular shape
1380 // function at all
1381 // q-points. if Dv is the
1382 // gradient of the shape
1383 // function on the unit
1384 // cell, then
1385 // (J^-T)Dv(J^-1) is the
1386 // value we want to have on
1387 // the real cell.
1388 for (unsigned int k = 0; k < n_q_points; ++k)
1389 fe_data.untransformed_shape_grads[k + offset] =
1390 fe_data.shape_grads[dof_index][k + offset];
1391
1392 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1393 make_array_view(fe_data.transformed_shape_grads,
1394 offset,
1395 n_q_points);
1396 mapping.transform(
1397 make_array_view(fe_data.untransformed_shape_grads,
1398 offset,
1399 n_q_points),
1400 mapping_covariant_gradient,
1401 mapping_internal,
1402 transformed_shape_grads);
1403
1404 for (unsigned int k = 0; k < n_q_points; ++k)
1405 for (unsigned int d = 0; d < spacedim; ++d)
1406 for (unsigned int n = 0; n < spacedim; ++n)
1407 transformed_shape_grads[k][d] -=
1408 output_data.shape_values(first + n, k) *
1409 mapping_data.jacobian_pushed_forward_grads[k][n][d];
1410
1411 for (unsigned int k = 0; k < n_q_points; ++k)
1412 for (unsigned int d = 0; d < dim; ++d)
1413 output_data.shape_gradients[first + d][k] =
1414 dof_sign * transformed_shape_grads[k][d];
1415
1416 break;
1417 }
1418
1419 default:
1420 Assert(false, ExcNotImplemented());
1421 }
1422 }
1423
1424 if (fe_data.update_each & update_hessians)
1425 {
1426 switch (mapping_kind)
1427 {
1428 case mapping_none:
1429 {
1430 const ArrayView<Tensor<3, spacedim>>
1431 transformed_shape_hessians =
1432 make_array_view(fe_data.transformed_shape_hessians,
1433 offset,
1434 n_q_points);
1435 mapping.transform(make_array_view(fe_data.shape_grad_grads,
1436 dof_index,
1437 offset,
1438 n_q_points),
1439 mapping_covariant_gradient,
1440 mapping_internal,
1441 transformed_shape_hessians);
1442
1443 for (unsigned int k = 0; k < n_q_points; ++k)
1444 for (unsigned int d = 0; d < spacedim; ++d)
1445 for (unsigned int n = 0; n < spacedim; ++n)
1446 transformed_shape_hessians[k][d] -=
1447 output_data.shape_gradients[first + d][k][n] *
1448 mapping_data.jacobian_pushed_forward_grads[k][n];
1449
1450 for (unsigned int k = 0; k < n_q_points; ++k)
1451 for (unsigned int d = 0; d < dim; ++d)
1452 output_data.shape_hessians[first + d][k] =
1453 transformed_shape_hessians[k][d];
1454
1455 break;
1456 }
1457 case mapping_covariant:
1458 {
1459 for (unsigned int k = 0; k < n_q_points; ++k)
1460 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1461 fe_data.shape_grad_grads[dof_index][k + offset];
1462
1463 const ArrayView<Tensor<3, spacedim>>
1464 transformed_shape_hessians =
1465 make_array_view(fe_data.transformed_shape_hessians,
1466 offset,
1467 n_q_points);
1468 mapping.transform(
1469 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1470 offset,
1471 n_q_points),
1472 mapping_covariant_hessian,
1473 mapping_internal,
1474 transformed_shape_hessians);
1475
1476 for (unsigned int k = 0; k < n_q_points; ++k)
1477 for (unsigned int d = 0; d < spacedim; ++d)
1478 for (unsigned int n = 0; n < spacedim; ++n)
1479 for (unsigned int i = 0; i < spacedim; ++i)
1480 for (unsigned int j = 0; j < spacedim; ++j)
1481 {
1482 transformed_shape_hessians[k][d][i][j] -=
1483 (output_data.shape_values(first + n, k) *
1484 mapping_data
1485 .jacobian_pushed_forward_2nd_derivatives
1486 [k][n][d][i][j]) +
1487 (output_data.shape_gradients[first + d][k][n] *
1488 mapping_data
1489 .jacobian_pushed_forward_grads[k][n][i][j]) +
1490 (output_data.shape_gradients[first + n][k][i] *
1491 mapping_data
1492 .jacobian_pushed_forward_grads[k][n][d][j]) +
1493 (output_data.shape_gradients[first + n][k][j] *
1494 mapping_data
1495 .jacobian_pushed_forward_grads[k][n][i][d]);
1496 }
1497
1498 for (unsigned int k = 0; k < n_q_points; ++k)
1499 for (unsigned int d = 0; d < dim; ++d)
1500 output_data.shape_hessians[first + d][k] =
1501 transformed_shape_hessians[k][d];
1502
1503 break;
1504 }
1505
1506 case mapping_contravariant:
1507 {
1508 for (unsigned int k = 0; k < n_q_points; ++k)
1509 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1510 fe_data.shape_grad_grads[dof_index][k + offset];
1511
1512 const ArrayView<Tensor<3, spacedim>>
1513 transformed_shape_hessians =
1514 make_array_view(fe_data.transformed_shape_hessians,
1515 offset,
1516 n_q_points);
1517 mapping.transform(
1518 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1519 offset,
1520 n_q_points),
1521 mapping_contravariant_hessian,
1522 mapping_internal,
1523 transformed_shape_hessians);
1524
1525 for (unsigned int k = 0; k < n_q_points; ++k)
1526 for (unsigned int d = 0; d < spacedim; ++d)
1527 for (unsigned int n = 0; n < spacedim; ++n)
1528 for (unsigned int i = 0; i < spacedim; ++i)
1529 for (unsigned int j = 0; j < spacedim; ++j)
1530 {
1531 transformed_shape_hessians[k][d][i][j] +=
1532 (output_data.shape_values(first + n, k) *
1533 mapping_data
1534 .jacobian_pushed_forward_2nd_derivatives
1535 [k][d][n][i][j]) +
1536 (output_data.shape_gradients[first + n][k][i] *
1537 mapping_data
1538 .jacobian_pushed_forward_grads[k][d][n][j]) +
1539 (output_data.shape_gradients[first + n][k][j] *
1540 mapping_data
1541 .jacobian_pushed_forward_grads[k][d][i][n]) -
1542 (output_data.shape_gradients[first + d][k][n] *
1543 mapping_data
1544 .jacobian_pushed_forward_grads[k][n][i][j]);
1545 for (unsigned int m = 0; m < spacedim; ++m)
1546 transformed_shape_hessians[k][d][i][j] -=
1547 (mapping_data
1548 .jacobian_pushed_forward_grads[k][d][i]
1549 [m] *
1550 mapping_data
1551 .jacobian_pushed_forward_grads[k][m][n]
1552 [j] *
1553 output_data.shape_values(first + n, k)) +
1554 (mapping_data
1555 .jacobian_pushed_forward_grads[k][d][m]
1556 [j] *
1557 mapping_data
1558 .jacobian_pushed_forward_grads[k][m][i]
1559 [n] *
1560 output_data.shape_values(first + n, k));
1561 }
1562
1563 for (unsigned int k = 0; k < n_q_points; ++k)
1564 for (unsigned int d = 0; d < dim; ++d)
1565 output_data.shape_hessians[first + d][k] =
1566 transformed_shape_hessians[k][d];
1567
1568 break;
1569 }
1570
1571 case mapping_raviart_thomas:
1572 case mapping_piola:
1573 {
1574 for (unsigned int k = 0; k < n_q_points; ++k)
1575 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1576 fe_data.shape_grad_grads[dof_index][k + offset];
1577
1578 const ArrayView<Tensor<3, spacedim>>
1579 transformed_shape_hessians =
1580 make_array_view(fe_data.transformed_shape_hessians,
1581 offset,
1582 n_q_points);
1583 mapping.transform(
1584 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1585 offset,
1586 n_q_points),
1587 mapping_piola_hessian,
1588 mapping_internal,
1589 transformed_shape_hessians);
1590
1591 for (unsigned int k = 0; k < n_q_points; ++k)
1592 for (unsigned int d = 0; d < spacedim; ++d)
1593 for (unsigned int n = 0; n < spacedim; ++n)
1594 for (unsigned int i = 0; i < spacedim; ++i)
1595 for (unsigned int j = 0; j < spacedim; ++j)
1596 {
1597 transformed_shape_hessians[k][d][i][j] +=
1598 (output_data.shape_values(first + n, k) *
1599 mapping_data
1600 .jacobian_pushed_forward_2nd_derivatives
1601 [k][d][n][i][j]) +
1602 (output_data.shape_gradients[first + n][k][i] *
1603 mapping_data
1604 .jacobian_pushed_forward_grads[k][d][n][j]) +
1605 (output_data.shape_gradients[first + n][k][j] *
1606 mapping_data
1607 .jacobian_pushed_forward_grads[k][d][i][n]) -
1608 (output_data.shape_gradients[first + d][k][n] *
1609 mapping_data
1610 .jacobian_pushed_forward_grads[k][n][i][j]);
1611
1612 transformed_shape_hessians[k][d][i][j] -=
1613 (output_data.shape_values(first + d, k) *
1614 mapping_data
1615 .jacobian_pushed_forward_2nd_derivatives
1616 [k][n][n][i][j]) +
1617 (output_data.shape_gradients[first + d][k][i] *
1618 mapping_data
1619 .jacobian_pushed_forward_grads[k][n][n][j]) +
1620 (output_data.shape_gradients[first + d][k][j] *
1621 mapping_data
1622 .jacobian_pushed_forward_grads[k][n][n][i]);
1623
1624 for (unsigned int m = 0; m < spacedim; ++m)
1625 {
1626 transformed_shape_hessians[k][d][i][j] -=
1627 (mapping_data
1628 .jacobian_pushed_forward_grads[k][d][i]
1629 [m] *
1630 mapping_data
1631 .jacobian_pushed_forward_grads[k][m][n]
1632 [j] *
1633 output_data.shape_values(first + n, k)) +
1634 (mapping_data
1635 .jacobian_pushed_forward_grads[k][d][m]
1636 [j] *
1637 mapping_data
1638 .jacobian_pushed_forward_grads[k][m][i]
1639 [n] *
1640 output_data.shape_values(first + n, k));
1641
1642 transformed_shape_hessians[k][d][i][j] +=
1643 (mapping_data
1644 .jacobian_pushed_forward_grads[k][n][i]
1645 [m] *
1646 mapping_data
1647 .jacobian_pushed_forward_grads[k][m][n]
1648 [j] *
1649 output_data.shape_values(first + d, k)) +
1650 (mapping_data
1651 .jacobian_pushed_forward_grads[k][n][m]
1652 [j] *
1653 mapping_data
1654 .jacobian_pushed_forward_grads[k][m][i]
1655 [n] *
1656 output_data.shape_values(first + d, k));
1657 }
1658 }
1659
1660 for (unsigned int k = 0; k < n_q_points; ++k)
1661 for (unsigned int d = 0; d < dim; ++d)
1662 output_data.shape_hessians[first + d][k] =
1663 dof_sign * transformed_shape_hessians[k][d];
1664
1665 break;
1666 }
1667
1668 case mapping_nedelec:
1669 {
1670 for (unsigned int k = 0; k < n_q_points; ++k)
1671 fe_data.untransformed_shape_hessian_tensors[k + offset] =
1672 fe_data.shape_grad_grads[dof_index][k + offset];
1673
1674 const ArrayView<Tensor<3, spacedim>>
1675 transformed_shape_hessians =
1676 make_array_view(fe_data.transformed_shape_hessians,
1677 offset,
1678 n_q_points);
1679 mapping.transform(
1680 make_array_view(fe_data.untransformed_shape_hessian_tensors,
1681 offset,
1682 n_q_points),
1683 mapping_covariant_hessian,
1684 mapping_internal,
1685 transformed_shape_hessians);
1686
1687 for (unsigned int k = 0; k < n_q_points; ++k)
1688 for (unsigned int d = 0; d < spacedim; ++d)
1689 for (unsigned int n = 0; n < spacedim; ++n)
1690 for (unsigned int i = 0; i < spacedim; ++i)
1691 for (unsigned int j = 0; j < spacedim; ++j)
1692 {
1693 transformed_shape_hessians[k][d][i][j] -=
1694 (output_data.shape_values(first + n, k) *
1695 mapping_data
1696 .jacobian_pushed_forward_2nd_derivatives
1697 [k][n][d][i][j]) +
1698 (output_data.shape_gradients[first + d][k][n] *
1699 mapping_data
1700 .jacobian_pushed_forward_grads[k][n][i][j]) +
1701 (output_data.shape_gradients[first + n][k][i] *
1702 mapping_data
1703 .jacobian_pushed_forward_grads[k][n][d][j]) +
1704 (output_data.shape_gradients[first + n][k][j] *
1705 mapping_data
1706 .jacobian_pushed_forward_grads[k][n][i][d]);
1707 }
1708
1709 for (unsigned int k = 0; k < n_q_points; ++k)
1710 for (unsigned int d = 0; d < dim; ++d)
1711 output_data.shape_hessians[first + d][k] =
1712 dof_sign * transformed_shape_hessians[k][d];
1713
1714 break;
1715 }
1716
1717 default:
1718 Assert(false, ExcNotImplemented());
1719 }
1720 }
1721
1722 // third derivatives are not implemented
1723 if (fe_data.update_each & update_3rd_derivatives)
1724 {
1725 Assert(false, ExcNotImplemented())
1726 }
1727 }
1728}
1729
1730
1731
1732template <int dim, int spacedim>
1733void
1734FE_PolyTensor<dim, spacedim>::fill_fe_subface_values(
1735 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1736 const unsigned int face_no,
1737 const unsigned int sub_no,
1738 const Quadrature<dim - 1> & quadrature,
1739 const Mapping<dim, spacedim> & mapping,
1740 const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
1741 const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1742 & mapping_data,
1743 const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
1744 ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
1745 spacedim>
1746 &output_data) const
1747{
1748 // convert data object to internal
1749 // data for this class. fails with
1750 // an exception if that is not
1751 // possible
1752 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
1753 ExcInternalError());
1754 const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
1755
1756 const unsigned int n_q_points = quadrature.size();
1757
1758 // offset determines which data set
1759 // to take (all data sets for all
1760 // sub-faces are stored contiguously)
1761 const auto offset =
1762 QProjector<dim>::DataSetDescriptor::subface(this->reference_cell(),
1763 face_no,
1764 sub_no,
1765 cell->face_orientation(face_no),
1766 cell->face_flip(face_no),
1767 cell->face_rotation(face_no),
1768 n_q_points,
1769 cell->subface_case(face_no));
1770
1771 // TODO: Size assertions
1772
1773 // TODO: The dof_sign_change only affects Nedelec elements and is not the
1774 // correct thing on complicated meshes for higher order Nedelec elements.
1775 // Something similar to FE_Q should be done to permute dofs and to change the
1776 // dof signs. A static way using tables (as done in the RaviartThomas<dim>
1777 // class) is preferable.
1778 std::fill(fe_data.dof_sign_change.begin(),
1779 fe_data.dof_sign_change.end(),
1780 1.0);
1781 internal::FE_PolyTensor::get_dof_sign_change_nedelec(cell,
1782 *this,
1783 this->mapping_kind,
1784 fe_data.dof_sign_change);
1785
1786 // TODO: This, similarly to the Nedelec case, is just a legacy function in 2d
1787 // and affects only face_dofs of H(div) conformal FEs. It does nothing in 1d.
1788 // Also nothing in 3d since we take care of it by using the
1789 // adjust_quad_dof_sign_for_face_orientation_table.
1790 internal::FE_PolyTensor::get_dof_sign_change_h_div(cell,
1791 *this,
1792 this->mapping_kind,
1793 fe_data.dof_sign_change);
1794
1795 // What is the first dof_index on a quad?
1796 const unsigned int first_quad_index = this->get_first_quad_index();
1797 // How many dofs per quad and how many quad dofs do we have at all?
1798 const unsigned int n_dofs_per_quad = this->n_dofs_per_quad();
1799 const unsigned int n_quad_dofs =
1800 n_dofs_per_quad * GeometryInfo<dim>::faces_per_cell;
1801
1802 for (unsigned int dof_index = 0; dof_index < this->n_dofs_per_cell();
1803 ++dof_index)
1804 {
1805 /*
1806 * This assumes that the dofs are ordered by first vertices, lines, quads
1807 * and volume dofs. Note that in 2d this always gives false.
1808 */
1809 const bool is_quad_dof =
1810 (dim == 2 ? false :
1811 (first_quad_index <= dof_index) &&
1812 (dof_index < first_quad_index + n_quad_dofs));
1813
1814 // TODO: This hack is not pretty and it is only here to handle the 2d
1815 // case and the Nedelec legacy case. In 2d dof_sign of a face_dof is never
1816 // handled by the
1817 // >>if(is_quad_dof){...}<< but still a possible dof sign change must be
1818 // handled, also for line_dofs in 3d such as in Nedelec. In these cases
1819 // this is encoded in the array fe_data.dof_sign_change[dof_index]. In 3d
1820 // it is handles with a table. This array is allocated in
1821 // fe_poly_tensor.h.
1822 double dof_sign = 1.0;
1823 // under some circumstances fe_data.dof_sign_change is not allocated
1824 if (fe_data.update_each & update_values)
1825 dof_sign = fe_data.dof_sign_change[dof_index];
1826
1827 if (is_quad_dof)
1828 {
1829 /*
1830 * Find the face belonging to this dof_index. This is integer
1831 * division.
1832 */
1833 unsigned int face_index_from_dof_index =
1834 (dof_index - first_quad_index) / (n_dofs_per_quad);
1835
1836 unsigned int local_quad_dof_index = dof_index % n_dofs_per_quad;
1837
1838 // Correct the dof_sign if necessary
1839 if (adjust_quad_dof_sign_for_face_orientation(
1840 local_quad_dof_index,
1841 face_index_from_dof_index,
1842 cell->face_orientation(face_index_from_dof_index),
1843 cell->face_flip(face_index_from_dof_index),
1844 cell->face_rotation(face_index_from_dof_index)))
1845 dof_sign = -1.0;
1846 }
1847
1848 const MappingKind mapping_kind = get_mapping_kind(dof_index);
1849
1850 const unsigned int first =
1851 output_data.shape_function_to_row_table
1852 [dof_index * this->n_components() +
1853 this->get_nonzero_components(dof_index).first_selected_component()];
1854
1855 if (fe_data.update_each & update_values)
1856 {
1857 switch (mapping_kind)
1858 {
1859 case mapping_none:
1860 {
1861 for (unsigned int k = 0; k < n_q_points; ++k)
1862 for (unsigned int d = 0; d < dim; ++d)
1863 output_data.shape_values(first + d, k) =
1864 fe_data.shape_values[dof_index][k + offset][d];
1865 break;
1866 }
1867
1868 case mapping_covariant:
1869 case mapping_contravariant:
1870 {
1871 const ArrayView<Tensor<1, spacedim>>
1872 transformed_shape_values =
1873 make_array_view(fe_data.transformed_shape_values,
1874 offset,
1875 n_q_points);
1876 mapping.transform(make_array_view(fe_data.shape_values,
1877 dof_index,
1878 offset,
1879 n_q_points),
1880 mapping_kind,
1881 mapping_internal,
1882 transformed_shape_values);
1883
1884 for (unsigned int k = 0; k < n_q_points; ++k)
1885 for (unsigned int d = 0; d < dim; ++d)
1886 output_data.shape_values(first + d, k) =
1887 transformed_shape_values[k][d];
1888
1889 break;
1890 }
1891
1892 case mapping_raviart_thomas:
1893 case mapping_piola:
1894 {
1895 const ArrayView<Tensor<1, spacedim>>
1896 transformed_shape_values =
1897 make_array_view(fe_data.transformed_shape_values,
1898 offset,
1899 n_q_points);
1900
1901 mapping.transform(make_array_view(fe_data.shape_values,
1902 dof_index,
1903 offset,
1904 n_q_points),
1905 mapping_piola,
1906 mapping_internal,
1907 transformed_shape_values);
1908 for (unsigned int k = 0; k < n_q_points; ++k)
1909 for (unsigned int d = 0; d < dim; ++d)
1910 output_data.shape_values(first + d, k) =
1911 dof_sign * transformed_shape_values[k][d];
1912 break;
1913 }
1914
1915 case mapping_nedelec:
1916 {
1917 const ArrayView<Tensor<1, spacedim>>
1918 transformed_shape_values =
1919 make_array_view(fe_data.transformed_shape_values,
1920 offset,
1921 n_q_points);
1922
1923 mapping.transform(make_array_view(fe_data.shape_values,
1924 dof_index,
1925 offset,
1926 n_q_points),
1927 mapping_covariant,
1928 mapping_internal,
1929 transformed_shape_values);
1930
1931 for (unsigned int k = 0; k < n_q_points; ++k)
1932 for (unsigned int d = 0; d < dim; ++d)
1933 output_data.shape_values(first + d, k) =
1934 dof_sign * transformed_shape_values[k][d];
1935
1936 break;
1937 }
1938
1939 default:
1940 Assert(false, ExcNotImplemented());
1941 }
1942 }
1943
1944 if (fe_data.update_each & update_gradients)
1945 {
1946 const ArrayView<Tensor<2, spacedim>> transformed_shape_grads =
1947 make_array_view(fe_data.transformed_shape_grads,
1948 offset,
1949 n_q_points);
1950 switch (mapping_kind)
1951 {
1952 case mapping_none:
1953 {
1954 mapping.transform(make_array_view(fe_data.shape_grads,
1955 dof_index,
1956 offset,
1957 n_q_points),
1958 mapping_covariant,
1959 mapping_internal,
1960 transformed_shape_grads);
1961 for (unsigned int k = 0; k < n_q_points; ++k)
1962 for (unsigned int d = 0; d < dim; ++d)
1963 output_data.shape_gradients[first + d][k] =
1964 transformed_shape_grads[k][d];
1965 break;
1966 }
1967
1968 case mapping_covariant:
1969 {
1970 mapping.transform(make_array_view(fe_data.shape_grads,
1971 dof_index,
1972 offset,
1973 n_q_points),
1974 mapping_covariant_gradient,
1975 mapping_internal,
1976 transformed_shape_grads);
1977
1978 for (unsigned int k = 0; k < n_q_points; ++k)
1979 for (unsigned int d = 0; d < spacedim; ++d)
1980 for (unsigned int n = 0; n < spacedim; ++n)
1981 transformed_shape_grads[k][d] -=
1982 output_data.shape_values(first + n, k) *
1983 mapping_data.jacobian_pushed_forward_grads[k][n][d];
1984
1985 for (unsigned int k = 0; k < n_q_points; ++k)
1986 for (unsigned int d = 0; d < dim; ++d)
1987 output_data.shape_gradients[first + d][k] =
1988 transformed_shape_grads[k][d];
1989
1990 break;
1991 }
1992
1993 case mapping_contravariant:
1994 {
1995 for (unsigned int k = 0; k < n_q_points; ++k)
1996 fe_data.untransformed_shape_grads[k + offset] =
1997 fe_data.shape_grads[dof_index][k + offset];
1998
1999 mapping.transform(
2000 make_array_view(fe_data.untransformed_shape_grads,
2001 offset,
2002 n_q_points),
2003 mapping_contravariant_gradient,
2004 mapping_internal,
2005 transformed_shape_grads);
2006
2007 for (unsigned int k = 0; k < n_q_points; ++k)
2008 for (unsigned int d = 0; d < spacedim; ++d)
2009 for (unsigned int n = 0; n < spacedim; ++n)
2010 transformed_shape_grads[k][d] +=
2011 output_data.shape_values(first + n, k) *
2012 mapping_data.jacobian_pushed_forward_grads[k][d][n];
2013
2014 for (unsigned int k = 0; k < n_q_points; ++k)
2015 for (unsigned int d = 0; d < dim; ++d)
2016 output_data.shape_gradients[first + d][k] =
2017 transformed_shape_grads[k][d];
2018
2019 break;
2020 }
2021
2022 case mapping_raviart_thomas:
2023 case mapping_piola:
2024 {
2025 for (unsigned int k = 0; k < n_q_points; ++k)
2026 fe_data.untransformed_shape_grads[k + offset] =
2027 fe_data.shape_grads[dof_index][k + offset];
2028
2029 mapping.transform(
2030 make_array_view(fe_data.untransformed_shape_grads,
2031 offset,
2032 n_q_points),
2033 mapping_piola_gradient,
2034 mapping_internal,
2035 transformed_shape_grads);
2036
2037 for (unsigned int k = 0; k < n_q_points; ++k)
2038 for (unsigned int d = 0; d < spacedim; ++d)
2039 for (unsigned int n = 0; n < spacedim; ++n)
2040 transformed_shape_grads[k][d] +=
2041 (output_data.shape_values(first + n, k) *
2042 mapping_data
2043 .jacobian_pushed_forward_grads[k][d][n]) -
2044 (output_data.shape_values(first + d, k) *
2045 mapping_data.jacobian_pushed_forward_grads[k][n][n]);
2046
2047 for (unsigned int k = 0; k < n_q_points; ++k)
2048 for (unsigned int d = 0; d < dim; ++d)
2049 output_data.shape_gradients[first + d][k] =
2050 dof_sign * transformed_shape_grads[k][d];
2051
2052 break;
2053 }
2054
2055 case mapping_nedelec:
2056 {
2057 // this particular shape
2058 // function at all
2059 // q-points. if Dv is the
2060 // gradient of the shape
2061 // function on the unit
2062 // cell, then
2063 // (J^-T)Dv(J^-1) is the
2064 // value we want to have on
2065 // the real cell.
2066 for (unsigned int k = 0; k < n_q_points; ++k)
2067 fe_data.untransformed_shape_grads[k + offset] =
2068 fe_data.shape_grads[dof_index][k + offset];
2069
2070 mapping.transform(
2071 make_array_view(fe_data.untransformed_shape_grads,
2072 offset,
2073 n_q_points),
2074 mapping_covariant_gradient,
2075 mapping_internal,
2076 transformed_shape_grads);
2077
2078 for (unsigned int k = 0; k < n_q_points; ++k)
2079 for (unsigned int d = 0; d < spacedim; ++d)
2080 for (unsigned int n = 0; n < spacedim; ++n)
2081 transformed_shape_grads[k][d] -=
2082 output_data.shape_values(first + n, k) *
2083 mapping_data.jacobian_pushed_forward_grads[k][n][d];
2084
2085 for (unsigned int k = 0; k < n_q_points; ++k)
2086 for (unsigned int d = 0; d < dim; ++d)
2087 output_data.shape_gradients[first + d][k] =
2088 dof_sign * transformed_shape_grads[k][d];
2089
2090 break;
2091 }
2092
2093 default:
2094 Assert(false, ExcNotImplemented());
2095 }
2096 }
2097
2098 if (fe_data.update_each & update_hessians)
2099 {
2100 switch (mapping_kind)
2101 {
2102 case mapping_none:
2103 {
2104 const ArrayView<Tensor<3, spacedim>>
2105 transformed_shape_hessians =
2106 make_array_view(fe_data.transformed_shape_hessians,
2107 offset,
2108 n_q_points);
2109 mapping.transform(make_array_view(fe_data.shape_grad_grads,
2110 dof_index,
2111 offset,
2112 n_q_points),
2113 mapping_covariant_gradient,
2114 mapping_internal,
2115 transformed_shape_hessians);
2116
2117 for (unsigned int k = 0; k < n_q_points; ++k)
2118 for (unsigned int d = 0; d < spacedim; ++d)
2119 for (unsigned int n = 0; n < spacedim; ++n)
2120 transformed_shape_hessians[k][d] -=
2121 output_data.shape_gradients[first + d][k][n] *
2122 mapping_data.jacobian_pushed_forward_grads[k][n];
2123
2124 for (unsigned int k = 0; k < n_q_points; ++k)
2125 for (unsigned int d = 0; d < dim; ++d)
2126 output_data.shape_hessians[first + d][k] =
2127 transformed_shape_hessians[k][d];
2128
2129 break;
2130 }
2131 case mapping_covariant:
2132 {
2133 for (unsigned int k = 0; k < n_q_points; ++k)
2134 fe_data.untransformed_shape_hessian_tensors[k + offset] =
2135 fe_data.shape_grad_grads[dof_index][k + offset];
2136
2137 const ArrayView<Tensor<3, spacedim>>
2138 transformed_shape_hessians =
2139 make_array_view(fe_data.transformed_shape_hessians,
2140 offset,
2141 n_q_points);
2142 mapping.transform(
2143 make_array_view(fe_data.untransformed_shape_hessian_tensors,
2144 offset,
2145 n_q_points),
2146 mapping_covariant_hessian,
2147 mapping_internal,
2148 transformed_shape_hessians);
2149
2150 for (unsigned int k = 0; k < n_q_points; ++k)
2151 for (unsigned int d = 0; d < spacedim; ++d)
2152 for (unsigned int n = 0; n < spacedim; ++n)
2153 for (unsigned int i = 0; i < spacedim; ++i)
2154 for (unsigned int j = 0; j < spacedim; ++j)
2155 {
2156 transformed_shape_hessians[k][d][i][j] -=
2157 (output_data.shape_values(first + n, k) *
2158 mapping_data
2159 .jacobian_pushed_forward_2nd_derivatives
2160 [k][n][d][i][j]) +
2161 (output_data.shape_gradients[first + d][k][n] *
2162 mapping_data
2163 .jacobian_pushed_forward_grads[k][n][i][j]) +
2164 (output_data.shape_gradients[first + n][k][i] *
2165 mapping_data
2166 .jacobian_pushed_forward_grads[k][n][d][j]) +
2167 (output_data.shape_gradients[first + n][k][j] *
2168 mapping_data
2169 .jacobian_pushed_forward_grads[k][n][i][d]);
2170 }
2171
2172 for (unsigned int k = 0; k < n_q_points; ++k)
2173 for (unsigned int d = 0; d < dim; ++d)
2174 output_data.shape_hessians[first + d][k] =
2175 transformed_shape_hessians[k][d];
2176
2177 break;
2178 }
2179
2180 case mapping_contravariant:
2181 {
2182 for (unsigned int k = 0; k < n_q_points; ++k)
2183 fe_data.untransformed_shape_hessian_tensors[k + offset] =
2184 fe_data.shape_grad_grads[dof_index][k + offset];
2185
2186 const ArrayView<Tensor<3, spacedim>>
2187 transformed_shape_hessians =
2188 make_array_view(fe_data.transformed_shape_hessians,
2189 offset,
2190 n_q_points);
2191 mapping.transform(
2192 make_array_view(fe_data.untransformed_shape_hessian_tensors,
2193 offset,
2194 n_q_points),
2195 mapping_contravariant_hessian,
2196 mapping_internal,
2197 transformed_shape_hessians);
2198
2199 for (unsigned int k = 0; k < n_q_points; ++k)
2200 for (unsigned int d = 0; d < spacedim; ++d)
2201 for (unsigned int n = 0; n < spacedim; ++n)
2202 for (unsigned int i = 0; i < spacedim; ++i)
2203 for (unsigned int j = 0; j < spacedim; ++j)
2204 {
2205 transformed_shape_hessians[k][d][i][j] +=
2206 (output_data.shape_values(first + n, k) *
2207 mapping_data
2208 .jacobian_pushed_forward_2nd_derivatives
2209 [k][d][n][i][j]) +
2210 (output_data.shape_gradients[first + n][k][i] *
2211 mapping_data
2212 .jacobian_pushed_forward_grads[k][d][n][j]) +
2213 (output_data.shape_gradients[first + n][k][j] *
2214 mapping_data
2215 .jacobian_pushed_forward_grads[k][d][i][n]) -
2216 (output_data.shape_gradients[first + d][k][n] *
2217 mapping_data
2218 .jacobian_pushed_forward_grads[k][n][i][j]);
2219 for (unsigned int m = 0; m < spacedim; ++m)
2220 transformed_shape_hessians[k][d][i][j] -=
2221 (mapping_data
2222 .jacobian_pushed_forward_grads[k][d][i]
2223 [m] *
2224 mapping_data
2225 .jacobian_pushed_forward_grads[k][m][n]
2226 [j] *
2227 output_data.shape_values(first + n, k)) +
2228 (mapping_data
2229 .jacobian_pushed_forward_grads[k][d][m]
2230 [j] *
2231 mapping_data
2232 .jacobian_pushed_forward_grads[k][m][i]
2233 [n] *
2234 output_data.shape_values(first + n, k));
2235 }
2236
2237 for (unsigned int k = 0; k < n_q_points; ++k)
2238 for (unsigned int d = 0; d < dim; ++d)
2239 output_data.shape_hessians[first + d][k] =
2240 transformed_shape_hessians[k][d];
2241
2242 break;
2243 }
2244
2245 case mapping_raviart_thomas:
2246 case mapping_piola:
2247 {
2248 for (unsigned int k = 0; k < n_q_points; ++k)
2249 fe_data.untransformed_shape_hessian_tensors[k + offset] =
2250 fe_data.shape_grad_grads[dof_index][k + offset];
2251
2252 const ArrayView<Tensor<3, spacedim>>
2253 transformed_shape_hessians =
2254 make_array_view(fe_data.transformed_shape_hessians,
2255 offset,
2256 n_q_points);
2257 mapping.transform(
2258 make_array_view(fe_data.untransformed_shape_hessian_tensors,
2259 offset,
2260 n_q_points),
2261 mapping_piola_hessian,
2262 mapping_internal,
2263 transformed_shape_hessians);
2264
2265 for (unsigned int k = 0; k < n_q_points; ++k)
2266 for (unsigned int d = 0; d < spacedim; ++d)
2267 for (unsigned int n = 0; n < spacedim; ++n)
2268 for (unsigned int i = 0; i < spacedim; ++i)
2269 for (unsigned int j = 0; j < spacedim; ++j)
2270 {
2271 transformed_shape_hessians[k][d][i][j] +=
2272 (output_data.shape_values(first + n, k) *
2273 mapping_data
2274 .jacobian_pushed_forward_2nd_derivatives
2275 [k][d][n][i][j]) +
2276 (output_data.shape_gradients[first + n][k][i] *
2277 mapping_data
2278 .jacobian_pushed_forward_grads[k][d][n][j]) +
2279 (output_data.shape_gradients[first + n][k][j] *
2280 mapping_data
2281 .jacobian_pushed_forward_grads[k][d][i][n]) -
2282 (output_data.shape_gradients[first + d][k][n] *
2283 mapping_data
2284 .jacobian_pushed_forward_grads[k][n][i][j]);
2285
2286 transformed_shape_hessians[k][d][i][j] -=
2287 (output_data.shape_values(first + d, k) *
2288 mapping_data
2289 .jacobian_pushed_forward_2nd_derivatives
2290 [k][n][n][i][j]) +
2291 (output_data.shape_gradients[first + d][k][i] *
2292 mapping_data
2293 .jacobian_pushed_forward_grads[k][n][n][j]) +
2294 (output_data.shape_gradients[first + d][k][j] *
2295 mapping_data
2296 .jacobian_pushed_forward_grads[k][n][n][i]);
2297 for (unsigned int m = 0; m < spacedim; ++m)
2298 {
2299 transformed_shape_hessians[k][d][i][j] -=
2300 (mapping_data
2301 .jacobian_pushed_forward_grads[k][d][i]
2302 [m] *
2303 mapping_data
2304 .jacobian_pushed_forward_grads[k][m][n]
2305 [j] *
2306 output_data.shape_values(first + n, k)) +
2307 (mapping_data
2308 .jacobian_pushed_forward_grads[k][d][m]
2309 [j] *
2310 mapping_data
2311 .jacobian_pushed_forward_grads[k][m][i]
2312 [n] *
2313 output_data.shape_values(first + n, k));
2314
2315 transformed_shape_hessians[k][d][i][j] +=
2316 (mapping_data
2317 .jacobian_pushed_forward_grads[k][n][i]
2318 [m] *
2319 mapping_data
2320 .jacobian_pushed_forward_grads[k][m][n]
2321 [j] *
2322 output_data.shape_values(first + d, k)) +
2323 (mapping_data
2324 .jacobian_pushed_forward_grads[k][n][m]
2325 [j] *
2326 mapping_data
2327 .jacobian_pushed_forward_grads[k][m][i]
2328 [n] *
2329 output_data.shape_values(first + d, k));
2330 }
2331 }
2332
2333 for (unsigned int k = 0; k < n_q_points; ++k)
2334 for (unsigned int d = 0; d < dim; ++d)
2335 output_data.shape_hessians[first + d][k] =
2336 dof_sign * transformed_shape_hessians[k][d];
2337
2338 break;
2339 }
2340
2341 case mapping_nedelec:
2342 {
2343 for (unsigned int k = 0; k < n_q_points; ++k)
2344 fe_data.untransformed_shape_hessian_tensors[k + offset] =
2345 fe_data.shape_grad_grads[dof_index][k + offset];
2346
2347 const ArrayView<Tensor<3, spacedim>>
2348 transformed_shape_hessians =
2349 make_array_view(fe_data.transformed_shape_hessians,
2350 offset,
2351 n_q_points);
2352 mapping.transform(
2353 make_array_view(fe_data.untransformed_shape_hessian_tensors,
2354 offset,
2355 n_q_points),
2356 mapping_covariant_hessian,
2357 mapping_internal,
2358 transformed_shape_hessians);
2359
2360 for (unsigned int k = 0; k < n_q_points; ++k)
2361 for (unsigned int d = 0; d < spacedim; ++d)
2362 for (unsigned int n = 0; n < spacedim; ++n)
2363 for (unsigned int i = 0; i < spacedim; ++i)
2364 for (unsigned int j = 0; j < spacedim; ++j)
2365 {
2366 transformed_shape_hessians[k][d][i][j] -=
2367 (output_data.shape_values(first + n, k) *
2368 mapping_data
2369 .jacobian_pushed_forward_2nd_derivatives
2370 [k][n][d][i][j]) +
2371 (output_data.shape_gradients[first + d][k][n] *
2372 mapping_data
2373 .jacobian_pushed_forward_grads[k][n][i][j]) +
2374 (output_data.shape_gradients[first + n][k][i] *
2375 mapping_data
2376 .jacobian_pushed_forward_grads[k][n][d][j]) +
2377 (output_data.shape_gradients[first + n][k][j] *
2378 mapping_data
2379 .jacobian_pushed_forward_grads[k][n][i][d]);
2380 }
2381
2382 for (unsigned int k = 0; k < n_q_points; ++k)
2383 for (unsigned int d = 0; d < dim; ++d)
2384 output_data.shape_hessians[first + d][k] =
2385 dof_sign * transformed_shape_hessians[k][d];
2386
2387 break;
2388 }
2389
2390 default:
2391 Assert(false, ExcNotImplemented());
2392 }
2393 }
2394
2395 // third derivatives are not implemented
2396 if (fe_data.update_each & update_3rd_derivatives)
2397 {
2398 Assert(false, ExcNotImplemented())
2399 }
2400 }
2401}
2402
2403
2404
2405template <int dim, int spacedim>
2406UpdateFlags
2407FE_PolyTensor<dim, spacedim>::requires_update_flags(
2408 const UpdateFlags flags) const
2409{
2410 UpdateFlags out = update_default;
2411
2412 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
2413 {
2414 const MappingKind mapping_kind = get_mapping_kind(i);
2415
2416 switch (mapping_kind)
2417 {
2418 case mapping_none:
2419 {
2420 if (flags & update_values)
2421 out |= update_values;
2422
2423 if (flags & update_gradients)
2424 out |= update_gradients | update_values |
2425 update_jacobian_pushed_forward_grads;
2426
2427 if (flags & update_hessians)
2428 out |= update_hessians | update_values | update_gradients |
2429 update_jacobian_pushed_forward_grads |
2430 update_jacobian_pushed_forward_2nd_derivatives;
2431 break;
2432 }
2433 case mapping_raviart_thomas:
2434 case mapping_piola:
2435 {
2436 if (flags & update_values)
2437 out |= update_values | update_piola;
2438
2439 if (flags & update_gradients)
2440 out |= update_gradients | update_values | update_piola |
2441 update_jacobian_pushed_forward_grads |
2442 update_covariant_transformation |
2443 update_contravariant_transformation;
2444
2445 if (flags & update_hessians)
2446 out |= update_hessians | update_piola | update_values |
2447 update_gradients | update_jacobian_pushed_forward_grads |
2448 update_jacobian_pushed_forward_2nd_derivatives |
2449 update_covariant_transformation;
2450
2451 break;
2452 }
2453
2454
2455 case mapping_contravariant:
2456 {
2457 if (flags & update_values)
2458 out |= update_values | update_piola;
2459
2460 if (flags & update_gradients)
2461 out |= update_gradients | update_values |
2462 update_jacobian_pushed_forward_grads |
2463 update_covariant_transformation |
2464 update_contravariant_transformation;
2465
2466 if (flags & update_hessians)
2467 out |= update_hessians | update_piola | update_values |
2468 update_gradients | update_jacobian_pushed_forward_grads |
2469 update_jacobian_pushed_forward_2nd_derivatives |
2470 update_covariant_transformation;
2471
2472 break;
2473 }
2474
2475 case mapping_nedelec:
2476 case mapping_covariant:
2477 {
2478 if (flags & update_values)
2479 out |= update_values | update_covariant_transformation;
2480
2481 if (flags & update_gradients)
2482 out |= update_gradients | update_values |
2483 update_jacobian_pushed_forward_grads |
2484 update_covariant_transformation;
2485
2486 if (flags & update_hessians)
2487 out |= update_hessians | update_values | update_gradients |
2488 update_jacobian_pushed_forward_grads |
2489 update_jacobian_pushed_forward_2nd_derivatives |
2490 update_covariant_transformation;
2491
2492 break;
2493 }
2494
2495 default:
2496 {
2497 Assert(false, ExcNotImplemented());
2498 }
2499 }
2500 }
2501
2502 return out;
2503}
2504
2505
2506// explicit instantiations
2507#include "fe_poly_tensor.inst"
2508
2509
2510DEAL_II_NAMESPACE_CLOSE
array_view.h
make_array_view
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:704
cell_id.h
FE_PolyTensor
Definition fe_poly_tensor.h:142
FE_PolyTensor::requires_update_flags
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
FE_PolyTensor::shape_grad
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const override
FE_PolyTensor::shape_grad_grad_component
virtual Tensor< 2, dim > shape_grad_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
FE_PolyTensor::shape_value_component
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
FE_PolyTensor::shape_grad_grad
virtual Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const override
FE_PolyTensor::fill_fe_face_values
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
FE_PolyTensor::shape_value
virtual double shape_value(const unsigned int i, const Point< dim > &p) const override
FE_PolyTensor::fill_fe_values
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
FE_PolyTensor::adjust_quad_dof_sign_for_face_orientation
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
FE_PolyTensor::fill_fe_subface_values
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
FE_PolyTensor::get_mapping_kind
MappingKind get_mapping_kind(const unsigned int i) const
FE_PolyTensor::FE_PolyTensor
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)
FE_PolyTensor::single_mapping_kind
bool single_mapping_kind() const
FE_PolyTensor::shape_grad_component
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
FiniteElementData::n_dofs_per_face
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
FiniteElement::face_to_cell_index
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false) const
Mapping
Abstract base class for mapping classes.
Definition mapping.h:317
Mapping::transform
virtual void transform(const ArrayView< const Tensor< 1, dim > > &input, const MappingKind kind, const typename Mapping< dim, spacedim >::InternalDataBase &internal, const ArrayView< Tensor< 1, spacedim > > &output) const =0
Point
Definition point.h:112
QProjector::DataSetDescriptor::subface
static DataSetDescriptor subface(const ReferenceCell &reference_cell, const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points, const internal::SubfaceCase< dim > ref_case=internal::SubfaceCase< dim >::case_isotropic)
QProjector::DataSetDescriptor::face
static DataSetDescriptor face(const ReferenceCell &reference_cell, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
Definition qprojector.cc:1318
Quadrature::size
unsigned int size() const
hp::Collection::size
unsigned int size() const
Definition collection.h:265
internal::FEValuesImplementation::MappingRelatedData::jacobian_pushed_forward_2nd_derivatives
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
Definition mapping_related_data.h:116
internal::FEValuesImplementation::MappingRelatedData::jacobian_pushed_forward_grads
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
Definition mapping_related_data.h:104
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
derivative_form.h
fe_values.h
fe_poly_tensor.h
first
Point< 2 > first
Definition grid_out.cc:4615
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1614
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:1787
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:1855
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
UpdateFlags
UpdateFlags
Definition fe_update_flags.h:65
update_jacobian_pushed_forward_2nd_derivatives
@ update_jacobian_pushed_forward_2nd_derivatives
Definition fe_update_flags.h:204
update_contravariant_transformation
@ update_contravariant_transformation
Contravariant transformation.
Definition fe_update_flags.h:173
update_jacobian_pushed_forward_grads
@ update_jacobian_pushed_forward_grads
Definition fe_update_flags.h:195
update_hessians
@ update_hessians
Second derivatives of shape functions.
Definition fe_update_flags.h:88
update_values
@ update_values
Shape function values.
Definition fe_update_flags.h:76
update_3rd_derivatives
@ update_3rd_derivatives
Third derivatives of shape functions.
Definition fe_update_flags.h:94
update_covariant_transformation
@ update_covariant_transformation
Covariant transformation.
Definition fe_update_flags.h:166
update_gradients
@ update_gradients
Shape function gradients.
Definition fe_update_flags.h:82
update_default
@ update_default
No update.
Definition fe_update_flags.h:67
update_piola
@ update_piola
Values needed for Piola transform.
Definition fe_update_flags.h:218
MappingKind
MappingKind
Definition mapping.h:78
mapping_piola
@ mapping_piola
Definition mapping.h:113
mapping_covariant_gradient
@ mapping_covariant_gradient
Definition mapping.h:99
mapping_covariant
@ mapping_covariant
Definition mapping.h:88
mapping_nedelec
@ mapping_nedelec
Definition mapping.h:128
mapping_contravariant
@ mapping_contravariant
Definition mapping.h:93
mapping_raviart_thomas
@ mapping_raviart_thomas
Definition mapping.h:133
mapping_contravariant_hessian
@ mapping_contravariant_hessian
Definition mapping.h:155
mapping_covariant_hessian
@ mapping_covariant_hessian
Definition mapping.h:149
mapping_contravariant_gradient
@ mapping_contravariant_gradient
Definition mapping.h:105
mapping_none
@ mapping_none
Definition mapping.h:83
mapping_piola_gradient
@ mapping_piola_gradient
Definition mapping.h:119
mapping_piola_hessian
@ mapping_piola_hessian
Definition mapping.h:161
mapping_cartesian.h
GridGenerator::reference_cell
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
polynomials_abf.h
polynomials_bdm.h
polynomials_bernardi_raugel.h
polynomials_nedelec.h
polynomials_raviart_thomas.h
polynomials_rt_bubbles.h
Generated by doxygen 1.9.7