deal.II version GIT relicensing-2167-g9622207b8f 2024-11-21 12:40: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 GeometryInfo<dim>::vertices_per_cell>
451MappingFEField<dim, spacedim, VectorType>::get_vertices(
452 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
453{
454 // we transform our tria iterator into a dof iterator so we can access
455 // data not associated with triangulations
456 const typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(
457 *cell, euler_dof_handler);
458
459 Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
460 AssertDimension(cell->n_vertices(), fe_values.n_quadrature_points);
461 AssertDimension(fe_to_real.size(),
462 euler_dof_handler->get_fe().n_components());
463 if (uses_level_dofs)
464 {
465 AssertIndexRange(cell->level(), euler_vector.size());
466 AssertDimension(euler_vector[cell->level()]->size(),
467 euler_dof_handler->n_dofs(cell->level()));
468 }
469 else
470 AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
471
472 {
473 std::lock_guard<std::mutex> lock(fe_values_mutex);
474 fe_values.reinit(dof_cell);
475 }
476 const unsigned int dofs_per_cell =
477 euler_dof_handler->get_fe().n_dofs_per_cell();
478 std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
479 if (uses_level_dofs)
480 dof_cell->get_mg_dof_indices(dof_indices);
481 else
482 dof_cell->get_dof_indices(dof_indices);
483
484 const VectorType &vector =
485 uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
486
487 boost::container::small_vector<Point<spacedim>,
488 GeometryInfo<dim>::vertices_per_cell>
489 vertices(cell->n_vertices());
490 for (unsigned int i = 0; i < dofs_per_cell; ++i)
491 {
492 const unsigned int comp = fe_to_real
493 [euler_dof_handler->get_fe().system_to_component_index(i).first];
494 if (comp != numbers::invalid_unsigned_int)
495 {
496 typename VectorType::value_type value =
497 internal::ElementAccess<VectorType>::get(vector, dof_indices[i]);
498 if (euler_dof_handler->get_fe().is_primitive(i))
499 for (const unsigned int v : cell->vertex_indices())
500 vertices[v][comp] += fe_values.shape_value(i, v) * value;
501 else
502 DEAL_II_NOT_IMPLEMENTED();
503 }
504 }
505
506 return vertices;
507}
508
509
510
511template <int dim, int spacedim, typename VectorType>
512UpdateFlags
513MappingFEField<dim, spacedim, VectorType>::requires_update_flags(
514 const UpdateFlags in) const
515{
516 // add flags if the respective quantities are necessary to compute
517 // what we need. note that some flags appear in both conditions and
518 // in subsequent set operations. this leads to some circular
519 // logic. the only way to treat this is to iterate. since there are
520 // 5 if-clauses in the loop, it will take at most 4 iterations to
521 // converge. do them:
522 UpdateFlags out = in;
523 for (unsigned int i = 0; i < 5; ++i)
524 {
525 // The following is a little incorrect:
526 // If not applied on a face,
527 // update_boundary_forms does not
528 // make sense. On the other hand,
529 // it is necessary on a
530 // face. Currently,
531 // update_boundary_forms is simply
532 // ignored for the interior of a
533 // cell.
534 if (out & (update_JxW_values | update_normal_vectors))
535 out |= update_boundary_forms;
536
537 if (out &
538 (update_covariant_transformation | update_jacobian_grads |
539 update_jacobians | update_boundary_forms | update_normal_vectors))
540 out |= update_contravariant_transformation;
541
542 if (out &
543 (update_inverse_jacobians | update_jacobian_pushed_forward_grads |
544 update_jacobian_pushed_forward_2nd_derivatives |
545 update_jacobian_pushed_forward_3rd_derivatives))
546 out |= update_covariant_transformation;
547
548 // The contravariant transformation is used in the Piola
549 // transformation, which requires the determinant of the Jacobi
550 // matrix of the transformation. Because we have no way of
551 // knowing here whether the finite element wants to use the
552 // contravariant or the Piola transforms, we add the volume elements
553 // to the list of flags to be updated for each cell.
554 if (out & update_contravariant_transformation)
555 out |= update_volume_elements;
556
557 if (out & update_normal_vectors)
558 out |= update_volume_elements;
559 }
560
561 return out;
562}
563
564
565template <int dim, int spacedim, typename VectorType>
566void
567MappingFEField<dim, spacedim, VectorType>::compute_face_data(
568 const unsigned int n_original_q_points,
569 InternalData &data) const
570{
571 // Set to the size of a single quadrature object for faces, as the size set
572 // in in reinit() is for all points
573 if (data.update_each & update_covariant_transformation)
574 data.covariant.resize(n_original_q_points);
575
576 if (data.update_each & update_contravariant_transformation)
577 data.contravariant.resize(n_original_q_points);
578
579 if (data.update_each & update_volume_elements)
580 data.volume_elements.resize(n_original_q_points);
581
582 if (dim > 1)
583 {
584 if (data.update_each & update_boundary_forms)
585 {
586 data.aux.resize(
587 dim - 1, std::vector<Tensor<1, spacedim>>(n_original_q_points));
588
589
590 // TODO: only a single reference cell type possible...
591 const auto n_faces = reference_cell.n_faces();
592
593 // Compute tangentials to the unit cell.
594 for (unsigned int i = 0; i < n_faces; ++i)
595 {
596 data.unit_tangentials[i].resize(n_original_q_points);
597 std::fill(
598 data.unit_tangentials[i].begin(),
599 data.unit_tangentials[i].end(),
600 reference_cell.template unit_tangential_vectors<dim>(i, 0));
601 if (dim > 2)
602 {
603 data.unit_tangentials[n_faces + i].resize(
604 n_original_q_points);
605 std::fill(
606 data.unit_tangentials[n_faces + i].begin(),
607 data.unit_tangentials[n_faces + i].end(),
608 reference_cell.template unit_tangential_vectors<dim>(i, 1));
609 }
610 }
611 }
612 }
613}
614
615
616
617template <int dim, int spacedim, typename VectorType>
618typename std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
619MappingFEField<dim, spacedim, VectorType>::get_data(
620 const UpdateFlags update_flags,
621 const Quadrature<dim> &quadrature) const
622{
623 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
624 std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
625 data_ptr->reinit(requires_update_flags(update_flags), quadrature);
626
627 return data_ptr;
628}
629
630
631
632template <int dim, int spacedim, typename VectorType>
633std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
634MappingFEField<dim, spacedim, VectorType>::get_face_data(
635 const UpdateFlags update_flags,
636 const hp::QCollection<dim - 1> &quadrature) const
637{
638 AssertDimension(quadrature.size(), 1);
639
640 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
641 std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
642 auto &data = dynamic_cast<InternalData &>(*data_ptr);
643
644 const Quadrature<dim> q(
645 QProjector<dim>::project_to_all_faces(reference_cell, quadrature[0]));
646 data.reinit(requires_update_flags(update_flags), q);
647 this->compute_face_data(quadrature[0].size(), data);
648
649 return data_ptr;
650}
651
652
653template <int dim, int spacedim, typename VectorType>
654std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
655MappingFEField<dim, spacedim, VectorType>::get_subface_data(
656 const UpdateFlags update_flags,
657 const Quadrature<dim - 1> &quadrature) const
658{
659 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
660 std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
661 auto &data = dynamic_cast<InternalData &>(*data_ptr);
662
663 const Quadrature<dim> q(
664 QProjector<dim>::project_to_all_subfaces(reference_cell, quadrature));
665 data.reinit(requires_update_flags(update_flags), q);
666 this->compute_face_data(quadrature.size(), data);
667
668 return data_ptr;
669}
670
671
672
673namespace internal
674{
675 namespace MappingFEFieldImplementation
676 {
677 namespace
678 {
685 template <int dim, int spacedim, typename VectorType>
686 void
687 maybe_compute_q_points(
688 const typename ::QProjector<dim>::DataSetDescriptor data_set,
689 const typename ::MappingFEField<dim, spacedim, VectorType>::
690 InternalData &data,
691 const FiniteElement<dim, spacedim> &fe,
692 const ComponentMask &fe_mask,
693 const std::vector<unsigned int> &fe_to_real,
694 std::vector<Point<spacedim>> &quadrature_points)
695 {
696 const UpdateFlags update_flags = data.update_each;
697
698 if (update_flags & update_quadrature_points)
699 {
700 for (unsigned int point = 0; point < quadrature_points.size();
701 ++point)
702 {
703 Point<spacedim> result;
704 const double *shape = &data.shape(point + data_set, 0);
705
706 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
707 {
708 const unsigned int comp_k =
709 fe.system_to_component_index(k).first;
710 if (fe_mask[comp_k])
711 result[fe_to_real[comp_k]] +=
712 data.local_dof_values[k] * shape[k];
713 }
714
715 quadrature_points[point] = result;
716 }
717 }
718 }
719
727 template <int dim, int spacedim, typename VectorType>
728 void
729 maybe_update_Jacobians(
730 const typename ::QProjector<dim>::DataSetDescriptor data_set,
731 const typename ::MappingFEField<dim, spacedim, VectorType>::
732 InternalData &data,
733 const FiniteElement<dim, spacedim> &fe,
734 const ComponentMask &fe_mask,
735 const std::vector<unsigned int> &fe_to_real)
736 {
737 const UpdateFlags update_flags = data.update_each;
738
739 // then Jacobians
740 if (update_flags & update_contravariant_transformation)
741 {
742 const unsigned int n_q_points = data.contravariant.size();
743
744 Assert(data.n_shape_functions > 0, ExcInternalError());
745
746 for (unsigned int point = 0; point < n_q_points; ++point)
747 {
748 const Tensor<1, dim> *data_derv =
749 &data.derivative(point + data_set, 0);
750
751 Tensor<1, dim> result[spacedim];
752
753 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
754 {
755 const unsigned int comp_k =
756 fe.system_to_component_index(k).first;
757 if (fe_mask[comp_k])
758 result[fe_to_real[comp_k]] +=
759 data.local_dof_values[k] * data_derv[k];
760 }
761
762 // write result into contravariant data
763 for (unsigned int i = 0; i < spacedim; ++i)
764 {
765 data.contravariant[point][i] = result[i];
766 }
767 }
768 }
769
770 if (update_flags & update_covariant_transformation)
771 {
772 AssertDimension(data.covariant.size(), data.contravariant.size());
773 for (unsigned int point = 0; point < data.contravariant.size();
774 ++point)
775 data.covariant[point] =
776 (data.contravariant[point]).covariant_form();
777 }
778
779 if (update_flags & update_volume_elements)
780 {
781 AssertDimension(data.contravariant.size(),
782 data.volume_elements.size());
783 for (unsigned int point = 0; point < data.contravariant.size();
784 ++point)
785 data.volume_elements[point] =
786 data.contravariant[point].determinant();
787 }
788 }
789
796 template <int dim, int spacedim, typename VectorType>
797 void
798 maybe_update_jacobian_grads(
799 const typename ::QProjector<dim>::DataSetDescriptor data_set,
800 const typename ::MappingFEField<dim, spacedim, VectorType>::
801 InternalData &data,
802 const FiniteElement<dim, spacedim> &fe,
803 const ComponentMask &fe_mask,
804 const std::vector<unsigned int> &fe_to_real,
805 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
806 {
807 const UpdateFlags update_flags = data.update_each;
808 if (update_flags & update_jacobian_grads)
809 {
810 const unsigned int n_q_points = jacobian_grads.size();
811
812 for (unsigned int point = 0; point < n_q_points; ++point)
813 {
814 const Tensor<2, dim> *second =
815 &data.second_derivative(point + data_set, 0);
816
817 DerivativeForm<2, dim, spacedim> result;
818
819 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
820 {
821 const unsigned int comp_k =
822 fe.system_to_component_index(k).first;
823 if (fe_mask[comp_k])
824 for (unsigned int j = 0; j < dim; ++j)
825 for (unsigned int l = 0; l < dim; ++l)
826 result[fe_to_real[comp_k]][j][l] +=
827 (second[k][j][l] * data.local_dof_values[k]);
828 }
829
830 // never touch any data for j=dim in case dim<spacedim, so
831 // it will always be zero as it was initialized
832 for (unsigned int i = 0; i < spacedim; ++i)
833 for (unsigned int j = 0; j < dim; ++j)
834 for (unsigned int l = 0; l < dim; ++l)
835 jacobian_grads[point][i][j][l] = result[i][j][l];
836 }
837 }
838 }
839
846 template <int dim, int spacedim, typename VectorType>
847 void
848 maybe_update_jacobian_pushed_forward_grads(
849 const typename ::QProjector<dim>::DataSetDescriptor data_set,
850 const typename ::MappingFEField<dim, spacedim, VectorType>::
851 InternalData &data,
852 const FiniteElement<dim, spacedim> &fe,
853 const ComponentMask &fe_mask,
854 const std::vector<unsigned int> &fe_to_real,
855 std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads)
856 {
857 const UpdateFlags update_flags = data.update_each;
858 if (update_flags & update_jacobian_pushed_forward_grads)
859 {
860 const unsigned int n_q_points =
861 jacobian_pushed_forward_grads.size();
862
863 double tmp[spacedim][spacedim][spacedim];
864 for (unsigned int point = 0; point < n_q_points; ++point)
865 {
866 const Tensor<2, dim> *second =
867 &data.second_derivative(point + data_set, 0);
868
869 DerivativeForm<2, dim, spacedim> result;
870
871 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
872 {
873 const unsigned int comp_k =
874 fe.system_to_component_index(k).first;
875 if (fe_mask[comp_k])
876 for (unsigned int j = 0; j < dim; ++j)
877 for (unsigned int l = 0; l < dim; ++l)
878 result[fe_to_real[comp_k]][j][l] +=
879 (second[k][j][l] * data.local_dof_values[k]);
880 }
881
882 // first push forward the j-components
883 for (unsigned int i = 0; i < spacedim; ++i)
884 for (unsigned int j = 0; j < spacedim; ++j)
885 for (unsigned int l = 0; l < dim; ++l)
886 {
887 tmp[i][j][l] =
888 result[i][0][l] * data.covariant[point][j][0];
889 for (unsigned int jr = 1; jr < dim; ++jr)
890 {
891 tmp[i][j][l] +=
892 result[i][jr][l] * data.covariant[point][j][jr];
893 }
894 }
895
896 // now, pushing forward the l-components
897 for (unsigned int i = 0; i < spacedim; ++i)
898 for (unsigned int j = 0; j < spacedim; ++j)
899 for (unsigned int l = 0; l < spacedim; ++l)
900 {
901 jacobian_pushed_forward_grads[point][i][j][l] =
902 tmp[i][j][0] * data.covariant[point][l][0];
903 for (unsigned int lr = 1; lr < dim; ++lr)
904 {
905 jacobian_pushed_forward_grads[point][i][j][l] +=
906 tmp[i][j][lr] * data.covariant[point][l][lr];
907 }
908 }
909 }
910 }
911 }
912
919 template <int dim, int spacedim, typename VectorType>
920 void
921 maybe_update_jacobian_2nd_derivatives(
922 const typename ::QProjector<dim>::DataSetDescriptor data_set,
923 const typename ::MappingFEField<dim, spacedim, VectorType>::
924 InternalData &data,
925 const FiniteElement<dim, spacedim> &fe,
926 const ComponentMask &fe_mask,
927 const std::vector<unsigned int> &fe_to_real,
928 std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
929 {
930 const UpdateFlags update_flags = data.update_each;
931 if (update_flags & update_jacobian_2nd_derivatives)
932 {
933 const unsigned int n_q_points = jacobian_2nd_derivatives.size();
934
935 for (unsigned int point = 0; point < n_q_points; ++point)
936 {
937 const Tensor<3, dim> *third =
938 &data.third_derivative(point + data_set, 0);
939
940 DerivativeForm<3, dim, spacedim> result;
941
942 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
943 {
944 const unsigned int comp_k =
945 fe.system_to_component_index(k).first;
946 if (fe_mask[comp_k])
947 for (unsigned int j = 0; j < dim; ++j)
948 for (unsigned int l = 0; l < dim; ++l)
949 for (unsigned int m = 0; m < dim; ++m)
950 result[fe_to_real[comp_k]][j][l][m] +=
951 (third[k][j][l][m] * data.local_dof_values[k]);
952 }
953
954 // never touch any data for j=dim in case dim<spacedim, so
955 // it will always be zero as it was initialized
956 for (unsigned int i = 0; i < spacedim; ++i)
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 jacobian_2nd_derivatives[point][i][j][l][m] =
961 result[i][j][l][m];
962 }
963 }
964 }
965
973 template <int dim, int spacedim, typename VectorType>
974 void
975 maybe_update_jacobian_pushed_forward_2nd_derivatives(
976 const typename ::QProjector<dim>::DataSetDescriptor data_set,
977 const typename ::MappingFEField<dim, spacedim, VectorType>::
978 InternalData &data,
979 const FiniteElement<dim, spacedim> &fe,
980 const ComponentMask &fe_mask,
981 const std::vector<unsigned int> &fe_to_real,
982 std::vector<Tensor<4, spacedim>>
983 &jacobian_pushed_forward_2nd_derivatives)
984 {
985 const UpdateFlags update_flags = data.update_each;
986 if (update_flags & update_jacobian_pushed_forward_2nd_derivatives)
987 {
988 const unsigned int n_q_points =
989 jacobian_pushed_forward_2nd_derivatives.size();
990
991 double tmp[spacedim][spacedim][spacedim][spacedim];
992 for (unsigned int point = 0; point < n_q_points; ++point)
993 {
994 const Tensor<3, dim> *third =
995 &data.third_derivative(point + data_set, 0);
996
997 DerivativeForm<3, dim, spacedim> result;
998
999 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1000 {
1001 const unsigned int comp_k =
1002 fe.system_to_component_index(k).first;
1003 if (fe_mask[comp_k])
1004 for (unsigned int j = 0; j < dim; ++j)
1005 for (unsigned int l = 0; l < dim; ++l)
1006 for (unsigned int m = 0; m < dim; ++m)
1007 result[fe_to_real[comp_k]][j][l][m] +=
1008 (third[k][j][l][m] * data.local_dof_values[k]);
1009 }
1010
1011 // push forward the j-coordinate
1012 for (unsigned int i = 0; i < spacedim; ++i)
1013 for (unsigned int j = 0; j < spacedim; ++j)
1014 for (unsigned int l = 0; l < dim; ++l)
1015 for (unsigned int m = 0; m < dim; ++m)
1016 {
1017 jacobian_pushed_forward_2nd_derivatives
1018 [point][i][j][l][m] =
1019 result[i][0][l][m] * data.covariant[point][j][0];
1020 for (unsigned int jr = 1; jr < dim; ++jr)
1021 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1022 [l][m] +=
1023 result[i][jr][l][m] *
1024 data.covariant[point][j][jr];
1025 }
1026
1027 // push forward the l-coordinate
1028 for (unsigned int i = 0; i < spacedim; ++i)
1029 for (unsigned int j = 0; j < spacedim; ++j)
1030 for (unsigned int l = 0; l < spacedim; ++l)
1031 for (unsigned int m = 0; m < dim; ++m)
1032 {
1033 tmp[i][j][l][m] =
1034 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1035 [0][m] *
1036 data.covariant[point][l][0];
1037 for (unsigned int lr = 1; lr < dim; ++lr)
1038 tmp[i][j][l][m] +=
1039 jacobian_pushed_forward_2nd_derivatives[point][i]
1040 [j][lr]
1041 [m] *
1042 data.covariant[point][l][lr];
1043 }
1044
1045 // push forward the m-coordinate
1046 for (unsigned int i = 0; i < spacedim; ++i)
1047 for (unsigned int j = 0; j < spacedim; ++j)
1048 for (unsigned int l = 0; l < spacedim; ++l)
1049 for (unsigned int m = 0; m < spacedim; ++m)
1050 {
1051 jacobian_pushed_forward_2nd_derivatives
1052 [point][i][j][l][m] =
1053 tmp[i][j][l][0] * data.covariant[point][m][0];
1054 for (unsigned int mr = 1; mr < dim; ++mr)
1055 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1056 [l][m] +=
1057 tmp[i][j][l][mr] * data.covariant[point][m][mr];
1058 }
1059 }
1060 }
1061 }
1062
1069 template <int dim, int spacedim, typename VectorType>
1070 void
1071 maybe_update_jacobian_3rd_derivatives(
1072 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1073 const typename ::MappingFEField<dim, spacedim, VectorType>::
1074 InternalData &data,
1075 const FiniteElement<dim, spacedim> &fe,
1076 const ComponentMask &fe_mask,
1077 const std::vector<unsigned int> &fe_to_real,
1078 std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1079 {
1080 const UpdateFlags update_flags = data.update_each;
1081 if (update_flags & update_jacobian_3rd_derivatives)
1082 {
1083 const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1084
1085 for (unsigned int point = 0; point < n_q_points; ++point)
1086 {
1087 const Tensor<4, dim> *fourth =
1088 &data.fourth_derivative(point + data_set, 0);
1089
1090 DerivativeForm<4, dim, spacedim> result;
1091
1092 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1093 {
1094 const unsigned int comp_k =
1095 fe.system_to_component_index(k).first;
1096 if (fe_mask[comp_k])
1097 for (unsigned int j = 0; j < dim; ++j)
1098 for (unsigned int l = 0; l < dim; ++l)
1099 for (unsigned int m = 0; m < dim; ++m)
1100 for (unsigned int n = 0; n < dim; ++n)
1101 result[fe_to_real[comp_k]][j][l][m][n] +=
1102 (fourth[k][j][l][m][n] *
1103 data.local_dof_values[k]);
1104 }
1105
1106 // never touch any data for j,l,m,n=dim in case
1107 // dim<spacedim, so it will always be zero as it was
1108 // initialized
1109 for (unsigned int i = 0; i < spacedim; ++i)
1110 for (unsigned int j = 0; j < dim; ++j)
1111 for (unsigned int l = 0; l < dim; ++l)
1112 for (unsigned int m = 0; m < dim; ++m)
1113 for (unsigned int n = 0; n < dim; ++n)
1114 jacobian_3rd_derivatives[point][i][j][l][m][n] =
1115 result[i][j][l][m][n];
1116 }
1117 }
1118 }
1119
1127 template <int dim, int spacedim, typename VectorType>
1128 void
1129 maybe_update_jacobian_pushed_forward_3rd_derivatives(
1130 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1131 const typename ::MappingFEField<dim, spacedim, VectorType>::
1132 InternalData &data,
1133 const FiniteElement<dim, spacedim> &fe,
1134 const ComponentMask &fe_mask,
1135 const std::vector<unsigned int> &fe_to_real,
1136 std::vector<Tensor<5, spacedim>>
1137 &jacobian_pushed_forward_3rd_derivatives)
1138 {
1139 const UpdateFlags update_flags = data.update_each;
1140 if (update_flags & update_jacobian_pushed_forward_3rd_derivatives)
1141 {
1142 const unsigned int n_q_points =
1143 jacobian_pushed_forward_3rd_derivatives.size();
1144
1145 double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
1146 for (unsigned int point = 0; point < n_q_points; ++point)
1147 {
1148 const Tensor<4, dim> *fourth =
1149 &data.fourth_derivative(point + data_set, 0);
1150
1151 DerivativeForm<4, dim, spacedim> result;
1152
1153 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1154 {
1155 const unsigned int comp_k =
1156 fe.system_to_component_index(k).first;
1157 if (fe_mask[comp_k])
1158 for (unsigned int j = 0; j < dim; ++j)
1159 for (unsigned int l = 0; l < dim; ++l)
1160 for (unsigned int m = 0; m < dim; ++m)
1161 for (unsigned int n = 0; n < dim; ++n)
1162 result[fe_to_real[comp_k]][j][l][m][n] +=
1163 (fourth[k][j][l][m][n] *
1164 data.local_dof_values[k]);
1165 }
1166
1167 // push-forward the j-coordinate
1168 for (unsigned int i = 0; i < spacedim; ++i)
1169 for (unsigned int j = 0; j < spacedim; ++j)
1170 for (unsigned int l = 0; l < dim; ++l)
1171 for (unsigned int m = 0; m < dim; ++m)
1172 for (unsigned int n = 0; n < dim; ++n)
1173 {
1174 tmp[i][j][l][m][n] = result[i][0][l][m][n] *
1175 data.covariant[point][j][0];
1176 for (unsigned int jr = 1; jr < dim; ++jr)
1177 tmp[i][j][l][m][n] +=
1178 result[i][jr][l][m][n] *
1179 data.covariant[point][j][jr];
1180 }
1181
1182 // push-forward the l-coordinate
1183 for (unsigned int i = 0; i < spacedim; ++i)
1184 for (unsigned int j = 0; j < spacedim; ++j)
1185 for (unsigned int l = 0; l < spacedim; ++l)
1186 for (unsigned int m = 0; m < dim; ++m)
1187 for (unsigned int n = 0; n < dim; ++n)
1188 {
1189 jacobian_pushed_forward_3rd_derivatives
1190 [point][i][j][l][m][n] =
1191 tmp[i][j][0][m][n] *
1192 data.covariant[point][l][0];
1193 for (unsigned int lr = 1; lr < dim; ++lr)
1194 jacobian_pushed_forward_3rd_derivatives[point][i]
1195 [j][l][m]
1196 [n] +=
1197 tmp[i][j][lr][m][n] *
1198 data.covariant[point][l][lr];
1199 }
1200
1201 // push-forward the m-coordinate
1202 for (unsigned int i = 0; i < spacedim; ++i)
1203 for (unsigned int j = 0; j < spacedim; ++j)
1204 for (unsigned int l = 0; l < spacedim; ++l)
1205 for (unsigned int m = 0; m < spacedim; ++m)
1206 for (unsigned int n = 0; n < dim; ++n)
1207 {
1208 tmp[i][j][l][m][n] =
1209 jacobian_pushed_forward_3rd_derivatives[point][i]
1210 [j][l][0]
1211 [n] *
1212 data.covariant[point][m][0];
1213 for (unsigned int mr = 1; mr < dim; ++mr)
1214 tmp[i][j][l][m][n] +=
1215 jacobian_pushed_forward_3rd_derivatives[point]
1216 [i][j][l]
1217 [mr][n] *
1218 data.covariant[point][m][mr];
1219 }
1220
1221 // push-forward the n-coordinate
1222 for (unsigned int i = 0; i < spacedim; ++i)
1223 for (unsigned int j = 0; j < spacedim; ++j)
1224 for (unsigned int l = 0; l < spacedim; ++l)
1225 for (unsigned int m = 0; m < spacedim; ++m)
1226 for (unsigned int n = 0; n < spacedim; ++n)
1227 {
1228 jacobian_pushed_forward_3rd_derivatives
1229 [point][i][j][l][m][n] =
1230 tmp[i][j][l][m][0] *
1231 data.covariant[point][n][0];
1232 for (unsigned int nr = 1; nr < dim; ++nr)
1233 jacobian_pushed_forward_3rd_derivatives[point][i]
1234 [j][l][m]
1235 [n] +=
1236 tmp[i][j][l][m][nr] *
1237 data.covariant[point][n][nr];
1238 }
1239 }
1240 }
1241 }
1242
1243
1253 template <int dim, int spacedim, typename VectorType>
1254 void
1255 maybe_compute_face_data(
1256 const ::Mapping<dim, spacedim> &mapping,
1257 const typename ::Triangulation<dim, spacedim>::cell_iterator
1258 &cell,
1259 const unsigned int face_no,
1260 const unsigned int subface_no,
1261 const typename QProjector<dim>::DataSetDescriptor data_set,
1262 const typename ::MappingFEField<dim, spacedim, VectorType>::
1263 InternalData &data,
1264 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1265 &output_data)
1266 {
1267 const UpdateFlags update_flags = data.update_each;
1268
1269 if (update_flags & update_boundary_forms)
1270 {
1271 const unsigned int n_q_points = output_data.boundary_forms.size();
1272 if (update_flags & update_normal_vectors)
1273 AssertDimension(output_data.normal_vectors.size(), n_q_points);
1274 if (update_flags & update_JxW_values)
1275 AssertDimension(output_data.JxW_values.size(), n_q_points);
1276
1277 // map the unit tangentials to the real cell. checking for d!=dim-1
1278 // eliminates compiler warnings regarding unsigned int expressions <
1279 // 0.
1280 for (unsigned int d = 0; d != dim - 1; ++d)
1281 {
1282 Assert(face_no + cell->n_faces() * d <
1283 data.unit_tangentials.size(),
1284 ExcInternalError());
1285 Assert(
1286 data.aux[d].size() <=
1287 data.unit_tangentials[face_no + cell->n_faces() * d].size(),
1288 ExcInternalError());
1289
1290 mapping.transform(
1291 make_array_view(
1292 data.unit_tangentials[face_no + cell->n_faces() * d]),
1293 mapping_contravariant,
1294 data,
1295 make_array_view(data.aux[d]));
1296 }
1297
1298 // if dim==spacedim, we can use the unit tangentials to compute the
1299 // boundary form by simply taking the cross product
1300 if (dim == spacedim)
1301 {
1302 for (unsigned int i = 0; i < n_q_points; ++i)
1303 switch (dim)
1304 {
1305 case 1:
1306 // in 1d, we don't have access to any of the data.aux
1307 // fields (because it has only dim-1 components), but we
1308 // can still compute the boundary form by simply looking
1309 // at the number of the face
1310 output_data.boundary_forms[i][0] =
1311 (face_no == 0 ? -1 : +1);
1312 break;
1313 case 2:
1314 output_data.boundary_forms[i] =
1315 cross_product_2d(data.aux[0][i]);
1316 break;
1317 case 3:
1318 output_data.boundary_forms[i] =
1319 cross_product_3d(data.aux[0][i], data.aux[1][i]);
1320 break;
1321 default:
1322 DEAL_II_NOT_IMPLEMENTED();
1323 }
1324 }
1325 else //(dim < spacedim)
1326 {
1327 // in the codim-one case, the boundary form results from the
1328 // cross product of all the face tangential vectors and the cell
1329 // normal vector
1330 //
1331 // to compute the cell normal, use the same method used in
1332 // fill_fe_values for cells above
1333 AssertDimension(data.contravariant.size(), n_q_points);
1334
1335 for (unsigned int point = 0; point < n_q_points; ++point)
1336 {
1337 if (dim == 1)
1338 {
1339 // J is a tangent vector
1340 output_data.boundary_forms[point] =
1341 data.contravariant[point].transpose()[0];
1342 output_data.boundary_forms[point] /=
1343 (face_no == 0 ? -1. : +1.) *
1344 output_data.boundary_forms[point].norm();
1345 }
1346
1347 if (dim == 2)
1348 {
1349 const DerivativeForm<1, spacedim, dim> DX_t =
1350 data.contravariant[point].transpose();
1351
1352 Tensor<1, spacedim> cell_normal =
1353 cross_product_3d(DX_t[0], DX_t[1]);
1354 cell_normal /= cell_normal.norm();
1355
1356 // then compute the face normal from the face tangent
1357 // and the cell normal:
1358 output_data.boundary_forms[point] =
1359 cross_product_3d(data.aux[0][point], cell_normal);
1360 }
1361 }
1362 }
1363
1364 if (update_flags & (update_normal_vectors | update_JxW_values))
1365 for (unsigned int i = 0; i < output_data.boundary_forms.size();
1366 ++i)
1367 {
1368 if (update_flags & update_JxW_values)
1369 {
1370 output_data.JxW_values[i] =
1371 output_data.boundary_forms[i].norm() *
1372 data.quadrature_weights[i + data_set];
1373
1374 if (subface_no != numbers::invalid_unsigned_int)
1375 {
1376 // TODO
1377 const double area_ratio =
1378 GeometryInfo<dim>::subface_ratio(
1379 cell->subface_case(face_no), subface_no);
1380 output_data.JxW_values[i] *= area_ratio;
1381 }
1382 }
1383
1384 if (update_flags & update_normal_vectors)
1385 output_data.normal_vectors[i] =
1386 Point<spacedim>(output_data.boundary_forms[i] /
1387 output_data.boundary_forms[i].norm());
1388 }
1389 }
1390 }
1391
1398 template <int dim, int spacedim, typename VectorType>
1399 void
1400 do_fill_fe_face_values(
1401 const ::Mapping<dim, spacedim> &mapping,
1402 const typename ::Triangulation<dim, spacedim>::cell_iterator
1403 &cell,
1404 const unsigned int face_no,
1405 const unsigned int subface_no,
1406 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1407 const typename ::MappingFEField<dim, spacedim, VectorType>::
1408 InternalData &data,
1409 const FiniteElement<dim, spacedim> &fe,
1410 const ComponentMask &fe_mask,
1411 const std::vector<unsigned int> &fe_to_real,
1412 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1413 &output_data)
1414 {
1415 maybe_compute_q_points<dim, spacedim, VectorType>(
1416 data_set,
1417 data,
1418 fe,
1419 fe_mask,
1420 fe_to_real,
1421 output_data.quadrature_points);
1422
1423 maybe_update_Jacobians<dim, spacedim, VectorType>(
1424 data_set, data, fe, fe_mask, fe_to_real);
1425
1426 const UpdateFlags update_flags = data.update_each;
1427 const unsigned int n_q_points = data.contravariant.size();
1428
1429 if (update_flags & update_jacobians)
1430 for (unsigned int point = 0; point < n_q_points; ++point)
1431 output_data.jacobians[point] = data.contravariant[point];
1432
1433 if (update_flags & update_inverse_jacobians)
1434 for (unsigned int point = 0; point < n_q_points; ++point)
1435 output_data.inverse_jacobians[point] =
1436 data.covariant[point].transpose();
1437
1438 maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1439 data_set, data, fe, fe_mask, fe_to_real, output_data.jacobian_grads);
1440
1441 maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1442 data_set,
1443 data,
1444 fe,
1445 fe_mask,
1446 fe_to_real,
1447 output_data.jacobian_pushed_forward_grads);
1448
1449 maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1450 data_set,
1451 data,
1452 fe,
1453 fe_mask,
1454 fe_to_real,
1455 output_data.jacobian_2nd_derivatives);
1456
1457 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1458 spacedim,
1459 VectorType>(
1460 data_set,
1461 data,
1462 fe,
1463 fe_mask,
1464 fe_to_real,
1465 output_data.jacobian_pushed_forward_2nd_derivatives);
1466
1467 maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1468 data_set,
1469 data,
1470 fe,
1471 fe_mask,
1472 fe_to_real,
1473 output_data.jacobian_3rd_derivatives);
1474
1475 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1476 spacedim,
1477 VectorType>(
1478 data_set,
1479 data,
1480 fe,
1481 fe_mask,
1482 fe_to_real,
1483 output_data.jacobian_pushed_forward_3rd_derivatives);
1484
1485 maybe_compute_face_data<dim, spacedim, VectorType>(
1486 mapping, cell, face_no, subface_no, data_set, data, output_data);
1487 }
1488 } // namespace
1489 } // namespace MappingFEFieldImplementation
1490} // namespace internal
1491
1492
1493// Note that the CellSimilarity flag is modifiable, since MappingFEField can
1494// need to recalculate data even when cells are similar.
1495template <int dim, int spacedim, typename VectorType>
1496CellSimilarity::Similarity
1497MappingFEField<dim, spacedim, VectorType>::fill_fe_values(
1498 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1499 const CellSimilarity::Similarity,
1500 const Quadrature<dim> &quadrature,
1501 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1502 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1503 &output_data) const
1504{
1505 // convert data object to internal data for this class. fails with an
1506 // exception if that is not possible
1507 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1508 ExcInternalError());
1509 const InternalData &data = static_cast<const InternalData &>(internal_data);
1510
1511 const unsigned int n_q_points = quadrature.size();
1512
1513 update_internal_dofs(cell, data);
1514
1515 internal::MappingFEFieldImplementation::
1516 maybe_compute_q_points<dim, spacedim, VectorType>(
1517 QProjector<dim>::DataSetDescriptor::cell(),
1518 data,
1519 euler_dof_handler->get_fe(),
1520 fe_mask,
1521 fe_to_real,
1522 output_data.quadrature_points);
1523
1524 internal::MappingFEFieldImplementation::
1525 maybe_update_Jacobians<dim, spacedim, VectorType>(
1526 QProjector<dim>::DataSetDescriptor::cell(),
1527 data,
1528 euler_dof_handler->get_fe(),
1529 fe_mask,
1530 fe_to_real);
1531
1532 const UpdateFlags update_flags = data.update_each;
1533 const std::vector<double> &weights = quadrature.get_weights();
1534
1535 // Multiply quadrature weights by absolute value of Jacobian determinants or
1536 // the area element g=sqrt(DX^t DX) in case of codim > 0
1537
1538 if (update_flags & (update_normal_vectors | update_JxW_values))
1539 {
1540 AssertDimension(output_data.JxW_values.size(), n_q_points);
1541
1542 Assert(!(update_flags & update_normal_vectors) ||
1543 (output_data.normal_vectors.size() == n_q_points),
1544 ExcDimensionMismatch(output_data.normal_vectors.size(),
1545 n_q_points));
1546
1547
1548 for (unsigned int point = 0; point < n_q_points; ++point)
1549 {
1550 if (dim == spacedim)
1551 {
1552 const double det = data.volume_elements[point];
1553
1554 // check for distorted cells.
1555
1556 // TODO: this allows for anisotropies of up to 1e6 in 3d and
1557 // 1e12 in 2d. might want to find a finer
1558 // (dimension-independent) criterion
1559 Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1560 cell->diameter() / std::sqrt(double(dim))),
1561 (typename Mapping<dim, spacedim>::ExcDistortedMappedCell(
1562 cell->center(), det, point)));
1563 output_data.JxW_values[point] = weights[point] * det;
1564 }
1565 // if dim==spacedim, then there is no cell normal to
1566 // compute. since this is for FEValues (and not FEFaceValues),
1567 // there are also no face normals to compute
1568 else // codim>0 case
1569 {
1570 Tensor<1, spacedim> DX_t[dim];
1571 for (unsigned int i = 0; i < spacedim; ++i)
1572 for (unsigned int j = 0; j < dim; ++j)
1573 DX_t[j][i] = data.contravariant[point][i][j];
1574
1575 Tensor<2, dim> G; // First fundamental form
1576 for (unsigned int i = 0; i < dim; ++i)
1577 for (unsigned int j = 0; j < dim; ++j)
1578 G[i][j] = DX_t[i] * DX_t[j];
1579
1580 output_data.JxW_values[point] =
1581 std::sqrt(determinant(G)) * weights[point];
1582
1583 if (update_flags & update_normal_vectors)
1584 {
1585 Assert(spacedim - dim == 1,
1586 ExcMessage("There is no cell normal in codim 2."));
1587
1588 if (dim == 1)
1589 output_data.normal_vectors[point] =
1590 cross_product_2d(-DX_t[0]);
1591 else
1592 {
1593 Assert(dim == 2, ExcInternalError());
1594
1595 // dim-1==1 for the second argument, but this
1596 // avoids a compiler warning about array bounds:
1597 output_data.normal_vectors[point] =
1598 cross_product_3d(DX_t[0], DX_t[dim - 1]);
1599 }
1600
1601 output_data.normal_vectors[point] /=
1602 output_data.normal_vectors[point].norm();
1603
1604 if (cell->direction_flag() == false)
1605 output_data.normal_vectors[point] *= -1.;
1606 }
1607 } // codim>0 case
1608 }
1609 }
1610
1611 // copy values from InternalData to vector given by reference
1612 if (update_flags & update_jacobians)
1613 {
1614 AssertDimension(output_data.jacobians.size(), n_q_points);
1615 for (unsigned int point = 0; point < n_q_points; ++point)
1616 output_data.jacobians[point] = data.contravariant[point];
1617 }
1618
1619 // copy values from InternalData to vector given by reference
1620 if (update_flags & update_inverse_jacobians)
1621 {
1622 AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1623 for (unsigned int point = 0; point < n_q_points; ++point)
1624 output_data.inverse_jacobians[point] =
1625 data.covariant[point].transpose();
1626 }
1627
1628 // calculate derivatives of the Jacobians
1629 internal::MappingFEFieldImplementation::
1630 maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1631 QProjector<dim>::DataSetDescriptor::cell(),
1632 data,
1633 euler_dof_handler->get_fe(),
1634 fe_mask,
1635 fe_to_real,
1636 output_data.jacobian_grads);
1637
1638 // calculate derivatives of the Jacobians pushed forward to real cell
1639 // coordinates
1640 internal::MappingFEFieldImplementation::
1641 maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1642 QProjector<dim>::DataSetDescriptor::cell(),
1643 data,
1644 euler_dof_handler->get_fe(),
1645 fe_mask,
1646 fe_to_real,
1647 output_data.jacobian_pushed_forward_grads);
1648
1649 // calculate hessians of the Jacobians
1650 internal::MappingFEFieldImplementation::
1651 maybe_update_jacobian_2nd_derivatives<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_2nd_derivatives);
1658
1659 // calculate hessians of the Jacobians pushed forward to real cell coordinates
1660 internal::MappingFEFieldImplementation::
1661 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1662 spacedim,
1663 VectorType>(
1664 QProjector<dim>::DataSetDescriptor::cell(),
1665 data,
1666 euler_dof_handler->get_fe(),
1667 fe_mask,
1668 fe_to_real,
1669 output_data.jacobian_pushed_forward_2nd_derivatives);
1670
1671 // calculate gradients of the hessians of the Jacobians
1672 internal::MappingFEFieldImplementation::
1673 maybe_update_jacobian_3rd_derivatives<dim, spacedim, 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_3rd_derivatives);
1680
1681 // calculate gradients of the hessians of the Jacobians pushed forward to real
1682 // cell coordinates
1683 internal::MappingFEFieldImplementation::
1684 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1685 spacedim,
1686 VectorType>(
1687 QProjector<dim>::DataSetDescriptor::cell(),
1688 data,
1689 euler_dof_handler->get_fe(),
1690 fe_mask,
1691 fe_to_real,
1692 output_data.jacobian_pushed_forward_3rd_derivatives);
1693
1694 return CellSimilarity::invalid_next_cell;
1695}
1696
1697
1698
1699template <int dim, int spacedim, typename VectorType>
1700void
1701MappingFEField<dim, spacedim, VectorType>::fill_fe_face_values(
1702 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1703 const unsigned int face_no,
1704 const hp::QCollection<dim - 1> &quadrature,
1705 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1706 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1707 &output_data) const
1708{
1709 AssertDimension(quadrature.size(), 1);
1710
1711 // convert data object to internal data for this class. fails with an
1712 // exception if that is not possible
1713 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1714 ExcInternalError());
1715 const InternalData &data = static_cast<const InternalData &>(internal_data);
1716
1717 update_internal_dofs(cell, data);
1718
1719 internal::MappingFEFieldImplementation::
1720 do_fill_fe_face_values<dim, spacedim, VectorType>(
1721 *this,
1722 cell,
1723 face_no,
1724 numbers::invalid_unsigned_int,
1725 QProjector<dim>::DataSetDescriptor::face(reference_cell,
1726 face_no,
1727 cell->combined_face_orientation(
1728 face_no),
1729 quadrature[0].size()),
1730 data,
1731 euler_dof_handler->get_fe(),
1732 fe_mask,
1733 fe_to_real,
1734 output_data);
1735}
1736
1737
1738template <int dim, int spacedim, typename VectorType>
1739void
1740MappingFEField<dim, spacedim, VectorType>::fill_fe_subface_values(
1741 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1742 const unsigned int face_no,
1743 const unsigned int subface_no,
1744 const Quadrature<dim - 1> &quadrature,
1745 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1746 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1747 &output_data) const
1748{
1749 // convert data object to internal data for this class. fails with an
1750 // exception if that is not possible
1751 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1752 ExcInternalError());
1753 const InternalData &data = static_cast<const InternalData &>(internal_data);
1754
1755 update_internal_dofs(cell, data);
1756
1757 internal::MappingFEFieldImplementation::do_fill_fe_face_values<dim,
1758 spacedim,
1759 VectorType>(
1760 *this,
1761 cell,
1762 face_no,
1763 numbers::invalid_unsigned_int,
1764 QProjector<dim>::DataSetDescriptor::subface(reference_cell,
1765 face_no,
1766 subface_no,
1767 cell->combined_face_orientation(
1768 face_no),
1769 quadrature.size(),
1770 cell->subface_case(face_no)),
1771 data,
1772 euler_dof_handler->get_fe(),
1773 fe_mask,
1774 fe_to_real,
1775 output_data);
1776}
1777
1778
1779
1780template <int dim, int spacedim, typename VectorType>
1781void
1782MappingFEField<dim, spacedim, VectorType>::fill_fe_immersed_surface_values(
1783 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1784 const NonMatching::ImmersedSurfaceQuadrature<dim> &quadrature,
1785 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1786 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1787 &output_data) const
1788{
1789 AssertDimension(dim, spacedim);
1790 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1791 ExcInternalError());
1792 const InternalData &data = static_cast<const InternalData &>(internal_data);
1793
1794 const unsigned int n_q_points = quadrature.size();
1795
1796 update_internal_dofs(cell, data);
1797
1798 internal::MappingFEFieldImplementation::
1799 maybe_compute_q_points<dim, spacedim, VectorType>(
1800 QProjector<dim>::DataSetDescriptor::cell(),
1801 data,
1802 euler_dof_handler->get_fe(),
1803 fe_mask,
1804 fe_to_real,
1805 output_data.quadrature_points);
1806
1807 internal::MappingFEFieldImplementation::
1808 maybe_update_Jacobians<dim, spacedim, VectorType>(
1809 QProjector<dim>::DataSetDescriptor::cell(),
1810 data,
1811 euler_dof_handler->get_fe(),
1812 fe_mask,
1813 fe_to_real);
1814
1815 const UpdateFlags update_flags = data.update_each;
1816 const std::vector<double> &weights = quadrature.get_weights();
1817
1818 if (update_flags & (update_normal_vectors | update_JxW_values))
1819 {
1820 AssertDimension(output_data.JxW_values.size(), n_q_points);
1821
1822 Assert(!(update_flags & update_normal_vectors) ||
1823 (output_data.normal_vectors.size() == n_q_points),
1824 ExcDimensionMismatch(output_data.normal_vectors.size(),
1825 n_q_points));
1826
1827
1828 for (unsigned int point = 0; point < n_q_points; ++point)
1829 {
1830 const double det = data.volume_elements[point];
1831
1832 // check for distorted cells.
1833
1834 // TODO: this allows for anisotropies of up to 1e6 in 3d and
1835 // 1e12 in 2d. might want to find a finer
1836 // (dimension-independent) criterion
1837 Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1838 cell->diameter() / std::sqrt(double(dim))),
1839 (typename Mapping<dim, spacedim>::ExcDistortedMappedCell(
1840 cell->center(), det, point)));
1841
1842 // The normals are n = J^{-T} * \hat{n} before normalizing.
1843 Tensor<1, spacedim> normal;
1844 for (unsigned int d = 0; d < spacedim; d++)
1845 normal[d] =
1846 data.covariant[point][d] * quadrature.normal_vector(point);
1847
1848 output_data.JxW_values[point] = weights[point] * det * normal.norm();
1849
1850 if ((update_flags & update_normal_vectors) != 0u)
1851 {
1852 normal /= normal.norm();
1853 output_data.normal_vectors[point] = normal;
1854 }
1855 }
1856
1857 // copy values from InternalData to vector given by reference
1858 if (update_flags & update_jacobians)
1859 {
1860 AssertDimension(output_data.jacobians.size(), n_q_points);
1861 for (unsigned int point = 0; point < n_q_points; ++point)
1862 output_data.jacobians[point] = data.contravariant[point];
1863 }
1864
1865 // copy values from InternalData to vector given by reference
1866 if (update_flags & update_inverse_jacobians)
1867 {
1868 AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1869 for (unsigned int point = 0; point < n_q_points; ++point)
1870 output_data.inverse_jacobians[point] =
1871 data.covariant[point].transpose();
1872 }
1873
1874 // calculate derivatives of the Jacobians
1875 internal::MappingFEFieldImplementation::
1876 maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1877 QProjector<dim>::DataSetDescriptor::cell(),
1878 data,
1879 euler_dof_handler->get_fe(),
1880 fe_mask,
1881 fe_to_real,
1882 output_data.jacobian_grads);
1883
1884 // calculate derivatives of the Jacobians pushed forward to real cell
1885 // coordinates
1886 internal::MappingFEFieldImplementation::
1887 maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1888 QProjector<dim>::DataSetDescriptor::cell(),
1889 data,
1890 euler_dof_handler->get_fe(),
1891 fe_mask,
1892 fe_to_real,
1893 output_data.jacobian_pushed_forward_grads);
1894
1895 // calculate hessians of the Jacobians
1896 internal::MappingFEFieldImplementation::
1897 maybe_update_jacobian_2nd_derivatives<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_2nd_derivatives);
1904
1905 // calculate hessians of the Jacobians pushed forward to real cell
1906 // coordinates
1907 internal::MappingFEFieldImplementation::
1908 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1909 spacedim,
1910 VectorType>(
1911 QProjector<dim>::DataSetDescriptor::cell(),
1912 data,
1913 euler_dof_handler->get_fe(),
1914 fe_mask,
1915 fe_to_real,
1916 output_data.jacobian_pushed_forward_2nd_derivatives);
1917
1918 // calculate gradients of the hessians of the Jacobians
1919 internal::MappingFEFieldImplementation::
1920 maybe_update_jacobian_3rd_derivatives<dim, spacedim, 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_3rd_derivatives);
1927
1928 // calculate gradients of the hessians of the Jacobians pushed forward to
1929 // real cell coordinates
1930 internal::MappingFEFieldImplementation::
1931 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1932 spacedim,
1933 VectorType>(
1934 QProjector<dim>::DataSetDescriptor::cell(),
1935 data,
1936 euler_dof_handler->get_fe(),
1937 fe_mask,
1938 fe_to_real,
1939 output_data.jacobian_pushed_forward_3rd_derivatives);
1940 }
1941}
1942
1943namespace internal
1944{
1945 namespace MappingFEFieldImplementation
1946 {
1947 namespace
1948 {
1949 template <int dim, int spacedim, int rank, typename VectorType>
1950 void
1951 transform_fields(
1952 const ArrayView<const Tensor<rank, dim>> &input,
1953 const MappingKind mapping_kind,
1954 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1955 const ArrayView<Tensor<rank, spacedim>> &output)
1956 {
1957 AssertDimension(input.size(), output.size());
1958 Assert((dynamic_cast<
1959 const typename ::
1960 MappingFEField<dim, spacedim, VectorType>::InternalData *>(
1961 &mapping_data) != nullptr),
1962 ExcInternalError());
1963 const typename ::MappingFEField<dim, spacedim, VectorType>::
1964 InternalData &data = static_cast<
1965 const typename ::MappingFEField<dim, spacedim, VectorType>::
1966 InternalData &>(mapping_data);
1967
1968 switch (mapping_kind)
1969 {
1970 case mapping_contravariant:
1971 {
1972 Assert(
1973 data.update_each & update_contravariant_transformation,
1974 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
1975 "update_contravariant_transformation"));
1976
1977 for (unsigned int i = 0; i < output.size(); ++i)
1978 output[i] =
1979 apply_transformation(data.contravariant[i], input[i]);
1980
1981 return;
1982 }
1983
1984 case mapping_piola:
1985 {
1986 Assert(
1987 data.update_each & update_contravariant_transformation,
1988 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
1989 "update_contravariant_transformation"));
1990 Assert(
1991 data.update_each & update_volume_elements,
1992 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
1993 "update_volume_elements"));
1994 Assert(rank == 1, ExcMessage("Only for rank 1"));
1995 for (unsigned int i = 0; i < output.size(); ++i)
1996 {
1997 output[i] =
1998 apply_transformation(data.contravariant[i], input[i]);
1999 output[i] /= data.volume_elements[i];
2000 }
2001 return;
2002 }
2003
2004
2005 // We still allow this operation as in the
2006 // reference cell Derivatives are Tensor
2007 // rather than DerivativeForm
2008 case mapping_covariant:
2009 {
2010 Assert(
2011 data.update_each & update_contravariant_transformation,
2012 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
2013 "update_contravariant_transformation"));
2014
2015 for (unsigned int i = 0; i < output.size(); ++i)
2016 output[i] = apply_transformation(data.covariant[i], input[i]);
2017
2018 return;
2019 }
2020
2021 default:
2022 DEAL_II_NOT_IMPLEMENTED();
2023 }
2024 }
2025
2026
2027 template <int dim, int spacedim, int rank, typename VectorType>
2028 void
2029 transform_differential_forms(
2030 const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
2031 const MappingKind mapping_kind,
2032 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2033 const ArrayView<Tensor<rank + 1, spacedim>> &output)
2034 {
2035 AssertDimension(input.size(), output.size());
2036 Assert((dynamic_cast<
2037 const typename ::
2038 MappingFEField<dim, spacedim, VectorType>::InternalData *>(
2039 &mapping_data) != nullptr),
2040 ExcInternalError());
2041 const typename ::MappingFEField<dim, spacedim, VectorType>::
2042 InternalData &data = static_cast<
2043 const typename ::MappingFEField<dim, spacedim, VectorType>::
2044 InternalData &>(mapping_data);
2045
2046 switch (mapping_kind)
2047 {
2048 case mapping_covariant:
2049 {
2050 Assert(
2051 data.update_each & update_contravariant_transformation,
2052 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
2053 "update_contravariant_transformation"));
2054
2055 for (unsigned int i = 0; i < output.size(); ++i)
2056 output[i] = apply_transformation(data.covariant[i], input[i]);
2057
2058 return;
2059 }
2060 default:
2061 DEAL_II_NOT_IMPLEMENTED();
2062 }
2063 }
2064 } // namespace
2065 } // namespace MappingFEFieldImplementation
2066} // namespace internal
2067
2068
2069
2070template <int dim, int spacedim, typename VectorType>
2071void
2072MappingFEField<dim, spacedim, VectorType>::transform(
2073 const ArrayView<const Tensor<1, dim>> &input,
2074 const MappingKind mapping_kind,
2075 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2076 const ArrayView<Tensor<1, spacedim>> &output) const
2077{
2078 AssertDimension(input.size(), output.size());
2079
2080 internal::MappingFEFieldImplementation::
2081 transform_fields<dim, spacedim, 1, VectorType>(input,
2082 mapping_kind,
2083 mapping_data,
2084 output);
2085}
2086
2087
2088
2089template <int dim, int spacedim, typename VectorType>
2090void
2091MappingFEField<dim, spacedim, VectorType>::transform(
2092 const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
2093 const MappingKind mapping_kind,
2094 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2095 const ArrayView<Tensor<2, spacedim>> &output) const
2096{
2097 AssertDimension(input.size(), output.size());
2098
2099 internal::MappingFEFieldImplementation::
2100 transform_differential_forms<dim, spacedim, 1, VectorType>(input,
2101 mapping_kind,
2102 mapping_data,
2103 output);
2104}
2105
2106
2107
2108template <int dim, int spacedim, typename VectorType>
2109void
2110MappingFEField<dim, spacedim, VectorType>::transform(
2111 const ArrayView<const Tensor<2, dim>> &input,
2112 const MappingKind,
2113 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2114 const ArrayView<Tensor<2, spacedim>> &output) const
2115{
2116 (void)input;
2117 (void)output;
2118 (void)mapping_data;
2119 AssertDimension(input.size(), output.size());
2120
2121 AssertThrow(false, ExcNotImplemented());
2122}
2123
2124
2125
2126template <int dim, int spacedim, typename VectorType>
2127void
2128MappingFEField<dim, spacedim, VectorType>::transform(
2129 const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
2130 const MappingKind mapping_kind,
2131 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2132 const ArrayView<Tensor<3, spacedim>> &output) const
2133{
2134 AssertDimension(input.size(), output.size());
2135 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
2136 ExcInternalError());
2137 const InternalData &data = static_cast<const InternalData &>(mapping_data);
2138
2139 switch (mapping_kind)
2140 {
2141 case mapping_covariant_gradient:
2142 {
2143 Assert(data.update_each & update_contravariant_transformation,
2144 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
2145 "update_covariant_transformation"));
2146
2147 for (unsigned int q = 0; q < output.size(); ++q)
2148 for (unsigned int i = 0; i < spacedim; ++i)
2149 for (unsigned int j = 0; j < spacedim; ++j)
2150 for (unsigned int k = 0; k < spacedim; ++k)
2151 {
2152 output[q][i][j][k] = data.covariant[q][j][0] *
2153 data.covariant[q][k][0] *
2154 input[q][i][0][0];
2155 for (unsigned int J = 0; J < dim; ++J)
2156 {
2157 const unsigned int K0 = (0 == J) ? 1 : 0;
2158 for (unsigned int K = K0; K < dim; ++K)
2159 output[q][i][j][k] += data.covariant[q][j][J] *
2160 data.covariant[q][k][K] *
2161 input[q][i][J][K];
2162 }
2163 }
2164 return;
2165 }
2166
2167 default:
2168 DEAL_II_NOT_IMPLEMENTED();
2169 }
2170}
2171
2172
2173
2174template <int dim, int spacedim, typename VectorType>
2175void
2176MappingFEField<dim, spacedim, VectorType>::transform(
2177 const ArrayView<const Tensor<3, dim>> &input,
2178 const MappingKind /*mapping_kind*/,
2179 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2180 const ArrayView<Tensor<3, spacedim>> &output) const
2181{
2182 (void)input;
2183 (void)output;
2184 (void)mapping_data;
2185 AssertDimension(input.size(), output.size());
2186
2187 AssertThrow(false, ExcNotImplemented());
2188}
2189
2190
2191
2192template <int dim, int spacedim, typename VectorType>
2193Point<spacedim>
2194MappingFEField<dim, spacedim, VectorType>::transform_unit_to_real_cell(
2195 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2196 const Point<dim> &p) const
2197{
2198 // Use the get_data function to create an InternalData with data vectors of
2199 // the right size and transformation shape values already computed at point
2200 // p.
2201 const Quadrature<dim> point_quadrature(p);
2202 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2203 get_data(update_quadrature_points | update_jacobians, point_quadrature));
2204 Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2205 ExcInternalError());
2206
2207 update_internal_dofs(cell, static_cast<InternalData &>(*mdata));
2208
2209 return do_transform_unit_to_real_cell(static_cast<InternalData &>(*mdata));
2210}
2211
2212
2213template <int dim, int spacedim, typename VectorType>
2214Point<spacedim>
2215MappingFEField<dim, spacedim, VectorType>::do_transform_unit_to_real_cell(
2216 const InternalData &data) const
2217{
2218 Point<spacedim> p_real;
2219
2220 for (unsigned int i = 0; i < data.n_shape_functions; ++i)
2221 {
2222 unsigned int comp_i =
2223 euler_dof_handler->get_fe().system_to_component_index(i).first;
2224 if (fe_mask[comp_i])
2225 p_real[fe_to_real[comp_i]] +=
2226 data.local_dof_values[i] * data.shape(0, i);
2227 }
2228
2229 return p_real;
2230}
2231
2232
2233
2234template <int dim, int spacedim, typename VectorType>
2235Point<dim>
2236MappingFEField<dim, spacedim, VectorType>::transform_real_to_unit_cell(
2237 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2238 const Point<spacedim> &p) const
2239{
2240 // first a Newton iteration based on the real mapping. It uses the center
2241 // point of the cell as a starting point
2242 Point<dim> initial_p_unit;
2243 try
2244 {
2245 initial_p_unit = get_default_linear_mapping(cell->get_triangulation())
2246 .transform_real_to_unit_cell(cell, p);
2247 }
2248 catch (const typename Mapping<dim, spacedim>::ExcTransformationFailed &)
2249 {
2250 // mirror the conditions of the code below to determine if we need to
2251 // use an arbitrary starting point or if we just need to rethrow the
2252 // exception
2253 for (unsigned int d = 0; d < dim; ++d)
2254 initial_p_unit[d] = 0.5;
2255 }
2256
2257 initial_p_unit = cell->reference_cell().closest_point(initial_p_unit);
2258
2259 UpdateFlags update_flags = update_quadrature_points | update_jacobians;
2260 if (spacedim > dim)
2261 update_flags |= update_jacobian_grads;
2262 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2263 get_data(update_flags, Quadrature<dim>(initial_p_unit)));
2264 Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2265 ExcInternalError());
2266
2267 update_internal_dofs(cell, static_cast<InternalData &>(*mdata));
2268
2269 return do_transform_real_to_unit_cell(cell,
2270 p,
2271 initial_p_unit,
2272 static_cast<InternalData &>(*mdata));
2273}
2274
2275
2276template <int dim, int spacedim, typename VectorType>
2277Point<dim>
2278MappingFEField<dim, spacedim, VectorType>::do_transform_real_to_unit_cell(
2279 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2280 const Point<spacedim> &p,
2281 const Point<dim> &starting_guess,
2282 InternalData &mdata) const
2283{
2284 const unsigned int n_shapes = mdata.shape_values.size();
2285 (void)n_shapes;
2286 Assert(n_shapes != 0, ExcInternalError());
2287 AssertDimension(mdata.shape_derivatives.size(), n_shapes);
2288
2289
2290 // Newton iteration to solve
2291 // f(x)=p(x)-p=0
2292 // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
2293 // The start value was set to be the
2294 // linear approximation to the cell
2295 // The shape values and derivatives
2296 // of the mapping at this point are
2297 // previously computed.
2298
2299 Point<dim> p_unit = starting_guess;
2300 Point<dim> f;
2301 mdata.reinit(mdata.update_each, Quadrature<dim>(starting_guess));
2302
2303 Point<spacedim> p_real(do_transform_unit_to_real_cell(mdata));
2304 Tensor<1, spacedim> p_minus_F = p - p_real;
2305 const double eps = 1.e-12 * cell->diameter();
2306 const unsigned int newton_iteration_limit = 20;
2307 unsigned int newton_iteration = 0;
2308 while (p_minus_F.norm_square() > eps * eps)
2309 {
2310 // f'(x)
2311 Point<spacedim> DF[dim];
2312 Tensor<2, dim> df;
2313 for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
2314 {
2315 const Tensor<1, dim> &grad_k = mdata.derivative(0, k);
2316 const unsigned int comp_k =
2317 euler_dof_handler->get_fe().system_to_component_index(k).first;
2318 if (fe_mask[comp_k])
2319 for (unsigned int j = 0; j < dim; ++j)
2320 DF[j][fe_to_real[comp_k]] +=
2321 mdata.local_dof_values[k] * grad_k[j];
2322 }
2323 for (unsigned int j = 0; j < dim; ++j)
2324 {
2325 f[j] = DF[j] * p_minus_F;
2326 for (unsigned int l = 0; l < dim; ++l)
2327 df[j][l] = -DF[j] * DF[l];
2328 }
2329 // Solve [f'(x)]d=f(x)
2330 const Tensor<1, dim> delta =
2331 invert(df) * static_cast<const Tensor<1, dim> &>(f);
2332 // do a line search
2333 double step_length = 1;
2334 do
2335 {
2336 // update of p_unit. The
2337 // spacedimth component of
2338 // transformed point is simply
2339 // ignored in codimension one
2340 // case. When this component is
2341 // not zero, then we are
2342 // projecting the point to the
2343 // surface or curve identified
2344 // by the cell.
2345 Point<dim> p_unit_trial = p_unit;
2346 for (unsigned int i = 0; i < dim; ++i)
2347 p_unit_trial[i] -= step_length * delta[i];
2348 // shape values and derivatives
2349 // at new p_unit point
2350 mdata.reinit(mdata.update_each, Quadrature<dim>(p_unit_trial));
2351 // f(x)
2352 const Point<spacedim> p_real_trial =
2353 do_transform_unit_to_real_cell(mdata);
2354 const Tensor<1, spacedim> f_trial = p - p_real_trial;
2355 // see if we are making progress with the current step length
2356 // and if not, reduce it by a factor of two and try again
2357 if (f_trial.norm() < p_minus_F.norm())
2358 {
2359 p_real = p_real_trial;
2360 p_unit = p_unit_trial;
2361 p_minus_F = f_trial;
2362 break;
2363 }
2364 else if (step_length > 0.05)
2365 step_length /= 2;
2366 else
2367 goto failure;
2368 }
2369 while (true);
2370 ++newton_iteration;
2371 if (newton_iteration > newton_iteration_limit)
2372 goto failure;
2373 }
2374 return p_unit;
2375 // if we get to the following label, then we have either run out
2376 // of Newton iterations, or the line search has not converged.
2377 // in either case, we need to give up, so throw an exception that
2378 // can then be caught
2379failure:
2380 AssertThrow(false,
2381 (typename Mapping<dim, spacedim>::ExcTransformationFailed()));
2382 // ...the compiler wants us to return something, though we can
2383 // of course never get here...
2384 return {};
2385}
2386
2387
2388template <int dim, int spacedim, typename VectorType>
2389unsigned int
2390MappingFEField<dim, spacedim, VectorType>::get_degree() const
2391{
2392 return euler_dof_handler->get_fe().degree;
2393}
2394
2395
2396
2397template <int dim, int spacedim, typename VectorType>
2398ComponentMask
2403
2404
2405template <int dim, int spacedim, typename VectorType>
2406std::unique_ptr<Mapping<dim, spacedim>>
2407MappingFEField<dim, spacedim, VectorType>::clone() const
2408{
2409 return std::make_unique<MappingFEField<dim, spacedim, VectorType>>(*this);
2410}
2411
2412
2413template <int dim, int spacedim, typename VectorType>
2414void
2415MappingFEField<dim, spacedim, VectorType>::update_internal_dofs(
2416 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2417 const typename MappingFEField<dim, spacedim, VectorType>::InternalData &data)
2418 const
2419{
2420 Assert(euler_dof_handler != nullptr,
2421 ExcMessage("euler_dof_handler is empty"));
2422
2423 typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(*cell,
2424 euler_dof_handler);
2425 Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
2426 if (uses_level_dofs)
2427 {
2428 AssertIndexRange(cell->level(), euler_vector.size());
2429 AssertDimension(euler_vector[cell->level()]->size(),
2430 euler_dof_handler->n_dofs(cell->level()));
2431 }
2432 else
2433 AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
2434
2435 if (uses_level_dofs)
2436 dof_cell->get_mg_dof_indices(data.local_dof_indices);
2437 else
2438 dof_cell->get_dof_indices(data.local_dof_indices);
2439
2440 const VectorType &vector =
2441 uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
2442
2443 for (unsigned int i = 0; i < data.local_dof_values.size(); ++i)
2444 data.local_dof_values[i] =
2445 internal::ElementAccess<VectorType>::get(vector,
2446 data.local_dof_indices[i]);
2447}
2448
2449// explicit instantiations
2450#define SPLIT_INSTANTIATIONS_COUNT 2
2451#ifndef SPLIT_INSTANTIATIONS_INDEX
2452# define SPLIT_INSTANTIATIONS_INDEX 0
2453#endif
2454#include "mapping_fe_field.inst"
2455
2456
2457DEAL_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:384
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:493
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:377
MappingFEField::get_vertices
virtual boost::container::small_vector< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
Definition mapping_fe_field.cc:451
MappingFEField::fe_to_real
std::vector< unsigned int > fe_to_real
Definition mapping_fe_field.h:647
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:655
MappingFEField::fe_mask
const ComponentMask fe_mask
Definition mapping_fe_field.h:636
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:2415
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:2072
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:2194
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:2278
MappingFEField::compute_face_data
void compute_face_data(const unsigned int n_original_q_points, InternalData &data) const
Definition mapping_fe_field.cc:567
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:2236
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:634
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:1782
MappingFEField::MappingFEField
friend class MappingFEField
Definition mapping_fe_field.h:666
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:619
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:575
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:1497
MappingFEField::reference_cell
ReferenceCell reference_cell
Definition mapping_fe_field.h:562
MappingFEField::get_degree
unsigned int get_degree() const
Definition mapping_fe_field.cc:2390
MappingFEField::requires_update_flags
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
Definition mapping_fe_field.cc:513
MappingFEField::do_transform_unit_to_real_cell
Point< spacedim > do_transform_unit_to_real_cell(const InternalData &mdata) const
Definition mapping_fe_field.cc:2215
MappingFEField::fe_values_mutex
Threads::Mutex fe_values_mutex
Definition mapping_fe_field.h:658
MappingFEField::uses_level_dofs
const bool uses_level_dofs
Definition mapping_fe_field.h:568
MappingFEField::get_component_mask
ComponentMask get_component_mask() const
Definition mapping_fe_field.cc:2399
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:1701
MappingFEField::euler_dof_handler
ObserverPointer< const DoFHandler< dim, spacedim >, MappingFEField< dim, spacedim, VectorType > > euler_dof_handler
Definition mapping_fe_field.h:582
MappingFEField::fe_values
FEValues< dim, spacedim > fe_values
Definition mapping_fe_field.h:653
MappingFEField::clone
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Definition mapping_fe_field.cc:2407
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:1740
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:1715
ReferenceCell::get_dimension
unsigned int get_dimension() const
Definition reference_cell.h:1499
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:4624
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:468
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:294
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:2503
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition types.h:220
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