Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
mapping_q_eulerian.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/quadrature_lib.h>
17 #include <deal.II/base/std_cxx14/memory.h>
18 #include <deal.II/base/utilities.h>
19 
20 #include <deal.II/dofs/dof_accessor.h>
21 #include <deal.II/dofs/dof_handler.h>
22 
23 #include <deal.II/fe/fe.h>
24 #include <deal.II/fe/fe_tools.h>
25 #include <deal.II/fe/mapping_q1_eulerian.h>
26 #include <deal.II/fe/mapping_q_eulerian.h>
27 
28 #include <deal.II/grid/tria_iterator.h>
29 
30 #include <deal.II/lac/block_vector.h>
31 #include <deal.II/lac/la_parallel_block_vector.h>
32 #include <deal.II/lac/la_parallel_vector.h>
33 #include <deal.II/lac/la_vector.h>
34 #include <deal.II/lac/petsc_block_vector.h>
35 #include <deal.II/lac/petsc_vector.h>
36 #include <deal.II/lac/trilinos_parallel_block_vector.h>
37 #include <deal.II/lac/trilinos_vector.h>
38 #include <deal.II/lac/vector.h>
39 
40 DEAL_II_NAMESPACE_OPEN
41 
42 
43 
44 // .... MAPPING Q EULERIAN CONSTRUCTOR
45 
46 template <int dim, class VectorType, int spacedim>
49  const unsigned int degree,
50  const MappingQEulerian<dim, VectorType, spacedim> &mapping_q_eulerian)
51  : MappingQGeneric<dim, spacedim>(degree)
52  , mapping_q_eulerian(mapping_q_eulerian)
53  , support_quadrature(degree)
54  , fe_values(mapping_q_eulerian.euler_dof_handler->get_fe(),
55  support_quadrature,
57 {}
58 
59 
60 
61 template <int dim, class VectorType, int spacedim>
63  const unsigned int degree,
65  const VectorType & euler_vector,
66  const unsigned int level)
67  : MappingQ<dim, spacedim>(degree, true)
70  , level(level)
71 {
72  // reset the q1 mapping we use for interior cells (and previously
73  // set by the MappingQ constructor) to a MappingQ1Eulerian with the
74  // current vector
75  this->q1_mapping =
76  std::make_shared<MappingQ1Eulerian<dim, VectorType, spacedim>>(
78 
79  // also reset the qp mapping pointer with our own class
80  this->qp_mapping = std::make_shared<MappingQEulerianGeneric>(degree, *this);
81 }
82 
83 
84 
85 template <int dim, class VectorType, int spacedim>
86 std::unique_ptr<Mapping<dim, spacedim>>
88 {
89  return std_cxx14::make_unique<MappingQEulerian<dim, VectorType, spacedim>>(
90  this->get_degree(), *euler_dof_handler, *euler_vector);
91 }
92 
93 
94 
95 // .... SUPPORT QUADRATURE CONSTRUCTOR
96 
97 template <int dim, class VectorType, int spacedim>
99  SupportQuadrature::SupportQuadrature(const unsigned int map_degree)
100  : Quadrature<dim>(Utilities::fixed_power<dim>(map_degree + 1))
101 {
102  // first we determine the support points on the unit cell in lexicographic
103  // order, which are (in accordance with MappingQ) the support points of
104  // QGaussLobatto.
105  const QGaussLobatto<dim> q_iterated(map_degree + 1);
106  const unsigned int n_q_points = q_iterated.size();
107 
108  // we then need to define a renumbering vector that allows us to go from a
109  // lexicographic numbering scheme to a hierarchic one. this fragment is
110  // taking almost verbatim from the MappingQ class.
111  std::vector<unsigned int> renumber(n_q_points);
112  std::vector<unsigned int> dpo(dim + 1, 1U);
113  for (unsigned int i = 1; i < dpo.size(); ++i)
114  dpo[i] = dpo[i - 1] * (map_degree - 1);
115 
117  FiniteElementData<dim>(dpo, 1, map_degree), renumber);
118 
119  // finally we assign the quadrature points in the required order.
120  for (unsigned int q = 0; q < n_q_points; ++q)
121  this->quadrature_points[renumber[q]] = q_iterated.point(q);
122 }
123 
124 
125 
126 // .... COMPUTE MAPPING SUPPORT POINTS
127 
128 template <int dim, class VectorType, int spacedim>
129 std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
131  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
132 {
133  // get the vertices as the first 2^dim mapping support points
134  const std::vector<Point<spacedim>> a =
135  dynamic_cast<const MappingQEulerianGeneric &>(*this->qp_mapping)
137 
138  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
139  vertex_locations;
140  std::copy(a.begin(),
142  vertex_locations.begin());
143 
144  return vertex_locations;
145 }
146 
147 
148 
149 template <int dim, class VectorType, int spacedim>
150 std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
153  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
154 {
155  return mapping_q_eulerian.get_vertices(cell);
156 }
157 
158 
159 
160 template <int dim, class VectorType, int spacedim>
161 bool
164 {
165  return false;
166 }
167 
168 
169 
170 template <int dim, class VectorType, int spacedim>
171 std::vector<Point<spacedim>>
174  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
175 {
176  const bool mg_vector =
178 
179  const types::global_dof_index n_dofs =
180  mg_vector ?
181  mapping_q_eulerian.euler_dof_handler->n_dofs(mapping_q_eulerian.level) :
182  mapping_q_eulerian.euler_dof_handler->n_dofs();
183  const types::global_dof_index vector_size =
184  mapping_q_eulerian.euler_vector->size();
185 
186  (void)n_dofs;
187  (void)vector_size;
188 
189  AssertDimension(vector_size, n_dofs);
190 
191  // we then transform our tria iterator into a dof iterator so we can access
192  // data not associated with triangulations
194  *cell, mapping_q_eulerian.euler_dof_handler);
195 
196  Assert(mg_vector || dof_cell->active() == true, ExcInactiveCell());
197 
198  // our quadrature rule is chosen so that each quadrature point corresponds
199  // to a support point in the undeformed configuration. We can then query
200  // the given displacement field at these points to determine the shift
201  // vector that maps the support points to the deformed configuration.
202 
203  // We assume that the given field contains dim displacement components, but
204  // that there may be other solution components as well (e.g. pressures).
205  // this class therefore assumes that the first dim components represent the
206  // actual shift vector we need, and simply ignores any components after
207  // that. This implies that the user should order components appropriately,
208  // or create a separate dof handler for the displacements.
209  const unsigned int n_support_pts = support_quadrature.size();
210  const unsigned int n_components =
211  mapping_q_eulerian.euler_dof_handler->get_fe(0).n_components();
212 
213  Assert(n_components >= spacedim,
214  ExcDimensionMismatch(n_components, spacedim));
215 
216  std::vector<Vector<typename VectorType::value_type>> shift_vector(
217  n_support_pts, Vector<typename VectorType::value_type>(n_components));
218 
219  std::vector<types::global_dof_index> dof_indices(
220  mapping_q_eulerian.euler_dof_handler->get_fe(0).dofs_per_cell);
221  // fill shift vector for each support point using an fe_values object. make
222  // sure that the fe_values variable isn't used simultaneously from different
223  // threads
224  std::lock_guard<std::mutex> lock(fe_values_mutex);
225  fe_values.reinit(dof_cell);
226  if (mg_vector)
227  {
228  dof_cell->get_mg_dof_indices(dof_indices);
229  fe_values.get_function_values(*mapping_q_eulerian.euler_vector,
230  dof_indices,
231  shift_vector);
232  }
233  else
234  fe_values.get_function_values(*mapping_q_eulerian.euler_vector,
235  shift_vector);
236 
237  // and finally compute the positions of the support points in the deformed
238  // configuration.
239  std::vector<Point<spacedim>> a(n_support_pts);
240  for (unsigned int q = 0; q < n_support_pts; ++q)
241  {
242  a[q] = fe_values.quadrature_point(q);
243  for (unsigned int d = 0; d < spacedim; ++d)
244  a[q](d) += shift_vector[q](d);
245  }
246 
247  return a;
248 }
249 
250 
251 
252 template <int dim, class VectorType, int spacedim>
255  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
257  const Quadrature<dim> & quadrature,
258  const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
260  &output_data) const
261 {
262  // call the function of the base class, but ignoring
263  // any potentially detected cell similarity between
264  // the current and the previous cell
267  quadrature,
268  internal_data,
269  output_data);
270  // also return the updated flag since any detected
271  // similarity wasn't based on the mapped field, but
272  // the original vertices which are meaningless
274 }
275 
276 
277 
278 // explicit instantiations
279 #include "mapping_q_eulerian.inst"
280 
281 
282 DEAL_II_NAMESPACE_CLOSE
Shape function values.
static const unsigned int invalid_unsigned_int
Definition: types.h:173
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
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:368
const unsigned int level
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
MappingQEulerianGeneric(const unsigned int degree, const MappingQEulerian< dim, VectorType, spacedim > &mapping_q_eulerian)
virtual std::array< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
#define Assert(cond, exc)
Definition: exceptions.h:1407
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
MappingQEulerian(const unsigned int degree, const DoFHandler< dim, spacedim > &euler_dof_handler, const VectorType &euler_vector, const unsigned int level=numbers::invalid_unsigned_int)
virtual std::array< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:284
unsigned int size() const
void lexicographic_to_hierarchic_numbering(const FiniteElementData< dim > &fe_data, std::vector< unsigned int > &l2h)
unsigned int global_dof_index
Definition: types.h:89
Definition: cuda.h:31
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
const MappingQEulerian< dim, VectorType, spacedim > & mapping_q_eulerian
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:352
virtual bool preserves_vertex_locations() const override
SmartPointer< const VectorType, MappingQEulerian< dim, VectorType, spacedim > > euler_vector
std::shared_ptr< const MappingQGeneric< dim, spacedim > > q1_mapping
Definition: mapping_q.h:352