Reference documentation for deal.II version 9.3.0
\(\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_q_eulerian.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
17 #include <deal.II/base/utilities.h>
18 
21 
22 #include <deal.II/fe/fe.h>
23 #include <deal.II/fe/fe_tools.h>
26 
28 
32 #include <deal.II/lac/la_vector.h>
37 #include <deal.II/lac/vector.h>
38 
39 #include <memory>
40 
42 
43 
44 
45 // .... MAPPING Q EULERIAN CONSTRUCTOR
46 
47 template <int dim, class VectorType, int spacedim>
50  const unsigned int degree,
51  const MappingQEulerian<dim, VectorType, spacedim> &mapping_q_eulerian)
52  : MappingQGeneric<dim, spacedim>(degree)
53  , mapping_q_eulerian(mapping_q_eulerian)
54  , support_quadrature(degree)
55  , fe_values(mapping_q_eulerian.euler_dof_handler->get_fe(),
56  support_quadrature,
58 {}
59 
60 
61 
62 template <int dim, class VectorType, int spacedim>
64  const unsigned int degree,
66  const VectorType & euler_vector,
67  const unsigned int level)
68  : MappingQ<dim, spacedim>(degree, true)
69  , euler_vector(&euler_vector)
70  , euler_dof_handler(&euler_dof_handler)
71  , level(level)
72 {
73  // reset the q1 mapping we use for interior cells (and previously
74  // set by the MappingQ constructor) to a MappingQ1Eulerian with the
75  // current vector
76  this->q1_mapping =
77  std::make_shared<MappingQ1Eulerian<dim, VectorType, spacedim>>(
79 
80  // also reset the qp mapping pointer with our own class
81  this->qp_mapping = std::make_shared<MappingQEulerianGeneric>(degree, *this);
82 }
83 
84 
85 
86 template <int dim, class VectorType, int spacedim>
87 std::unique_ptr<Mapping<dim, spacedim>>
89 {
90  return std::make_unique<MappingQEulerian<dim, VectorType, spacedim>>(
91  this->get_degree(), *euler_dof_handler, *euler_vector, this->level);
92 }
93 
94 
95 
96 // .... SUPPORT QUADRATURE CONSTRUCTOR
97 
98 template <int dim, class VectorType, int spacedim>
100  SupportQuadrature::SupportQuadrature(const unsigned int map_degree)
101  : Quadrature<dim>(Utilities::fixed_power<dim>(map_degree + 1))
102 {
103  // first we determine the support points on the unit cell in lexicographic
104  // order, which are (in accordance with MappingQ) the support points of
105  // QGaussLobatto.
106  const QGaussLobatto<dim> q_iterated(map_degree + 1);
107  const unsigned int n_q_points = q_iterated.size();
108 
109  // we then need to define a renumbering vector that allows us to go from a
110  // lexicographic numbering scheme to a hierarchic one.
111  const std::vector<unsigned int> renumber =
112  FETools::lexicographic_to_hierarchic_numbering<dim>(map_degree);
113 
114  // finally we assign the quadrature points in the required order.
115  for (unsigned int q = 0; q < n_q_points; ++q)
116  this->quadrature_points[renumber[q]] = q_iterated.point(q);
117 }
118 
119 
120 
121 // .... COMPUTE MAPPING SUPPORT POINTS
122 
123 template <int dim, class VectorType, int spacedim>
124 boost::container::small_vector<Point<spacedim>,
127  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
128 {
129  // get the vertices as the first 2^dim mapping support points
130  const std::vector<Point<spacedim>> a =
131  dynamic_cast<const MappingQEulerianGeneric &>(*this->qp_mapping)
133 
134  boost::container::small_vector<Point<spacedim>,
136  vertex_locations(a.begin(),
138 
139  return vertex_locations;
140 }
141 
142 
143 
144 template <int dim, class VectorType, int spacedim>
145 boost::container::small_vector<Point<spacedim>,
149  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
150 {
151  return mapping_q_eulerian.get_vertices(cell);
152 }
153 
154 
155 
156 template <int dim, class VectorType, int spacedim>
157 bool
160 {
161  return false;
162 }
163 
164 
165 
166 template <int dim, class VectorType, int spacedim>
167 std::vector<Point<spacedim>>
170  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
171 {
172  const bool mg_vector =
174 
175  const types::global_dof_index n_dofs =
176  mg_vector ?
177  mapping_q_eulerian.euler_dof_handler->n_dofs(mapping_q_eulerian.level) :
178  mapping_q_eulerian.euler_dof_handler->n_dofs();
179  const types::global_dof_index vector_size =
180  mapping_q_eulerian.euler_vector->size();
181 
182  (void)n_dofs;
183  (void)vector_size;
184 
185  AssertDimension(vector_size, n_dofs);
186 
187  // we then transform our tria iterator into a dof iterator so we can access
188  // data not associated with triangulations
190  *cell, mapping_q_eulerian.euler_dof_handler);
191 
192  Assert(mg_vector || dof_cell->is_active() == true, ExcInactiveCell());
193 
194  // our quadrature rule is chosen so that each quadrature point corresponds
195  // to a support point in the undeformed configuration. We can then query
196  // the given displacement field at these points to determine the shift
197  // vector that maps the support points to the deformed configuration.
198 
199  // We assume that the given field contains dim displacement components, but
200  // that there may be other solution components as well (e.g. pressures).
201  // this class therefore assumes that the first dim components represent the
202  // actual shift vector we need, and simply ignores any components after
203  // that. This implies that the user should order components appropriately,
204  // or create a separate dof handler for the displacements.
205  const unsigned int n_support_pts = support_quadrature.size();
206  const unsigned int n_components =
207  mapping_q_eulerian.euler_dof_handler->get_fe(0).n_components();
208 
209  Assert(n_components >= spacedim,
210  ExcDimensionMismatch(n_components, spacedim));
211 
212  std::vector<Vector<typename VectorType::value_type>> shift_vector(
213  n_support_pts, Vector<typename VectorType::value_type>(n_components));
214 
215  std::vector<types::global_dof_index> dof_indices(
216  mapping_q_eulerian.euler_dof_handler->get_fe(0).n_dofs_per_cell());
217  // fill shift vector for each support point using an fe_values object. make
218  // sure that the fe_values variable isn't used simultaneously from different
219  // threads
220  std::lock_guard<std::mutex> lock(fe_values_mutex);
221  fe_values.reinit(dof_cell);
222  if (mg_vector)
223  {
224  dof_cell->get_mg_dof_indices(dof_indices);
225  fe_values.get_function_values(*mapping_q_eulerian.euler_vector,
226  dof_indices,
227  shift_vector);
228  }
229  else
230  fe_values.get_function_values(*mapping_q_eulerian.euler_vector,
231  shift_vector);
232 
233  // and finally compute the positions of the support points in the deformed
234  // configuration.
235  std::vector<Point<spacedim>> a(n_support_pts);
236  for (unsigned int q = 0; q < n_support_pts; ++q)
237  {
238  a[q] = fe_values.quadrature_point(q);
239  for (unsigned int d = 0; d < spacedim; ++d)
240  a[q](d) += shift_vector[q](d);
241  }
242 
243  return a;
244 }
245 
246 
247 
248 template <int dim, class VectorType, int spacedim>
251  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
253  const Quadrature<dim> & quadrature,
254  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
256  &output_data) const
257 {
258  // call the function of the base class, but ignoring
259  // any potentially detected cell similarity between
260  // the current and the previous cell
263  quadrature,
264  internal_data,
265  output_data);
266  // also return the updated flag since any detected
267  // similarity wasn't based on the mapped field, but
268  // the original vertices which are meaningless
270 }
271 
272 
273 
274 template <int dim, class VectorType, int spacedim>
277  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
278 {
279  return BoundingBox<spacedim>(
280  dynamic_cast<const MappingQEulerianGeneric &>(*this->qp_mapping)
282 }
283 
284 
285 
286 // explicit instantiations
287 #include "mapping_q_eulerian.inst"
288 
289 
unsigned int get_degree() const
Definition: mapping_q.cc:122
Shape function values.
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
SmartPointer< const DoFHandler< dim, spacedim >, MappingQEulerian< dim, VectorType, spacedim > > euler_dof_handler
std::shared_ptr< const MappingQGeneric< dim, spacedim > > qp_mapping
Definition: mapping_q.h:377
virtual boost::container::small_vector< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
const unsigned int level
virtual boost::container::small_vector< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
static ::ExceptionBase & ExcInactiveCell()
const Point< dim > & point(const unsigned int i) const
Transformed quadrature points.
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
T fixed_power(const T t)
Definition: utilities.h:1081
MappingQEulerianGeneric(const unsigned int degree, const MappingQEulerian< dim, VectorType, spacedim > &mapping_q_eulerian)
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
virtual std::vector< Point< spacedim > > compute_mapping_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
MappingQEulerian(const unsigned int degree, const DoFHandler< dim, spacedim > &euler_dof_handler, const VectorType &euler_vector, const unsigned int level=numbers::invalid_unsigned_int)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:283
unsigned int size() const
virtual BoundingBox< spacedim > get_bounding_box(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
const MappingQEulerian< dim, VectorType, spacedim > & mapping_q_eulerian
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:466
virtual bool preserves_vertex_locations() const override
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
Definition: mapping_q.cc:252
SmartPointer< const VectorType, MappingQEulerian< dim, VectorType, spacedim > > euler_vector
std::shared_ptr< const MappingQGeneric< dim, spacedim > > q1_mapping
Definition: mapping_q.h:361