Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-2846-g6fb608615f 2025-03-15 04:10:00+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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
mapping_q_eulerian.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2008 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
17
20
21#include <deal.II/fe/fe.h>
22#include <deal.II/fe/fe_tools.h>
24
26
34#include <deal.II/lac/vector.h>
35
36#include <boost/container/small_vector.hpp>
37
38#include <memory>
39
40
42
43
44// .... MAPPING Q EULERIAN CONSTRUCTOR
45
46
47template <int dim, typename 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
65template <int dim, typename VectorType, int spacedim>
66std::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
77template <int dim, typename 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
102template <int dim, typename VectorType, int spacedim>
103boost::container::small_vector<Point<spacedim>,
104#ifndef _MSC_VER
105 ReferenceCells::max_n_vertices<dim>()
106#else
108#endif
109 >
111 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
112{
113 // get the vertices as the first 2^dim mapping support points
114 const std::vector<Point<spacedim>> a = compute_mapping_support_points(cell);
115
116 boost::container::small_vector<Point<spacedim>,
117#ifndef _MSC_VER
118 ReferenceCells::max_n_vertices<dim>()
119#else
121#endif
122 >
123 vertex_locations(a.begin(), a.begin() + cell->n_vertices());
124
125 return vertex_locations;
126}
127
128
129
130template <int dim, typename VectorType, int spacedim>
131std::vector<Point<spacedim>>
133 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
134{
135 const bool mg_vector = level != numbers::invalid_unsigned_int;
136
137 [[maybe_unused]] const types::global_dof_index n_dofs =
138 mg_vector ? euler_dof_handler->n_dofs(level) : euler_dof_handler->n_dofs();
139 [[maybe_unused]] const types::global_dof_index vector_size =
140 euler_vector->size();
141
142 AssertDimension(vector_size, n_dofs);
143
144 // we then transform our tria iterator into a dof iterator so we can access
145 // data not associated with triangulations
146 typename DoFHandler<dim, spacedim>::cell_iterator dof_cell(*cell,
148
149 Assert(mg_vector || dof_cell->is_active() == true, ExcInactiveCell());
150
151 // our quadrature rule is chosen so that each quadrature point corresponds
152 // to a support point in the undeformed configuration. We can then query
153 // the given displacement field at these points to determine the shift
154 // vector that maps the support points to the deformed configuration.
155
156 // We assume that the given field contains dim displacement components, but
157 // that there may be other solution components as well (e.g. pressures).
158 // this class therefore assumes that the first dim components represent the
159 // actual shift vector we need, and simply ignores any components after
160 // that. This implies that the user should order components appropriately,
161 // or create a separate dof handler for the displacements.
162 const unsigned int n_support_pts = support_quadrature.size();
163 const unsigned int n_components = euler_dof_handler->get_fe(0).n_components();
164
165 Assert(n_components >= spacedim,
166 ExcDimensionMismatch(n_components, spacedim));
167
168 std::vector<Vector<typename VectorType::value_type>> shift_vector(
169 n_support_pts, Vector<typename VectorType::value_type>(n_components));
170
171 std::vector<types::global_dof_index> dof_indices(
172 euler_dof_handler->get_fe(0).n_dofs_per_cell());
173 // fill shift vector for each support point using an fe_values object. make
174 // sure that the fe_values variable isn't used simultaneously from different
175 // threads
176 std::lock_guard<std::mutex> lock(fe_values_mutex);
177 fe_values.reinit(dof_cell);
178 if (mg_vector)
179 {
180 dof_cell->get_mg_dof_indices(dof_indices);
181 fe_values.get_function_values(*euler_vector, dof_indices, shift_vector);
182 }
183 else
184 fe_values.get_function_values(*euler_vector, shift_vector);
185
186 // and finally compute the positions of the support points in the deformed
187 // configuration.
188 std::vector<Point<spacedim>> a(n_support_pts);
189 for (unsigned int q = 0; q < n_support_pts; ++q)
190 {
191 a[q] = fe_values.quadrature_point(q);
192 for (unsigned int d = 0; d < spacedim; ++d)
193 a[q][d] += shift_vector[q][d];
194 }
195
196 return a;
197}
198
199
200
201template <int dim, typename VectorType, int spacedim>
206 const Quadrature<dim> &quadrature,
207 const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
209 &output_data) const
210{
211 // call the function of the base class, but ignoring
212 // any potentially detected cell similarity between
213 // the current and the previous cell
216 quadrature,
217 internal_data,
218 output_data);
219 // also return the updated flag since any detected
220 // similarity wasn't based on the mapped field, but
221 // the original vertices which are meaningless
223}
224
225
226// explicit instantiations
227#include "fe/mapping_q_eulerian.inst"
228
229
SupportQuadrature(const unsigned int map_degree)
FEValues< dim, spacedim > fe_values
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
MappingQEulerian(const unsigned int degree, const DoFHandler< dim, spacedim > &euler_dof_handler, const VectorType &euler_vector, const unsigned int level=numbers::invalid_unsigned_int)
ObserverPointer< const VectorType, MappingQEulerian< dim, VectorType, spacedim > > euler_vector
ObserverPointer< const DoFHandler< dim, spacedim >, MappingQEulerian< dim, VectorType, spacedim > > euler_dof_handler
Threads::Mutex fe_values_mutex
const SupportQuadrature support_quadrature
virtual boost::container::small_vector< Point< spacedim >, ReferenceCells::max_n_vertices< dim >() > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
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:833
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
unsigned int level
Definition grid_out.cc:4635
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInactiveCell()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
typename ActiveSelector::cell_iterator cell_iterator
@ update_values
Shape function values.
@ update_quadrature_points
Transformed quadrature points.
constexpr unsigned int invalid_unsigned_int
Definition types.h:232