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