Reference documentation for deal.II version Git 840703ecf3 2019-10-14 17:15:20 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
mapping_manifold.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 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 
17 #include <deal.II/base/array_view.h>
18 #include <deal.II/base/derivative_form.h>
19 #include <deal.II/base/memory_consumption.h>
20 #include <deal.II/base/qprojector.h>
21 #include <deal.II/base/quadrature.h>
22 #include <deal.II/base/quadrature_lib.h>
23 #include <deal.II/base/std_cxx14/memory.h>
24 #include <deal.II/base/tensor_product_polynomials.h>
25 
26 #include <deal.II/dofs/dof_accessor.h>
27 
28 #include <deal.II/fe/fe.h>
29 #include <deal.II/fe/fe_tools.h>
30 #include <deal.II/fe/fe_values.h>
31 #include <deal.II/fe/mapping_manifold.h>
32 #include <deal.II/fe/mapping_q1.h>
33 
34 #include <deal.II/grid/manifold.h>
35 #include <deal.II/grid/tria.h>
36 #include <deal.II/grid/tria_iterator.h>
37 
38 #include <deal.II/lac/full_matrix.h>
39 
40 #include <algorithm>
41 #include <array>
42 #include <cmath>
43 #include <memory>
44 #include <numeric>
45 
46 
47 DEAL_II_NAMESPACE_OPEN
48 
49 template <int dim, int spacedim>
50 std::size_t
52 {
53  return (
66 }
67 
68 
69 template <int dim, int spacedim>
70 void
72  const UpdateFlags update_flags,
73  const Quadrature<dim> &q,
74  const unsigned int n_original_q_points)
75 {
76  // store the flags in the internal data object so we can access them
77  // in fill_fe_*_values()
78  this->update_each = update_flags;
79 
80  // Store the quadrature
81  this->quad = q;
82 
83  // Resize the weights
85 
86  // see if we need the (transformation) shape function values
87  // and/or gradients and resize the necessary arrays
88  if (this->update_each &
91 
93  covariant.resize(n_original_q_points);
94 
96  contravariant.resize(n_original_q_points);
97 
99  volume_elements.resize(n_original_q_points);
100 }
101 
102 
103 template <int dim, int spacedim>
104 void
106  const UpdateFlags update_flags,
107  const Quadrature<dim> &q,
108  const unsigned int n_original_q_points)
109 {
110  initialize(update_flags, q, n_original_q_points);
111 
112  if (dim > 1)
113  {
115  {
116  aux.resize(dim - 1,
117  std::vector<Tensor<1, spacedim>>(n_original_q_points));
118 
119  // Compute tangentials to the unit cell.
120  for (unsigned int i = 0; i < unit_tangentials.size(); ++i)
121  unit_tangentials[i].resize(n_original_q_points);
122  switch (dim)
123  {
124  case 2:
125  {
126  // ensure a counterclockwise
127  // orientation of tangentials
128  static const int tangential_orientation[4] = {-1, 1, 1, -1};
129  for (unsigned int i = 0;
130  i < GeometryInfo<dim>::faces_per_cell;
131  ++i)
132  {
133  Tensor<1, dim> tang;
134  tang[1 - i / 2] = tangential_orientation[i];
135  std::fill(unit_tangentials[i].begin(),
136  unit_tangentials[i].end(),
137  tang);
138  }
139  break;
140  }
141  case 3:
142  {
143  for (unsigned int i = 0;
144  i < GeometryInfo<dim>::faces_per_cell;
145  ++i)
146  {
147  Tensor<1, dim> tang1, tang2;
148 
149  const unsigned int nd =
151 
152  // first tangential
153  // vector in direction
154  // of the (nd+1)%3 axis
155  // and inverted in case
156  // of unit inward normal
157  tang1[(nd + 1) % dim] =
159  // second tangential
160  // vector in direction
161  // of the (nd+2)%3 axis
162  tang2[(nd + 2) % dim] = 1.;
163 
164  // same unit tangents
165  // for all quadrature
166  // points on this face
167  std::fill(unit_tangentials[i].begin(),
168  unit_tangentials[i].end(),
169  tang1);
170  std::fill(
172  .begin(),
174  .end(),
175  tang2);
176  }
177  break;
178  }
179  default:
180  Assert(false, ExcNotImplemented());
181  }
182  }
183  }
184 }
185 
186 
187 
188 template <int dim, int spacedim>
191 {}
192 
193 
194 
195 template <int dim, int spacedim>
196 std::unique_ptr<Mapping<dim, spacedim>>
198 {
199  return std_cxx14::make_unique<MappingManifold<dim, spacedim>>(*this);
200 }
201 
202 
203 
204 template <int dim, int spacedim>
208  const Point<spacedim> &) const
209 {
210  Assert(false, ExcNotImplemented());
211  return Point<dim>();
212 }
213 
214 
215 
216 template <int dim, int spacedim>
220  const Point<dim> & p) const
221 {
222  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertices;
223  std::array<double, GeometryInfo<dim>::vertices_per_cell> weights;
224 
225  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
226  {
227  vertices[v] = cell->vertex(v);
229  }
230  return cell->get_manifold().get_new_point(
231  make_array_view(vertices.begin(), vertices.end()),
232  make_array_view(weights.begin(), weights.end()));
233 }
234 
235 
236 
237 // In the code below, GCC tries to instantiate MappingManifold<3,4> when
238 // seeing which of the overloaded versions of
239 // do_transform_real_to_unit_cell_internal() to call. This leads to bad
240 // error messages and, generally, nothing very good. Avoid this by ensuring
241 // that this class exists, but does not have an inner InternalData
242 // type, thereby ruling out the codim-1 version of the function
243 // below when doing overload resolution.
244 template <>
245 class MappingManifold<3, 4>
246 {};
247 
248 
249 
250 template <int dim, int spacedim>
253  const UpdateFlags in) const
254 {
255  // add flags if the respective quantities are necessary to compute
256  // what we need. note that some flags appear in both the conditions
257  // and in subsequent set operations. this leads to some circular
258  // logic. the only way to treat this is to iterate. since there are
259  // 5 if-clauses in the loop, it will take at most 5 iterations to
260  // converge. do them:
261  UpdateFlags out = in;
262  for (unsigned int i = 0; i < 5; ++i)
263  {
264  // The following is a little incorrect:
265  // If not applied on a face,
266  // update_boundary_forms does not
267  // make sense. On the other hand,
268  // it is necessary on a
269  // face. Currently,
270  // update_boundary_forms is simply
271  // ignored for the interior of a
272  // cell.
274  out |= update_boundary_forms;
275 
280 
281  if (out &
286 
287  // The contravariant transformation used in the Piola
288  // transformation, which requires the determinant of the Jacobi
289  // matrix of the transformation. Because we have no way of
290  // knowing here whether the finite elements wants to use the
291  // contravariant of the Piola transforms, we add the JxW values
292  // to the list of flags to be updated for each cell.
294  out |= update_JxW_values;
295 
296  if (out & update_normal_vectors)
297  out |= update_JxW_values;
298  }
299 
300  // Now throw an exception if we stumble upon something that was not
301  // implemented yet
306 
307  return out;
308 }
309 
310 
311 
312 template <int dim, int spacedim>
313 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
315  const Quadrature<dim> &q) const
316 {
317  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
318  std_cxx14::make_unique<InternalData>();
319  auto &data = dynamic_cast<InternalData &>(*data_ptr);
320  data.initialize(this->requires_update_flags(update_flags), q, q.size());
321 
322  return data_ptr;
323 }
324 
325 
326 
327 template <int dim, int spacedim>
328 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
330  const UpdateFlags update_flags,
331  const Quadrature<dim - 1> &quadrature) const
332 {
333  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
334  std_cxx14::make_unique<InternalData>();
335  auto &data = dynamic_cast<InternalData &>(*data_ptr);
336  data.initialize_face(this->requires_update_flags(update_flags),
338  quadrature.size());
339 
340  return data_ptr;
341 }
342 
343 
344 
345 template <int dim, int spacedim>
346 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
348  const UpdateFlags update_flags,
349  const Quadrature<dim - 1> &quadrature) const
350 {
351  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
352  std_cxx14::make_unique<InternalData>();
353  auto &data = dynamic_cast<InternalData &>(*data_ptr);
354  data.initialize_face(this->requires_update_flags(update_flags),
356  quadrature.size());
357 
358  return data_ptr;
359 }
360 
361 
362 
363 namespace internal
364 {
365  namespace MappingManifoldImplementation
366  {
367  namespace
368  {
375  template <int dim, int spacedim>
376  void
377  maybe_compute_q_points(
378  const typename QProjector<dim>::DataSetDescriptor data_set,
379  const typename ::MappingManifold<dim, spacedim>::InternalData
380  & data,
381  std::vector<Point<spacedim>> &quadrature_points)
382  {
383  const UpdateFlags update_flags = data.update_each;
384 
385  AssertDimension(data.vertices.size(),
387 
388  if (update_flags & update_quadrature_points)
389  {
390  for (unsigned int point = 0; point < quadrature_points.size();
391  ++point)
392  {
393  quadrature_points[point] = data.manifold->get_new_point(
394  make_array_view(data.vertices),
395  make_array_view(
396  data.cell_manifold_quadrature_weights[point + data_set]));
397  }
398  }
399  }
400 
401 
402 
409  template <int dim, int spacedim>
410  void
411  maybe_update_Jacobians(
412  const typename ::QProjector<dim>::DataSetDescriptor data_set,
413  const typename ::MappingManifold<dim, spacedim>::InternalData
414  &data)
415  {
416  const UpdateFlags update_flags = data.update_each;
417 
418  if (update_flags & update_contravariant_transformation)
419  {
420  const unsigned int n_q_points = data.contravariant.size();
421 
422  std::fill(data.contravariant.begin(),
423  data.contravariant.end(),
425 
427  data.vertices.size());
428  for (unsigned int point = 0; point < n_q_points; ++point)
429  {
430  // Start by figuring out how to compute the direction in
431  // the reference space:
432  const Point<dim> &p = data.quad.point(point + data_set);
433 
434  // And get its image on the manifold:
435  const Point<spacedim> P = data.manifold->get_new_point(
436  make_array_view(data.vertices),
437  make_array_view(
438  data.cell_manifold_quadrature_weights[point + data_set]));
439 
440  // To compute the Jacobian, we choose dim points aligned
441  // with the dim reference axes, which are still in the
442  // given cell, and ask for the tangent vector in these
443  // directions. Choosing the points is somewhat arbitrary,
444  // so we try to be smart and we pick points which are
445  // on the opposite quadrant w.r.t. the evaluation
446  // point.
447  for (unsigned int i = 0; i < dim; ++i)
448  {
449  const Point<dim> ei = Point<dim>::unit_vector(i);
450  const double pi = p[i];
451  Assert(pi >= 0 && pi <= 1.0,
453  "Was expecting a quadrature point "
454  "inside the unit reference element."));
455 
456  // In the length L, we store also the direction sign,
457  // which is positive, if the coordinate is < .5,
458  const double L = pi > .5 ? -pi : 1 - pi;
459 
460  const Point<dim> np(p + L * ei);
461 
462  // Get the weights to compute the np point in real space
463  for (unsigned int j = 0;
464  j < GeometryInfo<dim>::vertices_per_cell;
465  ++j)
466  data.vertex_weights[j] =
468 
469  const Point<spacedim> NP = data.manifold->get_new_point(
470  make_array_view(data.vertices),
471  make_array_view(data.vertex_weights));
472 
473  const Tensor<1, spacedim> T =
474  data.manifold->get_tangent_vector(P, NP);
475 
476  for (unsigned int d = 0; d < spacedim; ++d)
477  data.contravariant[point][d][i] = T[d] / L;
478  }
479  }
480 
481  if (update_flags & update_covariant_transformation)
482  {
483  const unsigned int n_q_points = data.contravariant.size();
484  for (unsigned int point = 0; point < n_q_points; ++point)
485  {
486  data.covariant[point] =
487  (data.contravariant[point]).covariant_form();
488  }
489  }
490 
491  if (update_flags & update_volume_elements)
492  {
493  const unsigned int n_q_points = data.contravariant.size();
494  for (unsigned int point = 0; point < n_q_points; ++point)
495  data.volume_elements[point] =
496  data.contravariant[point].determinant();
497  }
498  }
499  }
500  } // namespace
501  } // namespace MappingManifoldImplementation
502 } // namespace internal
503 
504 
505 
506 template <int dim, int spacedim>
510  const CellSimilarity::Similarity cell_similarity,
511  const Quadrature<dim> & quadrature,
512  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
514  &output_data) const
515 {
516  // ensure that the following static_cast is really correct:
517  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
518  ExcInternalError());
519  const InternalData &data = static_cast<const InternalData &>(internal_data);
520 
521  const unsigned int n_q_points = quadrature.size();
522 
523  data.store_vertices(cell);
524  data.manifold = &(cell->get_manifold());
525 
526  internal::MappingManifoldImplementation::maybe_compute_q_points<dim,
527  spacedim>(
529  data,
530  output_data.quadrature_points);
531 
532  internal::MappingManifoldImplementation::maybe_update_Jacobians<dim,
533  spacedim>(
535 
536  const UpdateFlags update_flags = data.update_each;
537  const std::vector<double> &weights = quadrature.get_weights();
538 
539  // Multiply quadrature weights by absolute value of Jacobian determinants or
540  // the area element g=sqrt(DX^t DX) in case of codim > 0
541 
542  if (update_flags & (update_normal_vectors | update_JxW_values))
543  {
544  AssertDimension(output_data.JxW_values.size(), n_q_points);
545 
546  Assert(!(update_flags & update_normal_vectors) ||
547  (output_data.normal_vectors.size() == n_q_points),
548  ExcDimensionMismatch(output_data.normal_vectors.size(),
549  n_q_points));
550 
551 
552  for (unsigned int point = 0; point < n_q_points; ++point)
553  {
554  if (dim == spacedim)
555  {
556  const double det = data.contravariant[point].determinant();
557 
558  // check for distorted cells.
559 
560  // TODO: this allows for anisotropies of up to 1e6 in 3D and
561  // 1e12 in 2D. might want to find a finer
562  // (dimension-independent) criterion
563  Assert(det > 1e-12 * Utilities::fixed_power<dim>(
564  cell->diameter() / std::sqrt(double(dim))),
566  cell->center(), det, point)));
567 
568  output_data.JxW_values[point] = weights[point] * det;
569  }
570  // if dim==spacedim, then there is no cell normal to
571  // compute. since this is for FEValues (and not FEFaceValues),
572  // there are also no face normals to compute
573  else // codim>0 case
574  {
575  Tensor<1, spacedim> DX_t[dim];
576  for (unsigned int i = 0; i < spacedim; ++i)
577  for (unsigned int j = 0; j < dim; ++j)
578  DX_t[j][i] = data.contravariant[point][i][j];
579 
580  Tensor<2, dim> G; // First fundamental form
581  for (unsigned int i = 0; i < dim; ++i)
582  for (unsigned int j = 0; j < dim; ++j)
583  G[i][j] = DX_t[i] * DX_t[j];
584 
585  output_data.JxW_values[point] =
586  std::sqrt(determinant(G)) * weights[point];
587 
588  if (cell_similarity == CellSimilarity::inverted_translation)
589  {
590  // we only need to flip the normal
591  if (update_flags & update_normal_vectors)
592  output_data.normal_vectors[point] *= -1.;
593  }
594  else
595  {
596  if (update_flags & update_normal_vectors)
597  {
598  Assert(spacedim == dim + 1,
599  ExcMessage(
600  "There is no (unique) cell normal for " +
602  "-dimensional cells in " +
603  Utilities::int_to_string(spacedim) +
604  "-dimensional space. This only works if the "
605  "space dimension is one greater than the "
606  "dimensionality of the mesh cells."));
607 
608  if (dim == 1)
609  output_data.normal_vectors[point] =
610  cross_product_2d(-DX_t[0]);
611  else // dim == 2
612  output_data.normal_vectors[point] =
613  cross_product_3d(DX_t[0], DX_t[1]);
614 
615  output_data.normal_vectors[point] /=
616  output_data.normal_vectors[point].norm();
617 
618  if (cell->direction_flag() == false)
619  output_data.normal_vectors[point] *= -1.;
620  }
621  }
622  } // codim>0 case
623  }
624  }
625 
626 
627 
628  // copy values from InternalData to vector given by reference
629  if (update_flags & update_jacobians)
630  {
631  AssertDimension(output_data.jacobians.size(), n_q_points);
632  if (cell_similarity != CellSimilarity::translation)
633  for (unsigned int point = 0; point < n_q_points; ++point)
634  output_data.jacobians[point] = data.contravariant[point];
635  }
636 
637  // copy values from InternalData to vector given by reference
638  if (update_flags & update_inverse_jacobians)
639  {
640  AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
641  if (cell_similarity != CellSimilarity::translation)
642  for (unsigned int point = 0; point < n_q_points; ++point)
643  output_data.inverse_jacobians[point] =
644  data.covariant[point].transpose();
645  }
646 
647  return cell_similarity;
648 }
649 
650 
651 
652 namespace internal
653 {
654  namespace MappingManifoldImplementation
655  {
656  namespace
657  {
667  template <int dim, int spacedim>
668  void
669  maybe_compute_face_data(
670  const ::MappingManifold<dim, spacedim> &mapping,
671  const typename ::Triangulation<dim, spacedim>::cell_iterator
672  & cell,
673  const unsigned int face_no,
674  const unsigned int subface_no,
675  const unsigned int n_q_points,
676  const std::vector<double> &weights,
677  const typename ::MappingManifold<dim, spacedim>::InternalData
678  &data,
680  &output_data)
681  {
682  const UpdateFlags update_flags = data.update_each;
683 
684  if (update_flags & update_boundary_forms)
685  {
686  AssertDimension(output_data.boundary_forms.size(), n_q_points);
687  if (update_flags & update_normal_vectors)
688  AssertDimension(output_data.normal_vectors.size(), n_q_points);
689  if (update_flags & update_JxW_values)
690  AssertDimension(output_data.JxW_values.size(), n_q_points);
691 
692  // map the unit tangentials to the real cell. checking for d!=dim-1
693  // eliminates compiler warnings regarding unsigned int expressions <
694  // 0.
695  for (unsigned int d = 0; d != dim - 1; ++d)
696  {
698  data.unit_tangentials.size(),
699  ExcInternalError());
700  Assert(
701  data.aux[d].size() <=
702  data
703  .unit_tangentials[face_no +
705  .size(),
706  ExcInternalError());
707 
708  mapping.transform(
709  make_array_view(
710  data
711  .unit_tangentials[face_no +
714  data,
715  make_array_view(data.aux[d]));
716  }
717 
718  // if dim==spacedim, we can use the unit tangentials to compute the
719  // boundary form by simply taking the cross product
720  if (dim == spacedim)
721  {
722  for (unsigned int i = 0; i < n_q_points; ++i)
723  switch (dim)
724  {
725  case 1:
726  // in 1d, we don't have access to any of the data.aux
727  // fields (because it has only dim-1 components), but we
728  // can still compute the boundary form by simply
729  // looking at the number of the face
730  output_data.boundary_forms[i][0] =
731  (face_no == 0 ? -1 : +1);
732  break;
733  case 2:
734  output_data.boundary_forms[i] =
735  cross_product_2d(data.aux[0][i]);
736  break;
737  case 3:
738  output_data.boundary_forms[i] =
739  cross_product_3d(data.aux[0][i], data.aux[1][i]);
740  break;
741  default:
742  Assert(false, ExcNotImplemented());
743  }
744  }
745  else //(dim < spacedim)
746  {
747  // in the codim-one case, the boundary form results from the
748  // cross product of all the face tangential vectors and the cell
749  // normal vector
750  //
751  // to compute the cell normal, use the same method used in
752  // fill_fe_values for cells above
753  AssertDimension(data.contravariant.size(), n_q_points);
754 
755  for (unsigned int point = 0; point < n_q_points; ++point)
756  {
757  switch (dim)
758  {
759  case 1:
760  {
761  // J is a tangent vector
762  output_data.boundary_forms[point] =
763  data.contravariant[point].transpose()[0];
764  output_data.boundary_forms[point] /=
765  (face_no == 0 ? -1. : +1.) *
766  output_data.boundary_forms[point].norm();
767 
768  break;
769  }
770 
771  case 2:
772  {
774  data.contravariant[point].transpose();
775 
776  Tensor<1, spacedim> cell_normal =
777  cross_product_3d(DX_t[0], DX_t[1]);
778  cell_normal /= cell_normal.norm();
779 
780  // then compute the face normal from the face
781  // tangent and the cell normal:
782  output_data.boundary_forms[point] =
783  cross_product_3d(data.aux[0][point], cell_normal);
784 
785  break;
786  }
787 
788  default:
789  Assert(false, ExcNotImplemented());
790  }
791  }
792  }
793 
794  if (update_flags & (update_normal_vectors | update_JxW_values))
795  for (unsigned int i = 0; i < output_data.boundary_forms.size();
796  ++i)
797  {
798  if (update_flags & update_JxW_values)
799  {
800  output_data.JxW_values[i] =
801  output_data.boundary_forms[i].norm() * weights[i];
802 
803  if (subface_no != numbers::invalid_unsigned_int)
804  {
805  const double area_ratio =
807  cell->subface_case(face_no), subface_no);
808  output_data.JxW_values[i] *= area_ratio;
809  }
810  }
811 
812  if (update_flags & update_normal_vectors)
813  output_data.normal_vectors[i] =
814  Point<spacedim>(output_data.boundary_forms[i] /
815  output_data.boundary_forms[i].norm());
816  }
817 
818  if (update_flags & update_jacobians)
819  for (unsigned int point = 0; point < n_q_points; ++point)
820  output_data.jacobians[point] = data.contravariant[point];
821 
822  if (update_flags & update_inverse_jacobians)
823  for (unsigned int point = 0; point < n_q_points; ++point)
824  output_data.inverse_jacobians[point] =
825  data.covariant[point].transpose();
826  }
827  }
828 
829 
836  template <int dim, int spacedim>
837  void
838  do_fill_fe_face_values(
839  const ::MappingManifold<dim, spacedim> &mapping,
840  const typename ::Triangulation<dim, spacedim>::cell_iterator
841  & cell,
842  const unsigned int face_no,
843  const unsigned int subface_no,
844  const typename QProjector<dim>::DataSetDescriptor data_set,
845  const Quadrature<dim - 1> & quadrature,
846  const typename ::MappingManifold<dim, spacedim>::InternalData
847  &data,
849  &output_data)
850  {
851  data.store_vertices(cell);
852 
853  data.manifold = &cell->face(face_no)->get_manifold();
854 
855  maybe_compute_q_points<dim, spacedim>(data_set,
856  data,
857  output_data.quadrature_points);
858  maybe_update_Jacobians<dim, spacedim>(data_set, data);
859 
860  maybe_compute_face_data(mapping,
861  cell,
862  face_no,
863  subface_no,
864  quadrature.size(),
865  quadrature.get_weights(),
866  data,
867  output_data);
868  }
869 
870  template <int dim, int spacedim, int rank>
871  void
872  transform_fields(
873  const ArrayView<const Tensor<rank, dim>> & input,
874  const MappingKind mapping_kind,
875  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
876  const ArrayView<Tensor<rank, spacedim>> & output)
877  {
878  AssertDimension(input.size(), output.size());
879  Assert((dynamic_cast<const typename ::
880  MappingManifold<dim, spacedim>::InternalData *>(
881  &mapping_data) != nullptr),
882  ExcInternalError());
883  const typename ::MappingManifold<dim, spacedim>::InternalData
884  &data =
885  static_cast<const typename ::MappingManifold<dim, spacedim>::
886  InternalData &>(mapping_data);
887 
888  switch (mapping_kind)
889  {
891  {
892  Assert(
893  data.update_each & update_contravariant_transformation,
895  "update_contravariant_transformation"));
896 
897  for (unsigned int i = 0; i < output.size(); ++i)
898  output[i] =
899  apply_transformation(data.contravariant[i], input[i]);
900 
901  return;
902  }
903 
904  case mapping_piola:
905  {
906  Assert(
907  data.update_each & update_contravariant_transformation,
909  "update_contravariant_transformation"));
910  Assert(
911  data.update_each & update_volume_elements,
913  "update_volume_elements"));
914  Assert(rank == 1, ExcMessage("Only for rank 1"));
915  if (rank != 1)
916  return;
917 
918  for (unsigned int i = 0; i < output.size(); ++i)
919  {
920  output[i] =
921  apply_transformation(data.contravariant[i], input[i]);
922  output[i] /= data.volume_elements[i];
923  }
924  return;
925  }
926  // We still allow this operation as in the
927  // reference cell Derivatives are Tensor
928  // rather than DerivativeForm
929  case mapping_covariant:
930  {
931  Assert(
932  data.update_each & update_contravariant_transformation,
934  "update_covariant_transformation"));
935 
936  for (unsigned int i = 0; i < output.size(); ++i)
937  output[i] = apply_transformation(data.covariant[i], input[i]);
938 
939  return;
940  }
941 
942  default:
943  Assert(false, ExcNotImplemented());
944  }
945  }
946 
947 
948  template <int dim, int spacedim, int rank>
949  void
950  transform_gradients(
951  const ArrayView<const Tensor<rank, dim>> & input,
952  const MappingKind mapping_kind,
953  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
954  const ArrayView<Tensor<rank, spacedim>> & output)
955  {
956  AssertDimension(input.size(), output.size());
957  Assert((dynamic_cast<const typename ::
958  MappingManifold<dim, spacedim>::InternalData *>(
959  &mapping_data) != nullptr),
960  ExcInternalError());
961  const typename ::MappingManifold<dim, spacedim>::InternalData
962  &data =
963  static_cast<const typename ::MappingManifold<dim, spacedim>::
964  InternalData &>(mapping_data);
965 
966  switch (mapping_kind)
967  {
969  {
970  Assert(
971  data.update_each & update_covariant_transformation,
973  "update_covariant_transformation"));
974  Assert(
975  data.update_each & update_contravariant_transformation,
977  "update_contravariant_transformation"));
978  Assert(rank == 2, ExcMessage("Only for rank 2"));
979 
980  for (unsigned int i = 0; i < output.size(); ++i)
981  {
983  apply_transformation(data.contravariant[i],
984  transpose(input[i]));
985  output[i] =
986  apply_transformation(data.covariant[i], A.transpose());
987  }
988 
989  return;
990  }
991 
993  {
994  Assert(
995  data.update_each & update_covariant_transformation,
997  "update_covariant_transformation"));
998  Assert(rank == 2, ExcMessage("Only for rank 2"));
999 
1000  for (unsigned int i = 0; i < output.size(); ++i)
1001  {
1003  apply_transformation(data.covariant[i],
1004  transpose(input[i]));
1005  output[i] =
1006  apply_transformation(data.covariant[i], A.transpose());
1007  }
1008 
1009  return;
1010  }
1011 
1013  {
1014  Assert(
1015  data.update_each & update_covariant_transformation,
1017  "update_covariant_transformation"));
1018  Assert(
1019  data.update_each & update_contravariant_transformation,
1021  "update_contravariant_transformation"));
1022  Assert(
1023  data.update_each & update_volume_elements,
1025  "update_volume_elements"));
1026  Assert(rank == 2, ExcMessage("Only for rank 2"));
1027 
1028  for (unsigned int i = 0; i < output.size(); ++i)
1029  {
1031  apply_transformation(data.covariant[i], input[i]);
1033  apply_transformation(data.contravariant[i],
1034  A.transpose());
1035 
1036  output[i] = transpose(T);
1037  output[i] /= data.volume_elements[i];
1038  }
1039 
1040  return;
1041  }
1042 
1043  default:
1044  Assert(false, ExcNotImplemented());
1045  }
1046  }
1047 
1048 
1049 
1050  template <int dim, int spacedim>
1051  void
1052  transform_hessians(
1053  const ArrayView<const Tensor<3, dim>> & input,
1054  const MappingKind mapping_kind,
1055  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1056  const ArrayView<Tensor<3, spacedim>> & output)
1057  {
1058  AssertDimension(input.size(), output.size());
1059  Assert((dynamic_cast<const typename ::
1060  MappingManifold<dim, spacedim>::InternalData *>(
1061  &mapping_data) != nullptr),
1062  ExcInternalError());
1063  const typename ::MappingManifold<dim, spacedim>::InternalData
1064  &data =
1065  static_cast<const typename ::MappingManifold<dim, spacedim>::
1066  InternalData &>(mapping_data);
1067 
1068  switch (mapping_kind)
1069  {
1071  {
1072  Assert(
1073  data.update_each & update_covariant_transformation,
1075  "update_covariant_transformation"));
1076  Assert(
1077  data.update_each & update_contravariant_transformation,
1079  "update_contravariant_transformation"));
1080 
1081  for (unsigned int q = 0; q < output.size(); ++q)
1082  for (unsigned int i = 0; i < spacedim; ++i)
1083  {
1084  double tmp1[dim][dim];
1085  for (unsigned int J = 0; J < dim; ++J)
1086  for (unsigned int K = 0; K < dim; ++K)
1087  {
1088  tmp1[J][K] =
1089  data.contravariant[q][i][0] * input[q][0][J][K];
1090  for (unsigned int I = 1; I < dim; ++I)
1091  tmp1[J][K] +=
1092  data.contravariant[q][i][I] * input[q][I][J][K];
1093  }
1094  for (unsigned int j = 0; j < spacedim; ++j)
1095  {
1096  double tmp2[dim];
1097  for (unsigned int K = 0; K < dim; ++K)
1098  {
1099  tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1100  for (unsigned int J = 1; J < dim; ++J)
1101  tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1102  }
1103  for (unsigned int k = 0; k < spacedim; ++k)
1104  {
1105  output[q][i][j][k] =
1106  data.covariant[q][k][0] * tmp2[0];
1107  for (unsigned int K = 1; K < dim; ++K)
1108  output[q][i][j][k] +=
1109  data.covariant[q][k][K] * tmp2[K];
1110  }
1111  }
1112  }
1113  return;
1114  }
1115 
1117  {
1118  Assert(
1119  data.update_each & update_covariant_transformation,
1121  "update_covariant_transformation"));
1122 
1123  for (unsigned int q = 0; q < output.size(); ++q)
1124  for (unsigned int i = 0; i < spacedim; ++i)
1125  {
1126  double tmp1[dim][dim];
1127  for (unsigned int J = 0; J < dim; ++J)
1128  for (unsigned int K = 0; K < dim; ++K)
1129  {
1130  tmp1[J][K] =
1131  data.covariant[q][i][0] * input[q][0][J][K];
1132  for (unsigned int I = 1; I < dim; ++I)
1133  tmp1[J][K] +=
1134  data.covariant[q][i][I] * input[q][I][J][K];
1135  }
1136  for (unsigned int j = 0; j < spacedim; ++j)
1137  {
1138  double tmp2[dim];
1139  for (unsigned int K = 0; K < dim; ++K)
1140  {
1141  tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1142  for (unsigned int J = 1; J < dim; ++J)
1143  tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1144  }
1145  for (unsigned int k = 0; k < spacedim; ++k)
1146  {
1147  output[q][i][j][k] =
1148  data.covariant[q][k][0] * tmp2[0];
1149  for (unsigned int K = 1; K < dim; ++K)
1150  output[q][i][j][k] +=
1151  data.covariant[q][k][K] * tmp2[K];
1152  }
1153  }
1154  }
1155 
1156  return;
1157  }
1158 
1159  case mapping_piola_hessian:
1160  {
1161  Assert(
1162  data.update_each & update_covariant_transformation,
1164  "update_covariant_transformation"));
1165  Assert(
1166  data.update_each & update_contravariant_transformation,
1168  "update_contravariant_transformation"));
1169  Assert(
1170  data.update_each & update_volume_elements,
1172  "update_volume_elements"));
1173 
1174  for (unsigned int q = 0; q < output.size(); ++q)
1175  for (unsigned int i = 0; i < spacedim; ++i)
1176  {
1177  double factor[dim];
1178  for (unsigned int I = 0; I < dim; ++I)
1179  factor[I] =
1180  data.contravariant[q][i][I] / data.volume_elements[q];
1181  double tmp1[dim][dim];
1182  for (unsigned int J = 0; J < dim; ++J)
1183  for (unsigned int K = 0; K < dim; ++K)
1184  {
1185  tmp1[J][K] = factor[0] * input[q][0][J][K];
1186  for (unsigned int I = 1; I < dim; ++I)
1187  tmp1[J][K] += factor[I] * input[q][I][J][K];
1188  }
1189  for (unsigned int j = 0; j < spacedim; ++j)
1190  {
1191  double tmp2[dim];
1192  for (unsigned int K = 0; K < dim; ++K)
1193  {
1194  tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
1195  for (unsigned int J = 1; J < dim; ++J)
1196  tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
1197  }
1198  for (unsigned int k = 0; k < spacedim; ++k)
1199  {
1200  output[q][i][j][k] =
1201  data.covariant[q][k][0] * tmp2[0];
1202  for (unsigned int K = 1; K < dim; ++K)
1203  output[q][i][j][k] +=
1204  data.covariant[q][k][K] * tmp2[K];
1205  }
1206  }
1207  }
1208 
1209  return;
1210  }
1211 
1212  default:
1213  Assert(false, ExcNotImplemented());
1214  }
1215  }
1216 
1217 
1218 
1219  template <int dim, int spacedim, int rank>
1220  void
1221  transform_differential_forms(
1222  const ArrayView<const DerivativeForm<rank, dim, spacedim>> &input,
1223  const MappingKind mapping_kind,
1224  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1225  const ArrayView<Tensor<rank + 1, spacedim>> & output)
1226  {
1227  AssertDimension(input.size(), output.size());
1228  Assert((dynamic_cast<const typename ::
1229  MappingManifold<dim, spacedim>::InternalData *>(
1230  &mapping_data) != nullptr),
1231  ExcInternalError());
1232  const typename ::MappingManifold<dim, spacedim>::InternalData
1233  &data =
1234  static_cast<const typename ::MappingManifold<dim, spacedim>::
1235  InternalData &>(mapping_data);
1236 
1237  switch (mapping_kind)
1238  {
1239  case mapping_covariant:
1240  {
1241  Assert(
1242  data.update_each & update_contravariant_transformation,
1244  "update_covariant_transformation"));
1245 
1246  for (unsigned int i = 0; i < output.size(); ++i)
1247  output[i] = apply_transformation(data.covariant[i], input[i]);
1248 
1249  return;
1250  }
1251  default:
1252  Assert(false, ExcNotImplemented());
1253  }
1254  }
1255  } // namespace
1256  } // namespace MappingManifoldImplementation
1257 } // namespace internal
1258 
1259 
1260 
1261 template <int dim, int spacedim>
1262 void
1265  const unsigned int face_no,
1266  const Quadrature<dim - 1> & quadrature,
1267  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1269  &output_data) const
1270 {
1271  // ensure that the following cast is really correct:
1272  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1273  ExcInternalError());
1274  const InternalData &data = static_cast<const InternalData &>(internal_data);
1275 
1276  internal::MappingManifoldImplementation::do_fill_fe_face_values(
1277  *this,
1278  cell,
1279  face_no,
1282  cell->face_orientation(face_no),
1283  cell->face_flip(face_no),
1284  cell->face_rotation(face_no),
1285  quadrature.size()),
1286  quadrature,
1287  data,
1288  output_data);
1289 }
1290 
1291 
1292 
1293 template <int dim, int spacedim>
1294 void
1297  const unsigned int face_no,
1298  const unsigned int subface_no,
1299  const Quadrature<dim - 1> & quadrature,
1300  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
1302  &output_data) const
1303 {
1304  // ensure that the following cast is really correct:
1305  Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
1306  ExcInternalError());
1307  const InternalData &data = static_cast<const InternalData &>(internal_data);
1308 
1309  internal::MappingManifoldImplementation::do_fill_fe_face_values(
1310  *this,
1311  cell,
1312  face_no,
1313  subface_no,
1315  subface_no,
1316  cell->face_orientation(face_no),
1317  cell->face_flip(face_no),
1318  cell->face_rotation(face_no),
1319  quadrature.size(),
1320  cell->subface_case(face_no)),
1321  quadrature,
1322  data,
1323  output_data);
1324 }
1325 
1326 
1327 
1328 template <int dim, int spacedim>
1329 void
1331  const ArrayView<const Tensor<1, dim>> & input,
1332  const MappingKind mapping_kind,
1333  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1334  const ArrayView<Tensor<1, spacedim>> & output) const
1335 {
1336  internal::MappingManifoldImplementation::transform_fields(input,
1337  mapping_kind,
1338  mapping_data,
1339  output);
1340 }
1341 
1342 
1343 
1344 template <int dim, int spacedim>
1345 void
1347  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
1348  const MappingKind mapping_kind,
1349  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1350  const ArrayView<Tensor<2, spacedim>> & output) const
1351 {
1352  internal::MappingManifoldImplementation::transform_differential_forms(
1353  input, mapping_kind, mapping_data, output);
1354 }
1355 
1356 
1357 
1358 template <int dim, int spacedim>
1359 void
1361  const ArrayView<const Tensor<2, dim>> & input,
1362  const MappingKind mapping_kind,
1363  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1364  const ArrayView<Tensor<2, spacedim>> & output) const
1365 {
1366  switch (mapping_kind)
1367  {
1368  case mapping_contravariant:
1369  internal::MappingManifoldImplementation::transform_fields(input,
1370  mapping_kind,
1371  mapping_data,
1372  output);
1373  return;
1374 
1378  internal::MappingManifoldImplementation::transform_gradients(
1379  input, mapping_kind, mapping_data, output);
1380  return;
1381  default:
1382  Assert(false, ExcNotImplemented());
1383  }
1384 }
1385 
1386 
1387 
1388 template <int dim, int spacedim>
1389 void
1391  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
1392  const MappingKind mapping_kind,
1393  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1394  const ArrayView<Tensor<3, spacedim>> & output) const
1395 {
1396  AssertDimension(input.size(), output.size());
1397  Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
1398  ExcInternalError());
1399  const InternalData &data = static_cast<const InternalData &>(mapping_data);
1400 
1401  switch (mapping_kind)
1402  {
1404  {
1407  "update_covariant_transformation"));
1408 
1409  for (unsigned int q = 0; q < output.size(); ++q)
1410  for (unsigned int i = 0; i < spacedim; ++i)
1411  for (unsigned int j = 0; j < spacedim; ++j)
1412  {
1413  double tmp[dim];
1414  for (unsigned int K = 0; K < dim; ++K)
1415  {
1416  tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
1417  for (unsigned int J = 1; J < dim; ++J)
1418  tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
1419  }
1420  for (unsigned int k = 0; k < spacedim; ++k)
1421  {
1422  output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
1423  for (unsigned int K = 1; K < dim; ++K)
1424  output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
1425  }
1426  }
1427  return;
1428  }
1429 
1430  default:
1431  Assert(false, ExcNotImplemented());
1432  }
1433 }
1434 
1435 
1436 
1437 template <int dim, int spacedim>
1438 void
1440  const ArrayView<const Tensor<3, dim>> & input,
1441  const MappingKind mapping_kind,
1442  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
1443  const ArrayView<Tensor<3, spacedim>> & output) const
1444 {
1445  switch (mapping_kind)
1446  {
1447  case mapping_piola_hessian:
1450  internal::MappingManifoldImplementation::transform_hessians(
1451  input, mapping_kind, mapping_data, output);
1452  return;
1453  default:
1454  Assert(false, ExcNotImplemented());
1455  }
1456 }
1457 
1458 //--------------------------- Explicit instantiations -----------------------
1459 #include "mapping_manifold.inst"
1460 
1461 
1462 DEAL_II_NAMESPACE_CLOSE
Transformed quadrature weights.
std::vector< DerivativeForm< 1, dim, spacedim > > covariant
static const unsigned int invalid_unsigned_int
Definition: types.h:187
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
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
SmartPointer< const Manifold< dim, spacedim > > manifold
Contravariant transformation.
std::vector< DerivativeForm< 1, dim, spacedim > > contravariant
const std::vector< double > & get_weights() const
std::array< std::vector< Tensor< 1, dim > >, GeometryInfo< dim >::faces_per_cell *(dim - 1)> unit_tangentials
std::vector< Tensor< 1, spacedim > > boundary_forms
Volume element.
static Point< dim, Number > unit_vector(const unsigned int i)
Outer normal vector, not normalized.
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
Determinant of the Jacobian.
Transformed quadrature points.
void store_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
virtual std::size_t memory_consumption() const override
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1416
MappingKind
Definition: mapping.h:62
static DataSetDescriptor cell()
Definition: qprojector.h:344
Triangulation< dim, spacedim >::cell_iterator cell
std::vector< Point< spacedim > > vertices
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
std::vector< std::vector< Tensor< 1, spacedim > > > aux
DEAL_II_CONSTEXPR Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2369
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
UpdateFlags
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Abstract base class for mapping classes.
Definition: mapping.h:302
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
void initialize(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
DerivativeForm< 1, spacedim, dim, Number > transpose() const
Gradient of volume element.
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:383
std::vector< double > volume_elements
unsigned int size() const
std::vector< double > vertex_weights
void compute_manifold_quadrature_weights(const Quadrature< dim > &quadrature)
std::vector< Point< spacedim > > quadrature_points
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
DEAL_II_CONSTEXPR Number determinant(const SymmetricTensor< 2, dim, Number > &)
Normal vectors.
std::vector< DerivativeForm< 1, dim, spacedim > > jacobians
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const override
virtual void fill_fe_subface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, ::internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
static ::ExceptionBase & ExcNotImplemented()
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const override
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
void initialize_face(const UpdateFlags update_flags, const Quadrature< dim > &quadrature, const unsigned int n_original_q_points)
std::vector< DerivativeForm< 1, spacedim, dim > > inverse_jacobians
static double subface_ratio(const internal::SubfaceCase< dim > &subface_case, const unsigned int subface_no)
Covariant transformation.
std::vector< std::vector< double > > cell_manifold_quadrature_weights
std::vector< Tensor< 1, spacedim > > normal_vectors
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
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
MappingManifold()=default
UpdateFlags update_each
Definition: mapping.h:629
static ::ExceptionBase & ExcInternalError()