Reference documentation for deal.II version GIT relicensing-16-gf4f69b3b42 2024-02-27 08:20: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\}}\)
Loading...
Searching...
No Matches
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#ifdef DEBUG
52 if (const auto *p_tria = dynamic_cast<
54 Assert(p_tria->n_global_active_cells() == tria.n_cells(0),
56#endif
57
58 // the algorithm used simply traverses all cells and picks out the
59 // boundary vertices. it may or may not be faster to simply get all
60 // vectors, don't mark boundary vertices, and compute the distances
61 // thereof, but at least as the mesh is refined, it seems better to
62 // first mark boundary nodes, as marking is O(N) in the number of
63 // cells/vertices, while computing the maximal distance is O(N*N)
64 const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
65 std::vector<bool> boundary_vertices(vertices.size(), false);
66
70 tria.end();
71 for (; cell != endc; ++cell)
72 for (const unsigned int face : cell->face_indices())
73 if (cell->face(face)->at_boundary())
74 for (unsigned int i = 0; i < cell->face(face)->n_vertices(); ++i)
75 boundary_vertices[cell->face(face)->vertex_index(i)] = true;
76
77 // now traverse the list of boundary vertices and check distances.
78 // since distances are symmetric, we only have to check one half
79 double max_distance_sqr = 0;
80 std::vector<bool>::const_iterator pi = boundary_vertices.begin();
81 const unsigned int N = boundary_vertices.size();
82 for (unsigned int i = 0; i < N; ++i, ++pi)
83 {
84 std::vector<bool>::const_iterator pj = pi + 1;
85 for (unsigned int j = i + 1; j < N; ++j, ++pj)
86 if ((*pi == true) && (*pj == true) &&
87 ((vertices[i] - vertices[j]).norm_square() > max_distance_sqr))
88 max_distance_sqr = (vertices[i] - vertices[j]).norm_square();
89 }
90
91 return std::sqrt(max_distance_sqr);
92 }
93
94
95
96 template <int dim, int spacedim>
97 double
99 {
102 const ReferenceCell reference_cell = triangulation.get_reference_cells()[0];
103 return volume(
105 reference_cell.template get_default_linear_mapping<dim, spacedim>());
106 }
107
108
109
110 template <int dim, int spacedim>
111 double
114 {
115 // get the degree of the mapping if possible. if not, just assume 1
116 unsigned int mapping_degree = 1;
117 if (const auto *p = dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
118 mapping_degree = p->get_degree();
119 else if (const auto *p =
120 dynamic_cast<const MappingFE<dim, spacedim> *>(&mapping))
121 mapping_degree = p->get_degree();
122
123 // then initialize an appropriate quadrature formula
126 const ReferenceCell reference_cell = triangulation.get_reference_cells()[0];
127 const Quadrature<dim> quadrature_formula =
128 reference_cell.template get_gauss_type_quadrature<dim>(mapping_degree +
129 1);
130 const unsigned int n_q_points = quadrature_formula.size();
131
132 // we really want the JxW values from the FEValues object, but it
133 // wants a finite element. create a cheap element as a dummy
134 // element
135 FE_Nothing<dim, spacedim> dummy_fe(reference_cell);
137 dummy_fe,
138 quadrature_formula,
140
141 double local_volume = 0;
142
143 // compute the integral quantities by quadrature
144 for (const auto &cell : triangulation.active_cell_iterators())
145 if (cell->is_locally_owned())
146 {
147 fe_values.reinit(cell);
148 for (unsigned int q = 0; q < n_q_points; ++q)
149 local_volume += fe_values.JxW(q);
150 }
151
152 const double global_volume =
154
155 return global_volume;
156 }
157
158
159
160 template <int dim, int spacedim>
161 std::pair<unsigned int, double>
164 {
165 double max_ratio = 1;
166 unsigned int index = 0;
167
168 for (unsigned int i = 0; i < dim; ++i)
169 for (unsigned int j = i + 1; j < dim; ++j)
170 {
171 unsigned int ax = i % dim;
172 unsigned int next_ax = j % dim;
173
174 double ratio =
175 cell->extent_in_direction(ax) / cell->extent_in_direction(next_ax);
176
177 if (ratio > max_ratio)
178 {
179 max_ratio = ratio;
180 index = ax;
181 }
182 else if (1.0 / ratio > max_ratio)
183 {
184 max_ratio = 1.0 / ratio;
185 index = next_ax;
186 }
187 }
188 return std::make_pair(index, max_ratio);
189 }
190
191
192
193 namespace
194 {
209 template <int dim>
210 struct TransformR2UAffine
211 {
212 static const double KA[GeometryInfo<dim>::vertices_per_cell][dim];
214 };
215
216
217 /*
218 Octave code:
219 M=[0 1; 1 1];
220 K1 = transpose(M) * inverse (M*transpose(M));
221 printf ("{%f, %f},\n", K1' );
222 */
223 template <>
224 const double TransformR2UAffine<1>::KA[GeometryInfo<1>::vertices_per_cell]
225 [1] = {{-1.000000}, {1.000000}};
226
227 template <>
228 const double TransformR2UAffine<1>::Kb[GeometryInfo<1>::vertices_per_cell] =
229 {1.000000, 0.000000};
230
231
232 /*
233 Octave code:
234 M=[0 1 0 1;0 0 1 1;1 1 1 1];
235 K2 = transpose(M) * inverse (M*transpose(M));
236 printf ("{%f, %f, %f},\n", K2' );
237 */
238 template <>
239 const double TransformR2UAffine<2>::KA[GeometryInfo<2>::vertices_per_cell]
240 [2] = {{-0.500000, -0.500000},
241 {0.500000, -0.500000},
242 {-0.500000, 0.500000},
243 {0.500000, 0.500000}};
244
245 /*
246 Octave code:
247 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];
248 K3 = transpose(M) * inverse (M*transpose(M))
249 printf ("{%f, %f, %f, %f},\n", K3' );
250 */
251 template <>
252 const double TransformR2UAffine<2>::Kb[GeometryInfo<2>::vertices_per_cell] =
253 {0.750000, 0.250000, 0.250000, -0.250000};
254
255
256 template <>
257 const double TransformR2UAffine<3>::KA[GeometryInfo<3>::vertices_per_cell]
258 [3] = {
259 {-0.250000, -0.250000, -0.250000},
260 {0.250000, -0.250000, -0.250000},
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
268 };
269
270
271 template <>
272 const double TransformR2UAffine<3>::Kb[GeometryInfo<3>::vertices_per_cell] =
273 {0.500000,
274 0.250000,
275 0.250000,
276 0.000000,
277 0.250000,
278 0.000000,
279 0.000000,
280 -0.250000};
281 } // namespace
282
283
284
285 template <int dim, int spacedim>
286 std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
288 {
290
291 // A = vertex * KA
293
294 for (unsigned int d = 0; d < spacedim; ++d)
295 for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
296 for (unsigned int e = 0; e < dim; ++e)
297 A[d][e] += vertices[v][d] * TransformR2UAffine<dim>::KA[v][e];
298
299 // b = vertex * Kb
301 for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
302 b += vertices[v] * TransformR2UAffine<dim>::Kb[v];
303
304 return std::make_pair(A, b);
305 }
306
307
308
309 template <int dim>
313 const Quadrature<dim> &quadrature)
314 {
316 FEValues<dim> fe_values(mapping, fe, quadrature, update_jacobians);
317
318 Vector<double> aspect_ratio_vector(triangulation.n_active_cells());
319
320 // loop over cells of processor
321 for (const auto &cell : triangulation.active_cell_iterators())
322 {
323 if (cell->is_locally_owned())
324 {
325 double aspect_ratio_cell = 0.0;
326
327 fe_values.reinit(cell);
328
329 // loop over quadrature points
330 for (unsigned int q = 0; q < quadrature.size(); ++q)
331 {
332 const Tensor<2, dim, double> jacobian =
333 Tensor<2, dim, double>(fe_values.jacobian(q));
334
335 // We intentionally do not want to throw an exception in case of
336 // inverted elements since this is not the task of this
337 // function. Instead, inf is written into the vector in case of
338 // inverted elements.
339 if (determinant(jacobian) <= 0)
340 {
341 aspect_ratio_cell = std::numeric_limits<double>::infinity();
342 }
343 else
344 {
346 for (unsigned int i = 0; i < dim; ++i)
347 for (unsigned int j = 0; j < dim; ++j)
348 J(i, j) = jacobian[i][j];
349
350 J.compute_svd();
351
352 const double max_sv = J.singular_value(0);
353 const double min_sv = J.singular_value(dim - 1);
354 const double ar = max_sv / min_sv;
355
356 // Take the max between the previous and the current
357 // aspect ratio value; if we had previously encountered
358 // an inverted cell, we will have placed an infinity
359 // in the aspect_ratio_cell variable, and that value
360 // will survive this max operation.
361 aspect_ratio_cell = std::max(aspect_ratio_cell, ar);
362 }
363 }
364
365 // fill vector
366 aspect_ratio_vector(cell->active_cell_index()) = aspect_ratio_cell;
367 }
368 }
369
370 return aspect_ratio_vector;
371 }
372
373
374
375 template <int dim>
376 double
379 const Quadrature<dim> &quadrature)
380 {
381 Vector<double> aspect_ratio_vector =
383
385 aspect_ratio_vector,
387 }
388
389
390
391 template <int dim, int spacedim>
394 {
395 using iterator =
397 const auto predicate = [](const iterator &) { return true; };
398
400 tria, std::function<bool(const iterator &)>(predicate));
401 }
402
403
404
405 template <int dim, int spacedim>
406 double
409 {
410 double min_diameter = std::numeric_limits<double>::max();
411 for (const auto &cell : triangulation.active_cell_iterators())
412 if (!cell->is_artificial())
413 min_diameter = std::min(min_diameter, cell->diameter(mapping));
414
415 const double global_min_diameter =
417 return global_min_diameter;
418 }
419
420
421
422 template <int dim, int spacedim>
423 double
426 {
427 double max_diameter = 0.;
428 for (const auto &cell : triangulation.active_cell_iterators())
429 if (!cell->is_artificial())
430 max_diameter = std::max(max_diameter, cell->diameter(mapping));
431
432 const double global_max_diameter =
434 return global_max_diameter;
435 }
436} /* namespace GridTools */
437
438
439// explicit instantiations
440#include "grid_tools_geometry.inst"
441
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:316
Definition point.h:111
unsigned int size() const
virtual MPI_Comm get_communicator() const
unsigned int n_active_cells() const
const std::vector< Point< spacedim > > & get_vertices() const
const std::vector< ReferenceCell > & get_reference_cells() 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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static const double KA[GeometryInfo< dim >::vertices_per_cell][dim]
static const double Kb[GeometryInfo< dim >::vertices_per_cell]
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
@ update_JxW_values
Transformed quadrature weights.
@ update_jacobians
Volume element.
spacedim & mapping
spacedim const Point< spacedim > & p
Definition grid_tools.h:981
const Triangulation< dim, spacedim > & tria
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 >()))
const std::vector< Point< spacedim > > & vertices
Triangulation< dim, spacedim > & triangulation
Definition grid_tools.h:227
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 > &)
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)