Reference documentation for deal.II version GIT relicensing-487-ge9eb5ab491 2024-04-25 07: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 : 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 unsigned int size = 0;
295 for (unsigned int i = 0; i < fe_mask.size(); ++i)
296 {
297 if (fe_mask[i])
298 fe_to_real[i] = size++;
299 }
300 AssertDimension(size, spacedim);
301}
302
303
304
305template <int dim, int spacedim, typename VectorType>
307 const DoFHandler<dim, spacedim> &euler_dof_handler,
308 const std::vector<VectorType> &euler_vector,
309 const ComponentMask &mask)
310 : reference_cell(euler_dof_handler.get_fe().reference_cell())
311 , uses_level_dofs(true)
312 , euler_dof_handler(&euler_dof_handler)
313 , fe_mask(mask.size() != 0u ?
314 mask :
316 euler_dof_handler.get_fe().get_nonzero_components(0).size(),
317 true))
318 , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
319 , fe_values(this->euler_dof_handler->get_fe(),
320 reference_cell.template get_nodal_type_quadrature<dim>(),
322{
323 unsigned int size = 0;
324 for (unsigned int i = 0; i < fe_mask.size(); ++i)
325 {
326 if (fe_mask[i])
327 fe_to_real[i] = size++;
328 }
329 AssertDimension(size, spacedim);
330
331 Assert(euler_dof_handler.has_level_dofs(),
332 ExcMessage("The underlying DoFHandler object did not call "
333 "distribute_mg_dofs(). In this case, the construction via "
334 "level vectors does not make sense."));
336 euler_dof_handler.get_triangulation().n_global_levels());
337 this->euler_vector.clear();
338 this->euler_vector.resize(euler_vector.size());
339 for (unsigned int i = 0; i < euler_vector.size(); ++i)
340 this->euler_vector[i] = &euler_vector[i];
341}
342
343
344
345template <int dim, int spacedim, typename VectorType>
347 const DoFHandler<dim, spacedim> &euler_dof_handler,
348 const MGLevelObject<VectorType> &euler_vector,
349 const ComponentMask &mask)
350 : reference_cell(euler_dof_handler.get_fe().reference_cell())
351 , uses_level_dofs(true)
352 , euler_dof_handler(&euler_dof_handler)
353 , fe_mask(mask.size() != 0u ?
354 mask :
356 euler_dof_handler.get_fe().get_nonzero_components(0).size(),
357 true))
358 , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
359 , fe_values(this->euler_dof_handler->get_fe(),
360 reference_cell.template get_nodal_type_quadrature<dim>(),
362{
363 unsigned int size = 0;
364 for (unsigned int i = 0; i < fe_mask.size(); ++i)
365 {
366 if (fe_mask[i])
367 fe_to_real[i] = size++;
368 }
369 AssertDimension(size, spacedim);
370
371 Assert(euler_dof_handler.has_level_dofs(),
372 ExcMessage("The underlying DoFHandler object did not call "
373 "distribute_mg_dofs(). In this case, the construction via "
374 "level vectors does not make sense."));
375 AssertDimension(euler_vector.max_level() + 1,
376 euler_dof_handler.get_triangulation().n_global_levels());
377 this->euler_vector.clear();
378 this->euler_vector.resize(
379 euler_dof_handler.get_triangulation().n_global_levels());
380 for (unsigned int i = euler_vector.min_level(); i <= euler_vector.max_level();
381 ++i)
382 this->euler_vector[i] = &euler_vector[i];
383}
384
385
386
387template <int dim, int spacedim, typename VectorType>
390 : reference_cell(mapping.reference_cell)
391 , uses_level_dofs(mapping.uses_level_dofs)
392 , euler_vector(mapping.euler_vector)
393 , euler_dof_handler(mapping.euler_dof_handler)
394 , fe_mask(mapping.fe_mask)
395 , fe_to_real(mapping.fe_to_real)
396 , fe_values(mapping.euler_dof_handler->get_fe(),
397 reference_cell.template get_nodal_type_quadrature<dim>(),
399{}
400
401
402
403template <int dim, int spacedim, typename VectorType>
404inline const double &
406 const unsigned int qpoint,
407 const unsigned int shape_nr) const
408{
409 AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
410 return shape_values[qpoint * n_shape_functions + shape_nr];
411}
412
413
414
415template <int dim, int spacedim, typename VectorType>
416bool
421
422
423
424template <int dim, int spacedim, typename VectorType>
425bool
427 const ReferenceCell &reference_cell) const
428{
430 ExcMessage("The dimension of your mapping (" +
432 ") and the reference cell cell_type (" +
434 " ) do not agree."));
435
436 return this->reference_cell == reference_cell;
437}
438
439
440
441template <int dim, int spacedim, typename VectorType>
442boost::container::small_vector<Point<spacedim>,
445 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
446{
447 // we transform our tria iterator into a dof iterator so we can access
448 // data not associated with triangulations
449 const typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(
450 *cell, euler_dof_handler);
451
452 Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
453 AssertDimension(cell->n_vertices(), fe_values.n_quadrature_points);
455 euler_dof_handler->get_fe().n_components());
456 if (uses_level_dofs)
457 {
458 AssertIndexRange(cell->level(), euler_vector.size());
459 AssertDimension(euler_vector[cell->level()]->size(),
460 euler_dof_handler->n_dofs(cell->level()));
461 }
462 else
463 AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
464
465 {
466 std::lock_guard<std::mutex> lock(fe_values_mutex);
467 fe_values.reinit(dof_cell);
468 }
469 const unsigned int dofs_per_cell =
470 euler_dof_handler->get_fe().n_dofs_per_cell();
471 std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
472 if (uses_level_dofs)
473 dof_cell->get_mg_dof_indices(dof_indices);
474 else
475 dof_cell->get_dof_indices(dof_indices);
476
477 const VectorType &vector =
478 uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
479
480 boost::container::small_vector<Point<spacedim>,
482 vertices(cell->n_vertices());
483 for (unsigned int i = 0; i < dofs_per_cell; ++i)
484 {
485 const unsigned int comp = fe_to_real
486 [euler_dof_handler->get_fe().system_to_component_index(i).first];
488 {
489 typename VectorType::value_type value =
490 internal::ElementAccess<VectorType>::get(vector, dof_indices[i]);
491 if (euler_dof_handler->get_fe().is_primitive(i))
492 for (const unsigned int v : cell->vertex_indices())
493 vertices[v][comp] += fe_values.shape_value(i, v) * value;
494 else
496 }
497 }
498
499 return vertices;
500}
501
502
503
504template <int dim, int spacedim, typename VectorType>
507 const UpdateFlags in) const
508{
509 // add flags if the respective quantities are necessary to compute
510 // what we need. note that some flags appear in both conditions and
511 // in subsequent set operations. this leads to some circular
512 // logic. the only way to treat this is to iterate. since there are
513 // 5 if-clauses in the loop, it will take at most 4 iterations to
514 // converge. do them:
515 UpdateFlags out = in;
516 for (unsigned int i = 0; i < 5; ++i)
517 {
518 // The following is a little incorrect:
519 // If not applied on a face,
520 // update_boundary_forms does not
521 // make sense. On the other hand,
522 // it is necessary on a
523 // face. Currently,
524 // update_boundary_forms is simply
525 // ignored for the interior of a
526 // cell.
529
530 if (out &
534
535 if (out &
540
541 // The contravariant transformation is used in the Piola
542 // transformation, which requires the determinant of the Jacobi
543 // matrix of the transformation. Because we have no way of
544 // knowing here whether the finite element wants to use the
545 // contravariant or the Piola transforms, we add the volume elements
546 // to the list of flags to be updated for each cell.
549
550 if (out & update_normal_vectors)
552 }
553
554 return out;
555}
556
557
558template <int dim, int spacedim, typename VectorType>
559void
561 const unsigned int n_original_q_points,
562 InternalData &data) const
563{
564 // Set to the size of a single quadrature object for faces, as the size set
565 // in in reinit() is for all points
566 if (data.update_each & update_covariant_transformation)
567 data.covariant.resize(n_original_q_points);
568
569 if (data.update_each & update_contravariant_transformation)
570 data.contravariant.resize(n_original_q_points);
571
572 if (data.update_each & update_volume_elements)
573 data.volume_elements.resize(n_original_q_points);
574
575 if (dim > 1)
576 {
577 if (data.update_each & update_boundary_forms)
578 {
579 data.aux.resize(
580 dim - 1, std::vector<Tensor<1, spacedim>>(n_original_q_points));
581
582
583 // TODO: only a single reference cell type possible...
584 const auto n_faces = reference_cell.n_faces();
585
586 // Compute tangentials to the unit cell.
587 for (unsigned int i = 0; i < n_faces; ++i)
588 {
589 data.unit_tangentials[i].resize(n_original_q_points);
590 std::fill(
591 data.unit_tangentials[i].begin(),
592 data.unit_tangentials[i].end(),
593 reference_cell.template unit_tangential_vectors<dim>(i, 0));
594 if (dim > 2)
595 {
596 data.unit_tangentials[n_faces + i].resize(
597 n_original_q_points);
598 std::fill(
599 data.unit_tangentials[n_faces + i].begin(),
600 data.unit_tangentials[n_faces + i].end(),
601 reference_cell.template unit_tangential_vectors<dim>(i, 1));
602 }
603 }
604 }
605 }
606}
607
608
609
610template <int dim, int spacedim, typename VectorType>
611typename std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
613 const UpdateFlags update_flags,
614 const Quadrature<dim> &quadrature) const
615{
616 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
617 std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
618 data_ptr->reinit(requires_update_flags(update_flags), quadrature);
619
620 return data_ptr;
621}
622
623
624
625template <int dim, int spacedim, typename VectorType>
626std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
628 const UpdateFlags update_flags,
629 const hp::QCollection<dim - 1> &quadrature) const
630{
631 AssertDimension(quadrature.size(), 1);
632
633 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
634 std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
635 auto &data = dynamic_cast<InternalData &>(*data_ptr);
636
637 const Quadrature<dim> q(
639 data.reinit(requires_update_flags(update_flags), q);
640 this->compute_face_data(quadrature[0].size(), data);
641
642 return data_ptr;
643}
644
645
646template <int dim, int spacedim, typename VectorType>
647std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
649 const UpdateFlags update_flags,
650 const Quadrature<dim - 1> &quadrature) const
651{
652 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
653 std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
654 auto &data = dynamic_cast<InternalData &>(*data_ptr);
655
656 const Quadrature<dim> q(
658 data.reinit(requires_update_flags(update_flags), q);
659 this->compute_face_data(quadrature.size(), data);
660
661 return data_ptr;
662}
663
664
665
666namespace internal
667{
668 namespace MappingFEFieldImplementation
669 {
670 namespace
671 {
678 template <int dim, int spacedim, typename VectorType>
679 void
680 maybe_compute_q_points(
681 const typename ::QProjector<dim>::DataSetDescriptor data_set,
682 const typename ::MappingFEField<dim, spacedim, VectorType>::
683 InternalData &data,
685 const ComponentMask &fe_mask,
686 const std::vector<unsigned int> &fe_to_real,
687 std::vector<Point<spacedim>> &quadrature_points)
688 {
689 const UpdateFlags update_flags = data.update_each;
690
691 if (update_flags & update_quadrature_points)
692 {
693 for (unsigned int point = 0; point < quadrature_points.size();
694 ++point)
695 {
696 Point<spacedim> result;
697 const double *shape = &data.shape(point + data_set, 0);
698
699 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
700 {
701 const unsigned int comp_k =
702 fe.system_to_component_index(k).first;
703 if (fe_mask[comp_k])
704 result[fe_to_real[comp_k]] +=
705 data.local_dof_values[k] * shape[k];
706 }
707
708 quadrature_points[point] = result;
709 }
710 }
711 }
712
720 template <int dim, int spacedim, typename VectorType>
721 void
722 maybe_update_Jacobians(
723 const typename ::QProjector<dim>::DataSetDescriptor data_set,
724 const typename ::MappingFEField<dim, spacedim, VectorType>::
725 InternalData &data,
727 const ComponentMask &fe_mask,
728 const std::vector<unsigned int> &fe_to_real)
729 {
730 const UpdateFlags update_flags = data.update_each;
731
732 // then Jacobians
733 if (update_flags & update_contravariant_transformation)
734 {
735 const unsigned int n_q_points = data.contravariant.size();
736
737 Assert(data.n_shape_functions > 0, ExcInternalError());
738
739 for (unsigned int point = 0; point < n_q_points; ++point)
740 {
741 const Tensor<1, dim> *data_derv =
742 &data.derivative(point + data_set, 0);
743
744 Tensor<1, dim> result[spacedim];
745
746 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
747 {
748 const unsigned int comp_k =
749 fe.system_to_component_index(k).first;
750 if (fe_mask[comp_k])
751 result[fe_to_real[comp_k]] +=
752 data.local_dof_values[k] * data_derv[k];
753 }
754
755 // write result into contravariant data
756 for (unsigned int i = 0; i < spacedim; ++i)
757 {
758 data.contravariant[point][i] = result[i];
759 }
760 }
761 }
762
763 if (update_flags & update_covariant_transformation)
764 {
765 AssertDimension(data.covariant.size(), data.contravariant.size());
766 for (unsigned int point = 0; point < data.contravariant.size();
767 ++point)
768 data.covariant[point] =
769 (data.contravariant[point]).covariant_form();
770 }
771
772 if (update_flags & update_volume_elements)
773 {
774 AssertDimension(data.contravariant.size(),
775 data.volume_elements.size());
776 for (unsigned int point = 0; point < data.contravariant.size();
777 ++point)
778 data.volume_elements[point] =
779 data.contravariant[point].determinant();
780 }
781 }
782
789 template <int dim, int spacedim, typename VectorType>
790 void
791 maybe_update_jacobian_grads(
792 const typename ::QProjector<dim>::DataSetDescriptor data_set,
793 const typename ::MappingFEField<dim, spacedim, VectorType>::
794 InternalData &data,
796 const ComponentMask &fe_mask,
797 const std::vector<unsigned int> &fe_to_real,
798 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
799 {
800 const UpdateFlags update_flags = data.update_each;
801 if (update_flags & update_jacobian_grads)
802 {
803 const unsigned int n_q_points = jacobian_grads.size();
804
805 for (unsigned int point = 0; point < n_q_points; ++point)
806 {
807 const Tensor<2, dim> *second =
808 &data.second_derivative(point + data_set, 0);
809
811
812 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
813 {
814 const unsigned int comp_k =
815 fe.system_to_component_index(k).first;
816 if (fe_mask[comp_k])
817 for (unsigned int j = 0; j < dim; ++j)
818 for (unsigned int l = 0; l < dim; ++l)
819 result[fe_to_real[comp_k]][j][l] +=
820 (second[k][j][l] * data.local_dof_values[k]);
821 }
822
823 // never touch any data for j=dim in case dim<spacedim, so
824 // it will always be zero as it was initialized
825 for (unsigned int i = 0; i < spacedim; ++i)
826 for (unsigned int j = 0; j < dim; ++j)
827 for (unsigned int l = 0; l < dim; ++l)
828 jacobian_grads[point][i][j][l] = result[i][j][l];
829 }
830 }
831 }
832
839 template <int dim, int spacedim, typename VectorType>
840 void
841 maybe_update_jacobian_pushed_forward_grads(
842 const typename ::QProjector<dim>::DataSetDescriptor data_set,
843 const typename ::MappingFEField<dim, spacedim, VectorType>::
844 InternalData &data,
846 const ComponentMask &fe_mask,
847 const std::vector<unsigned int> &fe_to_real,
848 std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads)
849 {
850 const UpdateFlags update_flags = data.update_each;
851 if (update_flags & update_jacobian_pushed_forward_grads)
852 {
853 const unsigned int n_q_points =
854 jacobian_pushed_forward_grads.size();
855
856 double tmp[spacedim][spacedim][spacedim];
857 for (unsigned int point = 0; point < n_q_points; ++point)
858 {
859 const Tensor<2, dim> *second =
860 &data.second_derivative(point + data_set, 0);
861
863
864 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
865 {
866 const unsigned int comp_k =
867 fe.system_to_component_index(k).first;
868 if (fe_mask[comp_k])
869 for (unsigned int j = 0; j < dim; ++j)
870 for (unsigned int l = 0; l < dim; ++l)
871 result[fe_to_real[comp_k]][j][l] +=
872 (second[k][j][l] * data.local_dof_values[k]);
873 }
874
875 // first push forward the j-components
876 for (unsigned int i = 0; i < spacedim; ++i)
877 for (unsigned int j = 0; j < spacedim; ++j)
878 for (unsigned int l = 0; l < dim; ++l)
879 {
880 tmp[i][j][l] =
881 result[i][0][l] * data.covariant[point][j][0];
882 for (unsigned int jr = 1; jr < dim; ++jr)
883 {
884 tmp[i][j][l] +=
885 result[i][jr][l] * data.covariant[point][j][jr];
886 }
887 }
888
889 // now, pushing forward the l-components
890 for (unsigned int i = 0; i < spacedim; ++i)
891 for (unsigned int j = 0; j < spacedim; ++j)
892 for (unsigned int l = 0; l < spacedim; ++l)
893 {
894 jacobian_pushed_forward_grads[point][i][j][l] =
895 tmp[i][j][0] * data.covariant[point][l][0];
896 for (unsigned int lr = 1; lr < dim; ++lr)
897 {
898 jacobian_pushed_forward_grads[point][i][j][l] +=
899 tmp[i][j][lr] * data.covariant[point][l][lr];
900 }
901 }
902 }
903 }
904 }
905
912 template <int dim, int spacedim, typename VectorType>
913 void
914 maybe_update_jacobian_2nd_derivatives(
915 const typename ::QProjector<dim>::DataSetDescriptor data_set,
916 const typename ::MappingFEField<dim, spacedim, VectorType>::
917 InternalData &data,
919 const ComponentMask &fe_mask,
920 const std::vector<unsigned int> &fe_to_real,
921 std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
922 {
923 const UpdateFlags update_flags = data.update_each;
924 if (update_flags & update_jacobian_2nd_derivatives)
925 {
926 const unsigned int n_q_points = jacobian_2nd_derivatives.size();
927
928 for (unsigned int point = 0; point < n_q_points; ++point)
929 {
930 const Tensor<3, dim> *third =
931 &data.third_derivative(point + data_set, 0);
932
934
935 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
936 {
937 const unsigned int comp_k =
938 fe.system_to_component_index(k).first;
939 if (fe_mask[comp_k])
940 for (unsigned int j = 0; j < dim; ++j)
941 for (unsigned int l = 0; l < dim; ++l)
942 for (unsigned int m = 0; m < dim; ++m)
943 result[fe_to_real[comp_k]][j][l][m] +=
944 (third[k][j][l][m] * data.local_dof_values[k]);
945 }
946
947 // never touch any data for j=dim in case dim<spacedim, so
948 // it will always be zero as it was initialized
949 for (unsigned int i = 0; i < spacedim; ++i)
950 for (unsigned int j = 0; j < dim; ++j)
951 for (unsigned int l = 0; l < dim; ++l)
952 for (unsigned int m = 0; m < dim; ++m)
953 jacobian_2nd_derivatives[point][i][j][l][m] =
954 result[i][j][l][m];
955 }
956 }
957 }
958
966 template <int dim, int spacedim, typename VectorType>
967 void
968 maybe_update_jacobian_pushed_forward_2nd_derivatives(
969 const typename ::QProjector<dim>::DataSetDescriptor data_set,
970 const typename ::MappingFEField<dim, spacedim, VectorType>::
971 InternalData &data,
973 const ComponentMask &fe_mask,
974 const std::vector<unsigned int> &fe_to_real,
975 std::vector<Tensor<4, spacedim>>
976 &jacobian_pushed_forward_2nd_derivatives)
977 {
978 const UpdateFlags update_flags = data.update_each;
980 {
981 const unsigned int n_q_points =
982 jacobian_pushed_forward_2nd_derivatives.size();
983
984 double tmp[spacedim][spacedim][spacedim][spacedim];
985 for (unsigned int point = 0; point < n_q_points; ++point)
986 {
987 const Tensor<3, dim> *third =
988 &data.third_derivative(point + data_set, 0);
989
991
992 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
993 {
994 const unsigned int comp_k =
995 fe.system_to_component_index(k).first;
996 if (fe_mask[comp_k])
997 for (unsigned int j = 0; j < dim; ++j)
998 for (unsigned int l = 0; l < dim; ++l)
999 for (unsigned int m = 0; m < dim; ++m)
1000 result[fe_to_real[comp_k]][j][l][m] +=
1001 (third[k][j][l][m] * data.local_dof_values[k]);
1002 }
1003
1004 // push forward the j-coordinate
1005 for (unsigned int i = 0; i < spacedim; ++i)
1006 for (unsigned int j = 0; j < spacedim; ++j)
1007 for (unsigned int l = 0; l < dim; ++l)
1008 for (unsigned int m = 0; m < dim; ++m)
1009 {
1010 jacobian_pushed_forward_2nd_derivatives
1011 [point][i][j][l][m] =
1012 result[i][0][l][m] * data.covariant[point][j][0];
1013 for (unsigned int jr = 1; jr < dim; ++jr)
1014 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1015 [l][m] +=
1016 result[i][jr][l][m] *
1017 data.covariant[point][j][jr];
1018 }
1019
1020 // push forward the l-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 < spacedim; ++l)
1024 for (unsigned int m = 0; m < dim; ++m)
1025 {
1026 tmp[i][j][l][m] =
1027 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1028 [0][m] *
1029 data.covariant[point][l][0];
1030 for (unsigned int lr = 1; lr < dim; ++lr)
1031 tmp[i][j][l][m] +=
1032 jacobian_pushed_forward_2nd_derivatives[point][i]
1033 [j][lr]
1034 [m] *
1035 data.covariant[point][l][lr];
1036 }
1037
1038 // push forward the m-coordinate
1039 for (unsigned int i = 0; i < spacedim; ++i)
1040 for (unsigned int j = 0; j < spacedim; ++j)
1041 for (unsigned int l = 0; l < spacedim; ++l)
1042 for (unsigned int m = 0; m < spacedim; ++m)
1043 {
1044 jacobian_pushed_forward_2nd_derivatives
1045 [point][i][j][l][m] =
1046 tmp[i][j][l][0] * data.covariant[point][m][0];
1047 for (unsigned int mr = 1; mr < dim; ++mr)
1048 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1049 [l][m] +=
1050 tmp[i][j][l][mr] * data.covariant[point][m][mr];
1051 }
1052 }
1053 }
1054 }
1055
1062 template <int dim, int spacedim, typename VectorType>
1063 void
1064 maybe_update_jacobian_3rd_derivatives(
1065 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1066 const typename ::MappingFEField<dim, spacedim, VectorType>::
1067 InternalData &data,
1069 const ComponentMask &fe_mask,
1070 const std::vector<unsigned int> &fe_to_real,
1071 std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1072 {
1073 const UpdateFlags update_flags = data.update_each;
1074 if (update_flags & update_jacobian_3rd_derivatives)
1075 {
1076 const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1077
1078 for (unsigned int point = 0; point < n_q_points; ++point)
1079 {
1080 const Tensor<4, dim> *fourth =
1081 &data.fourth_derivative(point + data_set, 0);
1082
1084
1085 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1086 {
1087 const unsigned int comp_k =
1088 fe.system_to_component_index(k).first;
1089 if (fe_mask[comp_k])
1090 for (unsigned int j = 0; j < dim; ++j)
1091 for (unsigned int l = 0; l < dim; ++l)
1092 for (unsigned int m = 0; m < dim; ++m)
1093 for (unsigned int n = 0; n < dim; ++n)
1094 result[fe_to_real[comp_k]][j][l][m][n] +=
1095 (fourth[k][j][l][m][n] *
1096 data.local_dof_values[k]);
1097 }
1098
1099 // never touch any data for j,l,m,n=dim in case
1100 // dim<spacedim, so it will always be zero as it was
1101 // initialized
1102 for (unsigned int i = 0; i < spacedim; ++i)
1103 for (unsigned int j = 0; j < dim; ++j)
1104 for (unsigned int l = 0; l < dim; ++l)
1105 for (unsigned int m = 0; m < dim; ++m)
1106 for (unsigned int n = 0; n < dim; ++n)
1107 jacobian_3rd_derivatives[point][i][j][l][m][n] =
1108 result[i][j][l][m][n];
1109 }
1110 }
1111 }
1112
1120 template <int dim, int spacedim, typename VectorType>
1121 void
1122 maybe_update_jacobian_pushed_forward_3rd_derivatives(
1123 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1124 const typename ::MappingFEField<dim, spacedim, VectorType>::
1125 InternalData &data,
1127 const ComponentMask &fe_mask,
1128 const std::vector<unsigned int> &fe_to_real,
1129 std::vector<Tensor<5, spacedim>>
1130 &jacobian_pushed_forward_3rd_derivatives)
1131 {
1132 const UpdateFlags update_flags = data.update_each;
1134 {
1135 const unsigned int n_q_points =
1136 jacobian_pushed_forward_3rd_derivatives.size();
1137
1138 double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
1139 for (unsigned int point = 0; point < n_q_points; ++point)
1140 {
1141 const Tensor<4, dim> *fourth =
1142 &data.fourth_derivative(point + data_set, 0);
1143
1145
1146 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1147 {
1148 const unsigned int comp_k =
1149 fe.system_to_component_index(k).first;
1150 if (fe_mask[comp_k])
1151 for (unsigned int j = 0; j < dim; ++j)
1152 for (unsigned int l = 0; l < dim; ++l)
1153 for (unsigned int m = 0; m < dim; ++m)
1154 for (unsigned int n = 0; n < dim; ++n)
1155 result[fe_to_real[comp_k]][j][l][m][n] +=
1156 (fourth[k][j][l][m][n] *
1157 data.local_dof_values[k]);
1158 }
1159
1160 // push-forward the j-coordinate
1161 for (unsigned int i = 0; i < spacedim; ++i)
1162 for (unsigned int j = 0; j < spacedim; ++j)
1163 for (unsigned int l = 0; l < dim; ++l)
1164 for (unsigned int m = 0; m < dim; ++m)
1165 for (unsigned int n = 0; n < dim; ++n)
1166 {
1167 tmp[i][j][l][m][n] = result[i][0][l][m][n] *
1168 data.covariant[point][j][0];
1169 for (unsigned int jr = 1; jr < dim; ++jr)
1170 tmp[i][j][l][m][n] +=
1171 result[i][jr][l][m][n] *
1172 data.covariant[point][j][jr];
1173 }
1174
1175 // push-forward the l-coordinate
1176 for (unsigned int i = 0; i < spacedim; ++i)
1177 for (unsigned int j = 0; j < spacedim; ++j)
1178 for (unsigned int l = 0; l < spacedim; ++l)
1179 for (unsigned int m = 0; m < dim; ++m)
1180 for (unsigned int n = 0; n < dim; ++n)
1181 {
1182 jacobian_pushed_forward_3rd_derivatives
1183 [point][i][j][l][m][n] =
1184 tmp[i][j][0][m][n] *
1185 data.covariant[point][l][0];
1186 for (unsigned int lr = 1; lr < dim; ++lr)
1187 jacobian_pushed_forward_3rd_derivatives[point][i]
1188 [j][l][m]
1189 [n] +=
1190 tmp[i][j][lr][m][n] *
1191 data.covariant[point][l][lr];
1192 }
1193
1194 // push-forward the m-coordinate
1195 for (unsigned int i = 0; i < spacedim; ++i)
1196 for (unsigned int j = 0; j < spacedim; ++j)
1197 for (unsigned int l = 0; l < spacedim; ++l)
1198 for (unsigned int m = 0; m < spacedim; ++m)
1199 for (unsigned int n = 0; n < dim; ++n)
1200 {
1201 tmp[i][j][l][m][n] =
1202 jacobian_pushed_forward_3rd_derivatives[point][i]
1203 [j][l][0]
1204 [n] *
1205 data.covariant[point][m][0];
1206 for (unsigned int mr = 1; mr < dim; ++mr)
1207 tmp[i][j][l][m][n] +=
1208 jacobian_pushed_forward_3rd_derivatives[point]
1209 [i][j][l]
1210 [mr][n] *
1211 data.covariant[point][m][mr];
1212 }
1213
1214 // push-forward the n-coordinate
1215 for (unsigned int i = 0; i < spacedim; ++i)
1216 for (unsigned int j = 0; j < spacedim; ++j)
1217 for (unsigned int l = 0; l < spacedim; ++l)
1218 for (unsigned int m = 0; m < spacedim; ++m)
1219 for (unsigned int n = 0; n < spacedim; ++n)
1220 {
1221 jacobian_pushed_forward_3rd_derivatives
1222 [point][i][j][l][m][n] =
1223 tmp[i][j][l][m][0] *
1224 data.covariant[point][n][0];
1225 for (unsigned int nr = 1; nr < dim; ++nr)
1226 jacobian_pushed_forward_3rd_derivatives[point][i]
1227 [j][l][m]
1228 [n] +=
1229 tmp[i][j][l][m][nr] *
1230 data.covariant[point][n][nr];
1231 }
1232 }
1233 }
1234 }
1235
1236
1246 template <int dim, int spacedim, typename VectorType>
1247 void
1248 maybe_compute_face_data(
1249 const ::Mapping<dim, spacedim> &mapping,
1250 const typename ::Triangulation<dim, spacedim>::cell_iterator
1251 &cell,
1252 const unsigned int face_no,
1253 const unsigned int subface_no,
1254 const typename QProjector<dim>::DataSetDescriptor data_set,
1255 const typename ::MappingFEField<dim, spacedim, VectorType>::
1256 InternalData &data,
1258 &output_data)
1259 {
1260 const UpdateFlags update_flags = data.update_each;
1261
1262 if (update_flags & update_boundary_forms)
1263 {
1264 const unsigned int n_q_points = output_data.boundary_forms.size();
1265 if (update_flags & update_normal_vectors)
1266 AssertDimension(output_data.normal_vectors.size(), n_q_points);
1267 if (update_flags & update_JxW_values)
1268 AssertDimension(output_data.JxW_values.size(), n_q_points);
1269
1270 // map the unit tangentials to the real cell. checking for d!=dim-1
1271 // eliminates compiler warnings regarding unsigned int expressions <
1272 // 0.
1273 for (unsigned int d = 0; d != dim - 1; ++d)
1274 {
1275 Assert(face_no + cell->n_faces() * d <
1276 data.unit_tangentials.size(),
1278 Assert(
1279 data.aux[d].size() <=
1280 data.unit_tangentials[face_no + cell->n_faces() * d].size(),
1282
1283 mapping.transform(
1285 data.unit_tangentials[face_no + cell->n_faces() * d]),
1287 data,
1288 make_array_view(data.aux[d]));
1289 }
1290
1291 // if dim==spacedim, we can use the unit tangentials to compute the
1292 // boundary form by simply taking the cross product
1293 if (dim == spacedim)
1294 {
1295 for (unsigned int i = 0; i < n_q_points; ++i)
1296 switch (dim)
1297 {
1298 case 1:
1299 // in 1d, we don't have access to any of the data.aux
1300 // fields (because it has only dim-1 components), but we
1301 // can still compute the boundary form by simply looking
1302 // at the number of the face
1303 output_data.boundary_forms[i][0] =
1304 (face_no == 0 ? -1 : +1);
1305 break;
1306 case 2:
1307 output_data.boundary_forms[i] =
1308 cross_product_2d(data.aux[0][i]);
1309 break;
1310 case 3:
1311 output_data.boundary_forms[i] =
1312 cross_product_3d(data.aux[0][i], data.aux[1][i]);
1313 break;
1314 default:
1316 }
1317 }
1318 else //(dim < spacedim)
1319 {
1320 // in the codim-one case, the boundary form results from the
1321 // cross product of all the face tangential vectors and the cell
1322 // normal vector
1323 //
1324 // to compute the cell normal, use the same method used in
1325 // fill_fe_values for cells above
1326 AssertDimension(data.contravariant.size(), n_q_points);
1327
1328 for (unsigned int point = 0; point < n_q_points; ++point)
1329 {
1330 if (dim == 1)
1331 {
1332 // J is a tangent vector
1333 output_data.boundary_forms[point] =
1334 data.contravariant[point].transpose()[0];
1335 output_data.boundary_forms[point] /=
1336 (face_no == 0 ? -1. : +1.) *
1337 output_data.boundary_forms[point].norm();
1338 }
1339
1340 if (dim == 2)
1341 {
1343 data.contravariant[point].transpose();
1344
1345 Tensor<1, spacedim> cell_normal =
1346 cross_product_3d(DX_t[0], DX_t[1]);
1347 cell_normal /= cell_normal.norm();
1348
1349 // then compute the face normal from the face tangent
1350 // and the cell normal:
1351 output_data.boundary_forms[point] =
1352 cross_product_3d(data.aux[0][point], cell_normal);
1353 }
1354 }
1355 }
1356
1357 if (update_flags & (update_normal_vectors | update_JxW_values))
1358 for (unsigned int i = 0; i < output_data.boundary_forms.size();
1359 ++i)
1360 {
1361 if (update_flags & update_JxW_values)
1362 {
1363 output_data.JxW_values[i] =
1364 output_data.boundary_forms[i].norm() *
1365 data.quadrature_weights[i + data_set];
1366
1367 if (subface_no != numbers::invalid_unsigned_int)
1368 {
1369 // TODO
1370 const double area_ratio =
1372 cell->subface_case(face_no), subface_no);
1373 output_data.JxW_values[i] *= area_ratio;
1374 }
1375 }
1376
1377 if (update_flags & update_normal_vectors)
1378 output_data.normal_vectors[i] =
1379 Point<spacedim>(output_data.boundary_forms[i] /
1380 output_data.boundary_forms[i].norm());
1381 }
1382 }
1383 }
1384
1391 template <int dim, int spacedim, typename VectorType>
1392 void
1393 do_fill_fe_face_values(
1394 const ::Mapping<dim, spacedim> &mapping,
1395 const typename ::Triangulation<dim, spacedim>::cell_iterator
1396 &cell,
1397 const unsigned int face_no,
1398 const unsigned int subface_no,
1399 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1400 const typename ::MappingFEField<dim, spacedim, VectorType>::
1401 InternalData &data,
1403 const ComponentMask &fe_mask,
1404 const std::vector<unsigned int> &fe_to_real,
1406 &output_data)
1407 {
1408 maybe_compute_q_points<dim, spacedim, VectorType>(
1409 data_set,
1410 data,
1411 fe,
1412 fe_mask,
1413 fe_to_real,
1414 output_data.quadrature_points);
1415
1416 maybe_update_Jacobians<dim, spacedim, VectorType>(
1417 data_set, data, fe, fe_mask, fe_to_real);
1418
1419 const UpdateFlags update_flags = data.update_each;
1420 const unsigned int n_q_points = data.contravariant.size();
1421
1422 if (update_flags & update_jacobians)
1423 for (unsigned int point = 0; point < n_q_points; ++point)
1424 output_data.jacobians[point] = data.contravariant[point];
1425
1426 if (update_flags & update_inverse_jacobians)
1427 for (unsigned int point = 0; point < n_q_points; ++point)
1428 output_data.inverse_jacobians[point] =
1429 data.covariant[point].transpose();
1430
1431 maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1432 data_set, data, fe, fe_mask, fe_to_real, output_data.jacobian_grads);
1433
1434 maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1435 data_set,
1436 data,
1437 fe,
1438 fe_mask,
1439 fe_to_real,
1440 output_data.jacobian_pushed_forward_grads);
1441
1442 maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1443 data_set,
1444 data,
1445 fe,
1446 fe_mask,
1447 fe_to_real,
1448 output_data.jacobian_2nd_derivatives);
1449
1450 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1451 spacedim,
1452 VectorType>(
1453 data_set,
1454 data,
1455 fe,
1456 fe_mask,
1457 fe_to_real,
1458 output_data.jacobian_pushed_forward_2nd_derivatives);
1459
1460 maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1461 data_set,
1462 data,
1463 fe,
1464 fe_mask,
1465 fe_to_real,
1466 output_data.jacobian_3rd_derivatives);
1467
1468 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1469 spacedim,
1470 VectorType>(
1471 data_set,
1472 data,
1473 fe,
1474 fe_mask,
1475 fe_to_real,
1476 output_data.jacobian_pushed_forward_3rd_derivatives);
1477
1478 maybe_compute_face_data<dim, spacedim, VectorType>(
1479 mapping, cell, face_no, subface_no, data_set, data, output_data);
1480 }
1481 } // namespace
1482 } // namespace MappingFEFieldImplementation
1483} // namespace internal
1484
1485
1486// Note that the CellSimilarity flag is modifiable, since MappingFEField can
1487// need to recalculate data even when cells are similar.
1488template <int dim, int spacedim, typename VectorType>
1493 const Quadrature<dim> &quadrature,
1494 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1496 &output_data) const
1497{
1498 // convert data object to internal data for this class. fails with an
1499 // exception if that is not possible
1500 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1502 const InternalData &data = static_cast<const InternalData &>(internal_data);
1503
1504 const unsigned int n_q_points = quadrature.size();
1505
1506 update_internal_dofs(cell, data);
1507
1508 internal::MappingFEFieldImplementation::
1509 maybe_compute_q_points<dim, spacedim, VectorType>(
1511 data,
1512 euler_dof_handler->get_fe(),
1513 fe_mask,
1514 fe_to_real,
1515 output_data.quadrature_points);
1516
1517 internal::MappingFEFieldImplementation::
1518 maybe_update_Jacobians<dim, spacedim, VectorType>(
1520 data,
1521 euler_dof_handler->get_fe(),
1522 fe_mask,
1523 fe_to_real);
1524
1525 const UpdateFlags update_flags = data.update_each;
1526 const std::vector<double> &weights = quadrature.get_weights();
1527
1528 // Multiply quadrature weights by absolute value of Jacobian determinants or
1529 // the area element g=sqrt(DX^t DX) in case of codim > 0
1530
1531 if (update_flags & (update_normal_vectors | update_JxW_values))
1532 {
1533 AssertDimension(output_data.JxW_values.size(), n_q_points);
1534
1535 Assert(!(update_flags & update_normal_vectors) ||
1536 (output_data.normal_vectors.size() == n_q_points),
1537 ExcDimensionMismatch(output_data.normal_vectors.size(),
1538 n_q_points));
1539
1540
1541 for (unsigned int point = 0; point < n_q_points; ++point)
1542 {
1543 if (dim == spacedim)
1544 {
1545 const double det = data.volume_elements[point];
1546
1547 // check for distorted cells.
1548
1549 // TODO: this allows for anisotropies of up to 1e6 in 3d and
1550 // 1e12 in 2d. might want to find a finer
1551 // (dimension-independent) criterion
1552 Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1553 cell->diameter() / std::sqrt(double(dim))),
1555 cell->center(), det, point)));
1556 output_data.JxW_values[point] = weights[point] * det;
1557 }
1558 // if dim==spacedim, then there is no cell normal to
1559 // compute. since this is for FEValues (and not FEFaceValues),
1560 // there are also no face normals to compute
1561 else // codim>0 case
1562 {
1563 Tensor<1, spacedim> DX_t[dim];
1564 for (unsigned int i = 0; i < spacedim; ++i)
1565 for (unsigned int j = 0; j < dim; ++j)
1566 DX_t[j][i] = data.contravariant[point][i][j];
1567
1568 Tensor<2, dim> G; // First fundamental form
1569 for (unsigned int i = 0; i < dim; ++i)
1570 for (unsigned int j = 0; j < dim; ++j)
1571 G[i][j] = DX_t[i] * DX_t[j];
1572
1573 output_data.JxW_values[point] =
1574 std::sqrt(determinant(G)) * weights[point];
1575
1576 if (update_flags & update_normal_vectors)
1577 {
1578 Assert(spacedim - dim == 1,
1579 ExcMessage("There is no cell normal in codim 2."));
1580
1581 if (dim == 1)
1582 output_data.normal_vectors[point] =
1583 cross_product_2d(-DX_t[0]);
1584 else
1585 {
1586 Assert(dim == 2, ExcInternalError());
1587
1588 // dim-1==1 for the second argument, but this
1589 // avoids a compiler warning about array bounds:
1590 output_data.normal_vectors[point] =
1591 cross_product_3d(DX_t[0], DX_t[dim - 1]);
1592 }
1593
1594 output_data.normal_vectors[point] /=
1595 output_data.normal_vectors[point].norm();
1596
1597 if (cell->direction_flag() == false)
1598 output_data.normal_vectors[point] *= -1.;
1599 }
1600 } // codim>0 case
1601 }
1602 }
1603
1604 // copy values from InternalData to vector given by reference
1605 if (update_flags & update_jacobians)
1606 {
1607 AssertDimension(output_data.jacobians.size(), n_q_points);
1608 for (unsigned int point = 0; point < n_q_points; ++point)
1609 output_data.jacobians[point] = data.contravariant[point];
1610 }
1611
1612 // copy values from InternalData to vector given by reference
1613 if (update_flags & update_inverse_jacobians)
1614 {
1615 AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1616 for (unsigned int point = 0; point < n_q_points; ++point)
1617 output_data.inverse_jacobians[point] =
1618 data.covariant[point].transpose();
1619 }
1620
1621 // calculate derivatives of the Jacobians
1622 internal::MappingFEFieldImplementation::
1623 maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1625 data,
1626 euler_dof_handler->get_fe(),
1627 fe_mask,
1628 fe_to_real,
1629 output_data.jacobian_grads);
1630
1631 // calculate derivatives of the Jacobians pushed forward to real cell
1632 // coordinates
1633 internal::MappingFEFieldImplementation::
1634 maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1636 data,
1637 euler_dof_handler->get_fe(),
1638 fe_mask,
1639 fe_to_real,
1640 output_data.jacobian_pushed_forward_grads);
1641
1642 // calculate hessians of the Jacobians
1643 internal::MappingFEFieldImplementation::
1644 maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1646 data,
1647 euler_dof_handler->get_fe(),
1648 fe_mask,
1649 fe_to_real,
1650 output_data.jacobian_2nd_derivatives);
1651
1652 // calculate hessians of the Jacobians pushed forward to real cell coordinates
1653 internal::MappingFEFieldImplementation::
1654 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1655 spacedim,
1656 VectorType>(
1658 data,
1659 euler_dof_handler->get_fe(),
1660 fe_mask,
1661 fe_to_real,
1663
1664 // calculate gradients of the hessians of the Jacobians
1665 internal::MappingFEFieldImplementation::
1666 maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1668 data,
1669 euler_dof_handler->get_fe(),
1670 fe_mask,
1671 fe_to_real,
1672 output_data.jacobian_3rd_derivatives);
1673
1674 // calculate gradients of the hessians of the Jacobians pushed forward to real
1675 // cell coordinates
1676 internal::MappingFEFieldImplementation::
1677 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1678 spacedim,
1679 VectorType>(
1681 data,
1682 euler_dof_handler->get_fe(),
1683 fe_mask,
1684 fe_to_real,
1686
1688}
1689
1690
1691
1692template <int dim, int spacedim, typename VectorType>
1693void
1696 const unsigned int face_no,
1697 const hp::QCollection<dim - 1> &quadrature,
1698 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1700 &output_data) const
1701{
1702 AssertDimension(quadrature.size(), 1);
1703
1704 // convert data object to internal data for this class. fails with an
1705 // exception if that is not possible
1706 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1708 const InternalData &data = static_cast<const InternalData &>(internal_data);
1709
1710 update_internal_dofs(cell, data);
1711
1712 internal::MappingFEFieldImplementation::
1713 do_fill_fe_face_values<dim, spacedim, VectorType>(
1714 *this,
1715 cell,
1716 face_no,
1719 face_no,
1720 cell->combined_face_orientation(
1721 face_no),
1722 quadrature[0].size()),
1723 data,
1724 euler_dof_handler->get_fe(),
1725 fe_mask,
1726 fe_to_real,
1727 output_data);
1728}
1729
1730
1731template <int dim, int spacedim, typename VectorType>
1732void
1735 const unsigned int face_no,
1736 const unsigned int subface_no,
1737 const Quadrature<dim - 1> &quadrature,
1738 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1740 &output_data) const
1741{
1742 // convert data object to internal data for this class. fails with an
1743 // exception if that is not possible
1744 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1746 const InternalData &data = static_cast<const InternalData &>(internal_data);
1747
1748 update_internal_dofs(cell, data);
1749
1750 internal::MappingFEFieldImplementation::do_fill_fe_face_values<dim,
1751 spacedim,
1752 VectorType>(
1753 *this,
1754 cell,
1755 face_no,
1758 face_no,
1759 subface_no,
1760 cell->combined_face_orientation(
1761 face_no),
1762 quadrature.size(),
1763 cell->subface_case(face_no)),
1764 data,
1765 euler_dof_handler->get_fe(),
1766 fe_mask,
1767 fe_to_real,
1768 output_data);
1769}
1770
1771
1772
1773template <int dim, int spacedim, typename VectorType>
1774void
1778 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1780 &output_data) const
1781{
1782 AssertDimension(dim, spacedim);
1783 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1785 const InternalData &data = static_cast<const InternalData &>(internal_data);
1786
1787 const unsigned int n_q_points = quadrature.size();
1788
1789 update_internal_dofs(cell, data);
1790
1791 internal::MappingFEFieldImplementation::
1792 maybe_compute_q_points<dim, spacedim, VectorType>(
1794 data,
1795 euler_dof_handler->get_fe(),
1796 fe_mask,
1797 fe_to_real,
1798 output_data.quadrature_points);
1799
1800 internal::MappingFEFieldImplementation::
1801 maybe_update_Jacobians<dim, spacedim, VectorType>(
1803 data,
1804 euler_dof_handler->get_fe(),
1805 fe_mask,
1806 fe_to_real);
1807
1808 const UpdateFlags update_flags = data.update_each;
1809 const std::vector<double> &weights = quadrature.get_weights();
1810
1811 if (update_flags & (update_normal_vectors | update_JxW_values))
1812 {
1813 AssertDimension(output_data.JxW_values.size(), n_q_points);
1814
1815 Assert(!(update_flags & update_normal_vectors) ||
1816 (output_data.normal_vectors.size() == n_q_points),
1817 ExcDimensionMismatch(output_data.normal_vectors.size(),
1818 n_q_points));
1819
1820
1821 for (unsigned int point = 0; point < n_q_points; ++point)
1822 {
1823 const double det = data.volume_elements[point];
1824
1825 // check for distorted cells.
1826
1827 // TODO: this allows for anisotropies of up to 1e6 in 3d and
1828 // 1e12 in 2d. might want to find a finer
1829 // (dimension-independent) criterion
1830 Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1831 cell->diameter() / std::sqrt(double(dim))),
1833 cell->center(), det, point)));
1834
1835 // The normals are n = J^{-T} * \hat{n} before normalizing.
1836 Tensor<1, spacedim> normal;
1837 for (unsigned int d = 0; d < spacedim; d++)
1838 normal[d] =
1839 data.covariant[point][d] * quadrature.normal_vector(point);
1840
1841 output_data.JxW_values[point] = weights[point] * det * normal.norm();
1842
1843 if ((update_flags & update_normal_vectors) != 0u)
1844 {
1845 normal /= normal.norm();
1846 output_data.normal_vectors[point] = normal;
1847 }
1848 }
1849
1850 // copy values from InternalData to vector given by reference
1851 if (update_flags & update_jacobians)
1852 {
1853 AssertDimension(output_data.jacobians.size(), n_q_points);
1854 for (unsigned int point = 0; point < n_q_points; ++point)
1855 output_data.jacobians[point] = data.contravariant[point];
1856 }
1857
1858 // copy values from InternalData to vector given by reference
1859 if (update_flags & update_inverse_jacobians)
1860 {
1861 AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1862 for (unsigned int point = 0; point < n_q_points; ++point)
1863 output_data.inverse_jacobians[point] =
1864 data.covariant[point].transpose();
1865 }
1866
1867 // calculate derivatives of the Jacobians
1868 internal::MappingFEFieldImplementation::
1869 maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1871 data,
1872 euler_dof_handler->get_fe(),
1873 fe_mask,
1874 fe_to_real,
1875 output_data.jacobian_grads);
1876
1877 // calculate derivatives of the Jacobians pushed forward to real cell
1878 // coordinates
1879 internal::MappingFEFieldImplementation::
1880 maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1882 data,
1883 euler_dof_handler->get_fe(),
1884 fe_mask,
1885 fe_to_real,
1886 output_data.jacobian_pushed_forward_grads);
1887
1888 // calculate hessians of the Jacobians
1889 internal::MappingFEFieldImplementation::
1890 maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1892 data,
1893 euler_dof_handler->get_fe(),
1894 fe_mask,
1895 fe_to_real,
1896 output_data.jacobian_2nd_derivatives);
1897
1898 // calculate hessians of the Jacobians pushed forward to real cell
1899 // coordinates
1900 internal::MappingFEFieldImplementation::
1901 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1902 spacedim,
1903 VectorType>(
1905 data,
1906 euler_dof_handler->get_fe(),
1907 fe_mask,
1908 fe_to_real,
1910
1911 // calculate gradients of the hessians of the Jacobians
1912 internal::MappingFEFieldImplementation::
1913 maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1915 data,
1916 euler_dof_handler->get_fe(),
1917 fe_mask,
1918 fe_to_real,
1919 output_data.jacobian_3rd_derivatives);
1920
1921 // calculate gradients of the hessians of the Jacobians pushed forward to
1922 // real cell coordinates
1923 internal::MappingFEFieldImplementation::
1924 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1925 spacedim,
1926 VectorType>(
1928 data,
1929 euler_dof_handler->get_fe(),
1930 fe_mask,
1931 fe_to_real,
1933 }
1934}
1935
1936namespace internal
1937{
1938 namespace MappingFEFieldImplementation
1939 {
1940 namespace
1941 {
1942 template <int dim, int spacedim, int rank, typename VectorType>
1943 void
1944 transform_fields(
1945 const ArrayView<const Tensor<rank, dim>> &input,
1946 const MappingKind mapping_kind,
1947 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1948 const ArrayView<Tensor<rank, spacedim>> &output)
1949 {
1950 AssertDimension(input.size(), output.size());
1951 Assert((dynamic_cast<
1952 const typename ::
1953 MappingFEField<dim, spacedim, VectorType>::InternalData *>(
1954 &mapping_data) != nullptr),
1956 const typename ::MappingFEField<dim, spacedim, VectorType>::
1957 InternalData &data = static_cast<
1958 const typename ::MappingFEField<dim, spacedim, VectorType>::
1959 InternalData &>(mapping_data);
1960
1961 switch (mapping_kind)
1962 {
1964 {
1965 Assert(
1966 data.update_each & update_contravariant_transformation,
1968 "update_contravariant_transformation"));
1969
1970 for (unsigned int i = 0; i < output.size(); ++i)
1971 output[i] =
1972 apply_transformation(data.contravariant[i], input[i]);
1973
1974 return;
1975 }
1976
1977 case mapping_piola:
1978 {
1979 Assert(
1980 data.update_each & update_contravariant_transformation,
1982 "update_contravariant_transformation"));
1983 Assert(
1984 data.update_each & update_volume_elements,
1986 "update_volume_elements"));
1987 Assert(rank == 1, ExcMessage("Only for rank 1"));
1988 for (unsigned int i = 0; i < output.size(); ++i)
1989 {
1990 output[i] =
1991 apply_transformation(data.contravariant[i], input[i]);
1992 output[i] /= data.volume_elements[i];
1993 }
1994 return;
1995 }
1996
1997
1998 // We still allow this operation as in the
1999 // reference cell Derivatives are Tensor
2000 // rather than DerivativeForm
2001 case mapping_covariant:
2002 {
2003 Assert(
2004 data.update_each & update_contravariant_transformation,
2006 "update_contravariant_transformation"));
2007
2008 for (unsigned int i = 0; i < output.size(); ++i)
2009 output[i] = apply_transformation(data.covariant[i], input[i]);
2010
2011 return;
2012 }
2013
2014 default:
2016 }
2017 }
2018
2019
2020 template <int dim, int spacedim, int rank, typename VectorType>
2021 void
2024 const MappingKind mapping_kind,
2025 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2027 {
2028 AssertDimension(input.size(), output.size());
2029 Assert((dynamic_cast<
2030 const typename ::
2031 MappingFEField<dim, spacedim, VectorType>::InternalData *>(
2032 &mapping_data) != nullptr),
2034 const typename ::MappingFEField<dim, spacedim, VectorType>::
2035 InternalData &data = static_cast<
2036 const typename ::MappingFEField<dim, spacedim, VectorType>::
2037 InternalData &>(mapping_data);
2038
2039 switch (mapping_kind)
2040 {
2041 case mapping_covariant:
2042 {
2043 Assert(
2044 data.update_each & update_contravariant_transformation,
2046 "update_contravariant_transformation"));
2047
2048 for (unsigned int i = 0; i < output.size(); ++i)
2049 output[i] = apply_transformation(data.covariant[i], input[i]);
2050
2051 return;
2052 }
2053 default:
2055 }
2056 }
2057 } // namespace
2058 } // namespace MappingFEFieldImplementation
2059} // namespace internal
2060
2061
2062
2063template <int dim, int spacedim, typename VectorType>
2064void
2066 const ArrayView<const Tensor<1, dim>> &input,
2067 const MappingKind mapping_kind,
2068 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2069 const ArrayView<Tensor<1, spacedim>> &output) const
2070{
2071 AssertDimension(input.size(), output.size());
2072
2073 internal::MappingFEFieldImplementation::
2074 transform_fields<dim, spacedim, 1, VectorType>(input,
2075 mapping_kind,
2076 mapping_data,
2077 output);
2078}
2079
2080
2081
2082template <int dim, int spacedim, typename VectorType>
2083void
2085 const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
2086 const MappingKind mapping_kind,
2087 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2088 const ArrayView<Tensor<2, spacedim>> &output) const
2089{
2090 AssertDimension(input.size(), output.size());
2091
2092 internal::MappingFEFieldImplementation::
2093 transform_differential_forms<dim, spacedim, 1, VectorType>(input,
2094 mapping_kind,
2095 mapping_data,
2096 output);
2097}
2098
2099
2100
2101template <int dim, int spacedim, typename VectorType>
2102void
2104 const ArrayView<const Tensor<2, dim>> &input,
2105 const MappingKind,
2106 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2107 const ArrayView<Tensor<2, spacedim>> &output) const
2108{
2109 (void)input;
2110 (void)output;
2111 (void)mapping_data;
2112 AssertDimension(input.size(), output.size());
2113
2115}
2116
2117
2118
2119template <int dim, int spacedim, typename VectorType>
2120void
2122 const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
2123 const MappingKind mapping_kind,
2124 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2125 const ArrayView<Tensor<3, spacedim>> &output) const
2126{
2127 AssertDimension(input.size(), output.size());
2128 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
2130 const InternalData &data = static_cast<const InternalData &>(mapping_data);
2131
2132 switch (mapping_kind)
2133 {
2135 {
2138 "update_covariant_transformation"));
2139
2140 for (unsigned int q = 0; q < output.size(); ++q)
2141 for (unsigned int i = 0; i < spacedim; ++i)
2142 for (unsigned int j = 0; j < spacedim; ++j)
2143 for (unsigned int k = 0; k < spacedim; ++k)
2144 {
2145 output[q][i][j][k] = data.covariant[q][j][0] *
2146 data.covariant[q][k][0] *
2147 input[q][i][0][0];
2148 for (unsigned int J = 0; J < dim; ++J)
2149 {
2150 const unsigned int K0 = (0 == J) ? 1 : 0;
2151 for (unsigned int K = K0; K < dim; ++K)
2152 output[q][i][j][k] += data.covariant[q][j][J] *
2153 data.covariant[q][k][K] *
2154 input[q][i][J][K];
2155 }
2156 }
2157 return;
2158 }
2159
2160 default:
2162 }
2163}
2164
2165
2166
2167template <int dim, int spacedim, typename VectorType>
2168void
2170 const ArrayView<const Tensor<3, dim>> &input,
2171 const MappingKind /*mapping_kind*/,
2172 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2173 const ArrayView<Tensor<3, spacedim>> &output) const
2174{
2175 (void)input;
2176 (void)output;
2177 (void)mapping_data;
2178 AssertDimension(input.size(), output.size());
2179
2181}
2182
2183
2184
2185template <int dim, int spacedim, typename VectorType>
2189 const Point<dim> &p) const
2190{
2191 // Use the get_data function to create an InternalData with data vectors of
2192 // the right size and transformation shape values already computed at point
2193 // p.
2194 const Quadrature<dim> point_quadrature(p);
2195 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2197 Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2199
2200 update_internal_dofs(cell, static_cast<InternalData &>(*mdata));
2201
2202 return do_transform_unit_to_real_cell(static_cast<InternalData &>(*mdata));
2203}
2204
2205
2206template <int dim, int spacedim, typename VectorType>
2209 const InternalData &data) const
2210{
2211 Point<spacedim> p_real;
2212
2213 for (unsigned int i = 0; i < data.n_shape_functions; ++i)
2214 {
2215 unsigned int comp_i =
2216 euler_dof_handler->get_fe().system_to_component_index(i).first;
2217 if (fe_mask[comp_i])
2218 p_real[fe_to_real[comp_i]] +=
2219 data.local_dof_values[i] * data.shape(0, i);
2220 }
2221
2222 return p_real;
2223}
2224
2225
2226
2227template <int dim, int spacedim, typename VectorType>
2231 const Point<spacedim> &p) const
2232{
2233 // first a Newton iteration based on the real mapping. It uses the center
2234 // point of the cell as a starting point
2235 Point<dim> initial_p_unit;
2236 try
2237 {
2238 initial_p_unit = get_default_linear_mapping(cell->get_triangulation())
2240 }
2242 {
2243 // mirror the conditions of the code below to determine if we need to
2244 // use an arbitrary starting point or if we just need to rethrow the
2245 // exception
2246 for (unsigned int d = 0; d < dim; ++d)
2247 initial_p_unit[d] = 0.5;
2248 }
2249
2250 initial_p_unit = cell->reference_cell().closest_point(initial_p_unit);
2251
2252 Quadrature<dim> point_quadrature(initial_p_unit);
2253
2255 if (spacedim > dim)
2256 update_flags |= update_jacobian_grads;
2257 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2258 get_data(update_flags, point_quadrature));
2259 Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2261
2262 update_internal_dofs(cell, static_cast<InternalData &>(*mdata));
2263
2265 p,
2266 point_quadrature,
2267 static_cast<InternalData &>(*mdata));
2268}
2269
2270
2271template <int dim, int spacedim, typename VectorType>
2275 const Point<spacedim> &p,
2276 Quadrature<dim> &point_quadrature,
2277 InternalData &mdata) const
2278{
2279 const unsigned int n_shapes = mdata.shape_values.size();
2280 (void)n_shapes;
2281 Assert(n_shapes != 0, ExcInternalError());
2282 AssertDimension(mdata.shape_derivatives.size(), n_shapes);
2283
2284
2285 // Newton iteration to solve
2286 // f(x)=p(x)-p=0
2287 // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
2288 // The start value was set to be the
2289 // linear approximation to the cell
2290 // The shape values and derivatives
2291 // of the mapping at this point are
2292 // previously computed.
2293
2294 AssertDimension(point_quadrature.size(), 1);
2295 Point<dim> p_unit = point_quadrature.point(0);
2296 Point<dim> f;
2297 mdata.reinit(mdata.update_each, point_quadrature);
2298
2300 Tensor<1, spacedim> p_minus_F = p - p_real;
2301 const double eps = 1.e-12 * cell->diameter();
2302 const unsigned int newton_iteration_limit = 20;
2303 unsigned int newton_iteration = 0;
2304 while (p_minus_F.norm_square() > eps * eps)
2305 {
2306 // f'(x)
2307 Point<spacedim> DF[dim];
2308 Tensor<2, dim> df;
2309 for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
2310 {
2311 const Tensor<1, dim> &grad_k = mdata.derivative(0, k);
2312 const unsigned int comp_k =
2313 euler_dof_handler->get_fe().system_to_component_index(k).first;
2314 if (fe_mask[comp_k])
2315 for (unsigned int j = 0; j < dim; ++j)
2316 DF[j][fe_to_real[comp_k]] +=
2317 mdata.local_dof_values[k] * grad_k[j];
2318 }
2319 for (unsigned int j = 0; j < dim; ++j)
2320 {
2321 f[j] = DF[j] * p_minus_F;
2322 for (unsigned int l = 0; l < dim; ++l)
2323 df[j][l] = -DF[j] * DF[l];
2324 }
2325 // Solve [f'(x)]d=f(x)
2326 const Tensor<1, dim> delta =
2327 invert(df) * static_cast<const Tensor<1, dim> &>(f);
2328 // do a line search
2329 double step_length = 1;
2330 do
2331 {
2332 // update of p_unit. The
2333 // spacedimth component of
2334 // transformed point is simply
2335 // ignored in codimension one
2336 // case. When this component is
2337 // not zero, then we are
2338 // projecting the point to the
2339 // surface or curve identified
2340 // by the cell.
2341 Point<dim> p_unit_trial = p_unit;
2342 for (unsigned int i = 0; i < dim; ++i)
2343 p_unit_trial[i] -= step_length * delta[i];
2344 // shape values and derivatives
2345 // at new p_unit point
2346 point_quadrature.initialize(
2347 ArrayView<const Point<dim>>(&p_unit_trial, 1));
2348 mdata.reinit(mdata.update_each, point_quadrature);
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< DerivativeForm< 1, dim, spacedim > > contravariant
std::vector< std::vector< Tensor< 1, spacedim > > > aux
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 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
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< 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
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
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
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
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
Point< dim > do_transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p, Quadrature< dim > &point_quadrature, 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
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:318
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 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:461
void initialize(const ArrayView< const Point< dim > > &points, const ArrayView< const double > &weights={})
Definition quadrature.cc:51
const Point< dim > & point(const unsigned int i) const
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:4624
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:79
@ mapping_piola
Definition mapping.h:114
@ mapping_covariant_gradient
Definition mapping.h:100
@ mapping_covariant
Definition mapping.h:89
@ mapping_contravariant
Definition mapping.h:94
#define DEAL_II_NOT_IMPLEMENTED()
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Definition mapping.cc:294
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 > &)