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