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