Reference documentation for deal.II version 9.4.1
\(\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 - 2022 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
16#include <deal.II/base/array_view.h>
17#include <deal.II/base/memory_consumption.h>
18#include <deal.II/base/polynomial.h>
19#include <deal.II/base/qprojector.h>
20#include <deal.II/base/quadrature.h>
21#include <deal.II/base/quadrature_lib.h>
22#include <deal.II/base/tensor_product_polynomials.h>
23#include <deal.II/base/utilities.h>
24
25#include <deal.II/dofs/dof_accessor.h>
26
27#include <deal.II/fe/fe_q.h>
28#include <deal.II/fe/fe_system.h>
29#include <deal.II/fe/fe_tools.h>
30#include <deal.II/fe/fe_values.h>
31#include <deal.II/fe/mapping.h>
32#include <deal.II/fe/mapping_fe_field.h>
33#include <deal.II/fe/mapping_q1.h>
34
35#include <deal.II/grid/tria_iterator.h>
36
37#include <deal.II/lac/block_vector.h>
38#include <deal.II/lac/full_matrix.h>
39#include <deal.II/lac/la_parallel_block_vector.h>
40#include <deal.II/lac/la_parallel_vector.h>
41#include <deal.II/lac/la_vector.h>
42#include <deal.II/lac/petsc_block_vector.h>
43#include <deal.II/lac/petsc_vector.h>
44#include <deal.II/lac/trilinos_epetra_vector.h>
45#include <deal.II/lac/trilinos_parallel_block_vector.h>
46#include <deal.II/lac/trilinos_tpetra_vector.h>
47#include <deal.II/lac/trilinos_vector.h>
48#include <deal.II/lac/vector.h>
49
50#include <deal.II/numerics/vector_tools.h>
51
52#include <fstream>
53#include <memory>
54#include <numeric>
55
56
57
58DEAL_II_NAMESPACE_OPEN
59
60
61template <int dim, int spacedim, typename VectorType>
62MappingFEField<dim, spacedim, VectorType, void>::InternalData::InternalData(
63 const FiniteElement<dim, spacedim> &fe,
64 const ComponentMask & mask)
65 : unit_tangentials()
66 , n_shape_functions(fe.n_dofs_per_cell())
67 , mask(mask)
68 , local_dof_indices(fe.n_dofs_per_cell())
69 , local_dof_values(fe.n_dofs_per_cell())
70{}
71
72
73
74template <int dim, int spacedim, typename VectorType>
75std::size_t
76MappingFEField<dim, spacedim, VectorType, void>::InternalData::
77 memory_consumption() const
78{
79 Assert(false, ExcNotImplemented());
80 return 0;
81}
82
83
84
85template <int dim, int spacedim, typename VectorType>
86double &
87MappingFEField<dim, spacedim, VectorType, void>::InternalData::shape(
88 const unsigned int qpoint,
89 const unsigned int shape_nr)
90{
91 AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
92 return shape_values[qpoint * n_shape_functions + shape_nr];
93}
94
95
96template <int dim, int spacedim, typename VectorType>
97const Tensor<1, dim> &
98MappingFEField<dim, spacedim, VectorType, void>::InternalData::derivative(
99 const unsigned int qpoint,
100 const unsigned int shape_nr) const
101{
102 AssertIndexRange(qpoint * n_shape_functions + shape_nr,
103 shape_derivatives.size());
104 return shape_derivatives[qpoint * n_shape_functions + shape_nr];
105}
106
107
108
109template <int dim, int spacedim, typename VectorType>
110Tensor<1, dim> &
111MappingFEField<dim, spacedim, VectorType, void>::InternalData::derivative(
112 const unsigned int qpoint,
113 const unsigned int shape_nr)
114{
115 AssertIndexRange(qpoint * n_shape_functions + shape_nr,
116 shape_derivatives.size());
117 return shape_derivatives[qpoint * n_shape_functions + shape_nr];
118}
119
120
121template <int dim, int spacedim, typename VectorType>
122const Tensor<2, dim> &
123MappingFEField<dim, spacedim, VectorType, void>::InternalData::
124 second_derivative(const unsigned int qpoint,
125 const unsigned int shape_nr) const
126{
127 AssertIndexRange(qpoint * n_shape_functions + shape_nr,
128 shape_second_derivatives.size());
129 return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
130}
131
132
133
134template <int dim, int spacedim, typename VectorType>
135Tensor<2, dim> &
136MappingFEField<dim, spacedim, VectorType, void>::InternalData::
137 second_derivative(const unsigned int qpoint, 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> &
147MappingFEField<dim, spacedim, VectorType, void>::InternalData::third_derivative(
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>
159Tensor<3, dim> &
160MappingFEField<dim, spacedim, VectorType, void>::InternalData::third_derivative(
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> &
172MappingFEField<dim, spacedim, VectorType, void>::InternalData::
173 fourth_derivative(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>
184Tensor<4, dim> &
185MappingFEField<dim, spacedim, VectorType, void>::InternalData::
186 fourth_derivative(const unsigned int qpoint, const unsigned int shape_nr)
187{
188 AssertIndexRange(qpoint * n_shape_functions + shape_nr,
189 shape_fourth_derivatives.size());
190 return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
191}
192
193
194
195template <int dim, int spacedim, typename VectorType>
196MappingFEField<dim, spacedim, VectorType, void>::MappingFEField(
197 const DoFHandler<dim, spacedim> &euler_dof_handler,
198 const VectorType & euler_vector,
199 const ComponentMask & mask)
200 : reference_cell(euler_dof_handler.get_fe().reference_cell())
201 , uses_level_dofs(false)
202 , euler_vector({&euler_vector})
203 , euler_dof_handler(&euler_dof_handler)
204 , fe_mask(mask.size() != 0u ?
205 mask :
206 ComponentMask(
207 euler_dof_handler.get_fe().get_nonzero_components(0).size(),
208 true))
209 , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
210 , fe_values(this->euler_dof_handler->get_fe(),
211 reference_cell.template get_nodal_type_quadrature<dim>(),
212 update_values)
213{
214 unsigned int size = 0;
215 for (unsigned int i = 0; i < fe_mask.size(); ++i)
216 {
217 if (fe_mask[i])
218 fe_to_real[i] = size++;
219 }
220 AssertDimension(size, spacedim);
221}
222
223
224
225template <int dim, int spacedim, typename VectorType>
226MappingFEField<dim, spacedim, VectorType, void>::MappingFEField(
227 const DoFHandler<dim, spacedim> &euler_dof_handler,
228 const std::vector<VectorType> & euler_vector,
229 const ComponentMask & mask)
230 : reference_cell(euler_dof_handler.get_fe().reference_cell())
231 , uses_level_dofs(true)
232 , euler_dof_handler(&euler_dof_handler)
233 , fe_mask(mask.size() != 0u ?
234 mask :
235 ComponentMask(
236 euler_dof_handler.get_fe().get_nonzero_components(0).size(),
237 true))
238 , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
239 , fe_values(this->euler_dof_handler->get_fe(),
240 reference_cell.template get_nodal_type_quadrature<dim>(),
241 update_values)
242{
243 unsigned int size = 0;
244 for (unsigned int i = 0; i < fe_mask.size(); ++i)
245 {
246 if (fe_mask[i])
247 fe_to_real[i] = size++;
248 }
249 AssertDimension(size, spacedim);
250
251 Assert(euler_dof_handler.has_level_dofs(),
252 ExcMessage("The underlying DoFHandler object did not call "
253 "distribute_mg_dofs(). In this case, the construction via "
254 "level vectors does not make sense."));
255 AssertDimension(euler_vector.size(),
256 euler_dof_handler.get_triangulation().n_global_levels());
257 this->euler_vector.clear();
258 this->euler_vector.resize(euler_vector.size());
259 for (unsigned int i = 0; i < euler_vector.size(); ++i)
260 this->euler_vector[i] = &euler_vector[i];
261}
262
263
264
265template <int dim, int spacedim, typename VectorType>
266MappingFEField<dim, spacedim, VectorType, void>::MappingFEField(
267 const DoFHandler<dim, spacedim> &euler_dof_handler,
268 const MGLevelObject<VectorType> &euler_vector,
269 const ComponentMask & mask)
270 : reference_cell(euler_dof_handler.get_fe().reference_cell())
271 , uses_level_dofs(true)
272 , euler_dof_handler(&euler_dof_handler)
273 , fe_mask(mask.size() != 0u ?
274 mask :
275 ComponentMask(
276 euler_dof_handler.get_fe().get_nonzero_components(0).size(),
277 true))
278 , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
279 , fe_values(this->euler_dof_handler->get_fe(),
280 reference_cell.template get_nodal_type_quadrature<dim>(),
281 update_values)
282{
283 unsigned int size = 0;
284 for (unsigned int i = 0; i < fe_mask.size(); ++i)
285 {
286 if (fe_mask[i])
287 fe_to_real[i] = size++;
288 }
289 AssertDimension(size, spacedim);
290
291 Assert(euler_dof_handler.has_level_dofs(),
292 ExcMessage("The underlying DoFHandler object did not call "
293 "distribute_mg_dofs(). In this case, the construction via "
294 "level vectors does not make sense."));
295 AssertDimension(euler_vector.max_level() + 1,
296 euler_dof_handler.get_triangulation().n_global_levels());
297 this->euler_vector.clear();
298 this->euler_vector.resize(
299 euler_dof_handler.get_triangulation().n_global_levels());
300 for (unsigned int i = euler_vector.min_level(); i <= euler_vector.max_level();
301 ++i)
302 this->euler_vector[i] = &euler_vector[i];
303}
304
305
306
307template <int dim, int spacedim, typename VectorType>
308MappingFEField<dim, spacedim, VectorType, void>::MappingFEField(
309 const MappingFEField<dim, spacedim, VectorType, void> &mapping)
310 : reference_cell(mapping.reference_cell)
311 , uses_level_dofs(mapping.uses_level_dofs)
312 , euler_vector(mapping.euler_vector)
313 , euler_dof_handler(mapping.euler_dof_handler)
314 , fe_mask(mapping.fe_mask)
315 , fe_to_real(mapping.fe_to_real)
316 , fe_values(mapping.euler_dof_handler->get_fe(),
317 reference_cell.template get_nodal_type_quadrature<dim>(),
318 update_values)
319{}
320
321
322
323template <int dim, int spacedim, typename VectorType>
324inline const double &
325MappingFEField<dim, spacedim, VectorType, void>::InternalData::shape(
326 const unsigned int qpoint,
327 const unsigned int shape_nr) const
328{
329 AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
330 return shape_values[qpoint * n_shape_functions + shape_nr];
331}
332
333
334
335template <int dim, int spacedim, typename VectorType>
336bool
337MappingFEField<dim, spacedim, VectorType, void>::preserves_vertex_locations()
338 const
339{
340 return false;
341}
342
343
344
345template <int dim, int spacedim, typename VectorType>
346bool
347MappingFEField<dim, spacedim, VectorType, void>::is_compatible_with(
348 const ReferenceCell &reference_cell) const
349{
350 Assert(dim == reference_cell.get_dimension(),
351 ExcMessage("The dimension of your mapping (" +
352 Utilities::to_string(dim) +
353 ") and the reference cell cell_type (" +
354 Utilities::to_string(reference_cell.get_dimension()) +
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>,
364 GeometryInfo<dim>::vertices_per_cell>
365MappingFEField<dim, spacedim, VectorType, void>::get_vertices(
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);
375 AssertDimension(fe_to_real.size(),
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>,
402 GeometryInfo<dim>::vertices_per_cell>
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];
408 if (comp != numbers::invalid_unsigned_int)
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
427MappingFEField<dim, spacedim, VectorType, void>::compute_shapes_virtual(
428 const std::vector<Point<dim>> &unit_points,
429 typename MappingFEField<dim, spacedim, VectorType, void>::InternalData &data)
430 const
431{
432 const auto fe = &euler_dof_handler->get_fe();
433 const unsigned int n_points = unit_points.size();
434
435 for (unsigned int point = 0; point < n_points; ++point)
436 {
437 if (data.shape_values.size() != 0)
438 for (unsigned int i = 0; i < data.n_shape_functions; ++i)
439 data.shape(point, i) = fe->shape_value(i, unit_points[point]);
440
441 if (data.shape_derivatives.size() != 0)
442 for (unsigned int i = 0; i < data.n_shape_functions; ++i)
443 data.derivative(point, i) = fe->shape_grad(i, unit_points[point]);
444
445 if (data.shape_second_derivatives.size() != 0)
446 for (unsigned int i = 0; i < data.n_shape_functions; ++i)
447 data.second_derivative(point, i) =
448 fe->shape_grad_grad(i, unit_points[point]);
449
450 if (data.shape_third_derivatives.size() != 0)
451 for (unsigned int i = 0; i < data.n_shape_functions; ++i)
452 data.third_derivative(point, i) =
453 fe->shape_3rd_derivative(i, unit_points[point]);
454
455 if (data.shape_fourth_derivatives.size() != 0)
456 for (unsigned int i = 0; i < data.n_shape_functions; ++i)
457 data.fourth_derivative(point, i) =
458 fe->shape_4th_derivative(i, unit_points[point]);
459 }
460}
461
462
463template <int dim, int spacedim, typename VectorType>
464UpdateFlags
465MappingFEField<dim, spacedim, VectorType, void>::requires_update_flags(
466 const UpdateFlags in) const
467{
468 // add flags if the respective quantities are necessary to compute
469 // what we need. note that some flags appear in both conditions and
470 // in subsequent set operations. this leads to some circular
471 // logic. the only way to treat this is to iterate. since there are
472 // 5 if-clauses in the loop, it will take at most 4 iterations to
473 // converge. do them:
474 UpdateFlags out = in;
475 for (unsigned int i = 0; i < 5; ++i)
476 {
477 // The following is a little incorrect:
478 // If not applied on a face,
479 // update_boundary_forms does not
480 // make sense. On the other hand,
481 // it is necessary on a
482 // face. Currently,
483 // update_boundary_forms is simply
484 // ignored for the interior of a
485 // cell.
486 if (out & (update_JxW_values | update_normal_vectors))
487 out |= update_boundary_forms;
488
489 if (out & (update_covariant_transformation | update_JxW_values |
490 update_jacobians | update_jacobian_grads |
491 update_boundary_forms | update_normal_vectors))
492 out |= update_contravariant_transformation;
493
494 if (out &
495 (update_inverse_jacobians | update_jacobian_pushed_forward_grads |
496 update_jacobian_pushed_forward_2nd_derivatives |
497 update_jacobian_pushed_forward_3rd_derivatives))
498 out |= update_covariant_transformation;
499
500 // The contravariant transformation
501 // is a Piola transformation, which
502 // requires the determinant of the
503 // Jacobi matrix of the transformation.
504 // Therefore these values have to be
505 // updated for each cell.
506 if (out & update_contravariant_transformation)
507 out |= update_JxW_values;
508
509 if (out & update_normal_vectors)
510 out |= update_JxW_values;
511 }
512
513 return out;
514}
515
516
517
518template <int dim, int spacedim, typename VectorType>
519void
520MappingFEField<dim, spacedim, VectorType, void>::compute_data(
521 const UpdateFlags update_flags,
522 const Quadrature<dim> &q,
523 const unsigned int n_original_q_points,
524 InternalData & data) const
525{
526 // store the flags in the internal data object so we can access them
527 // in fill_fe_*_values(). use the transitive hull of the required
528 // flags
529 data.update_each = requires_update_flags(update_flags);
530
531 const unsigned int n_q_points = q.size();
532
533 // see if we need the (transformation) shape function values
534 // and/or gradients and resize the necessary arrays
535 if (data.update_each & update_quadrature_points)
536 data.shape_values.resize(data.n_shape_functions * n_q_points);
537
538 if (data.update_each &
539 (update_covariant_transformation | update_contravariant_transformation |
540 update_JxW_values | update_boundary_forms | update_normal_vectors |
541 update_jacobians | update_jacobian_grads | update_inverse_jacobians))
542 data.shape_derivatives.resize(data.n_shape_functions * n_q_points);
543
544 if (data.update_each & update_covariant_transformation)
545 data.covariant.resize(n_original_q_points);
546
547 if (data.update_each & update_contravariant_transformation)
548 data.contravariant.resize(n_original_q_points);
549
550 if (data.update_each & update_volume_elements)
551 data.volume_elements.resize(n_original_q_points);
552
553 if (data.update_each &
554 (update_jacobian_grads | update_jacobian_pushed_forward_grads))
555 data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points);
556
557 if (data.update_each & (update_jacobian_2nd_derivatives |
558 update_jacobian_pushed_forward_2nd_derivatives))
559 data.shape_third_derivatives.resize(data.n_shape_functions * n_q_points);
560
561 if (data.update_each & (update_jacobian_3rd_derivatives |
562 update_jacobian_pushed_forward_3rd_derivatives))
563 data.shape_fourth_derivatives.resize(data.n_shape_functions * n_q_points);
564
565 compute_shapes_virtual(q.get_points(), data);
566
567 // This (for face values and simplices) can be different for different calls,
568 // so always copy
569 data.quadrature_weights = q.get_weights();
570}
571
572
573template <int dim, int spacedim, typename VectorType>
574void
575MappingFEField<dim, spacedim, VectorType, void>::compute_face_data(
576 const UpdateFlags update_flags,
577 const Quadrature<dim> &q,
578 const unsigned int n_original_q_points,
579 InternalData & data) const
580{
581 compute_data(update_flags, q, n_original_q_points, data);
582
583 if (dim > 1)
584 {
585 if (data.update_each & update_boundary_forms)
586 {
587 data.aux.resize(
588 dim - 1, std::vector<Tensor<1, spacedim>>(n_original_q_points));
589
590
591 // TODO: only a single reference cell type possible...
592 const auto n_faces = reference_cell.n_faces();
593
594 // Compute tangentials to the unit cell.
595 for (unsigned int i = 0; i < n_faces; ++i)
596 {
597 data.unit_tangentials[i].resize(n_original_q_points);
598 std::fill(
599 data.unit_tangentials[i].begin(),
600 data.unit_tangentials[i].end(),
601 reference_cell.template unit_tangential_vectors<dim>(i, 0));
602 if (dim > 2)
603 {
604 data.unit_tangentials[n_faces + i].resize(
605 n_original_q_points);
606 std::fill(
607 data.unit_tangentials[n_faces + i].begin(),
608 data.unit_tangentials[n_faces + i].end(),
609 reference_cell.template unit_tangential_vectors<dim>(i, 1));
610 }
611 }
612 }
613 }
614}
615
616
617template <int dim, int spacedim, typename VectorType>
618typename std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
619MappingFEField<dim, spacedim, VectorType, void>::get_data(
620 const UpdateFlags update_flags,
621 const Quadrature<dim> &quadrature) const
622{
623 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
624 std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
625 auto &data = dynamic_cast<InternalData &>(*data_ptr);
626 this->compute_data(update_flags, quadrature, quadrature.size(), data);
627
628 return data_ptr;
629}
630
631
632
633template <int dim, int spacedim, typename VectorType>
634std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
635MappingFEField<dim, spacedim, VectorType, void>::get_face_data(
636 const UpdateFlags update_flags,
637 const hp::QCollection<dim - 1> &quadrature) const
638{
639 AssertDimension(quadrature.size(), 1);
640
641 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
642 std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
643 auto & data = dynamic_cast<InternalData &>(*data_ptr);
644 const Quadrature<dim> q(
645 QProjector<dim>::project_to_all_faces(reference_cell, quadrature[0]));
646 this->compute_face_data(update_flags, q, quadrature[0].size(), data);
647
648 return data_ptr;
649}
650
651
652template <int dim, int spacedim, typename VectorType>
653std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
654MappingFEField<dim, spacedim, VectorType, void>::get_subface_data(
655 const UpdateFlags update_flags,
656 const Quadrature<dim - 1> &quadrature) const
657{
658 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
659 std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
660 auto & data = dynamic_cast<InternalData &>(*data_ptr);
661 const Quadrature<dim> q(
662 QProjector<dim>::project_to_all_subfaces(reference_cell, quadrature));
663 this->compute_face_data(update_flags, q, quadrature.size(), data);
664
665 return data_ptr;
666}
667
668
669
670namespace internal
671{
672 namespace MappingFEFieldImplementation
673 {
674 namespace
675 {
682 template <int dim, int spacedim, typename VectorType>
683 void
684 maybe_compute_q_points(
685 const typename ::QProjector<dim>::DataSetDescriptor data_set,
686 const typename ::MappingFEField<dim, spacedim, VectorType, void>::
687 InternalData & data,
688 const FiniteElement<dim, spacedim> &fe,
689 const ComponentMask & fe_mask,
690 const std::vector<unsigned int> & fe_to_real,
691 std::vector<Point<spacedim>> & quadrature_points)
692 {
693 const UpdateFlags update_flags = data.update_each;
694
695 if (update_flags & update_quadrature_points)
696 {
697 for (unsigned int point = 0; point < quadrature_points.size();
698 ++point)
699 {
700 Point<spacedim> result;
701 const double * shape = &data.shape(point + data_set, 0);
702
703 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
704 {
705 const unsigned int comp_k =
706 fe.system_to_component_index(k).first;
707 if (fe_mask[comp_k])
708 result[fe_to_real[comp_k]] +=
709 data.local_dof_values[k] * shape[k];
710 }
711
712 quadrature_points[point] = result;
713 }
714 }
715 }
716
724 template <int dim, int spacedim, typename VectorType>
725 void
726 maybe_update_Jacobians(
727 const typename ::QProjector<dim>::DataSetDescriptor data_set,
728 const typename ::MappingFEField<dim, spacedim, VectorType, void>::
729 InternalData & data,
730 const FiniteElement<dim, spacedim> &fe,
731 const ComponentMask & fe_mask,
732 const std::vector<unsigned int> & fe_to_real)
733 {
734 const UpdateFlags update_flags = data.update_each;
735
736 // then Jacobians
737 if (update_flags & update_contravariant_transformation)
738 {
739 const unsigned int n_q_points = data.contravariant.size();
740
741 Assert(data.n_shape_functions > 0, ExcInternalError());
742
743 for (unsigned int point = 0; point < n_q_points; ++point)
744 {
745 const Tensor<1, dim> *data_derv =
746 &data.derivative(point + data_set, 0);
747
748 Tensor<1, dim> result[spacedim];
749
750 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
751 {
752 const unsigned int comp_k =
753 fe.system_to_component_index(k).first;
754 if (fe_mask[comp_k])
755 result[fe_to_real[comp_k]] +=
756 data.local_dof_values[k] * data_derv[k];
757 }
758
759 // write result into contravariant data
760 for (unsigned int i = 0; i < spacedim; ++i)
761 {
762 data.contravariant[point][i] = result[i];
763 }
764 }
765 }
766
767 if (update_flags & update_covariant_transformation)
768 {
769 AssertDimension(data.covariant.size(), data.contravariant.size());
770 for (unsigned int point = 0; point < data.contravariant.size();
771 ++point)
772 data.covariant[point] =
773 (data.contravariant[point]).covariant_form();
774 }
775
776 if (update_flags & update_volume_elements)
777 {
778 AssertDimension(data.contravariant.size(),
779 data.volume_elements.size());
780 for (unsigned int point = 0; point < data.contravariant.size();
781 ++point)
782 data.volume_elements[point] =
783 data.contravariant[point].determinant();
784 }
785 }
786
793 template <int dim, int spacedim, typename VectorType>
794 void
795 maybe_update_jacobian_grads(
796 const typename ::QProjector<dim>::DataSetDescriptor data_set,
797 const typename ::MappingFEField<dim, spacedim, VectorType, void>::
798 InternalData & data,
799 const FiniteElement<dim, spacedim> & fe,
800 const ComponentMask & fe_mask,
801 const std::vector<unsigned int> & fe_to_real,
802 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
803 {
804 const UpdateFlags update_flags = data.update_each;
805 if (update_flags & update_jacobian_grads)
806 {
807 const unsigned int n_q_points = jacobian_grads.size();
808
809 for (unsigned int point = 0; point < n_q_points; ++point)
810 {
811 const Tensor<2, dim> *second =
812 &data.second_derivative(point + data_set, 0);
813
814 DerivativeForm<2, dim, spacedim> result;
815
816 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
817 {
818 const unsigned int comp_k =
819 fe.system_to_component_index(k).first;
820 if (fe_mask[comp_k])
821 for (unsigned int j = 0; j < dim; ++j)
822 for (unsigned int l = 0; l < dim; ++l)
823 result[fe_to_real[comp_k]][j][l] +=
824 (second[k][j][l] * data.local_dof_values[k]);
825 }
826
827 // never touch any data for j=dim in case dim<spacedim, so
828 // it will always be zero as it was initialized
829 for (unsigned int i = 0; i < spacedim; ++i)
830 for (unsigned int j = 0; j < dim; ++j)
831 for (unsigned int l = 0; l < dim; ++l)
832 jacobian_grads[point][i][j][l] = result[i][j][l];
833 }
834 }
835 }
836
843 template <int dim, int spacedim, typename VectorType>
844 void
845 maybe_update_jacobian_pushed_forward_grads(
846 const typename ::QProjector<dim>::DataSetDescriptor data_set,
847 const typename ::MappingFEField<dim, spacedim, VectorType, void>::
848 InternalData & data,
849 const FiniteElement<dim, spacedim> &fe,
850 const ComponentMask & fe_mask,
851 const std::vector<unsigned int> & fe_to_real,
852 std::vector<Tensor<3, spacedim>> & jacobian_pushed_forward_grads)
853 {
854 const UpdateFlags update_flags = data.update_each;
855 if (update_flags & update_jacobian_pushed_forward_grads)
856 {
857 const unsigned int n_q_points =
858 jacobian_pushed_forward_grads.size();
859
860 double tmp[spacedim][spacedim][spacedim];
861 for (unsigned int point = 0; point < n_q_points; ++point)
862 {
863 const Tensor<2, dim> *second =
864 &data.second_derivative(point + data_set, 0);
865
866 DerivativeForm<2, dim, spacedim> result;
867
868 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
869 {
870 const unsigned int comp_k =
871 fe.system_to_component_index(k).first;
872 if (fe_mask[comp_k])
873 for (unsigned int j = 0; j < dim; ++j)
874 for (unsigned int l = 0; l < dim; ++l)
875 result[fe_to_real[comp_k]][j][l] +=
876 (second[k][j][l] * data.local_dof_values[k]);
877 }
878
879 // first push forward the j-components
880 for (unsigned int i = 0; i < spacedim; ++i)
881 for (unsigned int j = 0; j < spacedim; ++j)
882 for (unsigned int l = 0; l < dim; ++l)
883 {
884 tmp[i][j][l] =
885 result[i][0][l] * data.covariant[point][j][0];
886 for (unsigned int jr = 1; jr < dim; ++jr)
887 {
888 tmp[i][j][l] +=
889 result[i][jr][l] * data.covariant[point][j][jr];
890 }
891 }
892
893 // now, pushing forward the l-components
894 for (unsigned int i = 0; i < spacedim; ++i)
895 for (unsigned int j = 0; j < spacedim; ++j)
896 for (unsigned int l = 0; l < spacedim; ++l)
897 {
898 jacobian_pushed_forward_grads[point][i][j][l] =
899 tmp[i][j][0] * data.covariant[point][l][0];
900 for (unsigned int lr = 1; lr < dim; ++lr)
901 {
902 jacobian_pushed_forward_grads[point][i][j][l] +=
903 tmp[i][j][lr] * data.covariant[point][l][lr];
904 }
905 }
906 }
907 }
908 }
909
916 template <int dim, int spacedim, typename VectorType>
917 void
918 maybe_update_jacobian_2nd_derivatives(
919 const typename ::QProjector<dim>::DataSetDescriptor data_set,
920 const typename ::MappingFEField<dim, spacedim, VectorType, void>::
921 InternalData & data,
922 const FiniteElement<dim, spacedim> & fe,
923 const ComponentMask & fe_mask,
924 const std::vector<unsigned int> & fe_to_real,
925 std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
926 {
927 const UpdateFlags update_flags = data.update_each;
928 if (update_flags & update_jacobian_2nd_derivatives)
929 {
930 const unsigned int n_q_points = jacobian_2nd_derivatives.size();
931
932 for (unsigned int point = 0; point < n_q_points; ++point)
933 {
934 const Tensor<3, dim> *third =
935 &data.third_derivative(point + data_set, 0);
936
937 DerivativeForm<3, dim, spacedim> result;
938
939 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
940 {
941 const unsigned int comp_k =
942 fe.system_to_component_index(k).first;
943 if (fe_mask[comp_k])
944 for (unsigned int j = 0; j < dim; ++j)
945 for (unsigned int l = 0; l < dim; ++l)
946 for (unsigned int m = 0; m < dim; ++m)
947 result[fe_to_real[comp_k]][j][l][m] +=
948 (third[k][j][l][m] * data.local_dof_values[k]);
949 }
950
951 // never touch any data for j=dim in case dim<spacedim, so
952 // it will always be zero as it was initialized
953 for (unsigned int i = 0; i < spacedim; ++i)
954 for (unsigned int j = 0; j < dim; ++j)
955 for (unsigned int l = 0; l < dim; ++l)
956 for (unsigned int m = 0; m < dim; ++m)
957 jacobian_2nd_derivatives[point][i][j][l][m] =
958 result[i][j][l][m];
959 }
960 }
961 }
962
970 template <int dim, int spacedim, typename VectorType>
971 void
972 maybe_update_jacobian_pushed_forward_2nd_derivatives(
973 const typename ::QProjector<dim>::DataSetDescriptor data_set,
974 const typename ::MappingFEField<dim, spacedim, VectorType, void>::
975 InternalData & data,
976 const FiniteElement<dim, spacedim> &fe,
977 const ComponentMask & fe_mask,
978 const std::vector<unsigned int> & fe_to_real,
979 std::vector<Tensor<4, spacedim>>
980 &jacobian_pushed_forward_2nd_derivatives)
981 {
982 const UpdateFlags update_flags = data.update_each;
983 if (update_flags & update_jacobian_pushed_forward_2nd_derivatives)
984 {
985 const unsigned int n_q_points =
986 jacobian_pushed_forward_2nd_derivatives.size();
987
988 double tmp[spacedim][spacedim][spacedim][spacedim];
989 for (unsigned int point = 0; point < n_q_points; ++point)
990 {
991 const Tensor<3, dim> *third =
992 &data.third_derivative(point + data_set, 0);
993
994 DerivativeForm<3, dim, spacedim> result;
995
996 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
997 {
998 const unsigned int comp_k =
999 fe.system_to_component_index(k).first;
1000 if (fe_mask[comp_k])
1001 for (unsigned int j = 0; j < dim; ++j)
1002 for (unsigned int l = 0; l < dim; ++l)
1003 for (unsigned int m = 0; m < dim; ++m)
1004 result[fe_to_real[comp_k]][j][l][m] +=
1005 (third[k][j][l][m] * data.local_dof_values[k]);
1006 }
1007
1008 // push forward the j-coordinate
1009 for (unsigned int i = 0; i < spacedim; ++i)
1010 for (unsigned int j = 0; j < spacedim; ++j)
1011 for (unsigned int l = 0; l < dim; ++l)
1012 for (unsigned int m = 0; m < dim; ++m)
1013 {
1014 jacobian_pushed_forward_2nd_derivatives
1015 [point][i][j][l][m] =
1016 result[i][0][l][m] * data.covariant[point][j][0];
1017 for (unsigned int jr = 1; jr < dim; ++jr)
1018 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1019 [l][m] +=
1020 result[i][jr][l][m] *
1021 data.covariant[point][j][jr];
1022 }
1023
1024 // push forward the l-coordinate
1025 for (unsigned int i = 0; i < spacedim; ++i)
1026 for (unsigned int j = 0; j < spacedim; ++j)
1027 for (unsigned int l = 0; l < spacedim; ++l)
1028 for (unsigned int m = 0; m < dim; ++m)
1029 {
1030 tmp[i][j][l][m] =
1031 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1032 [0][m] *
1033 data.covariant[point][l][0];
1034 for (unsigned int lr = 1; lr < dim; ++lr)
1035 tmp[i][j][l][m] +=
1036 jacobian_pushed_forward_2nd_derivatives[point][i]
1037 [j][lr]
1038 [m] *
1039 data.covariant[point][l][lr];
1040 }
1041
1042 // push forward the m-coordinate
1043 for (unsigned int i = 0; i < spacedim; ++i)
1044 for (unsigned int j = 0; j < spacedim; ++j)
1045 for (unsigned int l = 0; l < spacedim; ++l)
1046 for (unsigned int m = 0; m < spacedim; ++m)
1047 {
1048 jacobian_pushed_forward_2nd_derivatives
1049 [point][i][j][l][m] =
1050 tmp[i][j][l][0] * data.covariant[point][m][0];
1051 for (unsigned int mr = 1; mr < dim; ++mr)
1052 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1053 [l][m] +=
1054 tmp[i][j][l][mr] * data.covariant[point][m][mr];
1055 }
1056 }
1057 }
1058 }
1059
1066 template <int dim, int spacedim, typename VectorType>
1067 void
1068 maybe_update_jacobian_3rd_derivatives(
1069 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1070 const typename ::MappingFEField<dim, spacedim, VectorType, void>::
1071 InternalData & data,
1072 const FiniteElement<dim, spacedim> & fe,
1073 const ComponentMask & fe_mask,
1074 const std::vector<unsigned int> & fe_to_real,
1075 std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1076 {
1077 const UpdateFlags update_flags = data.update_each;
1078 if (update_flags & update_jacobian_3rd_derivatives)
1079 {
1080 const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1081
1082 for (unsigned int point = 0; point < n_q_points; ++point)
1083 {
1084 const Tensor<4, dim> *fourth =
1085 &data.fourth_derivative(point + data_set, 0);
1086
1087 DerivativeForm<4, dim, spacedim> result;
1088
1089 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1090 {
1091 const unsigned int comp_k =
1092 fe.system_to_component_index(k).first;
1093 if (fe_mask[comp_k])
1094 for (unsigned int j = 0; j < dim; ++j)
1095 for (unsigned int l = 0; l < dim; ++l)
1096 for (unsigned int m = 0; m < dim; ++m)
1097 for (unsigned int n = 0; n < dim; ++n)
1098 result[fe_to_real[comp_k]][j][l][m][n] +=
1099 (fourth[k][j][l][m][n] *
1100 data.local_dof_values[k]);
1101 }
1102
1103 // never touch any data for j,l,m,n=dim in case
1104 // dim<spacedim, so it will always be zero as it was
1105 // initialized
1106 for (unsigned int i = 0; i < spacedim; ++i)
1107 for (unsigned int j = 0; j < dim; ++j)
1108 for (unsigned int l = 0; l < dim; ++l)
1109 for (unsigned int m = 0; m < dim; ++m)
1110 for (unsigned int n = 0; n < dim; ++n)
1111 jacobian_3rd_derivatives[point][i][j][l][m][n] =
1112 result[i][j][l][m][n];
1113 }
1114 }
1115 }
1116
1124 template <int dim, int spacedim, typename VectorType>
1125 void
1126 maybe_update_jacobian_pushed_forward_3rd_derivatives(
1127 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1128 const typename ::MappingFEField<dim, spacedim, VectorType, void>::
1129 InternalData & data,
1130 const FiniteElement<dim, spacedim> &fe,
1131 const ComponentMask & fe_mask,
1132 const std::vector<unsigned int> & fe_to_real,
1133 std::vector<Tensor<5, spacedim>>
1134 &jacobian_pushed_forward_3rd_derivatives)
1135 {
1136 const UpdateFlags update_flags = data.update_each;
1137 if (update_flags & update_jacobian_pushed_forward_3rd_derivatives)
1138 {
1139 const unsigned int n_q_points =
1140 jacobian_pushed_forward_3rd_derivatives.size();
1141
1142 double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
1143 for (unsigned int point = 0; point < n_q_points; ++point)
1144 {
1145 const Tensor<4, dim> *fourth =
1146 &data.fourth_derivative(point + data_set, 0);
1147
1148 DerivativeForm<4, dim, spacedim> result;
1149
1150 for (unsigned int k = 0; k < data.n_shape_functions; ++k)
1151 {
1152 const unsigned int comp_k =
1153 fe.system_to_component_index(k).first;
1154 if (fe_mask[comp_k])
1155 for (unsigned int j = 0; j < dim; ++j)
1156 for (unsigned int l = 0; l < dim; ++l)
1157 for (unsigned int m = 0; m < dim; ++m)
1158 for (unsigned int n = 0; n < dim; ++n)
1159 result[fe_to_real[comp_k]][j][l][m][n] +=
1160 (fourth[k][j][l][m][n] *
1161 data.local_dof_values[k]);
1162 }
1163
1164 // push-forward the j-coordinate
1165 for (unsigned int i = 0; i < spacedim; ++i)
1166 for (unsigned int j = 0; j < spacedim; ++j)
1167 for (unsigned int l = 0; l < dim; ++l)
1168 for (unsigned int m = 0; m < dim; ++m)
1169 for (unsigned int n = 0; n < dim; ++n)
1170 {
1171 tmp[i][j][l][m][n] = result[i][0][l][m][n] *
1172 data.covariant[point][j][0];
1173 for (unsigned int jr = 1; jr < dim; ++jr)
1174 tmp[i][j][l][m][n] +=
1175 result[i][jr][l][m][n] *
1176 data.covariant[point][j][jr];
1177 }
1178
1179 // push-forward the l-coordinate
1180 for (unsigned int i = 0; i < spacedim; ++i)
1181 for (unsigned int j = 0; j < spacedim; ++j)
1182 for (unsigned int l = 0; l < spacedim; ++l)
1183 for (unsigned int m = 0; m < dim; ++m)
1184 for (unsigned int n = 0; n < dim; ++n)
1185 {
1186 jacobian_pushed_forward_3rd_derivatives
1187 [point][i][j][l][m][n] =
1188 tmp[i][j][0][m][n] *
1189 data.covariant[point][l][0];
1190 for (unsigned int lr = 1; lr < dim; ++lr)
1191 jacobian_pushed_forward_3rd_derivatives[point][i]
1192 [j][l][m]
1193 [n] +=
1194 tmp[i][j][lr][m][n] *
1195 data.covariant[point][l][lr];
1196 }
1197
1198 // push-forward the m-coordinate
1199 for (unsigned int i = 0; i < spacedim; ++i)
1200 for (unsigned int j = 0; j < spacedim; ++j)
1201 for (unsigned int l = 0; l < spacedim; ++l)
1202 for (unsigned int m = 0; m < spacedim; ++m)
1203 for (unsigned int n = 0; n < dim; ++n)
1204 {
1205 tmp[i][j][l][m][n] =
1206 jacobian_pushed_forward_3rd_derivatives[point][i]
1207 [j][l][0]
1208 [n] *
1209 data.covariant[point][m][0];
1210 for (unsigned int mr = 1; mr < dim; ++mr)
1211 tmp[i][j][l][m][n] +=
1212 jacobian_pushed_forward_3rd_derivatives[point]
1213 [i][j][l]
1214 [mr][n] *
1215 data.covariant[point][m][mr];
1216 }
1217
1218 // push-forward the n-coordinate
1219 for (unsigned int i = 0; i < spacedim; ++i)
1220 for (unsigned int j = 0; j < spacedim; ++j)
1221 for (unsigned int l = 0; l < spacedim; ++l)
1222 for (unsigned int m = 0; m < spacedim; ++m)
1223 for (unsigned int n = 0; n < spacedim; ++n)
1224 {
1225 jacobian_pushed_forward_3rd_derivatives
1226 [point][i][j][l][m][n] =
1227 tmp[i][j][l][m][0] *
1228 data.covariant[point][n][0];
1229 for (unsigned int nr = 1; nr < dim; ++nr)
1230 jacobian_pushed_forward_3rd_derivatives[point][i]
1231 [j][l][m]
1232 [n] +=
1233 tmp[i][j][l][m][nr] *
1234 data.covariant[point][n][nr];
1235 }
1236 }
1237 }
1238 }
1239
1240
1250 template <int dim, int spacedim, typename VectorType>
1251 void
1252 maybe_compute_face_data(
1253 const ::Mapping<dim, spacedim> &mapping,
1254 const typename ::Triangulation<dim, spacedim>::cell_iterator
1255 & cell,
1256 const unsigned int face_no,
1257 const unsigned int subface_no,
1258 const typename QProjector<dim>::DataSetDescriptor data_set,
1259 const typename ::MappingFEField<dim, spacedim, VectorType, void>::
1260 InternalData &data,
1261 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1262 &output_data)
1263 {
1264 const UpdateFlags update_flags = data.update_each;
1265
1266 if (update_flags & update_boundary_forms)
1267 {
1268 const unsigned int n_q_points = output_data.boundary_forms.size();
1269 if (update_flags & update_normal_vectors)
1270 AssertDimension(output_data.normal_vectors.size(), n_q_points);
1271 if (update_flags & update_JxW_values)
1272 AssertDimension(output_data.JxW_values.size(), n_q_points);
1273
1274 // map the unit tangentials to the real cell. checking for d!=dim-1
1275 // eliminates compiler warnings regarding unsigned int expressions <
1276 // 0.
1277 for (unsigned int d = 0; d != dim - 1; ++d)
1278 {
1279 Assert(face_no + cell->n_faces() * d <
1280 data.unit_tangentials.size(),
1281 ExcInternalError());
1282 Assert(
1283 data.aux[d].size() <=
1284 data.unit_tangentials[face_no + cell->n_faces() * d].size(),
1285 ExcInternalError());
1286
1287 mapping.transform(
1288 make_array_view(
1289 data.unit_tangentials[face_no + cell->n_faces() * d]),
1290 mapping_contravariant,
1291 data,
1292 make_array_view(data.aux[d]));
1293 }
1294
1295 // if dim==spacedim, we can use the unit tangentials to compute the
1296 // boundary form by simply taking the cross product
1297 if (dim == spacedim)
1298 {
1299 for (unsigned int i = 0; i < n_q_points; ++i)
1300 switch (dim)
1301 {
1302 case 1:
1303 // in 1d, we don't have access to any of the data.aux
1304 // fields (because it has only dim-1 components), but we
1305 // can still compute the boundary form by simply looking
1306 // at the number of the face
1307 output_data.boundary_forms[i][0] =
1308 (face_no == 0 ? -1 : +1);
1309 break;
1310 case 2:
1311 output_data.boundary_forms[i] =
1312 cross_product_2d(data.aux[0][i]);
1313 break;
1314 case 3:
1315 output_data.boundary_forms[i] =
1316 cross_product_3d(data.aux[0][i], data.aux[1][i]);
1317 break;
1318 default:
1319 Assert(false, ExcNotImplemented());
1320 }
1321 }
1322 else //(dim < spacedim)
1323 {
1324 // in the codim-one case, the boundary form results from the
1325 // cross product of all the face tangential vectors and the cell
1326 // normal vector
1327 //
1328 // to compute the cell normal, use the same method used in
1329 // fill_fe_values for cells above
1330 AssertDimension(data.contravariant.size(), n_q_points);
1331
1332 for (unsigned int point = 0; point < n_q_points; ++point)
1333 {
1334 if (dim == 1)
1335 {
1336 // J is a tangent vector
1337 output_data.boundary_forms[point] =
1338 data.contravariant[point].transpose()[0];
1339 output_data.boundary_forms[point] /=
1340 (face_no == 0 ? -1. : +1.) *
1341 output_data.boundary_forms[point].norm();
1342 }
1343
1344 if (dim == 2)
1345 {
1346 const DerivativeForm<1, spacedim, dim> DX_t =
1347 data.contravariant[point].transpose();
1348
1349 Tensor<1, spacedim> cell_normal =
1350 cross_product_3d(DX_t[0], DX_t[1]);
1351 cell_normal /= cell_normal.norm();
1352
1353 // then compute the face normal from the face tangent
1354 // and the cell normal:
1355 output_data.boundary_forms[point] =
1356 cross_product_3d(data.aux[0][point], cell_normal);
1357 }
1358 }
1359 }
1360
1361 if (update_flags & (update_normal_vectors | update_JxW_values))
1362 for (unsigned int i = 0; i < output_data.boundary_forms.size();
1363 ++i)
1364 {
1365 if (update_flags & update_JxW_values)
1366 {
1367 output_data.JxW_values[i] =
1368 output_data.boundary_forms[i].norm() *
1369 data.quadrature_weights[i + data_set];
1370
1371 if (subface_no != numbers::invalid_unsigned_int)
1372 {
1373 // TODO
1374 const double area_ratio =
1375 GeometryInfo<dim>::subface_ratio(
1376 cell->subface_case(face_no), subface_no);
1377 output_data.JxW_values[i] *= area_ratio;
1378 }
1379 }
1380
1381 if (update_flags & update_normal_vectors)
1382 output_data.normal_vectors[i] =
1383 Point<spacedim>(output_data.boundary_forms[i] /
1384 output_data.boundary_forms[i].norm());
1385 }
1386
1387 if (update_flags & update_jacobians)
1388 for (unsigned int point = 0; point < n_q_points; ++point)
1389 output_data.jacobians[point] = data.contravariant[point];
1390
1391 if (update_flags & update_inverse_jacobians)
1392 for (unsigned int point = 0; point < n_q_points; ++point)
1393 output_data.inverse_jacobians[point] =
1394 data.covariant[point].transpose();
1395 }
1396 }
1397
1404 template <int dim, int spacedim, typename VectorType>
1405 void
1406 do_fill_fe_face_values(
1407 const ::Mapping<dim, spacedim> &mapping,
1408 const typename ::Triangulation<dim, spacedim>::cell_iterator
1409 & cell,
1410 const unsigned int face_no,
1411 const unsigned int subface_no,
1412 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1413 const typename ::MappingFEField<dim, spacedim, VectorType, void>::
1414 InternalData & data,
1415 const FiniteElement<dim, spacedim> &fe,
1416 const ComponentMask & fe_mask,
1417 const std::vector<unsigned int> & fe_to_real,
1418 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1419 &output_data)
1420 {
1421 maybe_compute_q_points<dim, spacedim, VectorType>(
1422 data_set,
1423 data,
1424 fe,
1425 fe_mask,
1426 fe_to_real,
1427 output_data.quadrature_points);
1428
1429 maybe_update_Jacobians<dim, spacedim, VectorType>(
1430 data_set, data, fe, fe_mask, fe_to_real);
1431
1432 maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1433 data_set, data, fe, fe_mask, fe_to_real, output_data.jacobian_grads);
1434
1435 maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1436 data_set,
1437 data,
1438 fe,
1439 fe_mask,
1440 fe_to_real,
1441 output_data.jacobian_pushed_forward_grads);
1442
1443 maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1444 data_set,
1445 data,
1446 fe,
1447 fe_mask,
1448 fe_to_real,
1449 output_data.jacobian_2nd_derivatives);
1450
1451 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1452 spacedim,
1453 VectorType>(
1454 data_set,
1455 data,
1456 fe,
1457 fe_mask,
1458 fe_to_real,
1459 output_data.jacobian_pushed_forward_2nd_derivatives);
1460
1461 maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1462 data_set,
1463 data,
1464 fe,
1465 fe_mask,
1466 fe_to_real,
1467 output_data.jacobian_3rd_derivatives);
1468
1469 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1470 spacedim,
1471 VectorType>(
1472 data_set,
1473 data,
1474 fe,
1475 fe_mask,
1476 fe_to_real,
1477 output_data.jacobian_pushed_forward_3rd_derivatives);
1478
1479 maybe_compute_face_data<dim, spacedim, VectorType>(
1480 mapping, cell, face_no, subface_no, data_set, data, output_data);
1481 }
1482 } // namespace
1483 } // namespace MappingFEFieldImplementation
1484} // namespace internal
1485
1486
1487// Note that the CellSimilarity flag is modifiable, since MappingFEField can
1488// need to recalculate data even when cells are similar.
1489template <int dim, int spacedim, typename VectorType>
1490CellSimilarity::Similarity
1491MappingFEField<dim, spacedim, VectorType, void>::fill_fe_values(
1492 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1493 const CellSimilarity::Similarity,
1494 const Quadrature<dim> & quadrature,
1495 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
1496 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1497 &output_data) const
1498{
1499 // convert data object to internal data for this class. fails with an
1500 // exception if that is not possible
1501 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1502 ExcInternalError());
1503 const InternalData &data = static_cast<const InternalData &>(internal_data);
1504
1505 const unsigned int n_q_points = quadrature.size();
1506
1507 update_internal_dofs(cell, data);
1508
1509 internal::MappingFEFieldImplementation::
1510 maybe_compute_q_points<dim, spacedim, VectorType>(
1511 QProjector<dim>::DataSetDescriptor::cell(),
1512 data,
1513 euler_dof_handler->get_fe(),
1514 fe_mask,
1515 fe_to_real,
1516 output_data.quadrature_points);
1517
1518 internal::MappingFEFieldImplementation::
1519 maybe_update_Jacobians<dim, spacedim, VectorType>(
1520 QProjector<dim>::DataSetDescriptor::cell(),
1521 data,
1522 euler_dof_handler->get_fe(),
1523 fe_mask,
1524 fe_to_real);
1525
1526 const UpdateFlags update_flags = data.update_each;
1527 const std::vector<double> &weights = quadrature.get_weights();
1528
1529 // Multiply quadrature weights by absolute value of Jacobian determinants or
1530 // the area element g=sqrt(DX^t DX) in case of codim > 0
1531
1532 if (update_flags & (update_normal_vectors | update_JxW_values))
1533 {
1534 AssertDimension(output_data.JxW_values.size(), n_q_points);
1535
1536 Assert(!(update_flags & update_normal_vectors) ||
1537 (output_data.normal_vectors.size() == n_q_points),
1538 ExcDimensionMismatch(output_data.normal_vectors.size(),
1539 n_q_points));
1540
1541
1542 for (unsigned int point = 0; point < n_q_points; ++point)
1543 {
1544 if (dim == spacedim)
1545 {
1546 const double det = data.contravariant[point].determinant();
1547
1548 // check for distorted cells.
1549
1550 // TODO: this allows for anisotropies of up to 1e6 in 3D and
1551 // 1e12 in 2D. might want to find a finer
1552 // (dimension-independent) criterion
1553 Assert(det > 1e-12 * Utilities::fixed_power<dim>(
1554 cell->diameter() / std::sqrt(double(dim))),
1555 (typename Mapping<dim, spacedim>::ExcDistortedMappedCell(
1556 cell->center(), det, point)));
1557 output_data.JxW_values[point] = weights[point] * det;
1558 }
1559 // if dim==spacedim, then there is no cell normal to
1560 // compute. since this is for FEValues (and not FEFaceValues),
1561 // there are also no face normals to compute
1562 else // codim>0 case
1563 {
1564 Tensor<1, spacedim> DX_t[dim];
1565 for (unsigned int i = 0; i < spacedim; ++i)
1566 for (unsigned int j = 0; j < dim; ++j)
1567 DX_t[j][i] = data.contravariant[point][i][j];
1568
1569 Tensor<2, dim> G; // First fundamental form
1570 for (unsigned int i = 0; i < dim; ++i)
1571 for (unsigned int j = 0; j < dim; ++j)
1572 G[i][j] = DX_t[i] * DX_t[j];
1573
1574 output_data.JxW_values[point] =
1575 std::sqrt(determinant(G)) * weights[point];
1576
1577 if (update_flags & update_normal_vectors)
1578 {
1579 Assert(spacedim - dim == 1,
1580 ExcMessage("There is no cell normal in codim 2."));
1581
1582 if (dim == 1)
1583 output_data.normal_vectors[point] =
1584 cross_product_2d(-DX_t[0]);
1585 else
1586 {
1587 Assert(dim == 2, ExcInternalError());
1588
1589 // dim-1==1 for the second argument, but this
1590 // avoids a compiler warning about array bounds:
1591 output_data.normal_vectors[point] =
1592 cross_product_3d(DX_t[0], DX_t[dim - 1]);
1593 }
1594
1595 output_data.normal_vectors[point] /=
1596 output_data.normal_vectors[point].norm();
1597
1598 if (cell->direction_flag() == false)
1599 output_data.normal_vectors[point] *= -1.;
1600 }
1601 } // codim>0 case
1602 }
1603 }
1604
1605 // copy values from InternalData to vector given by reference
1606 if (update_flags & update_jacobians)
1607 {
1608 AssertDimension(output_data.jacobians.size(), n_q_points);
1609 for (unsigned int point = 0; point < n_q_points; ++point)
1610 output_data.jacobians[point] = data.contravariant[point];
1611 }
1612
1613 // copy values from InternalData to vector given by reference
1614 if (update_flags & update_inverse_jacobians)
1615 {
1616 AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
1617 for (unsigned int point = 0; point < n_q_points; ++point)
1618 output_data.inverse_jacobians[point] =
1619 data.covariant[point].transpose();
1620 }
1621
1622 // calculate derivatives of the Jacobians
1623 internal::MappingFEFieldImplementation::
1624 maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1625 QProjector<dim>::DataSetDescriptor::cell(),
1626 data,
1627 euler_dof_handler->get_fe(),
1628 fe_mask,
1629 fe_to_real,
1630 output_data.jacobian_grads);
1631
1632 // calculate derivatives of the Jacobians pushed forward to real cell
1633 // coordinates
1634 internal::MappingFEFieldImplementation::
1635 maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1636 QProjector<dim>::DataSetDescriptor::cell(),
1637 data,
1638 euler_dof_handler->get_fe(),
1639 fe_mask,
1640 fe_to_real,
1641 output_data.jacobian_pushed_forward_grads);
1642
1643 // calculate hessians of the Jacobians
1644 internal::MappingFEFieldImplementation::
1645 maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1646 QProjector<dim>::DataSetDescriptor::cell(),
1647 data,
1648 euler_dof_handler->get_fe(),
1649 fe_mask,
1650 fe_to_real,
1651 output_data.jacobian_2nd_derivatives);
1652
1653 // calculate hessians of the Jacobians pushed forward to real cell coordinates
1654 internal::MappingFEFieldImplementation::
1655 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1656 spacedim,
1657 VectorType>(
1658 QProjector<dim>::DataSetDescriptor::cell(),
1659 data,
1660 euler_dof_handler->get_fe(),
1661 fe_mask,
1662 fe_to_real,
1663 output_data.jacobian_pushed_forward_2nd_derivatives);
1664
1665 // calculate gradients of the hessians of the Jacobians
1666 internal::MappingFEFieldImplementation::
1667 maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1668 QProjector<dim>::DataSetDescriptor::cell(),
1669 data,
1670 euler_dof_handler->get_fe(),
1671 fe_mask,
1672 fe_to_real,
1673 output_data.jacobian_3rd_derivatives);
1674
1675 // calculate gradients of the hessians of the Jacobians pushed forward to real
1676 // cell coordinates
1677 internal::MappingFEFieldImplementation::
1678 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1679 spacedim,
1680 VectorType>(
1681 QProjector<dim>::DataSetDescriptor::cell(),
1682 data,
1683 euler_dof_handler->get_fe(),
1684 fe_mask,
1685 fe_to_real,
1686 output_data.jacobian_pushed_forward_3rd_derivatives);
1687
1688 return CellSimilarity::invalid_next_cell;
1689}
1690
1691
1692
1693template <int dim, int spacedim, typename VectorType>
1694void
1695MappingFEField<dim, spacedim, VectorType, void>::fill_fe_face_values(
1696 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1697 const unsigned int face_no,
1698 const hp::QCollection<dim - 1> & quadrature,
1699 const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1700 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1701 &output_data) const
1702{
1703 AssertDimension(quadrature.size(), 1);
1704
1705 // convert data object to internal data for this class. fails with an
1706 // exception if that is not possible
1707 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1708 ExcInternalError());
1709 const InternalData &data = static_cast<const InternalData &>(internal_data);
1710
1711 update_internal_dofs(cell, data);
1712
1713 internal::MappingFEFieldImplementation::
1714 do_fill_fe_face_values<dim, spacedim, VectorType>(
1715 *this,
1716 cell,
1717 face_no,
1718 numbers::invalid_unsigned_int,
1719 QProjector<dim>::DataSetDescriptor::face(reference_cell,
1720 face_no,
1721 cell->face_orientation(face_no),
1722 cell->face_flip(face_no),
1723 cell->face_rotation(face_no),
1724 quadrature[0].size()),
1725 data,
1726 euler_dof_handler->get_fe(),
1727 fe_mask,
1728 fe_to_real,
1729 output_data);
1730}
1731
1732
1733template <int dim, int spacedim, typename VectorType>
1734void
1735MappingFEField<dim, spacedim, VectorType, void>::fill_fe_subface_values(
1736 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1737 const unsigned int face_no,
1738 const unsigned int subface_no,
1739 const Quadrature<dim - 1> & quadrature,
1740 const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1741 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
1742 &output_data) const
1743{
1744 // convert data object to internal data for this class. fails with an
1745 // exception if that is not possible
1746 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
1747 ExcInternalError());
1748 const InternalData &data = static_cast<const InternalData &>(internal_data);
1749
1750 update_internal_dofs(cell, data);
1751
1752 internal::MappingFEFieldImplementation::do_fill_fe_face_values<dim,
1753 spacedim,
1754 VectorType>(
1755 *this,
1756 cell,
1757 face_no,
1758 numbers::invalid_unsigned_int,
1759 QProjector<dim>::DataSetDescriptor::subface(reference_cell,
1760 face_no,
1761 subface_no,
1762 cell->face_orientation(face_no),
1763 cell->face_flip(face_no),
1764 cell->face_rotation(face_no),
1765 quadrature.size(),
1766 cell->subface_case(face_no)),
1767 data,
1768 euler_dof_handler->get_fe(),
1769 fe_mask,
1770 fe_to_real,
1771 output_data);
1772}
1773
1774
1775namespace internal
1776{
1777 namespace MappingFEFieldImplementation
1778 {
1779 namespace
1780 {
1781 template <int dim, int spacedim, int rank, typename VectorType>
1782 void
1783 transform_fields(
1784 const ArrayView<const Tensor<rank, dim>> & input,
1785 const MappingKind mapping_kind,
1786 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1787 const ArrayView<Tensor<rank, spacedim>> & output)
1788 {
1789 AssertDimension(input.size(), output.size());
1790 Assert(
1791 (dynamic_cast<
1792 const typename ::
1793 MappingFEField<dim, spacedim, VectorType, void>::InternalData *>(
1794 &mapping_data) != nullptr),
1795 ExcInternalError());
1796 const typename ::MappingFEField<dim, spacedim, VectorType, void>::
1797 InternalData &data = static_cast<
1798 const typename ::
1799 MappingFEField<dim, spacedim, VectorType, void>::InternalData &>(
1800 mapping_data);
1801
1802 switch (mapping_kind)
1803 {
1804 case mapping_contravariant:
1805 {
1806 Assert(
1807 data.update_each & update_contravariant_transformation,
1808 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
1809 "update_contravariant_transformation"));
1810
1811 for (unsigned int i = 0; i < output.size(); ++i)
1812 output[i] =
1813 apply_transformation(data.contravariant[i], input[i]);
1814
1815 return;
1816 }
1817
1818 case mapping_piola:
1819 {
1820 Assert(
1821 data.update_each & update_contravariant_transformation,
1822 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
1823 "update_contravariant_transformation"));
1824 Assert(
1825 data.update_each & update_volume_elements,
1826 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
1827 "update_volume_elements"));
1828 Assert(rank == 1, ExcMessage("Only for rank 1"));
1829 for (unsigned int i = 0; i < output.size(); ++i)
1830 {
1831 output[i] =
1832 apply_transformation(data.contravariant[i], input[i]);
1833 output[i] /= data.volume_elements[i];
1834 }
1835 return;
1836 }
1837
1838
1839 // We still allow this operation as in the
1840 // reference cell Derivatives are Tensor
1841 // rather than DerivativeForm
1842 case mapping_covariant:
1843 {
1844 Assert(
1845 data.update_each & update_contravariant_transformation,
1846 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
1847 "update_contravariant_transformation"));
1848
1849 for (unsigned int i = 0; i < output.size(); ++i)
1850 output[i] = apply_transformation(data.covariant[i], input[i]);
1851
1852 return;
1853 }
1854
1855 default:
1856 Assert(false, ExcNotImplemented());
1857 }
1858 }
1859
1860
1861 template <int dim, int spacedim, int rank, typename VectorType>
1862 void
1863 transform_differential_forms(
1864 const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
1865 const MappingKind mapping_kind,
1866 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1867 const ArrayView<Tensor<rank + 1, spacedim>> & output)
1868 {
1869 AssertDimension(input.size(), output.size());
1870 Assert(
1871 (dynamic_cast<
1872 const typename ::
1873 MappingFEField<dim, spacedim, VectorType, void>::InternalData *>(
1874 &mapping_data) != nullptr),
1875 ExcInternalError());
1876 const typename ::MappingFEField<dim, spacedim, VectorType, void>::
1877 InternalData &data = static_cast<
1878 const typename ::
1879 MappingFEField<dim, spacedim, VectorType, void>::InternalData &>(
1880 mapping_data);
1881
1882 switch (mapping_kind)
1883 {
1884 case mapping_covariant:
1885 {
1886 Assert(
1887 data.update_each & update_contravariant_transformation,
1888 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
1889 "update_contravariant_transformation"));
1890
1891 for (unsigned int i = 0; i < output.size(); ++i)
1892 output[i] = apply_transformation(data.covariant[i], input[i]);
1893
1894 return;
1895 }
1896 default:
1897 Assert(false, ExcNotImplemented());
1898 }
1899 }
1900 } // namespace
1901 } // namespace MappingFEFieldImplementation
1902} // namespace internal
1903
1904
1905
1906template <int dim, int spacedim, typename VectorType>
1907void
1908MappingFEField<dim, spacedim, VectorType, void>::transform(
1909 const ArrayView<const Tensor<1, dim>> & input,
1910 const MappingKind mapping_kind,
1911 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1912 const ArrayView<Tensor<1, spacedim>> & output) const
1913{
1914 AssertDimension(input.size(), output.size());
1915
1916 internal::MappingFEFieldImplementation::
1917 transform_fields<dim, spacedim, 1, VectorType>(input,
1918 mapping_kind,
1919 mapping_data,
1920 output);
1921}
1922
1923
1924
1925template <int dim, int spacedim, typename VectorType>
1926void
1927MappingFEField<dim, spacedim, VectorType, void>::transform(
1928 const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1929 const MappingKind mapping_kind,
1930 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1931 const ArrayView<Tensor<2, spacedim>> & output) const
1932{
1933 AssertDimension(input.size(), output.size());
1934
1935 internal::MappingFEFieldImplementation::
1936 transform_differential_forms<dim, spacedim, 1, VectorType>(input,
1937 mapping_kind,
1938 mapping_data,
1939 output);
1940}
1941
1942
1943
1944template <int dim, int spacedim, typename VectorType>
1945void
1946MappingFEField<dim, spacedim, VectorType, void>::transform(
1947 const ArrayView<const Tensor<2, dim>> &input,
1948 const MappingKind,
1949 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1950 const ArrayView<Tensor<2, spacedim>> & output) const
1951{
1952 (void)input;
1953 (void)output;
1954 (void)mapping_data;
1955 AssertDimension(input.size(), output.size());
1956
1957 AssertThrow(false, ExcNotImplemented());
1958}
1959
1960
1961
1962template <int dim, int spacedim, typename VectorType>
1963void
1964MappingFEField<dim, spacedim, VectorType, void>::transform(
1965 const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1966 const MappingKind mapping_kind,
1967 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1968 const ArrayView<Tensor<3, spacedim>> & output) const
1969{
1970 AssertDimension(input.size(), output.size());
1971 Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1972 ExcInternalError());
1973 const InternalData &data = static_cast<const InternalData &>(mapping_data);
1974
1975 switch (mapping_kind)
1976 {
1977 case mapping_covariant_gradient:
1978 {
1979 Assert(data.update_each & update_contravariant_transformation,
1980 typename FEValuesBase<dim>::ExcAccessToUninitializedField(
1981 "update_covariant_transformation"));
1982
1983 for (unsigned int q = 0; q < output.size(); ++q)
1984 for (unsigned int i = 0; i < spacedim; ++i)
1985 for (unsigned int j = 0; j < spacedim; ++j)
1986 for (unsigned int k = 0; k < spacedim; ++k)
1987 {
1988 output[q][i][j][k] = data.covariant[q][j][0] *
1989 data.covariant[q][k][0] *
1990 input[q][i][0][0];
1991 for (unsigned int J = 0; J < dim; ++J)
1992 {
1993 const unsigned int K0 = (0 == J) ? 1 : 0;
1994 for (unsigned int K = K0; K < dim; ++K)
1995 output[q][i][j][k] += data.covariant[q][j][J] *
1996 data.covariant[q][k][K] *
1997 input[q][i][J][K];
1998 }
1999 }
2000 return;
2001 }
2002
2003 default:
2004 Assert(false, ExcNotImplemented());
2005 }
2006}
2007
2008
2009
2010template <int dim, int spacedim, typename VectorType>
2011void
2012MappingFEField<dim, spacedim, VectorType, void>::transform(
2013 const ArrayView<const Tensor<3, dim>> &input,
2014 const MappingKind /*mapping_kind*/,
2015 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
2016 const ArrayView<Tensor<3, spacedim>> & output) const
2017{
2018 (void)input;
2019 (void)output;
2020 (void)mapping_data;
2021 AssertDimension(input.size(), output.size());
2022
2023 AssertThrow(false, ExcNotImplemented());
2024}
2025
2026
2027
2028template <int dim, int spacedim, typename VectorType>
2029Point<spacedim>
2030MappingFEField<dim, spacedim, VectorType, void>::transform_unit_to_real_cell(
2031 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2032 const Point<dim> & p) const
2033{
2034 // Use the get_data function to create an InternalData with data vectors of
2035 // the right size and transformation shape values already computed at point
2036 // p.
2037 const Quadrature<dim> point_quadrature(p);
2038 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2039 get_data(update_quadrature_points | update_jacobians, point_quadrature));
2040 Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2041 ExcInternalError());
2042
2043 update_internal_dofs(cell, dynamic_cast<InternalData &>(*mdata));
2044
2045 return do_transform_unit_to_real_cell(dynamic_cast<InternalData &>(*mdata));
2046}
2047
2048
2049template <int dim, int spacedim, typename VectorType>
2050Point<spacedim>
2051MappingFEField<dim, spacedim, VectorType, void>::do_transform_unit_to_real_cell(
2052 const InternalData &data) const
2053{
2054 Point<spacedim> p_real;
2055
2056 for (unsigned int i = 0; i < data.n_shape_functions; ++i)
2057 {
2058 unsigned int comp_i =
2059 euler_dof_handler->get_fe().system_to_component_index(i).first;
2060 if (fe_mask[comp_i])
2061 p_real[fe_to_real[comp_i]] +=
2062 data.local_dof_values[i] * data.shape(0, i);
2063 }
2064
2065 return p_real;
2066}
2067
2068
2069
2070template <int dim, int spacedim, typename VectorType>
2071Point<dim>
2072MappingFEField<dim, spacedim, VectorType, void>::transform_real_to_unit_cell(
2073 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2074 const Point<spacedim> & p) const
2075{
2076 // first a Newton iteration based on the real mapping. It uses the center
2077 // point of the cell as a starting point
2078 Point<dim> initial_p_unit;
2079 try
2080 {
2081 initial_p_unit = get_default_linear_mapping(cell->get_triangulation())
2082 .transform_real_to_unit_cell(cell, p);
2083 }
2084 catch (const typename Mapping<dim, spacedim>::ExcTransformationFailed &)
2085 {
2086 // mirror the conditions of the code below to determine if we need to
2087 // use an arbitrary starting point or if we just need to rethrow the
2088 // exception
2089 for (unsigned int d = 0; d < dim; ++d)
2090 initial_p_unit[d] = 0.5;
2091 }
2092
2093 // TODO
2094 initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);
2095
2096 // for (unsigned int d=0; d<dim; ++d)
2097 // initial_p_unit[d] = 0.;
2098
2099 const Quadrature<dim> point_quadrature(initial_p_unit);
2100
2101 UpdateFlags update_flags = update_quadrature_points | update_jacobians;
2102 if (spacedim > dim)
2103 update_flags |= update_jacobian_grads;
2104 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> mdata(
2105 get_data(update_flags, point_quadrature));
2106 Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
2107 ExcInternalError());
2108
2109 update_internal_dofs(cell, dynamic_cast<InternalData &>(*mdata));
2110
2111 return do_transform_real_to_unit_cell(cell,
2112 p,
2113 initial_p_unit,
2114 dynamic_cast<InternalData &>(*mdata));
2115}
2116
2117
2118template <int dim, int spacedim, typename VectorType>
2119Point<dim>
2120MappingFEField<dim, spacedim, VectorType, void>::do_transform_real_to_unit_cell(
2121 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2122 const Point<spacedim> & p,
2123 const Point<dim> & initial_p_unit,
2124 InternalData & mdata) const
2125{
2126 const unsigned int n_shapes = mdata.shape_values.size();
2127 (void)n_shapes;
2128 Assert(n_shapes != 0, ExcInternalError());
2129 AssertDimension(mdata.shape_derivatives.size(), n_shapes);
2130
2131
2132 // Newton iteration to solve
2133 // f(x)=p(x)-p=0
2134 // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
2135 // The start value was set to be the
2136 // linear approximation to the cell
2137 // The shape values and derivatives
2138 // of the mapping at this point are
2139 // previously computed.
2140 // f(x)
2141 Point<dim> p_unit = initial_p_unit;
2142 Point<dim> f;
2143 compute_shapes_virtual(std::vector<Point<dim>>(1, p_unit), mdata);
2144 Point<spacedim> p_real(do_transform_unit_to_real_cell(mdata));
2145 Tensor<1, spacedim> p_minus_F = p - p_real;
2146 const double eps = 1.e-12 * cell->diameter();
2147 const unsigned int newton_iteration_limit = 20;
2148 unsigned int newton_iteration = 0;
2149 while (p_minus_F.norm_square() > eps * eps)
2150 {
2151 // f'(x)
2152 Point<spacedim> DF[dim];
2153 Tensor<2, dim> df;
2154 for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
2155 {
2156 const Tensor<1, dim> &grad_k = mdata.derivative(0, k);
2157 const unsigned int comp_k =
2158 euler_dof_handler->get_fe().system_to_component_index(k).first;
2159 if (fe_mask[comp_k])
2160 for (unsigned int j = 0; j < dim; ++j)
2161 DF[j][fe_to_real[comp_k]] +=
2162 mdata.local_dof_values[k] * grad_k[j];
2163 }
2164 for (unsigned int j = 0; j < dim; ++j)
2165 {
2166 f[j] = DF[j] * p_minus_F;
2167 for (unsigned int l = 0; l < dim; ++l)
2168 df[j][l] = -DF[j] * DF[l];
2169 }
2170 // Solve [f'(x)]d=f(x)
2171 const Tensor<1, dim> delta =
2172 invert(df) * static_cast<const Tensor<1, dim> &>(f);
2173 // do a line search
2174 double step_length = 1;
2175 do
2176 {
2177 // update of p_unit. The
2178 // spacedimth component of
2179 // transformed point is simply
2180 // ignored in codimension one
2181 // case. When this component is
2182 // not zero, then we are
2183 // projecting the point to the
2184 // surface or curve identified
2185 // by the cell.
2186 Point<dim> p_unit_trial = p_unit;
2187 for (unsigned int i = 0; i < dim; ++i)
2188 p_unit_trial[i] -= step_length * delta[i];
2189 // shape values and derivatives
2190 // at new p_unit point
2191 compute_shapes_virtual(std::vector<Point<dim>>(1, p_unit_trial),
2192 mdata);
2193 // f(x)
2194 Point<spacedim> p_real_trial = do_transform_unit_to_real_cell(mdata);
2195 const Tensor<1, spacedim> f_trial = p - p_real_trial;
2196 // see if we are making progress with the current step length
2197 // and if not, reduce it by a factor of two and try again
2198 if (f_trial.norm() < p_minus_F.norm())
2199 {
2200 p_real = p_real_trial;
2201 p_unit = p_unit_trial;
2202 p_minus_F = f_trial;
2203 break;
2204 }
2205 else if (step_length > 0.05)
2206 step_length /= 2;
2207 else
2208 goto failure;
2209 }
2210 while (true);
2211 ++newton_iteration;
2212 if (newton_iteration > newton_iteration_limit)
2213 goto failure;
2214 }
2215 return p_unit;
2216 // if we get to the following label, then we have either run out
2217 // of Newton iterations, or the line search has not converged.
2218 // in either case, we need to give up, so throw an exception that
2219 // can then be caught
2220failure:
2221 AssertThrow(false,
2222 (typename Mapping<dim, spacedim>::ExcTransformationFailed()));
2223 // ...the compiler wants us to return something, though we can
2224 // of course never get here...
2225 return {};
2226}
2227
2228
2229template <int dim, int spacedim, typename VectorType>
2230unsigned int
2231MappingFEField<dim, spacedim, VectorType, void>::get_degree() const
2232{
2233 return euler_dof_handler->get_fe().degree;
2234}
2235
2236
2237
2238template <int dim, int spacedim, typename VectorType>
2239ComponentMask
2240MappingFEField<dim, spacedim, VectorType, void>::get_component_mask() const
2241{
2242 return this->fe_mask;
2243}
2244
2245
2246template <int dim, int spacedim, typename VectorType>
2247std::unique_ptr<Mapping<dim, spacedim>>
2248MappingFEField<dim, spacedim, VectorType, void>::clone() const
2249{
2250 return std::make_unique<MappingFEField<dim, spacedim, VectorType, void>>(
2251 *this);
2252}
2253
2254
2255template <int dim, int spacedim, typename VectorType>
2256void
2257MappingFEField<dim, spacedim, VectorType, void>::update_internal_dofs(
2258 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2259 const typename MappingFEField<dim, spacedim, VectorType, void>::InternalData
2260 &data) const
2261{
2262 Assert(euler_dof_handler != nullptr,
2263 ExcMessage("euler_dof_handler is empty"));
2264
2265 typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(*cell,
2266 euler_dof_handler);
2267 Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
2268 if (uses_level_dofs)
2269 {
2270 AssertIndexRange(cell->level(), euler_vector.size());
2271 AssertDimension(euler_vector[cell->level()]->size(),
2272 euler_dof_handler->n_dofs(cell->level()));
2273 }
2274 else
2275 AssertDimension(euler_vector[0]->size(), euler_dof_handler->n_dofs());
2276
2277 if (uses_level_dofs)
2278 dof_cell->get_mg_dof_indices(data.local_dof_indices);
2279 else
2280 dof_cell->get_dof_indices(data.local_dof_indices);
2281
2282 const VectorType &vector =
2283 uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0];
2284
2285 for (unsigned int i = 0; i < data.local_dof_values.size(); ++i)
2286 data.local_dof_values[i] =
2287 internal::ElementAccess<VectorType>::get(vector,
2288 data.local_dof_indices[i]);
2289}
2290
2291// explicit instantiations
2292#define SPLIT_INSTANTIATIONS_COUNT 2
2293#ifndef SPLIT_INSTANTIATIONS_INDEX
2294# define SPLIT_INSTANTIATIONS_INDEX 0
2295#endif
2296#include "mapping_fe_field.inst"
2297
2298
2299DEAL_II_NAMESPACE_CLOSE
array_view.h
make_array_view
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:699
block_vector.h
DerivativeForm::transpose
DerivativeForm< 1, spacedim, dim, Number > transpose() const
DoFHandler::get_triangulation
const Triangulation< dim, spacedim > & get_triangulation() const
DoFHandler::has_level_dofs
bool has_level_dofs() const
FiniteElement::system_to_component_index
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
MGLevelObject::resize
void resize(const unsigned int new_minlevel, const unsigned int new_maxlevel, Args &&...args)
Definition: mg_level_object.h:250
MGLevelObject::max_level
unsigned int max_level() const
Definition: mg_level_object.h:298
MGLevelObject::min_level
unsigned int min_level() const
Definition: mg_level_object.h:290
Mapping
Abstract base class for mapping classes.
Definition: mapping.h:311
Mapping::transform_real_to_unit_cell
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
Point
Definition: point.h:111
QProjector::DataSetDescriptor::cell
static DataSetDescriptor cell()
Definition: qprojector.h:563
QProjector::DataSetDescriptor::subface
static DataSetDescriptor subface(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)
Quadrature::get_weights
const std::vector< double > & get_weights() const
Quadrature::get_points
const std::vector< Point< dim > > & get_points() const
Quadrature::size
unsigned int size() const
Tensor
Definition: tensor.h:503
Tensor::norm
numbers::NumberTraits< Number >::real_type norm() const
Tensor::norm_square
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
Triangulation::n_global_levels
virtual unsigned int n_global_levels() const
hp::Collection::size
unsigned int size() const
Definition: collection.h:264
internal::FEValuesImplementation::MappingRelatedData::inverse_jacobians
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
Definition: fe_update_flags.h:449
internal::FEValuesImplementation::MappingRelatedData::normal_vectors
std::vector< Tensor< 1, spacedim > > normal_vectors
Definition: fe_update_flags.h:491
internal::FEValuesImplementation::MappingRelatedData::jacobian_pushed_forward_3rd_derivatives
std::vector< Tensor< 5, spacedim > > jacobian_pushed_forward_3rd_derivatives
Definition: fe_update_flags.h:479
internal::FEValuesImplementation::MappingRelatedData::jacobian_3rd_derivatives
std::vector< DerivativeForm< 4, dim, spacedim > > jacobian_3rd_derivatives
Definition: fe_update_flags.h:473
internal::FEValuesImplementation::MappingRelatedData::jacobian_2nd_derivatives
std::vector< DerivativeForm< 3, dim, spacedim > > jacobian_2nd_derivatives
Definition: fe_update_flags.h:461
internal::FEValuesImplementation::MappingRelatedData::quadrature_points
std::vector< Point< spacedim > > quadrature_points
Definition: fe_update_flags.h:486
internal::FEValuesImplementation::MappingRelatedData::jacobian_pushed_forward_2nd_derivatives
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
Definition: fe_update_flags.h:467
internal::FEValuesImplementation::MappingRelatedData::jacobian_pushed_forward_grads
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
Definition: fe_update_flags.h:455
internal::FEValuesImplementation::MappingRelatedData::jacobian_grads
std::vector< DerivativeForm< 2, dim, spacedim > > jacobian_grads
Definition: fe_update_flags.h:444
internal::FEValuesImplementation::MappingRelatedData::boundary_forms
std::vector< Tensor< 1, spacedim > > boundary_forms
Definition: fe_update_flags.h:496
internal::FEValuesImplementation::MappingRelatedData::jacobians
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
Definition: fe_update_flags.h:438
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
vertices
Point< 3 > vertices[4]
Definition: data_out_base.cc:192
apply_transformation
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
Definition: derivative_form.h:431
dof_accessor.h
fe_values.h
fe_q.h
fe_system.h
fe_tools.h
UpdateFlags
UpdateFlags
Definition: fe_update_flags.h:67
update_jacobian_pushed_forward_2nd_derivatives
@ update_jacobian_pushed_forward_2nd_derivatives
Definition: fe_update_flags.h:198
update_volume_elements
@ update_volume_elements
Determinant of the Jacobian.
Definition: fe_update_flags.h:184
update_contravariant_transformation
@ update_contravariant_transformation
Contravariant transformation.
Definition: fe_update_flags.h:167
update_jacobian_pushed_forward_grads
@ update_jacobian_pushed_forward_grads
Definition: fe_update_flags.h:189
update_jacobian_3rd_derivatives
@ update_jacobian_3rd_derivatives
Definition: fe_update_flags.h:202
update_values
@ update_values
Shape function values.
Definition: fe_update_flags.h:78
update_jacobian_grads
@ update_jacobian_grads
Gradient of volume element.
Definition: fe_update_flags.h:147
update_normal_vectors
@ update_normal_vectors
Normal vectors.
Definition: fe_update_flags.h:136
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:129
update_covariant_transformation
@ update_covariant_transformation
Covariant transformation.
Definition: fe_update_flags.h:160
update_jacobians
@ update_jacobians
Volume element.
Definition: fe_update_flags.h:142
update_inverse_jacobians
@ update_inverse_jacobians
Volume element.
Definition: fe_update_flags.h:153
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition: fe_update_flags.h:122
update_jacobian_pushed_forward_3rd_derivatives
@ update_jacobian_pushed_forward_3rd_derivatives
Definition: fe_update_flags.h:207
update_boundary_forms
@ update_boundary_forms
Outer normal vector, not normalized.
Definition: fe_update_flags.h:103
update_jacobian_2nd_derivatives
@ update_jacobian_2nd_derivatives
Definition: fe_update_flags.h:193
full_matrix.h
second
Point< 2 > second
Definition: grid_out.cc:4604
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1473
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
DoFHandler::cell_iterator
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:466
MappingKind
MappingKind
Definition: mapping.h:72
mapping_piola
@ mapping_piola
Definition: mapping.h:107
mapping_covariant_gradient
@ mapping_covariant_gradient
Definition: mapping.h:93
mapping_covariant
@ mapping_covariant
Definition: mapping.h:82
mapping_contravariant
@ mapping_contravariant
Definition: mapping.h:87
la_parallel_block_vector.h
la_parallel_vector.h
la_vector.h
get_default_linear_mapping
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Definition: mapping.cc:260
mapping_fe_field.h
mapping_q1.h
memory_consumption.h
Utilities::to_string
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:482
internal::MappingQImplementation::transform_differential_forms
void transform_differential_forms(const ArrayView< const DerivativeForm< rank, dim, spacedim > > &input, const MappingKind mapping_kind, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_data, const ArrayView< Tensor< rank+1, spacedim > > &output)
Definition: mapping_q_internal.h:2451
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:201
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5510
petsc_block_vector.h
petsc_vector.h
quadrature_lib.h
GeometryInfo::subface_ratio
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
GeometryInfo::project_to_unit_cell
static Point< dim, Number > project_to_unit_cell(const Point< dim, Number > &p)
internal::ElementAccess::get
static VectorType::value_type get(const VectorType &V, const types::global_dof_index i)
Definition: vector_element_access.h:73
determinant
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
Definition: symmetric_tensor.h:2767
invert
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
Definition: symmetric_tensor.h:3422
cross_product_2d
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2635
cross_product_3d
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition: tensor.h:2660
tensor_product_polynomials.h
tria_iterator.h
trilinos_epetra_vector.h
trilinos_parallel_block_vector.h
trilinos_tpetra_vector.h
trilinos_vector.h
vector_tools.h
Generated by doxygen 1.9.5