Reference documentation for deal.II version Git fe553f7db3 2020-12-03 08:11:48 +0100
\(\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_cartesian.cc
Go to the documentation of this file.
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 
21 #include <deal.II/base/tensor.h>
22 
24 
25 #include <deal.II/fe/fe_values.h>
27 
28 #include <deal.II/grid/tria.h>
30 
32 
33 #include <algorithm>
34 #include <cmath>
35 #include <memory>
36 
37 
39 
40 
41 
42 template <int dim, int spacedim>
44  const Quadrature<dim> &q)
45  : cell_extents(numbers::signaling_nan<Tensor<1, dim>>())
46  , volume_element(numbers::signaling_nan<double>())
47  , quadrature_points(q.get_points())
48 {}
49 
50 
51 
52 template <int dim, int spacedim>
53 std::size_t
55 {
60 }
61 
62 
63 
64 template <int dim, int spacedim>
65 bool
67 {
68  return true;
69 }
70 
71 
72 
73 template <int dim, int spacedim>
76  const UpdateFlags in) const
77 {
78  // this mapping is pretty simple in that it can basically compute
79  // every piece of information wanted by FEValues without requiring
80  // computing any other quantities. boundary forms are one exception
81  // since they can be computed from the normal vectors without much
82  // further ado
83  UpdateFlags out = in;
84  if (out & update_boundary_forms)
85  out |= update_normal_vectors;
86 
87  return out;
88 }
89 
90 
91 
92 template <int dim, int spacedim>
93 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
95  const Quadrature<dim> &q) const
96 {
97  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
98  std::make_unique<InternalData>(q);
99  auto &data = dynamic_cast<InternalData &>(*data_ptr);
100 
101  // store the flags in the internal data object so we can access them
102  // in fill_fe_*_values(). use the transitive hull of the required
103  // flags
104  data.update_each = requires_update_flags(update_flags);
105 
106  return data_ptr;
107 }
108 
109 
110 
111 template <int dim, int spacedim>
112 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
114  const UpdateFlags update_flags,
115  const Quadrature<dim - 1> &quadrature) const
116 {
117  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
118  std::make_unique<InternalData>(
120  quadrature));
121  auto &data = dynamic_cast<InternalData &>(*data_ptr);
122 
123  // verify that we have computed the transitive hull of the required
124  // flags and that FEValues has faithfully passed them on to us
125  Assert(update_flags == requires_update_flags(update_flags),
126  ExcInternalError());
127 
128  // store the flags in the internal data object so we can access them
129  // in fill_fe_*_values()
130  data.update_each = update_flags;
131 
132  return data_ptr;
133 }
134 
135 
136 
137 template <int dim, int spacedim>
138 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
140  const UpdateFlags update_flags,
141  const Quadrature<dim - 1> &quadrature) const
142 {
143  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
144  std::make_unique<InternalData>(QProjector<dim>::project_to_all_subfaces(
145  ReferenceCell::get_hypercube(dim), quadrature));
146  auto &data = dynamic_cast<InternalData &>(*data_ptr);
147 
148  // verify that we have computed the transitive hull of the required
149  // flags and that FEValues has faithfully passed them on to us
150  Assert(update_flags == requires_update_flags(update_flags),
151  ExcInternalError());
152 
153  // store the flags in the internal data object so we can access them
154  // in fill_fe_*_values()
155  data.update_each = update_flags;
156 
157  return data_ptr;
158 }
159 
160 
161 
162 template <int dim, int spacedim>
163 void
165  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
166  const CellSimilarity::Similarity cell_similarity,
167  const InternalData & data) const
168 {
169  // Compute start point and sizes
170  // along axes. Strange vertex
171  // numbering makes this complicated
172  // again.
173  if (cell_similarity != CellSimilarity::translation)
174  {
175  const Point<dim> &start = cell->vertex(0);
176  switch (dim)
177  {
178  case 1:
179  data.cell_extents[0] = cell->vertex(1)(0) - start(0);
180  break;
181  case 2:
182  data.cell_extents[0] = cell->vertex(1)(0) - start(0);
183  data.cell_extents[1] = cell->vertex(2)(1) - start(1);
184  break;
185  case 3:
186  data.cell_extents[0] = cell->vertex(1)(0) - start(0);
187  data.cell_extents[1] = cell->vertex(2)(1) - start(1);
188  data.cell_extents[2] = cell->vertex(4)(2) - start(2);
189  break;
190  default:
191  Assert(false, ExcNotImplemented());
192  }
193  }
194 }
195 
196 
197 
198 template <int dim, int spacedim>
199 void
201  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
202  const InternalData & data,
203  std::vector<Point<dim>> &quadrature_points) const
204 {
206  {
207  const auto offset = QProjector<dim>::DataSetDescriptor::cell();
208 
209  transform_quadrature_points(cell, data, offset, quadrature_points);
210  }
211 }
212 
213 
214 
215 template <int dim, int spacedim>
216 void
218  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
219  const unsigned int face_no,
220  const InternalData & data,
221  std::vector<Point<dim>> &quadrature_points) const
222 {
224 
226  {
227  const auto offset = QProjector<dim>::DataSetDescriptor::face(
229  face_no,
230  cell->face_orientation(face_no),
231  cell->face_flip(face_no),
232  cell->face_rotation(face_no),
233  quadrature_points.size());
234 
235 
236  transform_quadrature_points(cell, data, offset, quadrature_points);
237  }
238 }
239 
240 
241 
242 template <int dim, int spacedim>
243 void
245  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
246  const unsigned int face_no,
247  const unsigned int sub_no,
248  const InternalData & data,
249  std::vector<Point<dim>> &quadrature_points) const
250 {
253  if (cell->face(face_no)->has_children())
254  {
255  AssertIndexRange(sub_no, cell->face(face_no)->n_children());
256  }
257 
259  {
262  face_no,
263  sub_no,
264  cell->face_orientation(face_no),
265  cell->face_flip(face_no),
266  cell->face_rotation(face_no),
267  quadrature_points.size(),
268  cell->subface_case(face_no));
269 
270  transform_quadrature_points(cell, data, offset, quadrature_points);
271  }
272 }
273 
274 
275 
276 template <int dim, int spacedim>
277 void
279  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
280  const InternalData & data,
281  const typename QProjector<dim>::DataSetDescriptor & offset,
282  std::vector<Point<dim>> &quadrature_points) const
283 {
284  // let @p{start} be the origin of a local coordinate system. it is chosen
285  // as the (lower) left vertex
286  const Point<dim> &start = cell->vertex(0);
287 
288  for (unsigned int i = 0; i < quadrature_points.size(); ++i)
289  {
290  quadrature_points[i] = start;
291  for (unsigned int d = 0; d < dim; ++d)
292  quadrature_points[i](d) +=
293  data.cell_extents[d] * data.quadrature_points[i + offset](d);
294  }
295 }
296 
297 
298 
299 template <int dim, int spacedim>
300 void
302  const unsigned int face_no,
303  const InternalData & data,
304  std::vector<Tensor<1, dim>> &normal_vectors) const
305 {
306  // compute normal vectors. All normals on a face have the same value.
308  {
310  std::fill(normal_vectors.begin(),
311  normal_vectors.end(),
313  }
314 }
315 
316 
317 
318 template <int dim, int spacedim>
319 void
321  const InternalData & data,
322  const CellSimilarity::Similarity cell_similarity,
324  &output_data) const
325 {
326  if (cell_similarity != CellSimilarity::translation)
327  {
329  for (unsigned int i = 0; i < output_data.jacobian_grads.size(); ++i)
331 
333  for (unsigned int i = 0;
334  i < output_data.jacobian_pushed_forward_grads.size();
335  ++i)
337 
339  for (unsigned int i = 0;
340  i < output_data.jacobian_2nd_derivatives.size();
341  ++i)
342  output_data.jacobian_2nd_derivatives[i] =
344 
346  for (unsigned int i = 0;
347  i < output_data.jacobian_pushed_forward_2nd_derivatives.size();
348  ++i)
351 
353  for (unsigned int i = 0;
354  i < output_data.jacobian_3rd_derivatives.size();
355  ++i)
356  output_data.jacobian_3rd_derivatives[i] =
358 
360  for (unsigned int i = 0;
361  i < output_data.jacobian_pushed_forward_3rd_derivatives.size();
362  ++i)
365  }
366 }
367 
368 
369 
370 template <int dim, int spacedim>
373  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
374  const CellSimilarity::Similarity cell_similarity,
375  const Quadrature<dim> & quadrature,
376  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
378  &output_data) const
379 {
380  // convert data object to internal data for this class. fails with
381  // an exception if that is not possible
382  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
383  ExcInternalError());
384  const InternalData &data = static_cast<const InternalData &>(internal_data);
385 
386 
387  update_cell_extents(cell, cell_similarity, data);
388 
390  data,
391  output_data.quadrature_points);
392 
393  // compute Jacobian determinant. all values are equal and are the
394  // product of the local lengths in each coordinate direction
396  if (cell_similarity != CellSimilarity::translation)
397  {
398  double J = data.cell_extents[0];
399  for (unsigned int d = 1; d < dim; ++d)
400  J *= data.cell_extents[d];
401  data.volume_element = J;
402  if (data.update_each & update_JxW_values)
403  for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
404  output_data.JxW_values[i] = J * quadrature.weight(i);
405  }
406  // "compute" Jacobian at the quadrature points, which are all the
407  // same
408  if (data.update_each & update_jacobians)
409  if (cell_similarity != CellSimilarity::translation)
410  for (unsigned int i = 0; i < output_data.jacobians.size(); ++i)
411  {
412  output_data.jacobians[i] = DerivativeForm<1, dim, spacedim>();
413  for (unsigned int j = 0; j < dim; ++j)
414  output_data.jacobians[i][j][j] = data.cell_extents[j];
415  }
416 
417  maybe_update_jacobian_derivatives(data, cell_similarity, output_data);
418 
419  // "compute" inverse Jacobian at the quadrature points, which are
420  // all the same
422  if (cell_similarity != CellSimilarity::translation)
423  for (unsigned int i = 0; i < output_data.inverse_jacobians.size(); ++i)
424  {
425  output_data.inverse_jacobians[i] = Tensor<2, dim>();
426  for (unsigned int j = 0; j < dim; ++j)
427  output_data.inverse_jacobians[i][j][j] = 1. / data.cell_extents[j];
428  }
429 
430  return cell_similarity;
431 }
432 
433 
434 
435 template <int dim, int spacedim>
436 void
438  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
439  const unsigned int face_no,
440  const Quadrature<dim - 1> & quadrature,
441  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
443  &output_data) const
444 {
445  // convert data object to internal
446  // data for this class. fails with
447  // an exception if that is not
448  // possible
449  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
450  ExcInternalError());
451  const InternalData &data = static_cast<const InternalData &>(internal_data);
452 
454 
456  face_no,
457  data,
458  output_data.quadrature_points);
459 
460  maybe_update_normal_vectors(face_no, data, output_data.normal_vectors);
461 
462  // first compute Jacobian determinant, which is simply the product
463  // of the local lengths since the jacobian is diagonal
464  double J = 1.;
465  for (unsigned int d = 0; d < dim; ++d)
467  J *= data.cell_extents[d];
468 
469  if (data.update_each & update_JxW_values)
470  for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
471  output_data.JxW_values[i] = J * quadrature.weight(i);
472 
474  for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
475  output_data.boundary_forms[i] = J * output_data.normal_vectors[i];
476 
478  {
479  J = data.cell_extents[0];
480  for (unsigned int d = 1; d < dim; ++d)
481  J *= data.cell_extents[d];
482  data.volume_element = J;
483  }
484 
485  if (data.update_each & update_jacobians)
486  for (unsigned int i = 0; i < output_data.jacobians.size(); ++i)
487  {
488  output_data.jacobians[i] = DerivativeForm<1, dim, spacedim>();
489  for (unsigned int d = 0; d < dim; ++d)
490  output_data.jacobians[i][d][d] = data.cell_extents[d];
491  }
492 
494 
496  for (unsigned int i = 0; i < output_data.inverse_jacobians.size(); ++i)
497  {
499  for (unsigned int d = 0; d < dim; ++d)
500  output_data.inverse_jacobians[i][d][d] = 1. / data.cell_extents[d];
501  }
502 }
503 
504 
505 
506 template <int dim, int spacedim>
507 void
509  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
510  const unsigned int face_no,
511  const unsigned int subface_no,
512  const Quadrature<dim - 1> & quadrature,
513  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
515  &output_data) const
516 {
517  // convert data object to internal data for this class. fails with
518  // an exception if that is not possible
519  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
520  ExcInternalError());
521  const InternalData &data = static_cast<const InternalData &>(internal_data);
522 
524 
526  cell, face_no, subface_no, data, output_data.quadrature_points);
527 
528  maybe_update_normal_vectors(face_no, data, output_data.normal_vectors);
529 
530  // first compute Jacobian determinant, which is simply the product
531  // of the local lengths since the jacobian is diagonal
532  double J = 1.;
533  for (unsigned int d = 0; d < dim; ++d)
535  J *= data.cell_extents[d];
536 
537  if (data.update_each & update_JxW_values)
538  {
539  // Here, cell->face(face_no)->n_children() would be the right
540  // choice, but unfortunately the current function is also called
541  // for faces without children (see tests/fe/mapping.cc). Add
542  // following switch to avoid diffs in tests/fe/mapping.OK
543  const unsigned int n_subfaces =
544  cell->face(face_no)->has_children() ?
545  cell->face(face_no)->n_children() :
547  for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
548  output_data.JxW_values[i] = J * quadrature.weight(i) / n_subfaces;
549  }
550 
552  for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i)
553  output_data.boundary_forms[i] = J * output_data.normal_vectors[i];
554 
556  {
557  J = data.cell_extents[0];
558  for (unsigned int d = 1; d < dim; ++d)
559  J *= data.cell_extents[d];
560  data.volume_element = J;
561  }
562 
563  if (data.update_each & update_jacobians)
564  for (unsigned int i = 0; i < output_data.jacobians.size(); ++i)
565  {
566  output_data.jacobians[i] = DerivativeForm<1, dim, spacedim>();
567  for (unsigned int d = 0; d < dim; ++d)
568  output_data.jacobians[i][d][d] = data.cell_extents[d];
569  }
570 
572 
574  for (unsigned int i = 0; i < output_data.inverse_jacobians.size(); ++i)
575  {
577  for (unsigned int d = 0; d < dim; ++d)
578  output_data.inverse_jacobians[i][d][d] = 1. / data.cell_extents[d];
579  }
580 }
581 
582 
583 
584 template <int dim, int spacedim>
585 void
587  const ArrayView<const Tensor<1, dim>> & input,
588  const MappingKind mapping_kind,
589  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
590  const ArrayView<Tensor<1, spacedim>> & output) const
591 {
592  AssertDimension(input.size(), output.size());
593  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
594  ExcInternalError());
595  const InternalData &data = static_cast<const InternalData &>(mapping_data);
596 
597  switch (mapping_kind)
598  {
599  case mapping_covariant:
600  {
603  "update_covariant_transformation"));
604 
605  for (unsigned int i = 0; i < output.size(); ++i)
606  for (unsigned int d = 0; d < dim; ++d)
607  output[i][d] = input[i][d] / data.cell_extents[d];
608  return;
609  }
610 
612  {
615  "update_contravariant_transformation"));
616 
617  for (unsigned int i = 0; i < output.size(); ++i)
618  for (unsigned int d = 0; d < dim; ++d)
619  output[i][d] = input[i][d] * data.cell_extents[d];
620  return;
621  }
622  case mapping_piola:
623  {
626  "update_contravariant_transformation"));
629  "update_volume_elements"));
630 
631  for (unsigned int i = 0; i < output.size(); ++i)
632  for (unsigned int d = 0; d < dim; ++d)
633  output[i][d] =
634  input[i][d] * data.cell_extents[d] / data.volume_element;
635  return;
636  }
637  default:
638  Assert(false, ExcNotImplemented());
639  }
640 }
641 
642 
643 
644 template <int dim, int spacedim>
645 void
647  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
648  const MappingKind mapping_kind,
649  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
650  const ArrayView<Tensor<2, spacedim>> & output) const
651 {
652  AssertDimension(input.size(), output.size());
653  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
654  ExcInternalError());
655  const InternalData &data = static_cast<const InternalData &>(mapping_data);
656 
657  switch (mapping_kind)
658  {
659  case mapping_covariant:
660  {
663  "update_covariant_transformation"));
664 
665  for (unsigned int i = 0; i < output.size(); ++i)
666  for (unsigned int d1 = 0; d1 < dim; ++d1)
667  for (unsigned int d2 = 0; d2 < dim; ++d2)
668  output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2];
669  return;
670  }
671 
673  {
676  "update_contravariant_transformation"));
677 
678  for (unsigned int i = 0; i < output.size(); ++i)
679  for (unsigned int d1 = 0; d1 < dim; ++d1)
680  for (unsigned int d2 = 0; d2 < dim; ++d2)
681  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2];
682  return;
683  }
684 
686  {
689  "update_covariant_transformation"));
690 
691  for (unsigned int i = 0; i < output.size(); ++i)
692  for (unsigned int d1 = 0; d1 < dim; ++d1)
693  for (unsigned int d2 = 0; d2 < dim; ++d2)
694  output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2] /
695  data.cell_extents[d1];
696  return;
697  }
698 
700  {
703  "update_contravariant_transformation"));
704 
705  for (unsigned int i = 0; i < output.size(); ++i)
706  for (unsigned int d1 = 0; d1 < dim; ++d1)
707  for (unsigned int d2 = 0; d2 < dim; ++d2)
708  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
709  data.cell_extents[d1];
710  return;
711  }
712 
713  case mapping_piola:
714  {
717  "update_contravariant_transformation"));
720  "update_volume_elements"));
721 
722  for (unsigned int i = 0; i < output.size(); ++i)
723  for (unsigned int d1 = 0; d1 < dim; ++d1)
724  for (unsigned int d2 = 0; d2 < dim; ++d2)
725  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
726  data.volume_element;
727  return;
728  }
729 
731  {
734  "update_contravariant_transformation"));
737  "update_volume_elements"));
738 
739  for (unsigned int i = 0; i < output.size(); ++i)
740  for (unsigned int d1 = 0; d1 < dim; ++d1)
741  for (unsigned int d2 = 0; d2 < dim; ++d2)
742  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
743  data.cell_extents[d1] / data.volume_element;
744  return;
745  }
746 
747  default:
748  Assert(false, ExcNotImplemented());
749  }
750 }
751 
752 
753 
754 template <int dim, int spacedim>
755 void
757  const ArrayView<const Tensor<2, dim>> & input,
758  const MappingKind mapping_kind,
759  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
760  const ArrayView<Tensor<2, spacedim>> & output) const
761 {
762  AssertDimension(input.size(), output.size());
763  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
764  ExcInternalError());
765  const InternalData &data = static_cast<const InternalData &>(mapping_data);
766 
767  switch (mapping_kind)
768  {
769  case mapping_covariant:
770  {
773  "update_covariant_transformation"));
774 
775  for (unsigned int i = 0; i < output.size(); ++i)
776  for (unsigned int d1 = 0; d1 < dim; ++d1)
777  for (unsigned int d2 = 0; d2 < dim; ++d2)
778  output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2];
779  return;
780  }
781 
783  {
786  "update_contravariant_transformation"));
787 
788  for (unsigned int i = 0; i < output.size(); ++i)
789  for (unsigned int d1 = 0; d1 < dim; ++d1)
790  for (unsigned int d2 = 0; d2 < dim; ++d2)
791  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2];
792  return;
793  }
794 
796  {
799  "update_covariant_transformation"));
800 
801  for (unsigned int i = 0; i < output.size(); ++i)
802  for (unsigned int d1 = 0; d1 < dim; ++d1)
803  for (unsigned int d2 = 0; d2 < dim; ++d2)
804  output[i][d1][d2] = input[i][d1][d2] / data.cell_extents[d2] /
805  data.cell_extents[d1];
806  return;
807  }
808 
810  {
813  "update_contravariant_transformation"));
814 
815  for (unsigned int i = 0; i < output.size(); ++i)
816  for (unsigned int d1 = 0; d1 < dim; ++d1)
817  for (unsigned int d2 = 0; d2 < dim; ++d2)
818  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
819  data.cell_extents[d1];
820  return;
821  }
822 
823  case mapping_piola:
824  {
827  "update_contravariant_transformation"));
830  "update_volume_elements"));
831 
832  for (unsigned int i = 0; i < output.size(); ++i)
833  for (unsigned int d1 = 0; d1 < dim; ++d1)
834  for (unsigned int d2 = 0; d2 < dim; ++d2)
835  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
836  data.volume_element;
837  return;
838  }
839 
841  {
844  "update_contravariant_transformation"));
847  "update_volume_elements"));
848 
849  for (unsigned int i = 0; i < output.size(); ++i)
850  for (unsigned int d1 = 0; d1 < dim; ++d1)
851  for (unsigned int d2 = 0; d2 < dim; ++d2)
852  output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] /
853  data.cell_extents[d1] / data.volume_element;
854  return;
855  }
856 
857  default:
858  Assert(false, ExcNotImplemented());
859  }
860 }
861 
862 
863 
864 template <int dim, int spacedim>
865 void
867  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
868  const MappingKind mapping_kind,
869  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
870  const ArrayView<Tensor<3, spacedim>> & output) const
871 {
872  AssertDimension(input.size(), output.size());
873  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
874  ExcInternalError());
875  const InternalData &data = static_cast<const InternalData &>(mapping_data);
876 
877  switch (mapping_kind)
878  {
880  {
883  "update_covariant_transformation"));
884 
885  for (unsigned int q = 0; q < output.size(); ++q)
886  for (unsigned int i = 0; i < spacedim; ++i)
887  for (unsigned int j = 0; j < spacedim; ++j)
888  for (unsigned int k = 0; k < spacedim; ++k)
889  {
890  output[q][i][j][k] = input[q][i][j][k] /
891  data.cell_extents[j] /
892  data.cell_extents[k];
893  }
894  return;
895  }
896  default:
897  Assert(false, ExcNotImplemented());
898  }
899 }
900 
901 
902 
903 template <int dim, int spacedim>
904 void
906  const ArrayView<const Tensor<3, dim>> & input,
907  const MappingKind mapping_kind,
908  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
909  const ArrayView<Tensor<3, spacedim>> & output) const
910 {
911  AssertDimension(input.size(), output.size());
912  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
913  ExcInternalError());
914  const InternalData &data = static_cast<const InternalData &>(mapping_data);
915 
916  switch (mapping_kind)
917  {
919  {
922  "update_covariant_transformation"));
925  "update_contravariant_transformation"));
926 
927  for (unsigned int q = 0; q < output.size(); ++q)
928  for (unsigned int i = 0; i < spacedim; ++i)
929  for (unsigned int j = 0; j < spacedim; ++j)
930  for (unsigned int k = 0; k < spacedim; ++k)
931  {
932  output[q][i][j][k] =
933  input[q][i][j][k] * data.cell_extents[i] /
934  data.cell_extents[j] / data.cell_extents[k];
935  }
936  return;
937  }
938 
940  {
943  "update_covariant_transformation"));
944 
945  for (unsigned int q = 0; q < output.size(); ++q)
946  for (unsigned int i = 0; i < spacedim; ++i)
947  for (unsigned int j = 0; j < spacedim; ++j)
948  for (unsigned int k = 0; k < spacedim; ++k)
949  {
950  output[q][i][j][k] =
951  input[q][i][j][k] / data.cell_extents[i] /
952  data.cell_extents[j] / data.cell_extents[k];
953  }
954 
955  return;
956  }
957 
959  {
962  "update_covariant_transformation"));
965  "update_contravariant_transformation"));
968  "update_volume_elements"));
969 
970  for (unsigned int q = 0; q < output.size(); ++q)
971  for (unsigned int i = 0; i < spacedim; ++i)
972  for (unsigned int j = 0; j < spacedim; ++j)
973  for (unsigned int k = 0; k < spacedim; ++k)
974  {
975  output[q][i][j][k] =
976  input[q][i][j][k] * data.cell_extents[i] /
977  data.volume_element / data.cell_extents[j] /
978  data.cell_extents[k];
979  }
980 
981  return;
982  }
983 
984  default:
985  Assert(false, ExcNotImplemented());
986  }
987 }
988 
989 
990 
991 template <int dim, int spacedim>
994  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
995  const Point<dim> & p) const
996 {
997  Tensor<1, dim> length;
998  const Point<dim> start = cell->vertex(0);
999  switch (dim)
1000  {
1001  case 1:
1002  length[0] = cell->vertex(1)(0) - start(0);
1003  break;
1004  case 2:
1005  length[0] = cell->vertex(1)(0) - start(0);
1006  length[1] = cell->vertex(2)(1) - start(1);
1007  break;
1008  case 3:
1009  length[0] = cell->vertex(1)(0) - start(0);
1010  length[1] = cell->vertex(2)(1) - start(1);
1011  length[2] = cell->vertex(4)(2) - start(2);
1012  break;
1013  default:
1014  Assert(false, ExcNotImplemented());
1015  }
1016 
1017  Point<dim> p_real = cell->vertex(0);
1018  for (unsigned int d = 0; d < dim; ++d)
1019  p_real(d) += length[d] * p(d);
1020 
1021  return p_real;
1022 }
1023 
1024 
1025 
1026 template <int dim, int spacedim>
1027 Point<dim>
1029  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1030  const Point<spacedim> & p) const
1031 {
1032  if (dim != spacedim)
1033  Assert(false, ExcNotImplemented());
1034  const Point<dim> &start = cell->vertex(0);
1035  Point<dim> real = p;
1036  real -= start;
1037 
1038  switch (dim)
1039  {
1040  case 1:
1041  real(0) /= cell->vertex(1)(0) - start(0);
1042  break;
1043  case 2:
1044  real(0) /= cell->vertex(1)(0) - start(0);
1045  real(1) /= cell->vertex(2)(1) - start(1);
1046  break;
1047  case 3:
1048  real(0) /= cell->vertex(1)(0) - start(0);
1049  real(1) /= cell->vertex(2)(1) - start(1);
1050  real(2) /= cell->vertex(4)(2) - start(2);
1051  break;
1052  default:
1053  Assert(false, ExcNotImplemented());
1054  }
1055  return real;
1056 }
1057 
1058 
1059 
1060 template <int dim, int spacedim>
1061 std::unique_ptr<Mapping<dim, spacedim>>
1063 {
1064  return std::make_unique<MappingCartesian<dim, spacedim>>(*this);
1065 }
1066 
1067 
1068 //---------------------------------------------------------------------------
1069 // explicit instantiations
1070 #include "mapping_cartesian.inst"
1071 
1072 
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
virtual void transform(const ArrayView< const Tensor< 1, dim >> &input, const MappingKind kind, const typename Mapping< dim, spacedim >::InternalDataBase &internal, const ArrayView< Tensor< 1, spacedim >> &output) const override
void maybe_update_normal_vectors(const unsigned int face_no, const InternalData &data, std::vector< Tensor< 1, dim >> &normal_vectors) const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1623
Contravariant transformation.
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=StaticMappingQ1< dim, spacedim >::mapping, const std::vector< std::vector< double >> &properties={})
Definition: generators.cc:444
virtual CellSimilarity::Similarity fill_fe_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const Quadrature< dim > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
void update_cell_extents(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const InternalData &data) const
Type get_hypercube(const unsigned int dim)
void transform_quadrature_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const InternalData &data, const typename QProjector< dim >::DataSetDescriptor &offset, std::vector< Point< dim >> &quadrature_points) const
std::vector< Tensor< 1, spacedim > > boundary_forms
Volume element.
#define AssertIndexRange(index, range)
Definition: exceptions.h:1691
std::vector< DerivativeForm< 3, dim, spacedim > > jacobian_2nd_derivatives
Outer normal vector, not normalized.
virtual void fill_fe_subface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
std::vector< Point< dim > > quadrature_points
Determinant of the Jacobian.
std::vector< DerivativeForm< 4, dim, spacedim > > jacobian_3rd_derivatives
Transformed quadrature points.
void maybe_update_face_quadrature_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const InternalData &data, std::vector< Point< dim >> &quadrature_points) const
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
MappingKind
Definition: mapping.h:62
static DataSetDescriptor cell()
Definition: qprojector.h:564
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:1466
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
UpdateFlags
Abstract base class for mapping classes.
Definition: mapping.h:301
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
void maybe_update_jacobian_derivatives(const InternalData &data, const CellSimilarity::Similarity cell_similarity, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
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
void maybe_update_subface_quadrature_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const InternalData &data, std::vector< Point< dim >> &quadrature_points) const
std::vector< Point< spacedim > > quadrature_points
static Quadrature< dim > project_to_all_faces(const SubQuadrature &quadrature)
void maybe_update_cell_quadrature_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const InternalData &data, std::vector< Point< dim >> &quadrature_points) const
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: tensor.h:448
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
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
T signaling_nan()
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
static DataSetDescriptor subface(const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points, const internal::SubfaceCase< dim > ref_case=internal::SubfaceCase< dim >::case_isotropic)
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:643
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: qprojector.cc:1184