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