Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
mapping_cartesian.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2019 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/qprojector.h>
19 #include <deal.II/base/quadrature.h>
20 #include <deal.II/base/signaling_nan.h>
21 #include <deal.II/base/std_cxx14/memory.h>
22 #include <deal.II/base/tensor.h>
23 
24 #include <deal.II/dofs/dof_accessor.h>
25 
26 #include <deal.II/fe/fe_values.h>
27 #include <deal.II/fe/mapping_cartesian.h>
28 
29 #include <deal.II/grid/tria.h>
30 #include <deal.II/grid/tria_iterator.h>
31 
32 #include <deal.II/lac/full_matrix.h>
33 
34 #include <algorithm>
35 #include <cmath>
36 
37 
38 DEAL_II_NAMESPACE_OPEN
39 
40 
41 template <int dim, int spacedim>
43 
44 
45 
46 template <int dim, int spacedim>
48  const Quadrature<dim> &q)
49  : cell_extents(numbers::signaling_nan<Tensor<1, dim>>())
50  , volume_element(numbers::signaling_nan<double>())
51  , quadrature_points(q.get_points())
52 {}
53 
54 
55 
56 template <int dim, int spacedim>
57 std::size_t
59 {
63  MemoryConsumption::memory_consumption(quadrature_points));
64 }
65 
66 
67 
68 template <int dim, int spacedim>
69 bool
71 {
72  return true;
73 }
74 
75 
76 
77 template <int dim, int spacedim>
80  const UpdateFlags in) const
81 {
82  // this mapping is pretty simple in that it can basically compute
83  // every piece of information wanted by FEValues without requiring
84  // computing any other quantities. boundary forms are one exception
85  // since they can be computed from the normal vectors without much
86  // further ado
87  UpdateFlags out = in;
88  if (out & update_boundary_forms)
89  out |= update_normal_vectors;
90 
91  return out;
92 }
93 
94 
95 
96 template <int dim, int spacedim>
97 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
99  const Quadrature<dim> &q) const
100 {
101  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
102  std_cxx14::make_unique<InternalData>(q);
103  auto &data = dynamic_cast<InternalData &>(*data_ptr);
104 
105  // store the flags in the internal data object so we can access them
106  // in fill_fe_*_values(). use the transitive hull of the required
107  // flags
108  data.update_each = requires_update_flags(update_flags);
109 
110  return data_ptr;
111 }
112 
113 
114 
115 template <int dim, int spacedim>
116 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
118  const UpdateFlags update_flags,
119  const Quadrature<dim - 1> &quadrature) const
120 {
121  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
122  std_cxx14::make_unique<InternalData>(
124  auto &data = dynamic_cast<InternalData &>(*data_ptr);
125 
126  // verify that we have computed the transitive hull of the required
127  // flags and that FEValues has faithfully passed them on to us
128  Assert(update_flags == requires_update_flags(update_flags),
129  ExcInternalError());
130 
131  // store the flags in the internal data object so we can access them
132  // in fill_fe_*_values()
133  data.update_each = update_flags;
134 
135  return data_ptr;
136 }
137 
138 
139 
140 template <int dim, int spacedim>
141 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
143  const UpdateFlags update_flags,
144  const Quadrature<dim - 1> &quadrature) const
145 {
146  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
147  std_cxx14::make_unique<InternalData>(
149  auto &data = dynamic_cast<InternalData &>(*data_ptr);
150 
151  // verify that we have computed the transitive hull of the required
152  // flags and that FEValues has faithfully passed them on to us
153  Assert(update_flags == requires_update_flags(update_flags),
154  ExcInternalError());
155 
156  // store the flags in the internal data object so we can access them
157  // in fill_fe_*_values()
158  data.update_each = update_flags;
159 
160  return data_ptr;
161 }
162 
163 
164 
165 template <int dim, int spacedim>
166 void
168  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
169  const unsigned int face_no,
170  const unsigned int sub_no,
171  const CellSimilarity::Similarity cell_similarity,
172  const InternalData & data,
173  std::vector<Point<dim>> & quadrature_points,
174  std::vector<Tensor<1, dim>> &normal_vectors) const
175 {
176  const UpdateFlags update_flags = data.update_each;
177 
178  // some more sanity checks
179  if (face_no != invalid_face_number)
180  {
181  // Add 1 on both sides of
182  // assertion to avoid compiler
183  // warning about testing
184  // unsigned int < 0 in 1d.
185  Assert(face_no + 1 < GeometryInfo<dim>::faces_per_cell + 1,
187 
188  // We would like to check for
189  // sub_no < cell->face(face_no)->n_children(),
190  // but unfortunately the current
191  // function is also called for
192  // faces without children (see
193  // tests/fe/mapping.cc). Therefore,
194  // we must use following workaround
195  // of two separate assertions
196  Assert((sub_no == invalid_face_number) ||
197  cell->face(face_no)->has_children() ||
198  (sub_no + 1 < GeometryInfo<dim>::max_children_per_face + 1),
199  ExcIndexRange(sub_no,
200  0,
202  Assert((sub_no == invalid_face_number) ||
203  !cell->face(face_no)->has_children() ||
204  (sub_no < cell->face(face_no)->n_children()),
205  ExcIndexRange(sub_no, 0, cell->face(face_no)->n_children()));
206  }
207  else
208  // invalid face number, so
209  // subface should be invalid as
210  // well
212 
213  // let @p{start} be the origin of a
214  // local coordinate system. it is
215  // chosen as the (lower) left
216  // vertex
217  const Point<dim> start = cell->vertex(0);
218 
219  // Compute start point and sizes
220  // along axes. Strange vertex
221  // numbering makes this complicated
222  // again.
223  if (cell_similarity != CellSimilarity::translation)
224  {
225  switch (dim)
226  {
227  case 1:
228  data.cell_extents[0] = cell->vertex(1)(0) - start(0);
229  break;
230  case 2:
231  data.cell_extents[0] = cell->vertex(1)(0) - start(0);
232  data.cell_extents[1] = cell->vertex(2)(1) - start(1);
233  break;
234  case 3:
235  data.cell_extents[0] = cell->vertex(1)(0) - start(0);
236  data.cell_extents[1] = cell->vertex(2)(1) - start(1);
237  data.cell_extents[2] = cell->vertex(4)(2) - start(2);
238  break;
239  default:
240  Assert(false, ExcNotImplemented());
241  }
242  }
243 
244 
245  // transform quadrature point. this
246  // is obtained simply by scaling
247  // unit coordinates with lengths in
248  // each direction
249  if (update_flags & update_quadrature_points)
250  {
251  const typename QProjector<dim>::DataSetDescriptor offset =
252  (face_no == invalid_face_number ?
254  (sub_no == invalid_face_number ?
255  // called from FEFaceValues
257  face_no,
258  cell->face_orientation(face_no),
259  cell->face_flip(face_no),
260  cell->face_rotation(face_no),
261  quadrature_points.size()) :
262  // called from FESubfaceValues
264  face_no,
265  sub_no,
266  cell->face_orientation(face_no),
267  cell->face_flip(face_no),
268  cell->face_rotation(face_no),
269  quadrature_points.size(),
270  cell->subface_case(face_no))));
271 
272  for (unsigned int i = 0; i < quadrature_points.size(); ++i)
273  {
274  quadrature_points[i] = start;
275  for (unsigned int d = 0; d < dim; ++d)
276  quadrature_points[i](d) +=
277  data.cell_extents[d] * data.quadrature_points[i + offset](d);
278  }
279  }
280 
281 
282  // compute normal vectors. since
283  // cells are aligned to coordinate
284  // axes, they are simply vectors
285  // with exactly one entry equal to
286  // 1 or -1. Furthermore, all
287  // normals on a face have the same
288  // value
289  if (update_flags & update_normal_vectors)
290  {
292 
293  switch (dim)
294  {
295  case 1:
296  {
297  static const Point<dim> normals[GeometryInfo<1>::faces_per_cell] =
298  {Point<dim>(-1.), Point<dim>(1.)};
299  std::fill(normal_vectors.begin(),
300  normal_vectors.end(),
301  normals[face_no]);
302  break;
303  }
304 
305  case 2:
306  {
307  static const Point<dim> normals[GeometryInfo<2>::faces_per_cell] =
308  {Point<dim>(-1, 0),
309  Point<dim>(1, 0),
310  Point<dim>(0, -1),
311  Point<dim>(0, 1)};
312  std::fill(normal_vectors.begin(),
313  normal_vectors.end(),
314  normals[face_no]);
315  break;
316  }
317 
318  case 3:
319  {
320  static const Point<dim> normals[GeometryInfo<3>::faces_per_cell] =
321  {Point<dim>(-1, 0, 0),
322  Point<dim>(1, 0, 0),
323  Point<dim>(0, -1, 0),
324  Point<dim>(0, 1, 0),
325  Point<dim>(0, 0, -1),
326  Point<dim>(0, 0, 1)};
327  std::fill(normal_vectors.begin(),
328  normal_vectors.end(),
329  normals[face_no]);
330  break;
331  }
332 
333  default:
334  Assert(false, ExcNotImplemented());
335  }
336  }
337 }
338 
339 
340 
341 template <int dim, int spacedim>
344  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
345  const CellSimilarity::Similarity cell_similarity,
346  const Quadrature<dim> & quadrature,
347  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
349  &output_data) const
350 {
351  // convert data object to internal data for this class. fails with
352  // an exception if that is not possible
353  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
354  ExcInternalError());
355  const InternalData &data = static_cast<const InternalData &>(internal_data);
356 
357  std::vector<Tensor<1, dim>> dummy;
358 
359  compute_fill(cell,
362  cell_similarity,
363  data,
364  output_data.quadrature_points,
365  dummy);
366 
367  // compute Jacobian determinant. all values are equal and are the
368  // product of the local lengths in each coordinate direction
369  if (data.update_each & (update_JxW_values | update_volume_elements))
370  if (cell_similarity != CellSimilarity::translation)
371  {
372  double J = data.cell_extents[0];
373  for (unsigned int d = 1; d < dim; ++d)
374  J *= data.cell_extents[d];
375  data.volume_element = J;
376  if (data.update_each & update_JxW_values)
377  for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
378  output_data.JxW_values[i] = J * quadrature.weight(i);
379  }
380  // "compute" Jacobian at the quadrature points, which are all the
381  // same
382  if (data.update_each & update_jacobians)
383  if (cell_similarity != CellSimilarity::translation)
384  for (unsigned int i = 0; i < output_data.jacobians.size(); ++i)
385  {
386  output_data.jacobians[i] = DerivativeForm<1, dim, spacedim>();
387  for (unsigned int j = 0; j < dim; ++j)
388  output_data.jacobians[i][j][j] = data.cell_extents[j];
389  }
390  // "compute" the derivative of the Jacobian at the quadrature
391  // points, which are all zero of course
392  if (data.update_each & update_jacobian_grads)
393  if (cell_similarity != CellSimilarity::translation)
394  for (unsigned int i = 0; i < output_data.jacobian_grads.size(); ++i)
396 
397  if (data.update_each & update_jacobian_pushed_forward_grads)
398  if (cell_similarity != CellSimilarity::translation)
399  for (unsigned int i = 0;
400  i < output_data.jacobian_pushed_forward_grads.size();
401  ++i)
403 
404  // "compute" the hessian of the Jacobian at the quadrature points,
405  // which are all also zero of course
406  if (data.update_each & update_jacobian_2nd_derivatives)
407  if (cell_similarity != CellSimilarity::translation)
408  for (unsigned int i = 0; i < output_data.jacobian_2nd_derivatives.size();
409  ++i)
410  output_data.jacobian_2nd_derivatives[i] =
412 
413  if (data.update_each & update_jacobian_pushed_forward_2nd_derivatives)
414  if (cell_similarity != CellSimilarity::translation)
415  for (unsigned int i = 0;
416  i < output_data.jacobian_pushed_forward_2nd_derivatives.size();
417  ++i)
420 
421  if (data.update_each & update_jacobian_3rd_derivatives)
422  if (cell_similarity != CellSimilarity::translation)
423  for (unsigned int i = 0; i < output_data.jacobian_3rd_derivatives.size();
424  ++i)
425  output_data.jacobian_3rd_derivatives[i] =
427 
428  if (data.update_each & update_jacobian_pushed_forward_3rd_derivatives)
429  if (cell_similarity != CellSimilarity::translation)
430  for (unsigned int i = 0;
431  i < output_data.jacobian_pushed_forward_3rd_derivatives.size();
432  ++i)
435 
436  // "compute" inverse Jacobian at the quadrature points, which are
437  // all the same
438  if (data.update_each & update_inverse_jacobians)
439  if (cell_similarity != CellSimilarity::translation)
440  for (unsigned int i = 0; i < output_data.inverse_jacobians.size(); ++i)
441  {
442  output_data.inverse_jacobians[i] = Tensor<2, dim>();
443  for (unsigned int j = 0; j < dim; ++j)
444  output_data.inverse_jacobians[i][j][j] = 1. / data.cell_extents[j];
445  }
446 
447  return cell_similarity;
448 }
449 
450 
451 
452 template <int dim, int spacedim>
453 void
455  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
456  const unsigned int face_no,
457  const Quadrature<dim - 1> & quadrature,
458  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
460  &output_data) const
461 {
462  // convert data object to internal
463  // data for this class. fails with
464  // an exception if that is not
465  // possible
466  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
467  ExcInternalError());
468  const InternalData &data = static_cast<const InternalData &>(internal_data);
469 
470  compute_fill(cell,
471  face_no,
474  data,
475  output_data.quadrature_points,
476  output_data.normal_vectors);
477 
478  // first compute Jacobian determinant, which is simply the product
479  // of the local lengths since the jacobian is diagonal
480  double J = 1.;
481  for (unsigned int d = 0; d < dim; ++d)
483  J *= data.cell_extents[d];
484 
485  if (data.update_each & update_JxW_values)
486  for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
487  output_data.JxW_values[i] = J * quadrature.weight(i);
488 
489  if (data.update_each & update_boundary_forms)
490  for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
491  output_data.boundary_forms[i] = J * output_data.normal_vectors[i];
492 
493  if (data.update_each & update_volume_elements)
494  {
495  J = data.cell_extents[0];
496  for (unsigned int d = 1; d < dim; ++d)
497  J *= data.cell_extents[d];
498  data.volume_element = J;
499  }
500 
501  if (data.update_each & update_jacobians)
502  for (unsigned int i = 0; i < output_data.jacobians.size(); ++i)
503  {
504  output_data.jacobians[i] = DerivativeForm<1, dim, spacedim>();
505  for (unsigned int d = 0; d < dim; ++d)
506  output_data.jacobians[i][d][d] = data.cell_extents[d];
507  }
508 
509  if (data.update_each & update_jacobian_grads)
510  for (unsigned int i = 0; i < output_data.jacobian_grads.size(); ++i)
512 
513  if (data.update_each & update_jacobian_pushed_forward_grads)
514  for (unsigned int i = 0;
515  i < output_data.jacobian_pushed_forward_grads.size();
516  ++i)
518 
519  if (data.update_each & update_jacobian_2nd_derivatives)
520  for (unsigned int i = 0; i < output_data.jacobian_2nd_derivatives.size();
521  ++i)
522  output_data.jacobian_2nd_derivatives[i] =
524 
525  if (data.update_each & update_jacobian_pushed_forward_2nd_derivatives)
526  for (unsigned int i = 0;
527  i < output_data.jacobian_pushed_forward_2nd_derivatives.size();
528  ++i)
531 
532  if (data.update_each & update_jacobian_3rd_derivatives)
533  for (unsigned int i = 0; i < output_data.jacobian_3rd_derivatives.size();
534  ++i)
535  output_data.jacobian_3rd_derivatives[i] =
537 
538  if (data.update_each & update_jacobian_pushed_forward_3rd_derivatives)
539  for (unsigned int i = 0;
540  i < output_data.jacobian_pushed_forward_3rd_derivatives.size();
541  ++i)
544 
545  if (data.update_each & update_inverse_jacobians)
546  for (unsigned int i = 0; i < output_data.inverse_jacobians.size(); ++i)
547  {
549  for (unsigned int d = 0; d < dim; ++d)
550  output_data.inverse_jacobians[i][d][d] = 1. / data.cell_extents[d];
551  }
552 }
553 
554 
555 
556 template <int dim, int spacedim>
557 void
559  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
560  const unsigned int face_no,
561  const unsigned int subface_no,
562  const Quadrature<dim - 1> & quadrature,
563  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
565  &output_data) const
566 {
567  // convert data object to internal data for this class. fails with
568  // an exception if that is not possible
569  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
570  ExcInternalError());
571  const InternalData &data = static_cast<const InternalData &>(internal_data);
572 
573  compute_fill(cell,
574  face_no,
575  subface_no,
577  data,
578  output_data.quadrature_points,
579  output_data.normal_vectors);
580 
581  // first compute Jacobian determinant, which is simply the product
582  // of the local lengths since the jacobian is diagonal
583  double J = 1.;
584  for (unsigned int d = 0; d < dim; ++d)
586  J *= data.cell_extents[d];
587 
588  if (data.update_each & update_JxW_values)
589  {
590  // Here, cell->face(face_no)->n_children() would be the right
591  // choice, but unfortunately the current function is also called
592  // for faces without children (see tests/fe/mapping.cc). Add
593  // following switch to avoid diffs in tests/fe/mapping.OK
594  const unsigned int n_subfaces =
595  cell->face(face_no)->has_children() ?
596  cell->face(face_no)->n_children() :
598  for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
599  output_data.JxW_values[i] = J * quadrature.weight(i) / n_subfaces;
600  }
601 
602  if (data.update_each & update_boundary_forms)
603  for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
604  output_data.boundary_forms[i] = J * output_data.normal_vectors[i];
605 
606  if (data.update_each & update_volume_elements)
607  {
608  J = data.cell_extents[0];
609  for (unsigned int d = 1; d < dim; ++d)
610  J *= data.cell_extents[d];
611  data.volume_element = J;
612  }
613 
614  if (data.update_each & update_jacobians)
615  for (unsigned int i = 0; i < output_data.jacobians.size(); ++i)
616  {
617  output_data.jacobians[i] = DerivativeForm<1, dim, spacedim>();
618  for (unsigned int d = 0; d < dim; ++d)
619  output_data.jacobians[i][d][d] = data.cell_extents[d];
620  }
621 
622  if (data.update_each & update_jacobian_grads)
623  for (unsigned int i = 0; i < output_data.jacobian_grads.size(); ++i)
625 
626  if (data.update_each & update_jacobian_pushed_forward_grads)
627  for (unsigned int i = 0;
628  i < output_data.jacobian_pushed_forward_grads.size();
629  ++i)
631 
632  if (data.update_each & update_jacobian_2nd_derivatives)
633  for (unsigned int i = 0; i < output_data.jacobian_2nd_derivatives.size();
634  ++i)
635  output_data.jacobian_2nd_derivatives[i] =
637 
638  if (data.update_each & update_jacobian_pushed_forward_2nd_derivatives)
639  for (unsigned int i = 0;
640  i < output_data.jacobian_pushed_forward_2nd_derivatives.size();
641  ++i)
644 
645  if (data.update_each & update_jacobian_3rd_derivatives)
646  for (unsigned int i = 0; i < output_data.jacobian_3rd_derivatives.size();
647  ++i)
648  output_data.jacobian_3rd_derivatives[i] =
650 
651  if (data.update_each & update_jacobian_pushed_forward_3rd_derivatives)
652  for (unsigned int i = 0;
653  i < output_data.jacobian_pushed_forward_3rd_derivatives.size();
654  ++i)
657 
658  if (data.update_each & update_inverse_jacobians)
659  for (unsigned int i = 0; i < output_data.inverse_jacobians.size(); ++i)
660  {
662  for (unsigned int d = 0; d < dim; ++d)
663  output_data.inverse_jacobians[i][d][d] = 1. / data.cell_extents[d];
664  }
665 }
666 
667 
668 
669 template <int dim, int spacedim>
670 void
672  const ArrayView<const Tensor<1, dim>> & input,
673  const MappingType mapping_type,
674  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
675  const ArrayView<Tensor<1, spacedim>> & output) const
676 {
677  AssertDimension(input.size(), output.size());
678  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
679  ExcInternalError());
680  const InternalData &data = static_cast<const InternalData &>(mapping_data);
681 
682  switch (mapping_type)
683  {
684  case mapping_covariant:
685  {
688  "update_covariant_transformation"));
689 
690  for (unsigned int i = 0; i < output.size(); ++i)
691  for (unsigned int d = 0; d < dim; ++d)
692  output[i][d] = input[i][d] / data.cell_extents[d];
693  return;
694  }
695 
697  {
700  "update_contravariant_transformation"));
701 
702  for (unsigned int i = 0; i < output.size(); ++i)
703  for (unsigned int d = 0; d < dim; ++d)
704  output[i][d] = input[i][d] * data.cell_extents[d];
705  return;
706  }
707  case mapping_piola:
708  {
711  "update_contravariant_transformation"));
714  "update_volume_elements"));
715 
716  for (unsigned int i = 0; i < output.size(); ++i)
717  for (unsigned int d = 0; d < dim; ++d)
718  output[i][d] =
719  input[i][d] * data.cell_extents[d] / data.volume_element;
720  return;
721  }
722  default:
723  Assert(false, ExcNotImplemented());
724  }
725 }
726 
727 
728 
729 template <int dim, int spacedim>
730 void
732  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
733  const MappingType mapping_type,
734  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
735  const ArrayView<Tensor<2, spacedim>> & output) const
736 {
737  AssertDimension(input.size(), output.size());
738  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
739  ExcInternalError());
740  const InternalData &data = static_cast<const InternalData &>(mapping_data);
741 
742  switch (mapping_type)
743  {
744  case mapping_covariant:
745  {
748  "update_covariant_transformation"));
749 
750  for (unsigned int i = 0; i < output.size(); ++i)
751  for (unsigned int d1 = 0; d1 < dim; ++d1)
752  for (unsigned int d2 = 0; d2 < dim; ++d2)
753  output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2];
754  return;
755  }
756 
758  {
761  "update_contravariant_transformation"));
762 
763  for (unsigned int i = 0; i < output.size(); ++i)
764  for (unsigned int d1 = 0; d1 < dim; ++d1)
765  for (unsigned int d2 = 0; d2 < dim; ++d2)
766  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2];
767  return;
768  }
769 
771  {
774  "update_covariant_transformation"));
775 
776  for (unsigned int i = 0; i < output.size(); ++i)
777  for (unsigned int d1 = 0; d1 < dim; ++d1)
778  for (unsigned int d2 = 0; d2 < dim; ++d2)
779  output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2] /
780  data.cell_extents[d1];
781  return;
782  }
783 
785  {
788  "update_contravariant_transformation"));
789 
790  for (unsigned int i = 0; i < output.size(); ++i)
791  for (unsigned int d1 = 0; d1 < dim; ++d1)
792  for (unsigned int d2 = 0; d2 < dim; ++d2)
793  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
794  data.cell_extents[d1];
795  return;
796  }
797 
798  case mapping_piola:
799  {
802  "update_contravariant_transformation"));
805  "update_volume_elements"));
806 
807  for (unsigned int i = 0; i < output.size(); ++i)
808  for (unsigned int d1 = 0; d1 < dim; ++d1)
809  for (unsigned int d2 = 0; d2 < dim; ++d2)
810  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
811  data.volume_element;
812  return;
813  }
814 
816  {
819  "update_contravariant_transformation"));
822  "update_volume_elements"));
823 
824  for (unsigned int i = 0; i < output.size(); ++i)
825  for (unsigned int d1 = 0; d1 < dim; ++d1)
826  for (unsigned int d2 = 0; d2 < dim; ++d2)
827  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
828  data.cell_extents[d1] / data.volume_element;
829  return;
830  }
831 
832  default:
833  Assert(false, ExcNotImplemented());
834  }
835 }
836 
837 
838 
839 template <int dim, int spacedim>
840 void
842  const ArrayView<const Tensor<2, dim>> & input,
843  const MappingType mapping_type,
844  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
845  const ArrayView<Tensor<2, spacedim>> & output) const
846 {
847  AssertDimension(input.size(), output.size());
848  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
849  ExcInternalError());
850  const InternalData &data = static_cast<const InternalData &>(mapping_data);
851 
852  switch (mapping_type)
853  {
854  case mapping_covariant:
855  {
858  "update_covariant_transformation"));
859 
860  for (unsigned int i = 0; i < output.size(); ++i)
861  for (unsigned int d1 = 0; d1 < dim; ++d1)
862  for (unsigned int d2 = 0; d2 < dim; ++d2)
863  output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2];
864  return;
865  }
866 
868  {
871  "update_contravariant_transformation"));
872 
873  for (unsigned int i = 0; i < output.size(); ++i)
874  for (unsigned int d1 = 0; d1 < dim; ++d1)
875  for (unsigned int d2 = 0; d2 < dim; ++d2)
876  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2];
877  return;
878  }
879 
881  {
884  "update_covariant_transformation"));
885 
886  for (unsigned int i = 0; i < output.size(); ++i)
887  for (unsigned int d1 = 0; d1 < dim; ++d1)
888  for (unsigned int d2 = 0; d2 < dim; ++d2)
889  output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2] /
890  data.cell_extents[d1];
891  return;
892  }
893 
895  {
898  "update_contravariant_transformation"));
899 
900  for (unsigned int i = 0; i < output.size(); ++i)
901  for (unsigned int d1 = 0; d1 < dim; ++d1)
902  for (unsigned int d2 = 0; d2 < dim; ++d2)
903  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
904  data.cell_extents[d1];
905  return;
906  }
907 
908  case mapping_piola:
909  {
912  "update_contravariant_transformation"));
915  "update_volume_elements"));
916 
917  for (unsigned int i = 0; i < output.size(); ++i)
918  for (unsigned int d1 = 0; d1 < dim; ++d1)
919  for (unsigned int d2 = 0; d2 < dim; ++d2)
920  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
921  data.volume_element;
922  return;
923  }
924 
926  {
929  "update_contravariant_transformation"));
932  "update_volume_elements"));
933 
934  for (unsigned int i = 0; i < output.size(); ++i)
935  for (unsigned int d1 = 0; d1 < dim; ++d1)
936  for (unsigned int d2 = 0; d2 < dim; ++d2)
937  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
938  data.cell_extents[d1] / data.volume_element;
939  return;
940  }
941 
942  default:
943  Assert(false, ExcNotImplemented());
944  }
945 }
946 
947 
948 template <int dim, int spacedim>
949 void
951  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
952  const MappingType mapping_type,
953  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
954  const ArrayView<Tensor<3, spacedim>> & output) const
955 {
956  AssertDimension(input.size(), output.size());
957  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
958  ExcInternalError());
959  const InternalData &data = static_cast<const InternalData &>(mapping_data);
960 
961  switch (mapping_type)
962  {
964  {
967  "update_covariant_transformation"));
968 
969  for (unsigned int q = 0; q < output.size(); ++q)
970  for (unsigned int i = 0; i < spacedim; ++i)
971  for (unsigned int j = 0; j < spacedim; ++j)
972  for (unsigned int k = 0; k < spacedim; ++k)
973  {
974  output[q][i][j][k] = input[q][i][j][k] /
975  data.cell_extents[j] /
976  data.cell_extents[k];
977  }
978  return;
979  }
980  default:
981  Assert(false, ExcNotImplemented());
982  }
983 }
984 
985 template <int dim, int spacedim>
986 void
988  const ArrayView<const Tensor<3, dim>> & input,
989  const MappingType mapping_type,
990  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
991  const ArrayView<Tensor<3, spacedim>> & output) const
992 {
993  AssertDimension(input.size(), output.size());
994  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
995  ExcInternalError());
996  const InternalData &data = static_cast<const InternalData &>(mapping_data);
997 
998  switch (mapping_type)
999  {
1001  {
1004  "update_covariant_transformation"));
1007  "update_contravariant_transformation"));
1008 
1009  for (unsigned int q = 0; q < output.size(); ++q)
1010  for (unsigned int i = 0; i < spacedim; ++i)
1011  for (unsigned int j = 0; j < spacedim; ++j)
1012  for (unsigned int k = 0; k < spacedim; ++k)
1013  {
1014  output[q][i][j][k] =
1015  input[q][i][j][k] * data.cell_extents[i] /
1016  data.cell_extents[j] / data.cell_extents[k];
1017  }
1018  return;
1019  }
1020 
1022  {
1025  "update_covariant_transformation"));
1026 
1027  for (unsigned int q = 0; q < output.size(); ++q)
1028  for (unsigned int i = 0; i < spacedim; ++i)
1029  for (unsigned int j = 0; j < spacedim; ++j)
1030  for (unsigned int k = 0; k < spacedim; ++k)
1031  {
1032  output[q][i][j][k] =
1033  input[q][i][j][k] / data.cell_extents[i] /
1034  data.cell_extents[j] / data.cell_extents[k];
1035  }
1036 
1037  return;
1038  }
1039 
1040  case mapping_piola_hessian:
1041  {
1044  "update_covariant_transformation"));
1047  "update_contravariant_transformation"));
1050  "update_volume_elements"));
1051 
1052  for (unsigned int q = 0; q < output.size(); ++q)
1053  for (unsigned int i = 0; i < spacedim; ++i)
1054  for (unsigned int j = 0; j < spacedim; ++j)
1055  for (unsigned int k = 0; k < spacedim; ++k)
1056  {
1057  output[q][i][j][k] =
1058  input[q][i][j][k] * data.cell_extents[i] /
1059  data.volume_element / data.cell_extents[j] /
1060  data.cell_extents[k];
1061  }
1062 
1063  return;
1064  }
1065 
1066  default:
1067  Assert(false, ExcNotImplemented());
1068  }
1069 }
1070 
1071 
1072 template <int dim, int spacedim>
1075  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1076  const Point<dim> & p) const
1077 {
1078  Tensor<1, dim> length;
1079  const Point<dim> start = cell->vertex(0);
1080  switch (dim)
1081  {
1082  case 1:
1083  length[0] = cell->vertex(1)(0) - start(0);
1084  break;
1085  case 2:
1086  length[0] = cell->vertex(1)(0) - start(0);
1087  length[1] = cell->vertex(2)(1) - start(1);
1088  break;
1089  case 3:
1090  length[0] = cell->vertex(1)(0) - start(0);
1091  length[1] = cell->vertex(2)(1) - start(1);
1092  length[2] = cell->vertex(4)(2) - start(2);
1093  break;
1094  default:
1095  Assert(false, ExcNotImplemented());
1096  }
1097 
1098  Point<dim> p_real = cell->vertex(0);
1099  for (unsigned int d = 0; d < dim; ++d)
1100  p_real(d) += length[d] * p(d);
1101 
1102  return p_real;
1103 }
1104 
1105 
1106 
1107 template <int dim, int spacedim>
1108 Point<dim>
1110  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1111  const Point<spacedim> & p) const
1112 {
1113  if (dim != spacedim)
1114  Assert(false, ExcNotImplemented());
1115  const Point<dim> &start = cell->vertex(0);
1116  Point<dim> real = p;
1117  real -= start;
1118 
1119  switch (dim)
1120  {
1121  case 1:
1122  real(0) /= cell->vertex(1)(0) - start(0);
1123  break;
1124  case 2:
1125  real(0) /= cell->vertex(1)(0) - start(0);
1126  real(1) /= cell->vertex(2)(1) - start(1);
1127  break;
1128  case 3:
1129  real(0) /= cell->vertex(1)(0) - start(0);
1130  real(1) /= cell->vertex(2)(1) - start(1);
1131  real(2) /= cell->vertex(4)(2) - start(2);
1132  break;
1133  default:
1134  Assert(false, ExcNotImplemented());
1135  }
1136  return real;
1137 }
1138 
1139 
1140 template <int dim, int spacedim>
1141 std::unique_ptr<Mapping<dim, spacedim>>
1143 {
1144  return std_cxx14::make_unique<MappingCartesian<dim, spacedim>>(*this);
1145 }
1146 
1147 
1148 //---------------------------------------------------------------------------
1149 // explicit instantiations
1150 #include "mapping_cartesian.inst"
1151 
1152 
1153 DEAL_II_NAMESPACE_CLOSE
Transformed quadrature weights.
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const override
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
Contravariant transformation.
MappingType
Definition: mapping.h:61
std::vector< Tensor< 1, spacedim > > boundary_forms
Volume element.
std::vector< DerivativeForm< 3, dim, spacedim > > jacobian_2nd_derivatives
Outer normal vector, not normalized.
std::vector< Point< dim > > quadrature_points
Determinant of the Jacobian.
std::vector< DerivativeForm< 4, dim, spacedim > > jacobian_3rd_derivatives
Transformed quadrature points.
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
static DataSetDescriptor cell()
Definition: qprojector.h:344
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
static Quadrature< dim > project_to_all_subfaces(const SubQuadrature &quadrature)
virtual std::size_t memory_consumption() const override
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
#define Assert(cond, exc)
Definition: exceptions.h:1407
UpdateFlags
virtual void transform(const ArrayView< const Tensor< 1, dim >> &input, const MappingType type, const typename Mapping< dim, spacedim >::InternalDataBase &internal, const ArrayView< Tensor< 1, spacedim >> &output) const override
Abstract base class for mapping classes.
Definition: dof_tools.h:57
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
void compute_fill(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const CellSimilarity::Similarity cell_similarity, const InternalData &data, std::vector< Point< dim >> &quadrature_points, std::vector< Tensor< 1, dim >> &normal_vectors) const
Gradient of volume element.
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
static const unsigned int invalid_face_number
std::vector< Point< spacedim > > quadrature_points
static Quadrature< dim > project_to_all_faces(const SubQuadrature &quadrature)
InternalData(const Quadrature< dim > &quadrature)
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
Definition: mpi.h:90
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
virtual bool preserves_vertex_locations() const override
Normal vectors.
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const override
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
static ::ExceptionBase & ExcNotImplemented()
std::vector< Tensor< 5, spacedim > > jacobian_pushed_forward_3rd_derivatives
std::vector< DerivativeForm< 2, dim, spacedim > > jacobian_grads
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
Covariant transformation.
std::vector< Tensor< 1, spacedim > > normal_vectors
double weight(const unsigned int i) const
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
UpdateFlags update_each
Definition: mapping.h:628
static ::ExceptionBase & ExcInternalError()
static DataSetDescriptor face(const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
Definition: quadrature.cc:1139