1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2020 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/utilities.h>
18
19 #include <deal.II/dofs/dof_accessor.h>
20 #include <deal.II/dofs/dof_handler.h>
21
22 #include <deal.II/fe/fe.h>
23 #include <deal.II/fe/fe_tools.h>
24 #include <deal.II/fe/mapping_q1_eulerian.h>
25 #include <deal.II/fe/mapping_q_eulerian.h>
26
27 #include <deal.II/grid/tria_iterator.h>
28
29 #include <deal.II/lac/block_vector.h>
30 #include <deal.II/lac/la_parallel_block_vector.h>
31 #include <deal.II/lac/la_parallel_vector.h>
32 #include <deal.II/lac/la_vector.h>
33 #include <deal.II/lac/petsc_block_vector.h>
34 #include <deal.II/lac/petsc_vector.h>
35 #include <deal.II/lac/trilinos_parallel_block_vector.h>
36 #include <deal.II/lac/trilinos_vector.h>
37 #include <deal.II/lac/vector.h>
38
39 #include <memory>
40
41 DEAL_II_NAMESPACE_OPEN
42
43
44
45 // .... MAPPING Q EULERIAN CONSTRUCTOR
46
47 template <int dim, class VectorType, int spacedim>
48 MappingQEulerian<dim, VectorType, spacedim>::MappingQEulerianGeneric::
MappingQEulerianGeneric(const unsigned int degree,const MappingQEulerian<dim,VectorType,spacedim> & mapping_q_eulerian)49 MappingQEulerianGeneric(
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,
57 update_values | update_quadrature_points)
58 {}
59
60
61
62 template <int dim, class VectorType, int spacedim>
MappingQEulerian(const unsigned int degree,const DoFHandler<dim,spacedim> & euler_dof_handler,const VectorType & euler_vector,const unsigned int level)63 MappingQEulerian<dim, VectorType, spacedim>::MappingQEulerian(
64 const unsigned int degree,
65 const DoFHandler<dim, spacedim> &euler_dof_handler,
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>>(
78 euler_dof_handler, euler_vector);
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>>
clone() const88 MappingQEulerian<dim, VectorType, spacedim>::clone() const
89 {
90 return std::make_unique<MappingQEulerian<dim, VectorType, spacedim>>(
91 this->get_degree(), *euler_dof_handler, *euler_vector);
92 }
93
94
95
96 // .... SUPPORT QUADRATURE CONSTRUCTOR
97
98 template <int dim, class VectorType, int spacedim>
99 MappingQEulerian<dim, VectorType, spacedim>::MappingQEulerianGeneric::
SupportQuadrature(const unsigned int map_degree)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 std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
get_vertices(const typename Triangulation<dim,spacedim>::cell_iterator & cell) const125 MappingQEulerian<dim, VectorType, spacedim>::get_vertices(
126 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
127 {
128 // get the vertices as the first 2^dim mapping support points
129 const std::vector<Point<spacedim>> a =
130 dynamic_cast<const MappingQEulerianGeneric &>(*this->qp_mapping)
131 .compute_mapping_support_points(cell);
132
133 std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
134 vertex_locations;
135 std::copy(a.begin(),
136 a.begin() + GeometryInfo<dim>::vertices_per_cell,
137 vertex_locations.begin());
138
139 return vertex_locations;
140 }
141
142
143
144 template <int dim, class VectorType, int spacedim>
145 std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
146 MappingQEulerian<dim, VectorType, spacedim>::MappingQEulerianGeneric::
get_vertices(const typename Triangulation<dim,spacedim>::cell_iterator & cell) const147 get_vertices(
148 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
149 {
150 return mapping_q_eulerian.get_vertices(cell);
151 }
152
153
154
155 template <int dim, class VectorType, int spacedim>
156 bool
157 MappingQEulerian<dim, VectorType, spacedim>::MappingQEulerianGeneric::
preserves_vertex_locations() const158 preserves_vertex_locations() const
159 {
160 return false;
161 }
162
163
164
165 template <int dim, class VectorType, int spacedim>
166 std::vector<Point<spacedim>>
167 MappingQEulerian<dim, VectorType, spacedim>::MappingQEulerianGeneric::
compute_mapping_support_points(const typename Triangulation<dim,spacedim>::cell_iterator & cell) const168 compute_mapping_support_points(
169 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
170 {
171 const bool mg_vector =
172 mapping_q_eulerian.level != numbers::invalid_unsigned_int;
173
174 const types::global_dof_index n_dofs =
175 mg_vector ?
176 mapping_q_eulerian.euler_dof_handler->n_dofs(mapping_q_eulerian.level) :
177 mapping_q_eulerian.euler_dof_handler->n_dofs();
178 const types::global_dof_index vector_size =
179 mapping_q_eulerian.euler_vector->size();
180
181 (void)n_dofs;
182 (void)vector_size;
183
184 AssertDimension(vector_size, n_dofs);
185
186 // we then transform our tria iterator into a dof iterator so we can access
187 // data not associated with triangulations
188 typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(
189 *cell, mapping_q_eulerian.euler_dof_handler);
190
191 Assert(mg_vector || dof_cell->is_active() == true, ExcInactiveCell());
192
193 // our quadrature rule is chosen so that each quadrature point corresponds
194 // to a support point in the undeformed configuration. We can then query
195 // the given displacement field at these points to determine the shift
196 // vector that maps the support points to the deformed configuration.
197
198 // We assume that the given field contains dim displacement components, but
199 // that there may be other solution components as well (e.g. pressures).
200 // this class therefore assumes that the first dim components represent the
201 // actual shift vector we need, and simply ignores any components after
202 // that. This implies that the user should order components appropriately,
203 // or create a separate dof handler for the displacements.
204 const unsigned int n_support_pts = support_quadrature.size();
205 const unsigned int n_components =
206 mapping_q_eulerian.euler_dof_handler->get_fe(0).n_components();
207
208 Assert(n_components >= spacedim,
209 ExcDimensionMismatch(n_components, spacedim));
210
211 std::vector<Vector<typename VectorType::value_type>> shift_vector(
212 n_support_pts, Vector<typename VectorType::value_type>(n_components));
213
214 std::vector<types::global_dof_index> dof_indices(
215 mapping_q_eulerian.euler_dof_handler->get_fe(0).n_dofs_per_cell());
216 // fill shift vector for each support point using an fe_values object. make
217 // sure that the fe_values variable isn't used simultaneously from different
218 // threads
219 std::lock_guard<std::mutex> lock(fe_values_mutex);
220 fe_values.reinit(dof_cell);
221 if (mg_vector)
222 {
223 dof_cell->get_mg_dof_indices(dof_indices);
224 fe_values.get_function_values(*mapping_q_eulerian.euler_vector,
225 dof_indices,
226 shift_vector);
227 }
228 else
229 fe_values.get_function_values(*mapping_q_eulerian.euler_vector,
230 shift_vector);
231
232 // and finally compute the positions of the support points in the deformed
233 // configuration.
234 std::vector<Point<spacedim>> a(n_support_pts);
235 for (unsigned int q = 0; q < n_support_pts; ++q)
236 {
237 a[q] = fe_values.quadrature_point(q);
238 for (unsigned int d = 0; d < spacedim; ++d)
239 a[q](d) += shift_vector[q](d);
240 }
241
242 return a;
243 }
244
245
246
247 template <int dim, class VectorType, int spacedim>
248 CellSimilarity::Similarity
fill_fe_values(const typename Triangulation<dim,spacedim>::cell_iterator & cell,const CellSimilarity::Similarity,const Quadrature<dim> & quadrature,const typename Mapping<dim,spacedim>::InternalDataBase & internal_data,internal::FEValuesImplementation::MappingRelatedData<dim,spacedim> & output_data) const249 MappingQEulerian<dim, VectorType, spacedim>::fill_fe_values(
250 const typename Triangulation<dim, spacedim>::cell_iterator &cell,
251 const CellSimilarity::Similarity,
252 const Quadrature<dim> & quadrature,
253 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
254 internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
255 &output_data) const
256 {
257 // call the function of the base class, but ignoring
258 // any potentially detected cell similarity between
259 // the current and the previous cell
260 MappingQ<dim, spacedim>::fill_fe_values(cell,
261 CellSimilarity::invalid_next_cell,
262 quadrature,
263 internal_data,
264 output_data);
265 // also return the updated flag since any detected
266 // similarity wasn't based on the mapped field, but
267 // the original vertices which are meaningless
268 return CellSimilarity::invalid_next_cell;
269 }
270
271
272
273 // explicit instantiations
274 #include "mapping_q_eulerian.inst"
275
276
277 DEAL_II_NAMESPACE_CLOSE
278