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