Reference documentation for deal.II version GIT f8d87ffe95 2022-08-04 17:30: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\}}\)
grid_tools_nontemplates.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2022 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/point.h>
17 
19 
21 
22 #include <vector>
23 
24 // GridTools functions that are template specializations (i.e., only compiled
25 // once without expand_instantiations)
26 
28 
29 
30 namespace GridTools
31 {
32  template <>
33  double
34  cell_measure<1>(const std::vector<Point<1>> & all_vertices,
36  {
38 
39  return all_vertices[vertex_indices[1]][0] -
40  all_vertices[vertex_indices[0]][0];
41  }
42 
43 
44 
45  template <>
46  double
47  cell_measure<2>(const std::vector<Point<2>> & all_vertices,
49  {
50  if (vertex_indices.size() == 3) // triangle
51  {
52  const double x[3] = {all_vertices[vertex_indices[0]](0),
53  all_vertices[vertex_indices[1]](0),
54  all_vertices[vertex_indices[2]](0)};
55 
56  const double y[3] = {all_vertices[vertex_indices[0]](1),
57  all_vertices[vertex_indices[1]](1),
58  all_vertices[vertex_indices[2]](1)};
59 
60  return 0.5 *
61  ((x[0] - x[2]) * (y[1] - y[0]) - (x[1] - x[0]) * (y[0] - y[2]));
62  }
63 
65 
66  /*
67  Get the computation of the measure by this little Maple script. We
68  use the blinear mapping of the unit quad to the real quad. However,
69  every transformation mapping the unit faces to straight lines should
70  do.
71 
72  Remember that the area of the quad is given by
73  \int_K 1 dx dy = \int_{\hat K} |det J| d(xi) d(eta)
74 
75  # x and y are arrays holding the x- and y-values of the four vertices
76  # of this cell in real space.
77  x := array(0..3);
78  y := array(0..3);
79  z := array(0..3);
80  tphi[0] := (1-xi)*(1-eta):
81  tphi[1] := xi*(1-eta):
82  tphi[2] := (1-xi)*eta:
83  tphi[3] := xi*eta:
84  x_real := sum(x[s]*tphi[s], s=0..3):
85  y_real := sum(y[s]*tphi[s], s=0..3):
86  z_real := sum(z[s]*tphi[s], s=0..3):
87 
88  Jxi := <diff(x_real,xi) | diff(y_real,xi) | diff(z_real,xi)>;
89  Jeta := <diff(x_real,eta)| diff(y_real,eta)| diff(z_real,eta)>;
90  with(VectorCalculus):
91  J := CrossProduct(Jxi, Jeta);
92  detJ := sqrt(J[1]^2 + J[2]^2 +J[3]^2);
93 
94  # measure := evalf (Int (Int (detJ, xi=0..1, method = _NCrule ) ,
95  eta=0..1, method = _NCrule ) ): # readlib(C):
96 
97  # C(measure, optimized);
98 
99  additional optimization: divide by 2 only one time
100  */
101 
102  const double x[4] = {all_vertices[vertex_indices[0]](0),
103  all_vertices[vertex_indices[1]](0),
104  all_vertices[vertex_indices[2]](0),
105  all_vertices[vertex_indices[3]](0)};
106 
107  const double y[4] = {all_vertices[vertex_indices[0]](1),
108  all_vertices[vertex_indices[1]](1),
109  all_vertices[vertex_indices[2]](1),
110  all_vertices[vertex_indices[3]](1)};
111 
112  return (-x[1] * y[0] + x[1] * y[3] + y[0] * x[2] + x[0] * y[1] -
113  x[0] * y[2] - y[1] * x[3] - x[2] * y[3] + x[3] * y[2]) /
114  2;
115  }
116 
117 
118 
119  template <>
120  double
121  cell_measure<3>(const std::vector<Point<3>> & all_vertices,
123  {
124  if (vertex_indices.size() == 4) // tetrahedron
125  {
126  const auto &a = all_vertices[vertex_indices[0]];
127  const auto &b = all_vertices[vertex_indices[1]];
128  const auto &c = all_vertices[vertex_indices[2]];
129  const auto &d = all_vertices[vertex_indices[3]];
130 
131  return (1.0 / 6.0) * (d - a) * cross_product_3d(b - a, c - a);
132  }
133  else if (vertex_indices.size() == 5) // pyramid
134  {
135  // This remarkably simple formula comes from Equation 4 of
136  // "Calculation of the volume of a general hexahedron for flow
137  // predictions", Davies and Salmond, AIAA Journal vol. 23 no. 6.
138  const auto &x0 = all_vertices[vertex_indices[0]];
139  const auto &x1 = all_vertices[vertex_indices[1]];
140  const auto &x2 = all_vertices[vertex_indices[2]];
141  const auto &x3 = all_vertices[vertex_indices[3]];
142  const auto &x4 = all_vertices[vertex_indices[4]];
143 
144  const auto v01 = x1 - x0;
145  const auto v02 = x2 - x0;
146  const auto v03 = x3 - x0;
147  const auto v04 = x4 - x0;
148  const auto v21 = x2 - x1;
149 
150  // doing high - low consistently puts us off by -1 from the original
151  // paper in the first term
152  return -v04 * cross_product_3d(v21, v03) / 6.0 +
153  v03 * cross_product_3d(v01, v02) / 12.0;
154  }
155  else if (vertex_indices.size() == 6) // wedge
156  {
157  /* Script used to generate volume code:
158 
159  #!/usr/bin/env python
160  # coding: utf-8
161  import sympy as sp
162  from sympy.simplify.cse_main import cse
163  n_vertices = 6
164  xs = list(sp.symbols(" ".join(["x{}".format(i)
165  for i in range(n_vertices)])))
166  ys = list(sp.symbols(" ".join(["y{}".format(i)
167  for i in range(n_vertices)])))
168  zs = list(sp.symbols(" ".join(["z{}".format(i)
169  for i in range(n_vertices)])))
170  xi, eta, zeta = sp.symbols("xi eta zeta")
171  tphi = [(1 - xi - eta)*(1 - zeta),
172  (xi)*(1 - zeta),
173  (eta)*(1 - zeta),
174  (1 - xi - eta)*(zeta),
175  (xi)*(zeta),
176  (eta)*(zeta)]
177  x_real = sum(xs[i]*tphi[i] for i in range(n_vertices))
178  y_real = sum(ys[i]*tphi[i] for i in range(n_vertices))
179  z_real = sum(zs[i]*tphi[i] for i in range(n_vertices))
180  J = sp.Matrix([[var.diff(v) for v in [xi, eta, zeta]]
181  for var in [x_real, y_real, z_real]])
182  detJ = J.det()
183  detJ2 = detJ.expand().collect(zeta).collect(eta).collect(xi)
184  for x in xs:
185  detJ2 = detJ2.collect(x)
186  for y in ys:
187  detJ2 = detJ2.collect(y)
188  for z in zs:
189  detJ2 = detJ2.collect(z)
190  measure = sp.integrate(sp.integrate(
191  sp.integrate(detJ2, (eta, 0, 1 - xi)),
192  (xi, 0, 1)), (zeta, 0, 1))
193  measure2 = measure
194  for vs in [xs, ys, zs]:
195  for v in vs:
196  measure2 = measure2.collect(v)
197 
198  pairs, expression = cse(measure2)
199  for vertex_no in range(n_vertices):
200  for (coordinate, index) in [('x', 0), ('y', 1), ('z', 2)]:
201  print(
202  "const double {}{} = all_vertices[vertex_indices[{}]][{}];"
203  .format(coordinate, vertex_no, vertex_no, index))
204 
205  for pair in pairs:
206  print("const double " + sp.ccode(pair[0]) + " = "
207  + sp.ccode(pair[1]) + ";")
208  print("const double result = " + sp.ccode(expression[0]) + ";")
209  print("return result;")
210  */
211  const double x0 = all_vertices[vertex_indices[0]][0];
212  const double y0 = all_vertices[vertex_indices[0]][1];
213  const double z0 = all_vertices[vertex_indices[0]][2];
214  const double x1 = all_vertices[vertex_indices[1]][0];
215  const double y1 = all_vertices[vertex_indices[1]][1];
216  const double z1 = all_vertices[vertex_indices[1]][2];
217  const double x2 = all_vertices[vertex_indices[2]][0];
218  const double y2 = all_vertices[vertex_indices[2]][1];
219  const double z2 = all_vertices[vertex_indices[2]][2];
220  const double x3 = all_vertices[vertex_indices[3]][0];
221  const double y3 = all_vertices[vertex_indices[3]][1];
222  const double z3 = all_vertices[vertex_indices[3]][2];
223  const double x4 = all_vertices[vertex_indices[4]][0];
224  const double y4 = all_vertices[vertex_indices[4]][1];
225  const double z4 = all_vertices[vertex_indices[4]][2];
226  const double x5 = all_vertices[vertex_indices[5]][0];
227  const double y5 = all_vertices[vertex_indices[5]][1];
228  const double z5 = all_vertices[vertex_indices[5]][2];
229  const double x6 = (1.0 / 12.0) * z1;
230  const double x7 = -x6;
231  const double x8 = (1.0 / 12.0) * z3;
232  const double x9 = x7 + x8;
233  const double x10 = (1.0 / 12.0) * z2;
234  const double x11 = -x8;
235  const double x12 = x10 + x11;
236  const double x13 = (1.0 / 6.0) * z2;
237  const double x14 = (1.0 / 12.0) * z4;
238  const double x15 = (1.0 / 6.0) * z1;
239  const double x16 = (1.0 / 12.0) * z5;
240  const double x17 = -x16;
241  const double x18 = x16 + x7;
242  const double x19 = -x14;
243  const double x20 = x10 + x19;
244  const double x21 = (1.0 / 12.0) * z0;
245  const double x22 = x19 + x21;
246  const double x23 = -x10;
247  const double x24 = x14 + x23;
248  const double x25 = (1.0 / 6.0) * z0;
249  const double x26 = x17 + x21;
250  const double x27 = x23 + x8;
251  const double x28 = -x21;
252  const double x29 = x16 + x28;
253  const double x30 = x17 + x6;
254  const double x31 = x14 + x28;
255  const double x32 = x11 + x6;
256  const double x33 = (1.0 / 6.0) * z5;
257  const double x34 = (1.0 / 6.0) * z4;
258  const double x35 = (1.0 / 6.0) * z3;
259  const double result =
260  x0 * (x12 * y5 + x9 * y4 + y1 * (-x13 + x14 + x8) +
261  y2 * (x11 + x15 + x17) + y3 * (x18 + x20)) +
262  x1 * (x22 * y3 + x24 * y5 + y0 * (x11 + x13 + x19) +
263  y2 * (x14 + x16 - x25) + y4 * (x26 + x27)) +
264  x2 * (x29 * y3 + x30 * y4 + y0 * (-x15 + x16 + x8) +
265  y1 * (x17 + x19 + x25) + y5 * (x31 + x32)) +
266  x3 * (x26 * y2 + x31 * y1 + y0 * (x24 + x30) + y4 * (x28 + x33 + x7) +
267  y5 * (x10 + x21 - x34)) +
268  x4 * (x18 * y2 + x32 * y0 + y1 * (x12 + x29) + y3 * (x21 - x33 + x6) +
269  y5 * (x23 + x35 + x7)) +
270  x5 * (x20 * y1 + x27 * y0 + y2 * (x22 + x9) + y3 * (x23 + x28 + x34) +
271  y4 * (x10 - x35 + x6));
272  return result;
273  }
274 
276 
277  const double x[8] = {all_vertices[vertex_indices[0]](0),
278  all_vertices[vertex_indices[1]](0),
279  all_vertices[vertex_indices[2]](0),
280  all_vertices[vertex_indices[3]](0),
281  all_vertices[vertex_indices[4]](0),
282  all_vertices[vertex_indices[5]](0),
283  all_vertices[vertex_indices[6]](0),
284  all_vertices[vertex_indices[7]](0)};
285  const double y[8] = {all_vertices[vertex_indices[0]](1),
286  all_vertices[vertex_indices[1]](1),
287  all_vertices[vertex_indices[2]](1),
288  all_vertices[vertex_indices[3]](1),
289  all_vertices[vertex_indices[4]](1),
290  all_vertices[vertex_indices[5]](1),
291  all_vertices[vertex_indices[6]](1),
292  all_vertices[vertex_indices[7]](1)};
293  const double z[8] = {all_vertices[vertex_indices[0]](2),
294  all_vertices[vertex_indices[1]](2),
295  all_vertices[vertex_indices[2]](2),
296  all_vertices[vertex_indices[3]](2),
297  all_vertices[vertex_indices[4]](2),
298  all_vertices[vertex_indices[5]](2),
299  all_vertices[vertex_indices[6]](2),
300  all_vertices[vertex_indices[7]](2)};
301 
302  /*
303  This is the same Maple script as in the barycenter method above
304  except of that here the shape functions tphi[0]-tphi[7] are ordered
305  according to the lexicographic numbering.
306 
307  x := array(0..7):
308  y := array(0..7):
309  z := array(0..7):
310  tphi[0] := (1-xi)*(1-eta)*(1-zeta):
311  tphi[1] := xi*(1-eta)*(1-zeta):
312  tphi[2] := (1-xi)* eta*(1-zeta):
313  tphi[3] := xi* eta*(1-zeta):
314  tphi[4] := (1-xi)*(1-eta)*zeta:
315  tphi[5] := xi*(1-eta)*zeta:
316  tphi[6] := (1-xi)* eta*zeta:
317  tphi[7] := xi* eta*zeta:
318  x_real := sum(x[s]*tphi[s], s=0..7):
319  y_real := sum(y[s]*tphi[s], s=0..7):
320  z_real := sum(z[s]*tphi[s], s=0..7):
321  with (linalg):
322  J := matrix(3,3, [[diff(x_real, xi), diff(x_real, eta), diff(x_real,
323  zeta)], [diff(y_real, xi), diff(y_real, eta), diff(y_real, zeta)],
324  [diff(z_real, xi), diff(z_real, eta), diff(z_real, zeta)]]):
325  detJ := det (J):
326 
327  measure := simplify ( int ( int ( int (detJ, xi=0..1), eta=0..1),
328  zeta=0..1)):
329 
330  readlib(C):
331 
332  C(measure, optimized);
333 
334  The C code produced by this maple script is further optimized by
335  hand. In particular, division by 12 is performed only once, not
336  hundred of times.
337  */
338 
339  const double t3 = y[3] * x[2];
340  const double t5 = z[1] * x[5];
341  const double t9 = z[3] * x[2];
342  const double t11 = x[1] * y[0];
343  const double t14 = x[4] * y[0];
344  const double t18 = x[5] * y[7];
345  const double t20 = y[1] * x[3];
346  const double t22 = y[5] * x[4];
347  const double t26 = z[7] * x[6];
348  const double t28 = x[0] * y[4];
349  const double t34 =
350  z[3] * x[1] * y[2] + t3 * z[1] - t5 * y[7] + y[7] * x[4] * z[6] +
351  t9 * y[6] - t11 * z[4] - t5 * y[3] - t14 * z[2] + z[1] * x[4] * y[0] -
352  t18 * z[3] + t20 * z[0] - t22 * z[0] - y[0] * x[5] * z[4] - t26 * y[3] +
353  t28 * z[2] - t9 * y[1] - y[1] * x[4] * z[0] - t11 * z[5];
354  const double t37 = y[1] * x[0];
355  const double t44 = x[1] * y[5];
356  const double t46 = z[1] * x[0];
357  const double t49 = x[0] * y[2];
358  const double t52 = y[5] * x[7];
359  const double t54 = x[3] * y[7];
360  const double t56 = x[2] * z[0];
361  const double t58 = x[3] * y[2];
362  const double t64 = -x[6] * y[4] * z[2] - t37 * z[2] + t18 * z[6] -
363  x[3] * y[6] * z[2] + t11 * z[2] + t5 * y[0] +
364  t44 * z[4] - t46 * y[4] - t20 * z[7] - t49 * z[6] -
365  t22 * z[1] + t52 * z[3] - t54 * z[2] - t56 * y[4] -
366  t58 * z[0] + y[1] * x[2] * z[0] + t9 * y[7] + t37 * z[4];
367  const double t66 = x[1] * y[7];
368  const double t68 = y[0] * x[6];
369  const double t70 = x[7] * y[6];
370  const double t73 = z[5] * x[4];
371  const double t76 = x[6] * y[7];
372  const double t90 = x[4] * z[0];
373  const double t92 = x[1] * y[3];
374  const double t95 = -t66 * z[3] - t68 * z[2] - t70 * z[2] + t26 * y[5] -
375  t73 * y[6] - t14 * z[6] + t76 * z[2] - t3 * z[6] +
376  x[6] * y[2] * z[4] - z[3] * x[6] * y[2] + t26 * y[4] -
377  t44 * z[3] - x[1] * y[2] * z[0] + x[5] * y[6] * z[4] +
378  t54 * z[5] + t90 * y[2] - t92 * z[2] + t46 * y[2];
379  const double t102 = x[2] * y[0];
380  const double t107 = y[3] * x[7];
381  const double t114 = x[0] * y[6];
382  const double t125 =
383  y[0] * x[3] * z[2] - z[7] * x[5] * y[6] - x[2] * y[6] * z[4] +
384  t102 * z[6] - t52 * z[6] + x[2] * y[4] * z[6] - t107 * z[5] - t54 * z[6] +
385  t58 * z[6] - x[7] * y[4] * z[6] + t37 * z[5] - t114 * z[4] + t102 * z[4] -
386  z[1] * x[2] * y[0] + t28 * z[6] - y[5] * x[6] * z[4] -
387  z[5] * x[1] * y[4] - t73 * y[7];
388  const double t129 = z[0] * x[6];
389  const double t133 = y[1] * x[7];
390  const double t145 = y[1] * x[5];
391  const double t156 = t90 * y[6] - t129 * y[4] + z[7] * x[2] * y[6] -
392  t133 * z[5] + x[5] * y[3] * z[7] - t26 * y[2] -
393  t70 * z[3] + t46 * y[3] + z[5] * x[7] * y[4] +
394  z[7] * x[3] * y[6] - t49 * z[4] + t145 * z[7] -
395  x[2] * y[7] * z[6] + t70 * z[5] + t66 * z[5] -
396  z[7] * x[4] * y[6] + t18 * z[4] + x[1] * y[4] * z[0];
397  const double t160 = x[5] * y[4];
398  const double t165 = z[1] * x[7];
399  const double t178 = z[1] * x[3];
400  const double t181 =
401  t107 * z[6] + t22 * z[7] + t76 * z[3] + t160 * z[1] - x[4] * y[2] * z[6] +
402  t70 * z[4] + t165 * y[5] + x[7] * y[2] * z[6] - t76 * z[5] - t76 * z[4] +
403  t133 * z[3] - t58 * z[1] + y[5] * x[0] * z[4] + t114 * z[2] - t3 * z[7] +
404  t20 * z[2] + t178 * y[7] + t129 * y[2];
405  const double t207 = t92 * z[7] + t22 * z[6] + z[3] * x[0] * y[2] -
406  x[0] * y[3] * z[2] - z[3] * x[7] * y[2] - t165 * y[3] -
407  t9 * y[0] + t58 * z[7] + y[3] * x[6] * z[2] +
408  t107 * z[2] + t73 * y[0] - x[3] * y[5] * z[7] +
409  t3 * z[0] - t56 * y[6] - z[5] * x[0] * y[4] +
410  t73 * y[1] - t160 * z[6] + t160 * z[0];
411  const double t228 = -t44 * z[7] + z[5] * x[6] * y[4] - t52 * z[4] -
412  t145 * z[4] + t68 * z[4] + t92 * z[5] - t92 * z[0] +
413  t11 * z[3] + t44 * z[0] + t178 * y[5] - t46 * y[5] -
414  t178 * y[0] - t145 * z[0] - t20 * z[5] - t37 * z[3] -
415  t160 * z[7] + t145 * z[3] + x[4] * y[6] * z[2];
416 
417  return (t34 + t64 + t95 + t125 + t156 + t181 + t207 + t228) / 12.;
418  }
419 
420 
421 
422  namespace
423  {
424  // the following class is only needed in 2d, so avoid trouble with compilers
425  // warning otherwise
426  class Rotate2d
427  {
428  public:
429  explicit Rotate2d(const double angle)
430  : rotation_matrix(
431  Physics::Transformations::Rotations::rotation_matrix_2d(angle))
432  {}
433  Point<2>
434  operator()(const Point<2> &p) const
435  {
436  return static_cast<Point<2>>(rotation_matrix * p);
437  }
438 
439  private:
441  };
442  } // namespace
443 
444 
445 
446  template <>
447  void
449  {
450  transform(Rotate2d(angle), triangulation);
451  }
452 
453 
454 
455  template <>
456  void
458  {
459  (void)angle;
460  (void)triangulation;
461 
462  AssertThrow(
463  false, ExcMessage("GridTools::rotate() is not available for dim = 3."));
464  }
465 } /* namespace GridTools */
466 
Definition: tensor.h:503
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
unsigned int vertex_indices[2]
Definition: grid_tools.cc:1296
const Tensor< 2, 2, double > rotation_matrix
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
double cell_measure< 1 >(const std::vector< Point< 1 >> &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)
double cell_measure< 3 >(const std::vector< Point< 3 >> &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
void rotate(const double angle, Triangulation< dim > &triangulation)
double cell_measure< 2 >(const std::vector< Point< 2 >> &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Tensor< 2, 2, Number > rotation_matrix_2d(const Number &angle)
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation