Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3083-g7b89508ac7 2025-04-18 12:50: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
grid_tools_geometry.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2023 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#include <deal.II/base/mpi.h>
19
21
26
28#include <deal.II/grid/tria.h>
29
31
33
34#include <functional>
35
37
38
39namespace GridTools
40{
41 template <int dim, int spacedim>
42 double
44 {
45 // we can't deal with distributed meshes since we don't have all
46 // vertices locally. there is one exception, however: if the mesh has
47 // never been refined. the way to test this is not to ask
48 // tria.n_levels()==1, since this is something that can happen on one
49 // processor without being true on all. however, we can ask for the
50 // global number of active cells and use that
51 if constexpr (running_in_debug_mode())
52 {
53 if (const auto *p_tria = dynamic_cast<
55 &tria))
56 Assert(p_tria->n_global_active_cells() == tria.n_cells(0),
58 }
59
60 // the algorithm used simply traverses all cells and picks out the
61 // boundary vertices. it may or may not be faster to simply get all
62 // vectors, don't mark boundary vertices, and compute the distances
63 // thereof, but at least as the mesh is refined, it seems better to
64 // first mark boundary nodes, as marking is O(N) in the number of
65 // cells/vertices, while computing the maximal distance is O(N*N)
66 const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
67 std::vector<bool> boundary_vertices(vertices.size(), false);
68
70 tria.begin_active();
72 tria.end();
73 for (; cell != endc; ++cell)
74 for (const unsigned int face : cell->face_indices())
75 if (cell->face(face)->at_boundary())
76 for (unsigned int i = 0; i < cell->face(face)->n_vertices(); ++i)
77 boundary_vertices[cell->face(face)->vertex_index(i)] = true;
78
79 // now traverse the list of boundary vertices and check distances.
80 // since distances are symmetric, we only have to check one half
81 double max_distance_sqr = 0;
82 std::vector<bool>::const_iterator pi = boundary_vertices.begin();
83 const unsigned int N = boundary_vertices.size();
84 for (unsigned int i = 0; i < N; ++i, ++pi)
85 {
86 std::vector<bool>::const_iterator pj = pi + 1;
87 for (unsigned int j = i + 1; j < N; ++j, ++pj)
88 if ((*pi == true) && (*pj == true) &&
89 ((vertices[i] - vertices[j]).norm_square() > max_distance_sqr))
90 max_distance_sqr = (vertices[i] - vertices[j]).norm_square();
91 }
92
93 return std::sqrt(max_distance_sqr);
94 }
95
96
97
98 template <int dim, int spacedim>
99 double
101 {
102 Assert(triangulation.get_reference_cells().size() == 1,
104 const ReferenceCell reference_cell = triangulation.get_reference_cells()[0];
105 return volume(
107 reference_cell.template get_default_linear_mapping<dim, spacedim>());
108 }
109
110
111
112 template <int dim, int spacedim>
113 double
115 const Mapping<dim, spacedim> &mapping)
116 {
117 // get the degree of the mapping if possible. if not, just assume 1
118 unsigned int mapping_degree = 1;
119 if (const auto *p = dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
120 mapping_degree = p->get_degree();
121 else if (const auto *p =
122 dynamic_cast<const MappingFE<dim, spacedim> *>(&mapping))
123 mapping_degree = p->get_degree();
124
125 // then initialize an appropriate quadrature formula
126 Assert(triangulation.get_reference_cells().size() == 1,
128 const ReferenceCell reference_cell = triangulation.get_reference_cells()[0];
129 const Quadrature<dim> quadrature_formula =
130 reference_cell.template get_gauss_type_quadrature<dim>(mapping_degree +
131 1);
132 const unsigned int n_q_points = quadrature_formula.size();
133
134 // we really want the JxW values from the FEValues object, but it
135 // wants a finite element. create a cheap element as a dummy
136 // element
137 FE_Nothing<dim, spacedim> dummy_fe(reference_cell);
138 FEValues<dim, spacedim> fe_values(mapping,
139 dummy_fe,
140 quadrature_formula,
142
143 double local_volume = 0;
144
145 // compute the integral quantities by quadrature
146 for (const auto &cell : triangulation.active_cell_iterators())
147 if (cell->is_locally_owned())
148 {
149 fe_values.reinit(cell);
150 for (unsigned int q = 0; q < n_q_points; ++q)
151 local_volume += fe_values.JxW(q);
152 }
153
154 const double global_volume =
155 Utilities::MPI::sum(local_volume, triangulation.get_mpi_communicator());
156
157 return global_volume;
158 }
159
160
161
162 template <int dim, int spacedim>
163 std::pair<unsigned int, double>
166 {
167 double max_ratio = 1;
168 unsigned int index = 0;
169
170 for (unsigned int i = 0; i < dim; ++i)
171 for (unsigned int j = i + 1; j < dim; ++j)
172 {
173 unsigned int ax = i % dim;
174 unsigned int next_ax = j % dim;
175
176 double ratio =
177 cell->extent_in_direction(ax) / cell->extent_in_direction(next_ax);
178
179 if (ratio > max_ratio)
180 {
181 max_ratio = ratio;
182 index = ax;
183 }
184 else if (1.0 / ratio > max_ratio)
185 {
186 max_ratio = 1.0 / ratio;
187 index = next_ax;
188 }
189 }
190 return std::make_pair(index, max_ratio);
191 }
192
193
194
195 namespace
196 {
211 template <int dim>
212 struct TransformR2UAffine
213 {
214 static const double KA[GeometryInfo<dim>::vertices_per_cell][dim];
216 };
217
218
219 /*
220 Octave code:
221 M=[0 1; 1 1];
222 K1 = transpose(M) * inverse (M*transpose(M));
223 printf ("{%f, %f},\n", K1' );
224 */
225 template <>
226 const double TransformR2UAffine<1>::KA[GeometryInfo<1>::vertices_per_cell]
227 [1] = {{-1.000000}, {1.000000}};
228
229 template <>
230 const double TransformR2UAffine<1>::Kb[GeometryInfo<1>::vertices_per_cell] =
231 {1.000000, 0.000000};
232
233
234 /*
235 Octave code:
236 M=[0 1 0 1;0 0 1 1;1 1 1 1];
237 K2 = transpose(M) * inverse (M*transpose(M));
238 printf ("{%f, %f, %f},\n", K2' );
239 */
240 template <>
241 const double TransformR2UAffine<2>::KA[GeometryInfo<2>::vertices_per_cell]
242 [2] = {{-0.500000, -0.500000},
243 {0.500000, -0.500000},
244 {-0.500000, 0.500000},
245 {0.500000, 0.500000}};
246
247 /*
248 Octave code:
249 M=[0 1 0 1 0 1 0 1;0 0 1 1 0 0 1 1; 0 0 0 0 1 1 1 1; 1 1 1 1 1 1 1 1];
250 K3 = transpose(M) * inverse (M*transpose(M))
251 printf ("{%f, %f, %f, %f},\n", K3' );
252 */
253 template <>
254 const double TransformR2UAffine<2>::Kb[GeometryInfo<2>::vertices_per_cell] =
255 {0.750000, 0.250000, 0.250000, -0.250000};
256
257
258 template <>
259 const double TransformR2UAffine<3>::KA[GeometryInfo<3>::vertices_per_cell]
260 [3] = {
261 {-0.250000, -0.250000, -0.250000},
262 {0.250000, -0.250000, -0.250000},
263 {-0.250000, 0.250000, -0.250000},
264 {0.250000, 0.250000, -0.250000},
265 {-0.250000, -0.250000, 0.250000},
266 {0.250000, -0.250000, 0.250000},
267 {-0.250000, 0.250000, 0.250000},
268 {0.250000, 0.250000, 0.250000}
269
270 };
271
272
273 template <>
274 const double TransformR2UAffine<3>::Kb[GeometryInfo<3>::vertices_per_cell] =
275 {0.500000,
276 0.250000,
277 0.250000,
278 0.000000,
279 0.250000,
280 0.000000,
281 0.000000,
282 -0.250000};
283 } // namespace
284
285
286
287 template <int dim, int spacedim>
288 std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
290 {
292
293 // A = vertex * KA
295
296 for (unsigned int d = 0; d < spacedim; ++d)
297 for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
298 for (unsigned int e = 0; e < dim; ++e)
299 A[d][e] += vertices[v][d] * TransformR2UAffine<dim>::KA[v][e];
300
301 // b = vertex * Kb
303 for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
304 b += vertices[v] * TransformR2UAffine<dim>::Kb[v];
305
306 return std::make_pair(A, b);
307 }
308
309
310
311 template <int dim>
315 const Quadrature<dim> &quadrature)
316 {
318 FEValues<dim> fe_values(mapping, fe, quadrature, update_jacobians);
319
320 Vector<double> aspect_ratio_vector(triangulation.n_active_cells());
321
322 // loop over cells of processor
323 for (const auto &cell : triangulation.active_cell_iterators())
324 {
325 if (cell->is_locally_owned())
326 {
327 double aspect_ratio_cell = 0.0;
328
329 fe_values.reinit(cell);
330
331 // loop over quadrature points
332 for (unsigned int q = 0; q < quadrature.size(); ++q)
333 {
334 const Tensor<2, dim, double> jacobian =
335 Tensor<2, dim, double>(fe_values.jacobian(q));
336
337 // We intentionally do not want to throw an exception in case of
338 // inverted elements since this is not the task of this
339 // function. Instead, inf is written into the vector in case of
340 // inverted elements.
341 if (determinant(jacobian) <= 0)
342 {
343 aspect_ratio_cell = std::numeric_limits<double>::infinity();
344 }
345 else
346 {
348 for (unsigned int i = 0; i < dim; ++i)
349 for (unsigned int j = 0; j < dim; ++j)
350 J(i, j) = jacobian[i][j];
351
352 J.compute_svd();
353
354 const double max_sv = J.singular_value(0);
355 const double min_sv = J.singular_value(dim - 1);
356 const double ar = max_sv / min_sv;
357
358 // Take the max between the previous and the current
359 // aspect ratio value; if we had previously encountered
360 // an inverted cell, we will have placed an infinity
361 // in the aspect_ratio_cell variable, and that value
362 // will survive this max operation.
363 aspect_ratio_cell = std::max(aspect_ratio_cell, ar);
364 }
365 }
366
367 // fill vector
368 aspect_ratio_vector(cell->active_cell_index()) = aspect_ratio_cell;
369 }
370 }
371
372 return aspect_ratio_vector;
373 }
374
375
376
377 template <int dim>
378 double
381 const Quadrature<dim> &quadrature)
382 {
383 Vector<double> aspect_ratio_vector =
384 compute_aspect_ratio_of_cells(mapping, triangulation, quadrature);
385
387 aspect_ratio_vector,
389 }
390
391
392
393 template <int dim, int spacedim>
396 {
397 using iterator =
399 const auto predicate = [](const iterator &) { return true; };
400
402 tria, std::function<bool(const iterator &)>(predicate));
403 }
404
405
406
407 template <int dim, int spacedim>
408 double
410 const Mapping<dim, spacedim> &mapping)
411 {
412 double min_diameter = std::numeric_limits<double>::max();
413 for (const auto &cell : triangulation.active_cell_iterators())
414 if (!cell->is_artificial())
415 min_diameter = std::min(min_diameter, cell->diameter(mapping));
416
417 const double global_min_diameter =
418 Utilities::MPI::min(min_diameter, triangulation.get_mpi_communicator());
419 return global_min_diameter;
420 }
421
422
423
424 template <int dim, int spacedim>
425 double
427 const Mapping<dim, spacedim> &mapping)
428 {
429 double max_diameter = 0.;
430 for (const auto &cell : triangulation.active_cell_iterators())
431 if (!cell->is_artificial())
432 max_diameter = std::max(max_diameter, cell->diameter(mapping));
433
434 const double global_max_diameter =
435 Utilities::MPI::max(max_diameter, triangulation.get_mpi_communicator());
436 return global_max_diameter;
437 }
438} /* namespace GridTools */
439
440
441// explicit instantiations
442#include "grid/grid_tools_geometry.inst"
443
const DerivativeForm< 1, dim, spacedim > & jacobian(const unsigned int q_point) const
double JxW(const unsigned int q_point) const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
number singular_value(const size_type i) const
Abstract base class for mapping classes.
Definition mapping.h:320
Definition point.h:113
unsigned int size() const
const std::vector< Point< spacedim > > & get_vertices() const
cell_iterator end() const
unsigned int n_cells() const
active_cell_iterator begin_active(const unsigned int level=0) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
constexpr bool running_in_debug_mode()
Definition config.h:73
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
static const double KA[GeometryInfo< dim >::vertices_per_cell][dim]
static const double Kb[GeometryInfo< dim >::vertices_per_cell]
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
@ update_JxW_values
Transformed quadrature weights.
@ update_jacobians
Volume element.
std::pair< DerivativeForm< 1, dim, spacedim >, Tensor< 1, spacedim > > affine_cell_approximation(const ArrayView< const Point< spacedim > > &vertices)
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Vector< double > compute_aspect_ratio_of_cells(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
double compute_maximum_aspect_ratio(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
double volume(const Triangulation< dim, spacedim > &tria)
double diameter(const Triangulation< dim, spacedim > &tria)
BoundingBox< spacedim > compute_bounding_box(const Triangulation< dim, spacedim > &triangulation)
double maximal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
std::pair< unsigned int, double > get_longest_direction(typename Triangulation< dim, spacedim >::active_cell_iterator cell)
T sum(const T &t, const MPI_Comm mpi_communicator)
T max(const T &t, const MPI_Comm mpi_communicator)
T min(const T &t, const MPI_Comm mpi_communicator)
double compute_global_error(const Triangulation< dim, spacedim > &tria, const InVector &cellwise_error, const NormType &norm, const double exponent=2.)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)