deal.II version GIT relicensing-2430-g7f5db50250 2025-01-24 01: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
mapping_fe_field.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2015 - 2024 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#include <deal.II/base/array_view.h>
16#include <deal.II/base/memory_consumption.h>
17#include <deal.II/base/polynomial.h>
18#include <deal.II/base/qprojector.h>
19#include <deal.II/base/quadrature.h>
20#include <deal.II/base/quadrature_lib.h>
21#include <deal.II/base/tensor_product_polynomials.h>
22#include <deal.II/base/utilities.h>
23
24#include <deal.II/dofs/dof_accessor.h>
25
26#include <deal.II/fe/fe_q.h>
27#include <deal.II/fe/fe_system.h>
28#include <deal.II/fe/fe_tools.h>
29#include <deal.II/fe/fe_values.h>
30#include <deal.II/fe/mapping.h>
31#include <deal.II/fe/mapping_fe_field.h>
32
33#include <deal.II/grid/tria_iterator.h>
34
35#include <deal.II/lac/block_vector.h>
36#include <deal.II/lac/full_matrix.h>
37#include <deal.II/lac/la_parallel_block_vector.h>
38#include <deal.II/lac/la_parallel_vector.h>
39#include <deal.II/lac/petsc_block_vector.h>
40#include <deal.II/lac/petsc_vector.h>
41#include <deal.II/lac/trilinos_epetra_vector.h>
42#include <deal.II/lac/trilinos_parallel_block_vector.h>
43#include <deal.II/lac/trilinos_tpetra_vector.h>
44#include <deal.II/lac/trilinos_vector.h>
45#include <deal.II/lac/vector.h>
46
47#include <deal.II/numerics/vector_tools.h>
48
49#include <fstream>
50#include <memory>
51#include <numeric>
52
53
54
55DEAL_II_NAMESPACE_OPEN
56
57
58template <int dim, int spacedim, typename VectorType>
59MappingFEField<dim, spacedim, VectorType>::InternalData::InternalData(
60 const FiniteElement<dim, spacedim> &fe,
61 const ComponentMask &mask)
62 : fe(&fe)
63 , unit_tangentials()
64 , n_shape_functions(fe.n_dofs_per_cell())
65 , mask(mask)
66 , local_dof_indices(fe.n_dofs_per_cell())
67 , local_dof_values(fe.n_dofs_per_cell())
68{}
69
70
71
72template <int dim, int spacedim, typename VectorType>
73void
74MappingFEField<dim, spacedim, VectorType>::InternalData::reinit(
75 const UpdateFlags update_flags,
76 const Quadrature<dim> &quadrature)
77{
78 // store the flags in the internal data object so we can access them
79 // in fill_fe_*_values(). use the transitive hull of the required
80 // flags
81 this->update_each = update_flags;
82
83 const unsigned int n_q_points = quadrature.size();
84 const std::vector<Point<dim>> &points = quadrature.get_points();
85
86 // see if we need the (transformation) shape function values
87 // and/or gradients and resize the necessary arrays
88 if (update_flags & update_quadrature_points)
89 {
90 shape_values.resize(n_shape_functions * n_q_points);
91 for (unsigned int point = 0; point < n_q_points; ++point)
92 for (unsigned int i = 0; i < n_shape_functions; ++i)
93 shape(point, i) = fe->shape_value(i, points[point]);
94 }
95
96 if (update_flags &
97 (update_covariant_transformation | update_contravariant_transformation |
98 update_JxW_values | update_boundary_forms | update_normal_vectors |
99 update_jacobians | update_jacobian_grads | update_inverse_jacobians))
100 {
101 shape_derivatives.resize(n_shape_functions * n_q_points);
102 for (unsigned int point = 0; point < n_q_points; ++point)
103 for (unsigned int i = 0; i < n_shape_functions; ++i)
104 derivative(point, i) = fe->shape_grad(i, points[point]);
105 }
106
107 if (update_flags & update_covariant_transformation)
108 covariant.resize(n_q_points);
109
110 if (update_flags & update_contravariant_transformation)
111 contravariant.resize(n_q_points);
112
113 if (update_flags & update_volume_elements)
114 volume_elements.resize(n_q_points);
115
116 if (update_flags &
117 (update_jacobian_grads | update_jacobian_pushed_forward_grads))
118 {
119 shape_second_derivatives.resize(n_shape_functions * n_q_points);
120 for (unsigned int point = 0; point < n_q_points; ++point)
121 for (unsigned int i = 0; i < n_shape_functions; ++i)
122 second_derivative(point, i) = fe->shape_grad_grad(i, points[point]);
123 }
124
125 if (update_flags & (update_jacobian_2nd_derivatives |
126 update_jacobian_pushed_forward_2nd_derivatives))
127 {
128 shape_third_derivatives.resize(n_shape_functions * n_q_points);
129 for (unsigned int point = 0; point < n_q_points; ++point)
130 for (unsigned int i = 0; i < n_shape_functions; ++i)
131 third_derivative(point, i) =
132 fe->shape_3rd_derivative(i, points[point]);
133 }
134
135 if (update_flags & (update_jacobian_3rd_derivatives |
136 update_jacobian_pushed_forward_3rd_derivatives))
137 {
138 shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
139 for (unsigned int point = 0; point < n_q_points; ++point)
140 for (unsigned int i = 0; i < n_shape_functions; ++i)
141 fourth_derivative(point, i) =
142 fe->shape_4th_derivative(i, points[point]);
143 }
144
145 // This (for face values and simplices) can be different for different
146 // calls, so always copy
147 quadrature_weights = quadrature.get_weights();
148}
149
150
151
152template <int dim, int spacedim, typename VectorType>
153std::size_t
160
161
162
163template <int dim, int spacedim, typename VectorType>
164double &
165MappingFEField<dim, spacedim, VectorType>::InternalData::shape(
166 const unsigned int qpoint,
167 const unsigned int shape_nr)
168{
169 AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
170 return shape_values[qpoint * n_shape_functions + shape_nr];
171}
172
173
174template <int dim, int spacedim, typename VectorType>
175const Tensor<1, dim> &
176MappingFEField<dim, spacedim, VectorType>::InternalData::derivative(
177 const unsigned int qpoint,
178 const unsigned int shape_nr) const
179{
180 AssertIndexRange(qpoint * n_shape_functions + shape_nr,
181 shape_derivatives.size());
182 return shape_derivatives[qpoint * n_shape_functions + shape_nr];
183}
184
185
186
187template <int dim, int spacedim, typename VectorType>
188Tensor<1, dim> &
189MappingFEField<dim, spacedim, VectorType>::InternalData::derivative(
190 const unsigned int qpoint,
191 const unsigned int shape_nr)
192{
193 AssertIndexRange(qpoint * n_shape_functions + shape_nr,
194 shape_derivatives.size());
195 return shape_derivatives[qpoint * n_shape_functions + shape_nr];
196}
197
198
199template <int dim, int spacedim, typename VectorType>
200const Tensor<2, dim> &
201MappingFEField<dim, spacedim, VectorType>::InternalData::second_derivative(
202 const unsigned int qpoint,
203 const unsigned int shape_nr) const
204{
205 AssertIndexRange(qpoint * n_shape_functions + shape_nr,
206 shape_second_derivatives.size());
207 return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
208}
209
210
211
212template <int dim, int spacedim, typename VectorType>
213Tensor<2, dim> &
214MappingFEField<dim, spacedim, VectorType>::InternalData::second_derivative(
215 const unsigned int qpoint,
216 const unsigned int shape_nr)
217{
218 AssertIndexRange(qpoint * n_shape_functions + shape_nr,
219 shape_second_derivatives.size());
220 return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
221}
222
223
224template <int dim, int spacedim, typename VectorType>
225const Tensor<3, dim> &
226MappingFEField<dim, spacedim, VectorType>::InternalData::third_derivative(
227 const unsigned int qpoint,
228 const unsigned int shape_nr) const
229{
230 AssertIndexRange(qpoint * n_shape_functions + shape_nr,
231 shape_third_derivatives.size());
232 return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
233}
234
235
236
237template <int dim, int spacedim, typename VectorType>
238Tensor<3, dim> &
239MappingFEField<dim, spacedim, VectorType>::InternalData::third_derivative(
240 const unsigned int qpoint,
241 const unsigned int shape_nr)
242{
243 AssertIndexRange(qpoint * n_shape_functions + shape_nr,
244 shape_third_derivatives.size());
245 return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
246}
247
248
249template <int dim, int spacedim, typename VectorType>
250const Tensor<4, dim> &
251MappingFEField<dim, spacedim, VectorType>::InternalData::fourth_derivative(
252 const unsigned int qpoint,
253 const unsigned int shape_nr) const
254{
255 AssertIndexRange(qpoint * n_shape_functions + shape_nr,
256 shape_fourth_derivatives.size());
257 return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
258}
259
260
261
262template <int dim, int spacedim, typename VectorType>
263Tensor<4, dim> &
264MappingFEField<dim, spacedim, VectorType>::InternalData::fourth_derivative(
265 const unsigned int qpoint,
266 const unsigned int shape_nr)
267{
268 AssertIndexRange(qpoint * n_shape_functions + shape_nr,
269 shape_fourth_derivatives.size());
270 return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
271}
272
273
274
275template <int dim, int spacedim, typename VectorType>
276MappingFEField<dim, spacedim, VectorType>::MappingFEField(
277 const DoFHandler<dim, spacedim> &euler_dof_handler,
278 const VectorType &euler_vector,
279 const ComponentMask &mask)
280 : reference_cell(euler_dof_handler.get_fe().reference_cell())
281 , uses_level_dofs(false)
282 , euler_vector({&euler_vector})
283 , euler_dof_handler(&euler_dof_handler)
284 , fe_mask(mask.size() != 0u ?
285 mask :
286 ComponentMask(
287 euler_dof_handler.get_fe().get_nonzero_components(0).size(),
288 true))
289 , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
290 , fe_values(this->euler_dof_handler->get_fe(),
291 reference_cell.template get_nodal_type_quadrature<dim>(),
292 update_values)
293{
294 AssertDimension(euler_dof_handler.n_dofs(), euler_vector.size());
295 unsigned int size = 0;
296 for (unsigned int i = 0; i < fe_mask.size(); ++i)
297 {
298 if (fe_mask[i])
299 fe_to_real[i] = size++;
300 }
301 AssertDimension(size, spacedim);
302}
303
304
305
306template <int dim, int spacedim, typename VectorType>
307MappingFEField<dim, spacedim, VectorType>::MappingFEField(
308 const DoFHandler<dim, spacedim> &euler_dof_handler,
309 const std::vector<VectorType> &euler_vector,
310 const ComponentMask &mask)
311 : reference_cell(euler_dof_handler.get_fe().reference_cell())
312 , uses_level_dofs(true)
313 , euler_dof_handler(&euler_dof_handler)
314 , fe_mask(mask.size() != 0u ?
315 mask :
316 ComponentMask(
317 euler_dof_handler.get_fe().get_nonzero_components(0).size(),
318 true))
319 , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
320 , fe_values(this->euler_dof_handler->get_fe(),
321 reference_cell.template get_nodal_type_quadrature<dim>(),
322 update_values)
323{
324 unsigned int size = 0;
325 for (unsigned int i = 0; i < fe_mask.size(); ++i)
326 {
327 if (fe_mask[i])
328 fe_to_real[i] = size++;
329 }
330 AssertDimension(size, spacedim);
331
332 Assert(euler_dof_handler.has_level_dofs(),
333 ExcMessage("The underlying DoFHandler object did not call "
334 "distribute_mg_dofs(). In this case, the construction via "
335 "level vectors does not make sense."));
336 AssertDimension(euler_vector.size(),
337 euler_dof_handler.get_triangulation().n_global_levels());
338 this->euler_vector.clear();
339 this->euler_vector.resize(euler_vector.size());
340 for (unsigned int i = 0; i < euler_vector.size(); ++i)
341 {
342 AssertDimension(euler_dof_handler.n_dofs(i), euler_vector[i].size());
343 this->euler_vector[i] = &euler_vector[i];
344 }
345}
346
347
348
349template <int dim, int spacedim, typename VectorType>
350MappingFEField<dim, spacedim, VectorType>::MappingFEField(
351 const DoFHandler<dim, spacedim> &euler_dof_handler,
352 const MGLevelObject<VectorType> &euler_vector,
353 const ComponentMask &mask)
354 : reference_cell(euler_dof_handler.get_fe().reference_cell())
355 , uses_level_dofs(true)
356 , euler_dof_handler(&euler_dof_handler)
357 , fe_mask(mask.size() != 0u ?
358 mask :
359 ComponentMask(
360 euler_dof_handler.get_fe().get_nonzero_components(0).size(),
361 true))
362 , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
363 , fe_values(this->euler_dof_handler->get_fe(),
364 reference_cell.template get_nodal_type_quadrature<dim>(),
365 update_values)
366{
367 unsigned int size = 0;
368 for (unsigned int i = 0; i < fe_mask.size(); ++i)
369 {
370 if (fe_mask[i])
371 fe_to_real[i] = size++;
372 }
373 AssertDimension(size, spacedim);
374
375 Assert(euler_dof_handler.has_level_dofs(),
376 ExcMessage("The underlying DoFHandler object did not call "
377 "distribute_mg_dofs(). In this case, the construction via "
378 "level vectors does not make sense."));
379 AssertDimension(euler_vector.max_level() + 1,
380 euler_dof_handler.get_triangulation().n_global_levels());
381 this->euler_vector.clear();
382 this->euler_vector.resize(
383 euler_dof_handler.get_triangulation().n_global_levels());
384 for (unsigned int i = euler_vector.min_level(); i <= euler_vector.max_level();
385 ++i)
386 {
387 AssertDimension(euler_dof_handler.n_dofs(i), euler_vector[i].size());
388 this->euler_vector[i] = &euler_vector[i];
389 }
390}
391
392
393
394template <int dim, int spacedim, typename VectorType>
395MappingFEField<dim, spacedim, VectorType>::MappingFEField(
396 const MappingFEField<dim, spacedim, VectorType> &mapping)
397 : reference_cell(mapping.reference_cell)
398 , uses_level_dofs(mapping.uses_level_dofs)
399 , euler_vector(mapping.euler_vector)
400 , euler_dof_handler(mapping.euler_dof_handler)
401 , fe_mask(mapping.fe_mask)
402 , fe_to_real(mapping.fe_to_real)
403 , fe_values(mapping.euler_dof_handler->get_fe(),
404 reference_cell.template get_nodal_type_quadrature<dim>(),
405 update_values)
406{}
407
408
409
410template <int dim, int spacedim, typename VectorType>
411inline const double &
412MappingFEField<dim, spacedim, VectorType>::InternalData::shape(
413 const unsigned int qpoint,
414 const unsigned int shape_nr) const
415{
416 AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
417 return shape_values[qpoint * n_shape_functions + shape_nr];
418}
419
420
421
422template <int dim, int spacedim, typename VectorType>
423bool
428
429
430
431template <int dim, int spacedim, typename VectorType>
432bool
433MappingFEField<dim, spacedim, VectorType>::is_compatible_with(
434 const ReferenceCell &reference_cell) const
435{
436 Assert(dim == reference_cell.get_dimension(),
437 ExcMessage("The dimension of your mapping (" +
438 Utilities::to_string(dim) +
439 ") and the reference cell cell_type (" +
440 Utilities::to_string(reference_cell.get_dimension()) +
441 " ) do not agree."));
442
443 return this->reference_cell == reference_cell;
444}
445
446
447
448template <int dim, int spacedim, typename VectorType>
449boost::container::small_vector<Point<spacedim>,
450#ifndef _MSC_VER
451 ReferenceCells::max_n_vertices<dim>()
452#else
453 GeometryInfo<dim>::vertices_per_cell
454#endif
455 >
456MappingFEField<dim, spacedim, VectorType>::get_vertices(
457 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
458{
459 // we transform our tria iterator into a dof iterator so we can access
460 // data not associated with triangulations
461 const typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(
462 *cell, euler_dof_handler);
463
464 Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
465 AssertDimension(cell->n_vertices(), fe_values.n_quadrature_points);
466 AssertDimension(fe_to_real.size(),
467 euler_dof_handler->get_fe().n_components());
468 if (uses_level_dofs)
469 {
470 AssertIndexRange(cell->level(), euler_vector.size());
471 AssertDimension(euler_vector[cell->level()]->size(),
472 euler_dof_handler->n_dofs(cell->level()));
473 }
474 else
475 AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
476
477 {
478 std::lock_guard<std::mutex> lock(fe_values_mutex);
479 fe_values.reinit(dof_cell);
480 }
481 const unsigned int dofs_per_cell =
482 euler_dof_handler->get_fe().n_dofs_per_cell();
483 std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
484 if (uses_level_dofs)
485 dof_cell->get_mg_dof_indices(dof_indices);
486 else
487 dof_cell->get_dof_indices(dof_indices);
488
489 const VectorType &vector =
490 uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
491
492 boost::container::small_vector<Point<spacedim>,
493#ifndef _MSC_VER
494 ReferenceCells::max_n_vertices<dim>()
495#else
496 GeometryInfo<dim>::vertices_per_cell
497#endif
498 >
499 vertices(cell->n_vertices());
500 for (unsigned int i = 0; i < dofs_per_cell; ++i)
501 {
502 const unsigned int comp = fe_to_real
503 [euler_dof_handler->get_fe().system_to_component_index(i).first];
504 if (comp != numbers::invalid_unsigned_int)
505 {
506 typename VectorType::value_type value =
507 internal::ElementAccess<VectorType>::get(vector, dof_indices[i]);
508 if (euler_dof_handler->get_fe().is_primitive(i))
509 for (const unsigned int v : cell->vertex_indices())
510 vertices[v][comp] += fe_values.shape_value(i, v) * value;
511 else
512 DEAL_II_NOT_IMPLEMENTED();
513 }
514 }
515
516 return vertices;
517}
518
519
520
521template <int dim, int spacedim, typename VectorType>
522UpdateFlags
523MappingFEField<dim, spacedim, VectorType>::requires_update_flags(
524 const UpdateFlags in) const
525{
526 // add flags if the respective quantities are necessary to compute
527 // what we need. note that some flags appear in both conditions and
528 // in subsequent set operations. this leads to some circular
529 // logic. the only way to treat this is to iterate. since there are
530 // 5 if-clauses in the loop, it will take at most 4 iterations to
531 // converge. do them:
532 UpdateFlags out = in;
533 for (unsigned int i = 0; i < 5; ++i)
534 {
535 // The following is a little incorrect:
536 // If not applied on a face,
537 // update_boundary_forms does not
538 // make sense. On the other hand,
539 // it is necessary on a
540 // face. Currently,
541 // update_boundary_forms is simply
542 // ignored for the interior of a
543 // cell.
544 if (out & (update_JxW_values | update_normal_vectors))
545 out |= update_boundary_forms;
546
547 if (out &
548 (update_covariant_transformation | update_jacobian_grads |
549 update_jacobians | update_boundary_forms | update_normal_vectors))
550 out |= update_contravariant_transformation;
551
552 if (out &
553 (update_inverse_jacobians | update_jacobian_pushed_forward_grads |
554 update_jacobian_pushed_forward_2nd_derivatives |
555 update_jacobian_pushed_forward_3rd_derivatives))
556 out |= update_covariant_transformation;
557
558 // The contravariant transformation is used in the Piola
559 // transformation, which requires the determinant of the Jacobi
560 // matrix of the transformation. Because we have no way of
561 // knowing here whether the finite element wants to use the
562 // contravariant or the Piola transforms, we add the volume elements
563 // to the list of flags to be updated for each cell.
564 if (out & update_contravariant_transformation)
565 out |= update_volume_elements;
566
567 if (out & update_normal_vectors)
568 out |= update_volume_elements;
569 }
570
571 return out;
572}
573
574
575template <int dim, int spacedim, typename VectorType>
576void
577MappingFEField<dim, spacedim, VectorType>::compute_face_data(
578 const unsigned int n_original_q_points,
579 InternalData &data) const
580{
581 // Set to the size of a single quadrature object for faces, as the size set
582 // in in reinit() is for all points
583 if (data.update_each & update_covariant_transformation)
584 data.covariant.resize(n_original_q_points);
585
586 if (data.update_each & update_contravariant_transformation)
587 data.contravariant.resize(n_original_q_points);
588
589 if (data.update_each & update_volume_elements)
590 data.volume_elements.resize(n_original_q_points);
591
592 if (dim > 1)
593 {
594 if (data.update_each & update_boundary_forms)
595 {
596 data.aux.resize(
597 dim - 1, std::vector<Tensor<1, spacedim>>(n_original_q_points));
598
599
600 // TODO: only a single reference cell type possible...
601 const auto n_faces = reference_cell.n_faces();
602
603 // Compute tangentials to the unit cell.
604 for (unsigned int i = 0; i < n_faces; ++i)
605 {
606 data.unit_tangentials[i].resize(n_original_q_points);
607 std::fill(
608 data.unit_tangentials[i].begin(),
609 data.unit_tangentials[i].end(),
610 reference_cell.template unit_tangential_vectors<dim>(i, 0));
611 if (dim > 2)
612 {
613 data.unit_tangentials[n_faces + i].resize(
614 n_original_q_points);
615 std::fill(
616 data.unit_tangentials[n_faces + i].begin(),
617 data.unit_tangentials[n_faces + i].end(),
618 reference_cell.template unit_tangential_vectors<dim>(i, 1));
619 }
620 }
621 }
622 }
623}
624
625
626
627template <int dim, int spacedim, typename VectorType>
628typename std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
629MappingFEField<dim, spacedim, VectorType>::get_data(
630 const UpdateFlags update_flags,
631 const Quadrature<dim> &quadrature) const
632{
633 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
634 std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
635 data_ptr->reinit(requires_update_flags(update_flags), quadrature);
636
637 return data_ptr;
638}
639
640
641
642template <int dim, int spacedim, typename VectorType>
643std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
644MappingFEField<dim, spacedim, VectorType>::get_face_data(
645 const UpdateFlags update_flags,
646 const hp::QCollection<dim - 1> &quadrature) const
647{
648 AssertDimension(quadrature.size(), 1);
649
650 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
651 std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
652 auto &data = dynamic_cast<InternalData &>(*data_ptr);
653
654 const Quadrature<dim> q(
655 QProjector<dim>::project_to_all_faces(reference_cell, quadrature[0]));
656 data.reinit(requires_update_flags(update_flags), q);
657 this->compute_face_data(quadrature[0].size(), data);
658
659 return data_ptr;
660}
661
662
663template <int dim, int spacedim, typename VectorType>
664std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
665MappingFEField<dim, spacedim, VectorType>::get_subface_data(
666 const UpdateFlags update_flags,
667 const Quadrature<dim - 1> &quadrature) const
668{
669 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
670 std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
671 auto &data = dynamic_cast<InternalData &>(*data_ptr);
672
673 const Quadrature<dim> q(
674 QProjector<dim>::project_to_all_subfaces(reference_cell, quadrature));
675 data.reinit(requires_update_flags(update_flags), q);
676 this->compute_face_data(quadrature.size(), data);
677
678 return data_ptr;
679}
680
681
682
683namespace internal
684{
685 namespace MappingFEFieldImplementation
686 {
687 namespace
688 {
695 template <int dim, int spacedim, typename VectorType>
696 void
697 maybe_compute_q_points(
698 const typename ::QProjector<dim>::DataSetDescriptor data_set,
699 const typename ::MappingFEField<dim, spacedim, VectorType>::
700 InternalData &data,
701 const FiniteElement<dim, spacedim> &fe,
702 const ComponentMask &fe_mask,
703 const std::vector<unsigned int> &fe_to_real,
704 std::vector<Point<spacedim>> &quadrature_points)
705 {
706 const UpdateFlags update_flags = data.update_each;
707
708 if (update_flags & update_quadrature_points)
709 {
710 for (unsigned int point = 0; point < quadrature_points.size();
711 ++point)
712 {
713 Point<spacedim> result;
714 const double *shape = &data.shape(point + data_set, 0);
715
716 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
717 {
718 const unsigned int comp_k =
719 fe.system_to_component_index(k).first;
720 if (fe_mask[comp_k])
721 result[fe_to_real[comp_k]] +=
722 data.local_dof_values[k] * shape[k];
723 }
724
725 quadrature_points[point] = result;
726 }
727 }
728 }
729
737 template <int dim, int spacedim, typename VectorType>
738 void
739 maybe_update_Jacobians(
740 const typename ::QProjector<dim>::DataSetDescriptor data_set,
741 const typename ::MappingFEField<dim, spacedim, VectorType>::
742 InternalData &data,
743 const FiniteElement<dim, spacedim> &fe,
744 const ComponentMask &fe_mask,
745 const std::vector<unsigned int> &fe_to_real)
746 {
747 const UpdateFlags update_flags = data.update_each;
748
749 // then Jacobians
750 if (update_flags & update_contravariant_transformation)
751 {
752 const unsigned int n_q_points = data.contravariant.size();
753
754 Assert(data.n_shape_functions > 0, ExcInternalError());
755
756 for (unsigned int point = 0; point < n_q_points; ++point)
757 {
758 const Tensor<1, dim> *data_derv =
759 &data.derivative(point + data_set, 0);
760
761 Tensor<1, dim> result[spacedim];
762
763 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
764 {
765 const unsigned int comp_k =
766 fe.system_to_component_index(k).first;
767 if (fe_mask[comp_k])
768 result[fe_to_real[comp_k]] +=
769 data.local_dof_values[k] * data_derv[k];
770 }
771
772 // write result into contravariant data
773 for (unsigned int i = 0; i < spacedim; ++i)
774 {
775 data.contravariant[point][i] = result[i];
776 }
777 }
778 }
779
780 if (update_flags & update_covariant_transformation)
781 {
782 AssertDimension(data.covariant.size(), data.contravariant.size());
783 for (unsigned int point = 0; point < data.contravariant.size();
784 ++point)
785 data.covariant[point] =
786 (data.contravariant[point]).covariant_form();
787 }
788
789 if (update_flags & update_volume_elements)
790 {
791 AssertDimension(data.contravariant.size(),
792 data.volume_elements.size());
793 for (unsigned int point = 0; point < data.contravariant.size();
794 ++point)
795 data.volume_elements[point] =
796 data.contravariant[point].determinant();
797 }
798 }
799
806 template <int dim, int spacedim, typename VectorType>
807 void
808 maybe_update_jacobian_grads(
809 const typename ::QProjector<dim>::DataSetDescriptor data_set,
810 const typename ::MappingFEField<dim, spacedim, VectorType>::
811 InternalData &data,
812 const FiniteElement<dim, spacedim> &fe,
813 const ComponentMask &fe_mask,
814 const std::vector<unsigned int> &fe_to_real,
815 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
816 {
817 const UpdateFlags update_flags = data.update_each;
818 if (update_flags & update_jacobian_grads)
819 {
820 const unsigned int n_q_points = jacobian_grads.size();
821
822 for (unsigned int point = 0; point < n_q_points; ++point)
823 {
824 const Tensor<2, dim> *second =
825 &data.second_derivative(point + data_set, 0);
826
827 DerivativeForm<2, dim, spacedim> result;
828
829 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
830 {
831 const unsigned int comp_k =
832 fe.system_to_component_index(k).first;
833 if (fe_mask[comp_k])
834 for (unsigned int j = 0; j < dim; ++j)
835 for (unsigned int l = 0; l < dim; ++l)
836 result[fe_to_real[comp_k]][j][l] +=
837 (second[k][j][l] * data.local_dof_values[k]);
838 }
839
840 // never touch any data for j=dim in case dim<spacedim, so
841 // it will always be zero as it was initialized
842 for (unsigned int i = 0; i < spacedim; ++i)
843 for (unsigned int j = 0; j < dim; ++j)
844 for (unsigned int l = 0; l < dim; ++l)
845 jacobian_grads[point][i][j][l] = result[i][j][l];
846 }
847 }
848 }
849
856 template <int dim, int spacedim, typename VectorType>
857 void
858 maybe_update_jacobian_pushed_forward_grads(
859 const typename ::QProjector<dim>::DataSetDescriptor data_set,
860 const typename ::MappingFEField<dim, spacedim, VectorType>::
861 InternalData &data,
862 const FiniteElement<dim, spacedim> &fe,
863 const ComponentMask &fe_mask,
864 const std::vector<unsigned int> &fe_to_real,
865 std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads)
866 {
867 const UpdateFlags update_flags = data.update_each;
868 if (update_flags & update_jacobian_pushed_forward_grads)
869 {
870 const unsigned int n_q_points =
871 jacobian_pushed_forward_grads.size();
872
873 double tmp[spacedim][spacedim][spacedim];
874 for (unsigned int point = 0; point < n_q_points; ++point)
875 {
876 const Tensor<2, dim> *second =
877 &data.second_derivative(point + data_set, 0);
878
879 DerivativeForm<2, dim, spacedim> result;
880
881 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
882 {
883 const unsigned int comp_k =
884 fe.system_to_component_index(k).first;
885 if (fe_mask[comp_k])
886 for (unsigned int j = 0; j < dim; ++j)
887 for (unsigned int l = 0; l < dim; ++l)
888 result[fe_to_real[comp_k]][j][l] +=
889 (second[k][j][l] * data.local_dof_values[k]);
890 }
891
892 // first push forward the j-components
893 for (unsigned int i = 0; i < spacedim; ++i)
894 for (unsigned int j = 0; j < spacedim; ++j)
895 for (unsigned int l = 0; l < dim; ++l)
896 {
897 tmp[i][j][l] =
898 result[i][0][l] * data.covariant[point][j][0];
899 for (unsigned int jr = 1; jr < dim; ++jr)
900 {
901 tmp[i][j][l] +=
902 result[i][jr][l] * data.covariant[point][j][jr];
903 }
904 }
905
906 // now, pushing forward the l-components
907 for (unsigned int i = 0; i < spacedim; ++i)
908 for (unsigned int j = 0; j < spacedim; ++j)
909 for (unsigned int l = 0; l < spacedim; ++l)
910 {
911 jacobian_pushed_forward_grads[point][i][j][l] =
912 tmp[i][j][0] * data.covariant[point][l][0];
913 for (unsigned int lr = 1; lr < dim; ++lr)
914 {
915 jacobian_pushed_forward_grads[point][i][j][l] +=
916 tmp[i][j][lr] * data.covariant[point][l][lr];
917 }
918 }
919 }
920 }
921 }
922
929 template <int dim, int spacedim, typename VectorType>
930 void
931 maybe_update_jacobian_2nd_derivatives(
932 const typename ::QProjector<dim>::DataSetDescriptor data_set,
933 const typename ::MappingFEField<dim, spacedim, VectorType>::
934 InternalData &data,
935 const FiniteElement<dim, spacedim> &fe,
936 const ComponentMask &fe_mask,
937 const std::vector<unsigned int> &fe_to_real,
938 std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
939 {
940 const UpdateFlags update_flags = data.update_each;
941 if (update_flags & update_jacobian_2nd_derivatives)
942 {
943 const unsigned int n_q_points = jacobian_2nd_derivatives.size();
944
945 for (unsigned int point = 0; point < n_q_points; ++point)
946 {
947 const Tensor<3, dim> *third =
948 &data.third_derivative(point + data_set, 0);
949
950 DerivativeForm<3, dim, spacedim> result;
951
952 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
953 {
954 const unsigned int comp_k =
955 fe.system_to_component_index(k).first;
956 if (fe_mask[comp_k])
957 for (unsigned int j = 0; j < dim; ++j)
958 for (unsigned int l = 0; l < dim; ++l)
959 for (unsigned int m = 0; m < dim; ++m)
960 result[fe_to_real[comp_k]][j][l][m] +=
961 (third[k][j][l][m] * data.local_dof_values[k]);
962 }
963
964 // never touch any data for j=dim in case dim<spacedim, so
965 // it will always be zero as it was initialized
966 for (unsigned int i = 0; i < spacedim; ++i)
967 for (unsigned int j = 0; j < dim; ++j)
968 for (unsigned int l = 0; l < dim; ++l)
969 for (unsigned int m = 0; m < dim; ++m)
970 jacobian_2nd_derivatives[point][i][j][l][m] =
971 result[i][j][l][m];
972 }
973 }
974 }
975
983 template <int dim, int spacedim, typename VectorType>
984 void
985 maybe_update_jacobian_pushed_forward_2nd_derivatives(
986 const typename ::QProjector<dim>::DataSetDescriptor data_set,
987 const typename ::MappingFEField<dim, spacedim, VectorType>::
988 InternalData &data,
989 const FiniteElement<dim, spacedim> &fe,
990 const ComponentMask &fe_mask,
991 const std::vector<unsigned int> &fe_to_real,
992 std::vector<Tensor<4, spacedim>>
993 &jacobian_pushed_forward_2nd_derivatives)
994 {
995 const UpdateFlags update_flags = data.update_each;
996 if (update_flags & update_jacobian_pushed_forward_2nd_derivatives)
997 {
998 const unsigned int n_q_points =
999 jacobian_pushed_forward_2nd_derivatives.size();
1000
1001 double tmp[spacedim][spacedim][spacedim][spacedim];
1002 for (unsigned int point = 0; point < n_q_points; ++point)
1003 {
1004 const Tensor<3, dim> *third =
1005 &data.third_derivative(point + data_set, 0);
1006
1007 DerivativeForm<3, dim, spacedim> result;
1008
1009 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1010 {
1011 const unsigned int comp_k =
1012 fe.system_to_component_index(k).first;
1013 if (fe_mask[comp_k])
1014 for (unsigned int j = 0; j < dim; ++j)
1015 for (unsigned int l = 0; l < dim; ++l)
1016 for (unsigned int m = 0; m < dim; ++m)
1017 result[fe_to_real[comp_k]][j][l][m] +=
1018 (third[k][j][l][m] * data.local_dof_values[k]);
1019 }
1020
1021 // push forward the j-coordinate
1022 for (unsigned int i = 0; i < spacedim; ++i)
1023 for (unsigned int j = 0; j < spacedim; ++j)
1024 for (unsigned int l = 0; l < dim; ++l)
1025 for (unsigned int m = 0; m < dim; ++m)
1026 {
1027 jacobian_pushed_forward_2nd_derivatives
1028 [point][i][j][l][m] =
1029 result[i][0][l][m] * data.covariant[point][j][0];
1030 for (unsigned int jr = 1; jr < dim; ++jr)
1031 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1032 [l][m] +=
1033 result[i][jr][l][m] *
1034 data.covariant[point][j][jr];
1035 }
1036
1037 // push forward the l-coordinate
1038 for (unsigned int i = 0; i < spacedim; ++i)
1039 for (unsigned int j = 0; j < spacedim; ++j)
1040 for (unsigned int l = 0; l < spacedim; ++l)
1041 for (unsigned int m = 0; m < dim; ++m)
1042 {
1043 tmp[i][j][l][m] =
1044 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1045 [0][m] *
1046 data.covariant[point][l][0];
1047 for (unsigned int lr = 1; lr < dim; ++lr)
1048 tmp[i][j][l][m] +=
1049 jacobian_pushed_forward_2nd_derivatives[point][i]
1050 [j][lr]
1051 [m] *
1052 data.covariant[point][l][lr];
1053 }
1054
1055 // push forward the m-coordinate
1056 for (unsigned int i = 0; i < spacedim; ++i)
1057 for (unsigned int j = 0; j < spacedim; ++j)
1058 for (unsigned int l = 0; l < spacedim; ++l)
1059 for (unsigned int m = 0; m < spacedim; ++m)
1060 {
1061 jacobian_pushed_forward_2nd_derivatives
1062 [point][i][j][l][m] =
1063 tmp[i][j][l][0] * data.covariant[point][m][0];
1064 for (unsigned int mr = 1; mr < dim; ++mr)
1065 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1066 [l][m] +=
1067 tmp[i][j][l][mr] * data.covariant[point][m][mr];
1068 }
1069 }
1070 }
1071 }
1072
1079 template <int dim, int spacedim, typename VectorType>
1080 void
1081 maybe_update_jacobian_3rd_derivatives(
1082 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1083 const typename ::MappingFEField<dim, spacedim, VectorType>::
1084 InternalData &data,
1085 const FiniteElement<dim, spacedim> &fe,
1086 const ComponentMask &fe_mask,
1087 const std::vector<unsigned int> &fe_to_real,
1088 std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1089 {
1090 const UpdateFlags update_flags = data.update_each;
1091 if (update_flags & update_jacobian_3rd_derivatives)
1092 {
1093 const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1094
1095 for (unsigned int point = 0; point < n_q_points; ++point)
1096 {
1097 const Tensor<4, dim> *fourth =
1098 &data.fourth_derivative(point + data_set, 0);
1099
1100 DerivativeForm<4, dim, spacedim> result;
1101
1102 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1103 {
1104 const unsigned int comp_k =
1105 fe.system_to_component_index(k).first;
1106 if (fe_mask[comp_k])
1107 for (unsigned int j = 0; j < dim; ++j)
1108 for (unsigned int l = 0; l < dim; ++l)
1109 for (unsigned int m = 0; m < dim; ++m)
1110 for (unsigned int n = 0; n < dim; ++n)
1111 result[fe_to_real[comp_k]][j][l][m][n] +=
1112 (fourth[k][j][l][m][n] *
1113 data.local_dof_values[k]);
1114 }
1115
1116 // never touch any data for j,l,m,n=dim in case
1117 // dim<spacedim, so it will always be zero as it was
1118 // initialized
1119 for (unsigned int i = 0; i < spacedim; ++i)
1120 for (unsigned int j = 0; j < dim; ++j)
1121 for (unsigned int l = 0; l < dim; ++l)
1122 for (unsigned int m = 0; m < dim; ++m)
1123 for (unsigned int n = 0; n < dim; ++n)
1124 jacobian_3rd_derivatives[point][i][j][l][m][n] =
1125 result[i][j][l][m][n];
1126 }
1127 }
1128 }
1129
1137 template <int dim, int spacedim, typename VectorType>
1138 void
1139 maybe_update_jacobian_pushed_forward_3rd_derivatives(
1140 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1141 const typename ::MappingFEField<dim, spacedim, VectorType>::
1142 InternalData &data,
1143 const FiniteElement<dim, spacedim> &fe,
1144 const ComponentMask &fe_mask,
1145 const std::vector<unsigned int> &fe_to_real,
1146 std::vector<Tensor<5, spacedim>>
1147 &jacobian_pushed_forward_3rd_derivatives)
1148 {
1149 const UpdateFlags update_flags = data.update_each;
1150 if (update_flags & update_jacobian_pushed_forward_3rd_derivatives)
1151 {
1152 const unsigned int n_q_points =
1153 jacobian_pushed_forward_3rd_derivatives.size();
1154
1155 double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
1156 for (unsigned int point = 0; point < n_q_points; ++point)
1157 {
1158 const Tensor<4, dim> *fourth =
1159 &data.fourth_derivative(point + data_set, 0);
1160
1161 DerivativeForm<4, dim, spacedim> result;
1162
1163 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1164 {
1165 const unsigned int comp_k =
1166 fe.system_to_component_index(k).first;
1167 if (fe_mask[comp_k])
1168 for (unsigned int j = 0; j < dim; ++j)
1169 for (unsigned int l = 0; l < dim; ++l)
1170 for (unsigned int m = 0; m < dim; ++m)
1171 for (unsigned int n = 0; n < dim; ++n)
1172 result[fe_to_real[comp_k]][j][l][m][n] +=
1173 (fourth[k][j][l][m][n] *
1174 data.local_dof_values[k]);
1175 }
1176
1177 // push-forward the j-coordinate
1178 for (unsigned int i = 0; i < spacedim; ++i)
1179 for (unsigned int j = 0; j < spacedim; ++j)
1180 for (unsigned int l = 0; l < dim; ++l)
1181 for (unsigned int m = 0; m < dim; ++m)
1182 for (unsigned int n = 0; n < dim; ++n)
1183 {
1184 tmp[i][j][l][m][n] = result[i][0][l][m][n] *
1185 data.covariant[point][j][0];
1186 for (unsigned int jr = 1; jr < dim; ++jr)
1187 tmp[i][j][l][m][n] +=
1188 result[i][jr][l][m][n] *
1189 data.covariant[point][j][jr];
1190 }
1191
1192 // push-forward the l-coordinate
1193 for (unsigned int i = 0; i < spacedim; ++i)
1194 for (unsigned int j = 0; j < spacedim; ++j)
1195 for (unsigned int l = 0; l < spacedim; ++l)
1196 for (unsigned int m = 0; m < dim; ++m)
1197 for (unsigned int n = 0; n < dim; ++n)
1198 {
1199 jacobian_pushed_forward_3rd_derivatives
1200 [point][i][j][l][m][n] =
1201 tmp[i][j][0][m][n] *
1202 data.covariant[point][l][0];
1203 for (unsigned int lr = 1; lr < dim; ++lr)
1204 jacobian_pushed_forward_3rd_derivatives[point][i]
1205 [j][l][m]
1206 [n] +=
1207 tmp[i][j][lr][m][n] *
1208 data.covariant[point][l][lr];
1209 }
1210
1211 // push-forward the m-coordinate
1212 for (unsigned int i = 0; i < spacedim; ++i)
1213 for (unsigned int j = 0; j < spacedim; ++j)
1214 for (unsigned int l = 0; l < spacedim; ++l)
1215 for (unsigned int m = 0; m < spacedim; ++m)
1216 for (unsigned int n = 0; n < dim; ++n)
1217 {
1218 tmp[i][j][l][m][n] =
1219 jacobian_pushed_forward_3rd_derivatives[point][i]
1220 [j][l][0]
1221 [n] *
1222 data.covariant[point][m][0];
1223 for (unsigned int mr = 1; mr < dim; ++mr)
1224 tmp[i][j][l][m][n] +=
1225 jacobian_pushed_forward_3rd_derivatives[point]
1226 [i][j][l]
1227 [mr][n] *
1228 data.covariant[point][m][mr];
1229 }
1230
1231 // push-forward the n-coordinate
1232 for (unsigned int i = 0; i < spacedim; ++i)
1233 for (unsigned int j = 0; j < spacedim; ++j)
1234 for (unsigned int l = 0; l < spacedim; ++l)
1235 for (unsigned int m = 0; m < spacedim; ++m)
1236 for (unsigned int n = 0; n < spacedim; ++n)
1237 {
1238 jacobian_pushed_forward_3rd_derivatives
1239 [point][i][j][l][m][n] =
1240 tmp[i][j][l][m][0] *
1241 data.covariant[point][n][0];
1242 for (unsigned int nr = 1; nr < dim; ++nr)
1243 jacobian_pushed_forward_3rd_derivatives[point][i]
1244 [j][l][m]
1245 [n] +=
1246 tmp[i][j][l][m][nr] *
1247 data.covariant[point][n][nr];
1248 }
1249 }
1250 }
1251 }
1252
1253
1263 template <int dim, int spacedim, typename VectorType>
1264 void
1265 maybe_compute_face_data(
1266 const ::Mapping<dim, spacedim> &mapping,
1267 const typename ::Triangulation<dim, spacedim>::cell_iterator
1268 &cell,
1269 const unsigned int face_no,
1270 const unsigned int subface_no,
1271 const typename QProjector<dim>::DataSetDescriptor data_set,
1272 const typename ::MappingFEField<dim, spacedim, VectorType>::
1273 InternalData &data,
1274 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1275 &output_data)
1276 {
1277 const UpdateFlags update_flags = data.update_each;
1278
1279 if (update_flags & update_boundary_forms)
1280 {
1281 const unsigned int n_q_points = output_data.boundary_forms.size();
1282 if (update_flags & update_normal_vectors)
1283 AssertDimension(output_data.normal_vectors.size(), n_q_points);
1284 if (update_flags & update_JxW_values)
1285 AssertDimension(output_data.JxW_values.size(), n_q_points);
1286
1287 // map the unit tangentials to the real cell. checking for d!=dim-1
1288 // eliminates compiler warnings regarding unsigned int expressions <
1289 // 0.
1290 for (unsigned int d = 0; d != dim - 1; ++d)
1291 {
1292 Assert(face_no + cell->n_faces() * d <
1293 data.unit_tangentials.size(),
1294 ExcInternalError());
1295 Assert(
1296 data.aux[d].size() <=
1297 data.unit_tangentials[face_no + cell->n_faces() * d].size(),
1298 ExcInternalError());
1299
1300 mapping.transform(
1301 make_array_view(
1302 data.unit_tangentials[face_no + cell->n_faces() * d]),
1303 mapping_contravariant,
1304 data,
1305 make_array_view(data.aux[d]));
1306 }
1307
1308 // if dim==spacedim, we can use the unit tangentials to compute the
1309 // boundary form by simply taking the cross product
1310 if (dim == spacedim)
1311 {
1312 for (unsigned int i = 0; i < n_q_points; ++i)
1313 switch (dim)
1314 {
1315 case 1:
1316 // in 1d, we don't have access to any of the data.aux
1317 // fields (because it has only dim-1 components), but we
1318 // can still compute the boundary form by simply looking
1319 // at the number of the face
1320 output_data.boundary_forms[i][0] =
1321 (face_no == 0 ? -1 : +1);
1322 break;
1323 case 2:
1324 output_data.boundary_forms[i] =
1325 cross_product_2d(data.aux[0][i]);
1326 break;
1327 case 3:
1328 output_data.boundary_forms[i] =
1329 cross_product_3d(data.aux[0][i], data.aux[1][i]);
1330 break;
1331 default:
1332 DEAL_II_NOT_IMPLEMENTED();
1333 }
1334 }
1335 else //(dim < spacedim)
1336 {
1337 // in the codim-one case, the boundary form results from the
1338 // cross product of all the face tangential vectors and the cell
1339 // normal vector
1340 //
1341 // to compute the cell normal, use the same method used in
1342 // fill_fe_values for cells above
1343 AssertDimension(data.contravariant.size(), n_q_points);
1344
1345 for (unsigned int point = 0; point < n_q_points; ++point)
1346 {
1347 if (dim == 1)
1348 {
1349 // J is a tangent vector
1350 output_data.boundary_forms[point] =
1351 data.contravariant[point].transpose()[0];
1352 output_data.boundary_forms[point] /=
1353 (face_no == 0 ? -1. : +1.) *
1354 output_data.boundary_forms[point].norm();
1355 }
1356
1357 if (dim == 2)
1358 {
1359 const DerivativeForm<1, spacedim, dim> DX_t =
1360 data.contravariant[point].transpose();
1361
1362 Tensor<1, spacedim> cell_normal =
1363 cross_product_3d(DX_t[0], DX_t[1]);
1364 cell_normal /= cell_normal.norm();
1365
1366 // then compute the face normal from the face tangent
1367 // and the cell normal:
1368 output_data.boundary_forms[point] =
1369 cross_product_3d(data.aux[0][point], cell_normal);
1370 }
1371 }
1372 }
1373
1374 if (update_flags & (update_normal_vectors | update_JxW_values))
1375 for (unsigned int i = 0; i < output_data.boundary_forms.size();
1376 ++i)
1377 {
1378 if (update_flags & update_JxW_values)
1379 {
1380 output_data.JxW_values[i] =
1381 output_data.boundary_forms[i].norm() *
1382 data.quadrature_weights[i + data_set];
1383
1384 if (subface_no != numbers::invalid_unsigned_int)
1385 {
1386 // TODO
1387 const double area_ratio =
1388 GeometryInfo<dim>::subface_ratio(
1389 cell->subface_case(face_no), subface_no);
1390 output_data.JxW_values[i] *= area_ratio;
1391 }
1392 }
1393
1394 if (update_flags & update_normal_vectors)
1395 output_data.normal_vectors[i] =
1396 Point<spacedim>(output_data.boundary_forms[i] /
1397 output_data.boundary_forms[i].norm());
1398 }
1399 }
1400 }
1401
1408 template <int dim, int spacedim, typename VectorType>
1409 void
1410 do_fill_fe_face_values(
1411 const ::Mapping<dim, spacedim> &mapping,
1412 const typename ::Triangulation<dim, spacedim>::cell_iterator
1413 &cell,
1414 const unsigned int face_no,
1415 const unsigned int subface_no,
1416 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1417 const typename ::MappingFEField<dim, spacedim, VectorType>::
1418 InternalData &data,
1419 const FiniteElement<dim, spacedim> &fe,
1420 const ComponentMask &fe_mask,
1421 const std::vector<unsigned int> &fe_to_real,
1422 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1423 &output_data)
1424 {
1425 maybe_compute_q_points<dim, spacedim, VectorType>(
1426 data_set,
1427 data,
1428 fe,
1429 fe_mask,
1430 fe_to_real,
1431 output_data.quadrature_points);
1432
1433 maybe_update_Jacobians<dim, spacedim, VectorType>(
1434 data_set, data, fe, fe_mask, fe_to_real);
1435
1436 const UpdateFlags update_flags = data.update_each;
1437 const unsigned int n_q_points = data.contravariant.size();
1438
1439 if (update_flags & update_jacobians)
1440 for (unsigned int point = 0; point < n_q_points; ++point)
1441 output_data.jacobians[point] = data.contravariant[point];
1442
1443 if (update_flags & update_inverse_jacobians)
1444 for (unsigned int point = 0; point < n_q_points; ++point)
1445 output_data.inverse_jacobians[point] =
1446 data.covariant[point].transpose();
1447
1448 maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1449 data_set, data, fe, fe_mask, fe_to_real, output_data.jacobian_grads);
1450
1451 maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1452 data_set,
1453 data,
1454 fe,
1455 fe_mask,
1456 fe_to_real,
1457 output_data.jacobian_pushed_forward_grads);
1458
1459 maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1460 data_set,
1461 data,
1462 fe,
1463 fe_mask,
1464 fe_to_real,
1465 output_data.jacobian_2nd_derivatives);
1466
1467 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1468 spacedim,
1469 VectorType>(
1470 data_set,
1471 data,
1472 fe,
1473 fe_mask,
1474 fe_to_real,
1475 output_data.jacobian_pushed_forward_2nd_derivatives);
1476
1477 maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1478 data_set,
1479 data,
1480 fe,
1481 fe_mask,
1482 fe_to_real,
1483 output_data.jacobian_3rd_derivatives);
1484
1485 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1486 spacedim,
1487 VectorType>(
1488 data_set,
1489 data,
1490 fe,
1491 fe_mask,
1492 fe_to_real,
1493 output_data.jacobian_pushed_forward_3rd_derivatives);
1494
1495 maybe_compute_face_data<dim, spacedim, VectorType>(
1496 mapping, cell, face_no, subface_no, data_set, data, output_data);
1497 }
1498 } // namespace
1499 } // namespace MappingFEFieldImplementation
1500} // namespace internal
1501
1502
1503// Note that the CellSimilarity flag is modifiable, since MappingFEField can
1504// need to recalculate data even when cells are similar.
1505template <int dim, int spacedim, typename VectorType>
1506CellSimilarity::Similarity
1507MappingFEField<dim, spacedim, VectorType>::fill_fe_values(
1508 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1509 const CellSimilarity::Similarity,
1510 const Quadrature<dim> &quadrature,
1511 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1512 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1513 &output_data) const
1514{
1515 // convert data object to internal data for this class. fails with an
1516 // exception if that is not possible
1517 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1518 ExcInternalError());
1519 const InternalData &data = static_cast<const InternalData &>(internal_data);
1520
1521 const unsigned int n_q_points = quadrature.size();
1522
1523 update_internal_dofs(cell, data);
1524
1525 internal::MappingFEFieldImplementation::
1526 maybe_compute_q_points<dim, spacedim, VectorType>(
1527 QProjector<dim>::DataSetDescriptor::cell(),
1528 data,
1529 euler_dof_handler->get_fe(),
1530 fe_mask,
1531 fe_to_real,
1532 output_data.quadrature_points);
1533
1534 internal::MappingFEFieldImplementation::
1535 maybe_update_Jacobians<dim, spacedim, VectorType>(
1536 QProjector<dim>::DataSetDescriptor::cell(),
1537 data,
1538 euler_dof_handler->get_fe(),
1539 fe_mask,
1540 fe_to_real);
1541
1542 const UpdateFlags update_flags = data.update_each;
1543 const std::vector<double> &weights = quadrature.get_weights();
1544
1545 // Multiply quadrature weights by absolute value of Jacobian determinants or
1546 // the area element g=sqrt(DX^t DX) in case of codim > 0
1547
1548 if (update_flags & (update_normal_vectors | update_JxW_values))
1549 {
1550 AssertDimension(output_data.JxW_values.size(), n_q_points);
1551
1552 Assert(!(update_flags & update_normal_vectors) ||
1553 (output_data.normal_vectors.size() == n_q_points),
1554 ExcDimensionMismatch(output_data.normal_vectors.size(),
1555 n_q_points));
1556
1557
1558 for (unsigned int point = 0; point < n_q_points; ++point)
1559 {
1560 if (dim == spacedim)
1561 {
1562 const double det = data.volume_elements[point];
1563
1564 // check for distorted cells.
1565
1566 // TODO: this allows for anisotropies of up to 1e6 in 3d and
1567 // 1e12 in 2d. might want to find a finer
1568 // (dimension-independent) criterion
1569 Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1570 cell->diameter() / std::sqrt(double(dim))),
1571 (typename Mapping<dim, spacedim>::ExcDistortedMappedCell(
1572 cell->center(), det, point)));
1573 output_data.JxW_values[point] = weights[point] * det;
1574 }
1575 // if dim==spacedim, then there is no cell normal to
1576 // compute. since this is for FEValues (and not FEFaceValues),
1577 // there are also no face normals to compute
1578 else // codim>0 case
1579 {
1580 Tensor<1, spacedim> DX_t[dim];
1581 for (unsigned int i = 0; i < spacedim; ++i)
1582 for (unsigned int j = 0; j < dim; ++j)
1583 DX_t[j][i] = data.contravariant[point][i][j];
1584
1585 Tensor<2, dim> G; // First fundamental form
1586 for (unsigned int i = 0; i < dim; ++i)
1587 for (unsigned int j = 0; j < dim; ++j)
1588 G[i][j] = DX_t[i] * DX_t[j];
1589
1590 output_data.JxW_values[point] =
1591 std::sqrt(determinant(G)) * weights[point];
1592
1593 if (update_flags & update_normal_vectors)
1594 {
1595 Assert(spacedim - dim == 1,
1596 ExcMessage("There is no cell normal in codim 2."));
1597
1598 if (dim == 1)
1599 output_data.normal_vectors[point] =
1600 cross_product_2d(-DX_t[0]);
1601 else
1602 {
1603 Assert(dim == 2, ExcInternalError());
1604
1605 // dim-1==1 for the second argument, but this
1606 // avoids a compiler warning about array bounds:
1607 output_data.normal_vectors[point] =
1608 cross_product_3d(DX_t[0], DX_t[dim - 1]);
1609 }
1610
1611 output_data.normal_vectors[point] /=
1612 output_data.normal_vectors[point].norm();
1613
1614 if (cell->direction_flag() == false)
1615 output_data.normal_vectors[point] *= -1.;
1616 }
1617 } // codim>0 case
1618 }
1619 }
1620
1621 // copy values from InternalData to vector given by reference
1622 if (update_flags & update_jacobians)
1623 {
1624 AssertDimension(output_data.jacobians.size(), n_q_points);
1625 for (unsigned int point = 0; point < n_q_points; ++point)
1626 output_data.jacobians[point] = data.contravariant[point];
1627 }
1628
1629 // copy values from InternalData to vector given by reference
1630 if (update_flags & update_inverse_jacobians)
1631 {
1632 AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1633 for (unsigned int point = 0; point < n_q_points; ++point)
1634 output_data.inverse_jacobians[point] =
1635 data.covariant[point].transpose();
1636 }
1637
1638 // calculate derivatives of the Jacobians
1639 internal::MappingFEFieldImplementation::
1640 maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1641 QProjector<dim>::DataSetDescriptor::cell(),
1642 data,
1643 euler_dof_handler->get_fe(),
1644 fe_mask,
1645 fe_to_real,
1646 output_data.jacobian_grads);
1647
1648 // calculate derivatives of the Jacobians pushed forward to real cell
1649 // coordinates
1650 internal::MappingFEFieldImplementation::
1651 maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1652 QProjector<dim>::DataSetDescriptor::cell(),
1653 data,
1654 euler_dof_handler->get_fe(),
1655 fe_mask,
1656 fe_to_real,
1657 output_data.jacobian_pushed_forward_grads);
1658
1659 // calculate hessians of the Jacobians
1660 internal::MappingFEFieldImplementation::
1661 maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1662 QProjector<dim>::DataSetDescriptor::cell(),
1663 data,
1664 euler_dof_handler->get_fe(),
1665 fe_mask,
1666 fe_to_real,
1667 output_data.jacobian_2nd_derivatives);
1668
1669 // calculate hessians of the Jacobians pushed forward to real cell coordinates
1670 internal::MappingFEFieldImplementation::
1671 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1672 spacedim,
1673 VectorType>(
1674 QProjector<dim>::DataSetDescriptor::cell(),
1675 data,
1676 euler_dof_handler->get_fe(),
1677 fe_mask,
1678 fe_to_real,
1679 output_data.jacobian_pushed_forward_2nd_derivatives);
1680
1681 // calculate gradients of the hessians of the Jacobians
1682 internal::MappingFEFieldImplementation::
1683 maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1684 QProjector<dim>::DataSetDescriptor::cell(),
1685 data,
1686 euler_dof_handler->get_fe(),
1687 fe_mask,
1688 fe_to_real,
1689 output_data.jacobian_3rd_derivatives);
1690
1691 // calculate gradients of the hessians of the Jacobians pushed forward to real
1692 // cell coordinates
1693 internal::MappingFEFieldImplementation::
1694 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1695 spacedim,
1696 VectorType>(
1697 QProjector<dim>::DataSetDescriptor::cell(),
1698 data,
1699 euler_dof_handler->get_fe(),
1700 fe_mask,
1701 fe_to_real,
1702 output_data.jacobian_pushed_forward_3rd_derivatives);
1703
1704 return CellSimilarity::invalid_next_cell;
1705}
1706
1707
1708
1709template <int dim, int spacedim, typename VectorType>
1710void
1711MappingFEField<dim, spacedim, VectorType>::fill_fe_face_values(
1712 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1713 const unsigned int face_no,
1714 const hp::QCollection<dim - 1> &quadrature,
1715 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1716 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1717 &output_data) const
1718{
1719 AssertDimension(quadrature.size(), 1);
1720
1721 // convert data object to internal data for this class. fails with an
1722 // exception if that is not possible
1723 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1724 ExcInternalError());
1725 const InternalData &data = static_cast<const InternalData &>(internal_data);
1726
1727 update_internal_dofs(cell, data);
1728
1729 internal::MappingFEFieldImplementation::
1730 do_fill_fe_face_values<dim, spacedim, VectorType>(
1731 *this,
1732 cell,
1733 face_no,
1734 numbers::invalid_unsigned_int,
1735 QProjector<dim>::DataSetDescriptor::face(reference_cell,
1736 face_no,
1737 cell->combined_face_orientation(
1738 face_no),
1739 quadrature[0].size()),
1740 data,
1741 euler_dof_handler->get_fe(),
1742 fe_mask,
1743 fe_to_real,
1744 output_data);
1745}
1746
1747
1748template <int dim, int spacedim, typename VectorType>
1749void
1750MappingFEField<dim, spacedim, VectorType>::fill_fe_subface_values(
1751 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1752 const unsigned int face_no,
1753 const unsigned int subface_no,
1754 const Quadrature<dim - 1> &quadrature,
1755 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1756 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1757 &output_data) const
1758{
1759 // convert data object to internal data for this class. fails with an
1760 // exception if that is not possible
1761 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1762 ExcInternalError());
1763 const InternalData &data = static_cast<const InternalData &>(internal_data);
1764
1765 update_internal_dofs(cell, data);
1766
1767 internal::MappingFEFieldImplementation::do_fill_fe_face_values<dim,
1768 spacedim,
1769 VectorType>(
1770 *this,
1771 cell,
1772 face_no,
1773 numbers::invalid_unsigned_int,
1774 QProjector<dim>::DataSetDescriptor::subface(reference_cell,
1775 face_no,
1776 subface_no,
1777 cell->combined_face_orientation(
1778 face_no),
1779 quadrature.size(),
1780 cell->subface_case(face_no)),
1781 data,
1782 euler_dof_handler->get_fe(),
1783 fe_mask,
1784 fe_to_real,
1785 output_data);
1786}
1787
1788
1789
1790template <int dim, int spacedim, typename VectorType>
1791void
1792MappingFEField<dim, spacedim, VectorType>::fill_fe_immersed_surface_values(
1793 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1794 const NonMatching::ImmersedSurfaceQuadrature<dim> &quadrature,
1795 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1796 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1797 &output_data) const
1798{
1799 AssertDimension(dim, spacedim);
1800 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1801 ExcInternalError());
1802 const InternalData &data = static_cast<const InternalData &>(internal_data);
1803
1804 const unsigned int n_q_points = quadrature.size();
1805
1806 update_internal_dofs(cell, data);
1807
1808 internal::MappingFEFieldImplementation::
1809 maybe_compute_q_points<dim, spacedim, VectorType>(
1810 QProjector<dim>::DataSetDescriptor::cell(),
1811 data,
1812 euler_dof_handler->get_fe(),
1813 fe_mask,
1814 fe_to_real,
1815 output_data.quadrature_points);
1816
1817 internal::MappingFEFieldImplementation::
1818 maybe_update_Jacobians<dim, spacedim, VectorType>(
1819 QProjector<dim>::DataSetDescriptor::cell(),
1820 data,
1821 euler_dof_handler->get_fe(),
1822 fe_mask,
1823 fe_to_real);
1824
1825 const UpdateFlags update_flags = data.update_each;
1826 const std::vector<double> &weights = quadrature.get_weights();
1827
1828 if (update_flags & (update_normal_vectors | update_JxW_values))
1829 {
1830 AssertDimension(output_data.JxW_values.size(), n_q_points);
1831
1832 Assert(!(update_flags & update_normal_vectors) ||
1833 (output_data.normal_vectors.size() == n_q_points),
1834 ExcDimensionMismatch(output_data.normal_vectors.size(),
1835 n_q_points));
1836
1837
1838 for (unsigned int point = 0; point < n_q_points; ++point)
1839 {
1840 const double det = data.volume_elements[point];
1841
1842 // check for distorted cells.
1843
1844 // TODO: this allows for anisotropies of up to 1e6 in 3d and
1845 // 1e12 in 2d. might want to find a finer
1846 // (dimension-independent) criterion
1847 Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1848 cell->diameter() / std::sqrt(double(dim))),
1849 (typename Mapping<dim, spacedim>::ExcDistortedMappedCell(
1850 cell->center(), det, point)));
1851
1852 // The normals are n = J^{-T} * \hat{n} before normalizing.
1853 Tensor<1, spacedim> normal;
1854 for (unsigned int d = 0; d < spacedim; d++)
1855 normal[d] =
1856 data.covariant[point][d] * quadrature.normal_vector(point);
1857
1858 output_data.JxW_values[point] = weights[point] * det * normal.norm();
1859
1860 if ((update_flags & update_normal_vectors) != 0u)
1861 {
1862 normal /= normal.norm();
1863 output_data.normal_vectors[point] = normal;
1864 }
1865 }
1866
1867 // copy values from InternalData to vector given by reference
1868 if (update_flags & update_jacobians)
1869 {
1870 AssertDimension(output_data.jacobians.size(), n_q_points);
1871 for (unsigned int point = 0; point < n_q_points; ++point)
1872 output_data.jacobians[point] = data.contravariant[point];
1873 }
1874
1875 // copy values from InternalData to vector given by reference
1876 if (update_flags & update_inverse_jacobians)
1877 {
1878 AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1879 for (unsigned int point = 0; point < n_q_points; ++point)
1880 output_data.inverse_jacobians[point] =
1881 data.covariant[point].transpose();
1882 }
1883
1884 // calculate derivatives of the Jacobians
1885 internal::MappingFEFieldImplementation::
1886 maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1887 QProjector<dim>::DataSetDescriptor::cell(),
1888 data,
1889 euler_dof_handler->get_fe(),
1890 fe_mask,
1891 fe_to_real,
1892 output_data.jacobian_grads);
1893
1894 // calculate derivatives of the Jacobians pushed forward to real cell
1895 // coordinates
1896 internal::MappingFEFieldImplementation::
1897 maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1898 QProjector<dim>::DataSetDescriptor::cell(),
1899 data,
1900 euler_dof_handler->get_fe(),
1901 fe_mask,
1902 fe_to_real,
1903 output_data.jacobian_pushed_forward_grads);
1904
1905 // calculate hessians of the Jacobians
1906 internal::MappingFEFieldImplementation::
1907 maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1908 QProjector<dim>::DataSetDescriptor::cell(),
1909 data,
1910 euler_dof_handler->get_fe(),
1911 fe_mask,
1912 fe_to_real,
1913 output_data.jacobian_2nd_derivatives);
1914
1915 // calculate hessians of the Jacobians pushed forward to real cell
1916 // coordinates
1917 internal::MappingFEFieldImplementation::
1918 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1919 spacedim,
1920 VectorType>(
1921 QProjector<dim>::DataSetDescriptor::cell(),
1922 data,
1923 euler_dof_handler->get_fe(),
1924 fe_mask,
1925 fe_to_real,
1926 output_data.jacobian_pushed_forward_2nd_derivatives);
1927
1928 // calculate gradients of the hessians of the Jacobians
1929 internal::MappingFEFieldImplementation::
1930 maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1931 QProjector<dim>::DataSetDescriptor::cell(),
1932 data,
1933 euler_dof_handler->get_fe(),
1934 fe_mask,
1935 fe_to_real,
1936 output_data.jacobian_3rd_derivatives);
1937
1938 // calculate gradients of the hessians of the Jacobians pushed forward to
1939 // real cell coordinates
1940 internal::MappingFEFieldImplementation::
1941 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1942 spacedim,
1943 VectorType>(
1944 QProjector<dim>::DataSetDescriptor::cell(),
1945 data,
1946 euler_dof_handler->get_fe(),
1947 fe_mask,
1948 fe_to_real,
1949 output_data.jacobian_pushed_forward_3rd_derivatives);
1950 }
1951}
1952
1953namespace internal
1954{
1955 namespace MappingFEFieldImplementation
1956 {
1957 namespace
1958 {
1959 template <int dim, int spacedim, int rank, typename VectorType>
1960 void
1961 transform_fields(
1962 const ArrayView<const Tensor<rank, dim>> &input,
1963 const MappingKind mapping_kind,
1964 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1965 const ArrayView<Tensor<rank, spacedim>> &output)
1966 {
1967 AssertDimension(input.size(), output.size());
1968 Assert((dynamic_cast<
1969 const typename ::
1970 MappingFEField<dim, spacedim, VectorType>::InternalData *>(
1971 &mapping_data) != nullptr),
1972 ExcInternalError());
1973 const typename ::MappingFEField<dim, spacedim, VectorType>::
1974 InternalData &data = static_cast<
1975 const typename ::MappingFEField<dim, spacedim, VectorType>::
1976 InternalData &>(mapping_data);
1977
1978 switch (mapping_kind)
1979 {
1980 case mapping_contravariant:
1981 {
1982 Assert(
1983 data.update_each & update_contravariant_transformation,
1984 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
1985 "update_contravariant_transformation"));
1986
1987 for (unsigned int i = 0; i < output.size(); ++i)
1988 output[i] =
1989 apply_transformation(data.contravariant[i], input[i]);
1990
1991 return;
1992 }
1993
1994 case mapping_piola:
1995 {
1996 Assert(
1997 data.update_each & update_contravariant_transformation,
1998 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
1999 "update_contravariant_transformation"));
2000 Assert(
2001 data.update_each & update_volume_elements,
2002 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
2003 "update_volume_elements"));
2004 Assert(rank == 1, ExcMessage("Only for rank 1"));
2005 for (unsigned int i = 0; i < output.size(); ++i)
2006 {
2007 output[i] =
2008 apply_transformation(data.contravariant[i], input[i]);
2009 output[i] /= data.volume_elements[i];
2010 }
2011 return;
2012 }
2013
2014
2015 // We still allow this operation as in the
2016 // reference cell Derivatives are Tensor
2017 // rather than DerivativeForm
2018 case mapping_covariant:
2019 {
2020 Assert(
2021 data.update_each & update_contravariant_transformation,
2022 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
2023 "update_contravariant_transformation"));
2024
2025 for (unsigned int i = 0; i < output.size(); ++i)
2026 output[i] = apply_transformation(data.covariant[i], input[i]);
2027
2028 return;
2029 }
2030
2031 default:
2032 DEAL_II_NOT_IMPLEMENTED();
2033 }
2034 }
2035
2036
2037 template <int dim, int spacedim, int rank, typename VectorType>
2038 void
2039 transform_differential_forms(
2040 const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
2041 const MappingKind mapping_kind,
2042 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2043 const ArrayView<Tensor<rank + 1, spacedim>> &output)
2044 {
2045 AssertDimension(input.size(), output.size());
2046 Assert((dynamic_cast<
2047 const typename ::
2048 MappingFEField<dim, spacedim, VectorType>::InternalData *>(
2049 &mapping_data) != nullptr),
2050 ExcInternalError());
2051 const typename ::MappingFEField<dim, spacedim, VectorType>::
2052 InternalData &data = static_cast<
2053 const typename ::MappingFEField<dim, spacedim, VectorType>::
2054 InternalData &>(mapping_data);
2055
2056 switch (mapping_kind)
2057 {
2058 case mapping_covariant:
2059 {
2060 Assert(
2061 data.update_each & update_contravariant_transformation,
2062 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
2063 "update_contravariant_transformation"));
2064
2065 for (unsigned int i = 0; i < output.size(); ++i)
2066 output[i] = apply_transformation(data.covariant[i], input[i]);
2067
2068 return;
2069 }
2070 default:
2071 DEAL_II_NOT_IMPLEMENTED();
2072 }
2073 }
2074 } // namespace
2075 } // namespace MappingFEFieldImplementation
2076} // namespace internal
2077
2078
2079
2080template <int dim, int spacedim, typename VectorType>
2081void
2082MappingFEField<dim, spacedim, VectorType>::transform(
2083 const ArrayView<const Tensor<1, dim>> &input,
2084 const MappingKind mapping_kind,
2085 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2086 const ArrayView<Tensor<1, spacedim>> &output) const
2087{
2088 AssertDimension(input.size(), output.size());
2089
2090 internal::MappingFEFieldImplementation::
2091 transform_fields<dim, spacedim, 1, VectorType>(input,
2092 mapping_kind,
2093 mapping_data,
2094 output);
2095}
2096
2097
2098
2099template <int dim, int spacedim, typename VectorType>
2100void
2101MappingFEField<dim, spacedim, VectorType>::transform(
2102 const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
2103 const MappingKind mapping_kind,
2104 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2105 const ArrayView<Tensor<2, spacedim>> &output) const
2106{
2107 AssertDimension(input.size(), output.size());
2108
2109 internal::MappingFEFieldImplementation::
2110 transform_differential_forms<dim, spacedim, 1, VectorType>(input,
2111 mapping_kind,
2112 mapping_data,
2113 output);
2114}
2115
2116
2117
2118template <int dim, int spacedim, typename VectorType>
2119void
2120MappingFEField<dim, spacedim, VectorType>::transform(
2121 const ArrayView<const Tensor<2, dim>> &input,
2122 const MappingKind,
2123 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2124 const ArrayView<Tensor<2, spacedim>> &output) const
2125{
2126 (void)input;
2127 (void)output;
2128 (void)mapping_data;
2129 AssertDimension(input.size(), output.size());
2130
2131 AssertThrow(false, ExcNotImplemented());
2132}
2133
2134
2135
2136template <int dim, int spacedim, typename VectorType>
2137void
2138MappingFEField<dim, spacedim, VectorType>::transform(
2139 const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
2140 const MappingKind mapping_kind,
2141 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2142 const ArrayView<Tensor<3, spacedim>> &output) const
2143{
2144 AssertDimension(input.size(), output.size());
2145 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
2146 ExcInternalError());
2147 const InternalData &data = static_cast<const InternalData &>(mapping_data);
2148
2149 switch (mapping_kind)
2150 {
2151 case mapping_covariant_gradient:
2152 {
2153 Assert(data.update_each & update_contravariant_transformation,
2154 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
2155 "update_covariant_transformation"));
2156
2157 for (unsigned int q = 0; q < output.size(); ++q)
2158 for (unsigned int i = 0; i < spacedim; ++i)
2159 for (unsigned int j = 0; j < spacedim; ++j)
2160 for (unsigned int k = 0; k < spacedim; ++k)
2161 {
2162 output[q][i][j][k] = data.covariant[q][j][0] *
2163 data.covariant[q][k][0] *
2164 input[q][i][0][0];
2165 for (unsigned int J = 0; J < dim; ++J)
2166 {
2167 const unsigned int K0 = (0 == J) ? 1 : 0;
2168 for (unsigned int K = K0; K < dim; ++K)
2169 output[q][i][j][k] += data.covariant[q][j][J] *
2170 data.covariant[q][k][K] *
2171 input[q][i][J][K];
2172 }
2173 }
2174 return;
2175 }
2176
2177 default:
2178 DEAL_II_NOT_IMPLEMENTED();
2179 }
2180}
2181
2182
2183
2184template <int dim, int spacedim, typename VectorType>
2185void
2186MappingFEField<dim, spacedim, VectorType>::transform(
2187 const ArrayView<const Tensor<3, dim>> &input,
2188 const MappingKind /*mapping_kind*/,
2189 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2190 const ArrayView<Tensor<3, spacedim>> &output) const
2191{
2192 (void)input;
2193 (void)output;
2194 (void)mapping_data;
2195 AssertDimension(input.size(), output.size());
2196
2197 AssertThrow(false, ExcNotImplemented());
2198}
2199
2200
2201
2202template <int dim, int spacedim, typename VectorType>
2203Point<spacedim>
2204MappingFEField<dim, spacedim, VectorType>::transform_unit_to_real_cell(
2205 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2206 const Point<dim> &p) const
2207{
2208 // Use the get_data function to create an InternalData with data vectors of
2209 // the right size and transformation shape values already computed at point
2210 // p.
2211 const Quadrature<dim> point_quadrature(p);
2212 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2213 get_data(update_quadrature_points | update_jacobians, point_quadrature));
2214 Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2215 ExcInternalError());
2216
2217 update_internal_dofs(cell, static_cast<InternalData &>(*mdata));
2218
2219 return do_transform_unit_to_real_cell(static_cast<InternalData &>(*mdata));
2220}
2221
2222
2223template <int dim, int spacedim, typename VectorType>
2224Point<spacedim>
2225MappingFEField<dim, spacedim, VectorType>::do_transform_unit_to_real_cell(
2226 const InternalData &data) const
2227{
2228 Point<spacedim> p_real;
2229
2230 for (unsigned int i = 0; i < data.n_shape_functions; ++i)
2231 {
2232 unsigned int comp_i =
2233 euler_dof_handler->get_fe().system_to_component_index(i).first;
2234 if (fe_mask[comp_i])
2235 p_real[fe_to_real[comp_i]] +=
2236 data.local_dof_values[i] * data.shape(0, i);
2237 }
2238
2239 return p_real;
2240}
2241
2242
2243
2244template <int dim, int spacedim, typename VectorType>
2245Point<dim>
2246MappingFEField<dim, spacedim, VectorType>::transform_real_to_unit_cell(
2247 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2248 const Point<spacedim> &p) const
2249{
2250 // first a Newton iteration based on the real mapping. It uses the center
2251 // point of the cell as a starting point
2252 Point<dim> initial_p_unit;
2253 try
2254 {
2255 initial_p_unit = get_default_linear_mapping(cell->get_triangulation())
2256 .transform_real_to_unit_cell(cell, p);
2257 }
2258 catch (const typename Mapping<dim, spacedim>::ExcTransformationFailed &)
2259 {
2260 // mirror the conditions of the code below to determine if we need to
2261 // use an arbitrary starting point or if we just need to rethrow the
2262 // exception
2263 for (unsigned int d = 0; d < dim; ++d)
2264 initial_p_unit[d] = 0.5;
2265 }
2266
2267 initial_p_unit = cell->reference_cell().closest_point(initial_p_unit);
2268
2269 UpdateFlags update_flags = update_quadrature_points | update_jacobians;
2270 if (spacedim > dim)
2271 update_flags |= update_jacobian_grads;
2272 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2273 get_data(update_flags, Quadrature<dim>(initial_p_unit)));
2274 Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2275 ExcInternalError());
2276
2277 update_internal_dofs(cell, static_cast<InternalData &>(*mdata));
2278
2279 return do_transform_real_to_unit_cell(cell,
2280 p,
2281 initial_p_unit,
2282 static_cast<InternalData &>(*mdata));
2283}
2284
2285
2286template <int dim, int spacedim, typename VectorType>
2287Point<dim>
2288MappingFEField<dim, spacedim, VectorType>::do_transform_real_to_unit_cell(
2289 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2290 const Point<spacedim> &p,
2291 const Point<dim> &starting_guess,
2292 InternalData &mdata) const
2293{
2294 const unsigned int n_shapes = mdata.shape_values.size();
2295 (void)n_shapes;
2296 Assert(n_shapes != 0, ExcInternalError());
2297 AssertDimension(mdata.shape_derivatives.size(), n_shapes);
2298
2299
2300 // Newton iteration to solve
2301 // f(x)=p(x)-p=0
2302 // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
2303 // The start value was set to be the
2304 // linear approximation to the cell
2305 // The shape values and derivatives
2306 // of the mapping at this point are
2307 // previously computed.
2308
2309 Point<dim> p_unit = starting_guess;
2310 Point<dim> f;
2311 mdata.reinit(mdata.update_each, Quadrature<dim>(starting_guess));
2312
2313 Point<spacedim> p_real(do_transform_unit_to_real_cell(mdata));
2314 Tensor<1, spacedim> p_minus_F = p - p_real;
2315 const double eps = 1.e-12 * cell->diameter();
2316 const unsigned int newton_iteration_limit = 20;
2317 unsigned int newton_iteration = 0;
2318 while (p_minus_F.norm_square() > eps * eps)
2319 {
2320 // f'(x)
2321 Point<spacedim> DF[dim];
2322 Tensor<2, dim> df;
2323 for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
2324 {
2325 const Tensor<1, dim> &grad_k = mdata.derivative(0, k);
2326 const unsigned int comp_k =
2327 euler_dof_handler->get_fe().system_to_component_index(k).first;
2328 if (fe_mask[comp_k])
2329 for (unsigned int j = 0; j < dim; ++j)
2330 DF[j][fe_to_real[comp_k]] +=
2331 mdata.local_dof_values[k] * grad_k[j];
2332 }
2333 for (unsigned int j = 0; j < dim; ++j)
2334 {
2335 f[j] = DF[j] * p_minus_F;
2336 for (unsigned int l = 0; l < dim; ++l)
2337 df[j][l] = -DF[j] * DF[l];
2338 }
2339 // Solve [f'(x)]d=f(x)
2340 const Tensor<1, dim> delta =
2341 invert(df) * static_cast<const Tensor<1, dim> &>(f);
2342 // do a line search
2343 double step_length = 1;
2344 do
2345 {
2346 // update of p_unit. The
2347 // spacedimth component of
2348 // transformed point is simply
2349 // ignored in codimension one
2350 // case. When this component is
2351 // not zero, then we are
2352 // projecting the point to the
2353 // surface or curve identified
2354 // by the cell.
2355 Point<dim> p_unit_trial = p_unit;
2356 for (unsigned int i = 0; i < dim; ++i)
2357 p_unit_trial[i] -= step_length * delta[i];
2358 // shape values and derivatives
2359 // at new p_unit point
2360 mdata.reinit(mdata.update_each, Quadrature<dim>(p_unit_trial));
2361 // f(x)
2362 const Point<spacedim> p_real_trial =
2363 do_transform_unit_to_real_cell(mdata);
2364 const Tensor<1, spacedim> f_trial = p - p_real_trial;
2365 // see if we are making progress with the current step length
2366 // and if not, reduce it by a factor of two and try again
2367 if (f_trial.norm() < p_minus_F.norm())
2368 {
2369 p_real = p_real_trial;
2370 p_unit = p_unit_trial;
2371 p_minus_F = f_trial;
2372 break;
2373 }
2374 else if (step_length > 0.05)
2375 step_length /= 2;
2376 else
2377 goto failure;
2378 }
2379 while (true);
2380 ++newton_iteration;
2381 if (newton_iteration > newton_iteration_limit)
2382 goto failure;
2383 }
2384 return p_unit;
2385 // if we get to the following label, then we have either run out
2386 // of Newton iterations, or the line search has not converged.
2387 // in either case, we need to give up, so throw an exception that
2388 // can then be caught
2389failure:
2390 AssertThrow(false,
2391 (typename Mapping<dim, spacedim>::ExcTransformationFailed()));
2392 // ...the compiler wants us to return something, though we can
2393 // of course never get here...
2394 return {};
2395}
2396
2397
2398template <int dim, int spacedim, typename VectorType>
2399unsigned int
2400MappingFEField<dim, spacedim, VectorType>::get_degree() const
2401{
2402 return euler_dof_handler->get_fe().degree;
2403}
2404
2405
2406
2407template <int dim, int spacedim, typename VectorType>
2408ComponentMask
2413
2414
2415template <int dim, int spacedim, typename VectorType>
2416std::unique_ptr<Mapping<dim, spacedim>>
2417MappingFEField<dim, spacedim, VectorType>::clone() const
2418{
2419 return std::make_unique<MappingFEField<dim, spacedim, VectorType>>(*this);
2420}
2421
2422
2423template <int dim, int spacedim, typename VectorType>
2424void
2425MappingFEField<dim, spacedim, VectorType>::update_internal_dofs(
2426 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2427 const typename MappingFEField<dim, spacedim, VectorType>::InternalData &data)
2428 const
2429{
2430 Assert(euler_dof_handler != nullptr,
2431 ExcMessage("euler_dof_handler is empty"));
2432
2433 typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(*cell,
2434 euler_dof_handler);
2435 Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
2436 if (uses_level_dofs)
2437 {
2438 AssertIndexRange(cell->level(), euler_vector.size());
2439 AssertDimension(euler_vector[cell->level()]->size(),
2440 euler_dof_handler->n_dofs(cell->level()));
2441 }
2442 else
2443 AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
2444
2445 if (uses_level_dofs)
2446 dof_cell->get_mg_dof_indices(data.local_dof_indices);
2447 else
2448 dof_cell->get_dof_indices(data.local_dof_indices);
2449
2450 const VectorType &vector =
2451 uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
2452
2453 for (unsigned int i = 0; i < data.local_dof_values.size(); ++i)
2454 data.local_dof_values[i] =
2455 internal::ElementAccess<VectorType>::get(vector,
2456 data.local_dof_indices[i]);
2457}
2458
2459// explicit instantiations
2460#define SPLIT_INSTANTIATIONS_COUNT 2
2461#ifndef SPLIT_INSTANTIATIONS_INDEX
2462# define SPLIT_INSTANTIATIONS_INDEX 0
2463#endif
2464#include "mapping_fe_field.inst"
2465
2466
2467DEAL_II_NAMESPACE_CLOSE
array_view.h
make_array_view
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
block_vector.h
ComponentMask::size
unsigned int size() const
DerivativeForm::transpose
DerivativeForm< 1, spacedim, dim, Number > transpose() const
MGLevelObject::resize
void resize(const unsigned int new_minlevel, const unsigned int new_maxlevel, Args &&...args)
Definition mg_level_object.h:256
MappingFEField::InternalData::third_derivative
const Tensor< 3, dim > & third_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
Definition mapping_fe_field.cc:226
MappingFEField::InternalData::InternalData
InternalData(const FiniteElement< dim, spacedim > &fe, const ComponentMask &mask)
Definition mapping_fe_field.cc:59
MappingFEField::InternalData::second_derivative
const Tensor< 2, dim > & second_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
Definition mapping_fe_field.cc:201
MappingFEField::InternalData::reinit
virtual void reinit(const UpdateFlags update_flags, const Quadrature< dim > &quadrature) override
Definition mapping_fe_field.cc:74
MappingFEField::InternalData::memory_consumption
virtual std::size_t memory_consumption() const override
Definition mapping_fe_field.cc:154
MappingFEField::InternalData::fourth_derivative
const Tensor< 4, dim > & fourth_derivative(const unsigned int qpoint, const unsigned int shape_nr) const
Definition mapping_fe_field.cc:251
MappingFEField::InternalData::shape_derivatives
std::vector< Tensor< 1, dim > > shape_derivatives
Definition mapping_fe_field.h:389
MappingFEField::InternalData::derivative
const Tensor< 1, dim > & derivative(const unsigned int qpoint, const unsigned int shape_nr) const
Definition mapping_fe_field.cc:176
MappingFEField::InternalData::local_dof_values
std::vector< double > local_dof_values
Definition mapping_fe_field.h:498
MappingFEField::InternalData::shape
const double & shape(const unsigned int qpoint, const unsigned int shape_nr) const
Definition mapping_fe_field.cc:412
MappingFEField::InternalData::shape_values
std::vector< double > shape_values
Definition mapping_fe_field.h:382
MappingFEField::fe_to_real
std::vector< unsigned int > fe_to_real
Definition mapping_fe_field.h:652
MappingFEField::get_subface_data
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
Definition mapping_fe_field.cc:665
MappingFEField::fe_mask
const ComponentMask fe_mask
Definition mapping_fe_field.h:641
MappingFEField::update_internal_dofs
void update_internal_dofs(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const typename MappingFEField< dim, spacedim, VectorType >::InternalData &data) const
Definition mapping_fe_field.cc:2425
MappingFEField::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 override
Definition mapping_fe_field.cc:2082
MappingFEField::transform_unit_to_real_cell
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const override
Definition mapping_fe_field.cc:2204
MappingFEField::preserves_vertex_locations
virtual bool preserves_vertex_locations() const override
Definition mapping_fe_field.cc:424
MappingFEField::do_transform_real_to_unit_cell
Point< dim > do_transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p, const Point< dim > &starting_guess, InternalData &mdata) const
Definition mapping_fe_field.cc:2288
MappingFEField::compute_face_data
void compute_face_data(const unsigned int n_original_q_points, InternalData &data) const
Definition mapping_fe_field.cc:577
MappingFEField::transform_real_to_unit_cell
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const override
Definition mapping_fe_field.cc:2246
MappingFEField::get_face_data
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const hp::QCollection< dim - 1 > &quadrature) const override
Definition mapping_fe_field.cc:644
MappingFEField::fill_fe_immersed_surface_values
virtual void fill_fe_immersed_surface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const NonMatching::ImmersedSurfaceQuadrature< dim > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
Definition mapping_fe_field.cc:1792
MappingFEField::MappingFEField
friend class MappingFEField
Definition mapping_fe_field.h:671
MappingFEField::get_data
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
Definition mapping_fe_field.cc:629
MappingFEField::is_compatible_with
virtual bool is_compatible_with(const ReferenceCell &reference_cell) const override
Definition mapping_fe_field.cc:433
MappingFEField::euler_vector
std::vector< ObserverPointer< const VectorType, MappingFEField< dim, spacedim, VectorType > > > euler_vector
Definition mapping_fe_field.h:580
MappingFEField::fill_fe_values
virtual CellSimilarity::Similarity fill_fe_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const Quadrature< dim > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
Definition mapping_fe_field.cc:1507
MappingFEField::reference_cell
ReferenceCell reference_cell
Definition mapping_fe_field.h:567
MappingFEField::get_degree
unsigned int get_degree() const
Definition mapping_fe_field.cc:2400
MappingFEField::requires_update_flags
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
Definition mapping_fe_field.cc:523
MappingFEField::do_transform_unit_to_real_cell
Point< spacedim > do_transform_unit_to_real_cell(const InternalData &mdata) const
Definition mapping_fe_field.cc:2225
MappingFEField::fe_values_mutex
Threads::Mutex fe_values_mutex
Definition mapping_fe_field.h:663
MappingFEField::uses_level_dofs
const bool uses_level_dofs
Definition mapping_fe_field.h:573
MappingFEField::get_component_mask
ComponentMask get_component_mask() const
Definition mapping_fe_field.cc:2409
MappingFEField::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 typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
Definition mapping_fe_field.cc:1711
MappingFEField::euler_dof_handler
ObserverPointer< const DoFHandler< dim, spacedim >, MappingFEField< dim, spacedim, VectorType > > euler_dof_handler
Definition mapping_fe_field.h:587
MappingFEField::fe_values
FEValues< dim, spacedim > fe_values
Definition mapping_fe_field.h:658
MappingFEField::clone
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Definition mapping_fe_field.cc:2417
MappingFEField::get_vertices
virtual boost::container::small_vector< Point< spacedim >, ReferenceCells::max_n_vertices< dim >() > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
Definition mapping_fe_field.cc:456
MappingFEField::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 subface_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
Definition mapping_fe_field.cc:1750
Mapping
Abstract base class for mapping classes.
Definition mapping.h:318
Mapping::transform_real_to_unit_cell
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
NonMatching::ImmersedSurfaceQuadrature::normal_vector
const Tensor< 1, spacedim > & normal_vector(const unsigned int i) const
Definition immersed_surface_quadrature.cc:69
Point
Definition point.h:111
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)
Definition qprojector.cc:1541
QProjector::DataSetDescriptor::cell
static DataSetDescriptor cell()
Definition qprojector.h:461
Quadrature::get_weights
const std::vector< double > & get_weights() const
Quadrature::get_points
const std::vector< Point< dim > > & get_points() const
Quadrature::size
unsigned int size() const
ReferenceCell::n_faces
unsigned int n_faces() const
Definition reference_cell.h:1863
ReferenceCell::get_dimension
unsigned int get_dimension() const
Definition reference_cell.h:1635
Tensor::norm
numbers::NumberTraits< Number >::real_type norm() const
Tensor::norm_square
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
hp::Collection::size
unsigned int size() const
Definition collection.h:315
internal::FEValuesImplementation::MappingRelatedData::inverse_jacobians
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
Definition mapping_related_data.h:97
internal::FEValuesImplementation::MappingRelatedData::jacobian_pushed_forward_3rd_derivatives
std::vector< Tensor< 5, spacedim > > jacobian_pushed_forward_3rd_derivatives
Definition mapping_related_data.h:127
internal::FEValuesImplementation::MappingRelatedData::jacobian_3rd_derivatives
std::vector< DerivativeForm< 4, dim, spacedim > > jacobian_3rd_derivatives
Definition mapping_related_data.h:121
internal::FEValuesImplementation::MappingRelatedData::jacobian_2nd_derivatives
std::vector< DerivativeForm< 3, dim, spacedim > > jacobian_2nd_derivatives
Definition mapping_related_data.h:109
internal::FEValuesImplementation::MappingRelatedData::jacobian_pushed_forward_2nd_derivatives
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
Definition mapping_related_data.h:115
internal::FEValuesImplementation::MappingRelatedData::jacobian_pushed_forward_grads
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
Definition mapping_related_data.h:103
internal::FEValuesImplementation::MappingRelatedData::jacobian_grads
std::vector< DerivativeForm< 2, dim, spacedim > > jacobian_grads
Definition mapping_related_data.h:92
internal::FEValuesImplementation::MappingRelatedData::jacobians
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
Definition mapping_related_data.h:86
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
apply_transformation
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
Definition derivative_form.h:478
dof_accessor.h
fe_values.h
fe_q.h
fe_system.h
fe_tools.h
full_matrix.h
second
Point< 2 > second
Definition grid_out.cc:4630
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1659
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:2026
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:2115
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
MappingFEField::ExcInactiveCell
static ::ExceptionBase & ExcInactiveCell()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition exceptions.h:1780
DoFHandler::cell_iterator
typename ActiveSelector::cell_iterator cell_iterator
Definition dof_handler.h:475
UpdateFlags
UpdateFlags
Definition fe_update_flags.h:64
update_jacobian_pushed_forward_2nd_derivatives
@ update_jacobian_pushed_forward_2nd_derivatives
Definition fe_update_flags.h:203
update_volume_elements
@ update_volume_elements
Determinant of the Jacobian.
Definition fe_update_flags.h:189
update_contravariant_transformation
@ update_contravariant_transformation
Contravariant transformation.
Definition fe_update_flags.h:172
update_jacobian_pushed_forward_grads
@ update_jacobian_pushed_forward_grads
Definition fe_update_flags.h:194
update_jacobian_3rd_derivatives
@ update_jacobian_3rd_derivatives
Definition fe_update_flags.h:207
update_values
@ update_values
Shape function values.
Definition fe_update_flags.h:75
update_jacobian_grads
@ update_jacobian_grads
Gradient of volume element.
Definition fe_update_flags.h:152
update_normal_vectors
@ update_normal_vectors
Normal vectors.
Definition fe_update_flags.h:141
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition fe_update_flags.h:134
update_covariant_transformation
@ update_covariant_transformation
Covariant transformation.
Definition fe_update_flags.h:165
update_jacobians
@ update_jacobians
Volume element.
Definition fe_update_flags.h:147
update_inverse_jacobians
@ update_inverse_jacobians
Volume element.
Definition fe_update_flags.h:158
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition fe_update_flags.h:127
update_jacobian_pushed_forward_3rd_derivatives
@ update_jacobian_pushed_forward_3rd_derivatives
Definition fe_update_flags.h:212
update_boundary_forms
@ update_boundary_forms
Outer normal vector, not normalized.
Definition fe_update_flags.h:100
update_jacobian_2nd_derivatives
@ update_jacobian_2nd_derivatives
Definition fe_update_flags.h:198
MappingKind
MappingKind
Definition mapping.h:79
mapping_piola
@ mapping_piola
Definition mapping.h:114
mapping_covariant_gradient
@ mapping_covariant_gradient
Definition mapping.h:100
mapping_covariant
@ mapping_covariant
Definition mapping.h:89
mapping_contravariant
@ mapping_contravariant
Definition mapping.h:94
DEAL_II_NOT_IMPLEMENTED
#define DEAL_II_NOT_IMPLEMENTED()
Definition exceptions.h:1855
la_parallel_block_vector.h
la_parallel_vector.h
get_default_linear_mapping
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Definition mapping.cc:316
mapping_fe_field.h
memory_consumption.h
data
std::vector< index_type > data
Definition mpi.cc:735
size
std::size_t size
Definition mpi.cc:734
Utilities::to_string
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:479
internal::MappingQImplementation::transform_differential_forms
void transform_differential_forms(const ArrayView< const DerivativeForm< rank, dim, spacedim > > &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank+1, spacedim > > &output)
Definition mapping_q_internal.h:2509
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition types.h:232
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6869
petsc_block_vector.h
petsc_vector.h
quadrature_lib.h
GeometryInfo::subface_ratio
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
internal::ElementAccess::get
static VectorType::value_type get(const VectorType &V, const types::global_dof_index i)
Definition vector_element_access.h:72
determinant
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
Definition symmetric_tensor.h:2639
invert
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
Definition symmetric_tensor.h:3164
tensor_product_polynomials.h
tria_iterator.h
trilinos_epetra_vector.h
trilinos_parallel_block_vector.h
trilinos_tpetra_vector.h
trilinos_vector.h
vector_tools.h
Generated by doxygen 1.9.8