Reference documentation for deal.II version GIT f0f8c7fe18 2023-03-21 21:25:02+00:00
\(\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>
25 
27 
31 #include <deal.II/lac/la_vector.h>
36 #include <deal.II/lac/vector.h>
37 
38 #include <memory>
39 
41 
42 
43 
44 // .... MAPPING Q EULERIAN CONSTRUCTOR
45 
46 
47 template <int dim, class VectorType, int spacedim>
49  const unsigned int degree,
50  const DoFHandler<dim, spacedim> &euler_dof_handler,
51  const VectorType & euler_vector,
52  const unsigned int level)
53  : MappingQ<dim, spacedim>(degree)
54  , euler_vector(&euler_vector)
55  , euler_dof_handler(&euler_dof_handler)
56  , level(level)
57  , support_quadrature(degree)
58  , fe_values(euler_dof_handler.get_fe(),
59  support_quadrature,
61 {}
62 
63 
64 
65 template <int dim, class VectorType, int spacedim>
66 std::unique_ptr<Mapping<dim, spacedim>>
68 {
69  return std::make_unique<MappingQEulerian<dim, VectorType, spacedim>>(
70  this->get_degree(), *euler_dof_handler, *euler_vector, this->level);
71 }
72 
73 
74 
75 // .... SUPPORT QUADRATURE CONSTRUCTOR
76 
77 template <int dim, class VectorType, int spacedim>
79  SupportQuadrature(const unsigned int map_degree)
80  : Quadrature<dim>(Utilities::fixed_power<dim>(map_degree + 1))
81 {
82  // first we determine the support points on the unit cell in lexicographic
83  // order, which are (in accordance with MappingQ) the support points of
84  // QGaussLobatto.
85  const QGaussLobatto<dim> q_iterated(map_degree + 1);
86  const unsigned int n_q_points = q_iterated.size();
87 
88  // we then need to define a renumbering vector that allows us to go from a
89  // lexicographic numbering scheme to a hierarchic one.
90  const std::vector<unsigned int> renumber =
91  FETools::lexicographic_to_hierarchic_numbering<dim>(map_degree);
92 
93  // finally we assign the quadrature points in the required order.
94  for (unsigned int q = 0; q < n_q_points; ++q)
95  this->quadrature_points[renumber[q]] = q_iterated.point(q);
96 }
97 
98 
99 
100 // .... COMPUTE MAPPING SUPPORT POINTS
101 
102 template <int dim, class VectorType, int spacedim>
103 boost::container::small_vector<Point<spacedim>,
106  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
107 {
108  // get the vertices as the first 2^dim mapping support points
109  const std::vector<Point<spacedim>> a = compute_mapping_support_points(cell);
110 
111  boost::container::small_vector<Point<spacedim>,
113  vertex_locations(a.begin(),
115 
116  return vertex_locations;
117 }
118 
119 
120 
121 template <int dim, class VectorType, int spacedim>
122 std::vector<Point<spacedim>>
124  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
125 {
126  const bool mg_vector = level != numbers::invalid_unsigned_int;
127 
128  const types::global_dof_index n_dofs =
129  mg_vector ? euler_dof_handler->n_dofs(level) : euler_dof_handler->n_dofs();
130  const types::global_dof_index vector_size = euler_vector->size();
131 
132  (void)n_dofs;
133  (void)vector_size;
134 
135  AssertDimension(vector_size, n_dofs);
136 
137  // we then transform our tria iterator into a dof iterator so we can access
138  // data not associated with triangulations
139  typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(*cell,
141 
142  Assert(mg_vector || dof_cell->is_active() == true, ExcInactiveCell());
143 
144  // our quadrature rule is chosen so that each quadrature point corresponds
145  // to a support point in the undeformed configuration. We can then query
146  // the given displacement field at these points to determine the shift
147  // vector that maps the support points to the deformed configuration.
148 
149  // We assume that the given field contains dim displacement components, but
150  // that there may be other solution components as well (e.g. pressures).
151  // this class therefore assumes that the first dim components represent the
152  // actual shift vector we need, and simply ignores any components after
153  // that. This implies that the user should order components appropriately,
154  // or create a separate dof handler for the displacements.
155  const unsigned int n_support_pts = support_quadrature.size();
156  const unsigned int n_components = euler_dof_handler->get_fe(0).n_components();
157 
158  Assert(n_components >= spacedim,
159  ExcDimensionMismatch(n_components, spacedim));
160 
161  std::vector<Vector<typename VectorType::value_type>> shift_vector(
162  n_support_pts, Vector<typename VectorType::value_type>(n_components));
163 
164  std::vector<types::global_dof_index> dof_indices(
165  euler_dof_handler->get_fe(0).n_dofs_per_cell());
166  // fill shift vector for each support point using an fe_values object. make
167  // sure that the fe_values variable isn't used simultaneously from different
168  // threads
169  std::lock_guard<std::mutex> lock(fe_values_mutex);
170  fe_values.reinit(dof_cell);
171  if (mg_vector)
172  {
173  dof_cell->get_mg_dof_indices(dof_indices);
174  fe_values.get_function_values(*euler_vector, dof_indices, shift_vector);
175  }
176  else
177  fe_values.get_function_values(*euler_vector, shift_vector);
178 
179  // and finally compute the positions of the support points in the deformed
180  // configuration.
181  std::vector<Point<spacedim>> a(n_support_pts);
182  for (unsigned int q = 0; q < n_support_pts; ++q)
183  {
184  a[q] = fe_values.quadrature_point(q);
185  for (unsigned int d = 0; d < spacedim; ++d)
186  a[q](d) += shift_vector[q](d);
187  }
188 
189  return a;
190 }
191 
192 
193 
194 template <int dim, class VectorType, int spacedim>
197  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
199  const Quadrature<dim> & quadrature,
200  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
202  &output_data) const
203 {
204  // call the function of the base class, but ignoring
205  // any potentially detected cell similarity between
206  // the current and the previous cell
209  quadrature,
210  internal_data,
211  output_data);
212  // also return the updated flag since any detected
213  // similarity wasn't based on the mapped field, but
214  // the original vertices which are meaningless
216 }
217 
218 
219 // explicit instantiations
220 #include "mapping_q_eulerian.inst"
221 
222 
SupportQuadrature(const unsigned int map_degree)
FEValues< dim, spacedim > fe_values
SmartPointer< const DoFHandler< dim, spacedim >, MappingQEulerian< dim, VectorType, spacedim > > euler_dof_handler
virtual std::vector< Point< spacedim > > compute_mapping_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell) 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
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
virtual boost::container::small_vector< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
MappingQEulerian(const unsigned int degree, const DoFHandler< dim, spacedim > &euler_dof_handler, const VectorType &euler_vector, const unsigned int level=numbers::invalid_unsigned_int)
SmartPointer< const VectorType, MappingQEulerian< dim, VectorType, spacedim > > euler_vector
Threads::Mutex fe_values_mutex
const SupportQuadrature support_quadrature
const unsigned int level
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:944
Abstract base class for mapping classes.
Definition: mapping.h:314
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:326
const Point< dim > & point(const unsigned int i) const
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
unsigned int level
Definition: grid_out.cc:4618
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1586
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
static ::ExceptionBase & ExcInactiveCell()
@ update_values
Shape function values.
@ update_quadrature_points
Transformed quadrature points.
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T fixed_power(const T t)
Definition: utilities.h:966
unsigned int get_degree(const std::vector< typename BarycentricPolynomials< dim >::PolyType > &polys)
static const unsigned int invalid_unsigned_int
Definition: types.h:213
unsigned int global_dof_index
Definition: types.h:82