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