Reference documentation for deal.II version Git 99efdf013c 2021-02-27 20:52:41 -0500
\(\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 - 2020 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,
33  const ArrayView<const unsigned int> &vertex_indices)
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,
46  const ArrayView<const unsigned int> &vertex_indices)
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,
120  const ArrayView<const unsigned int> &vertex_indices)
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 
253  AssertDimension(vertex_indices.size(), GeometryInfo<3>::vertices_per_cell);
254  // note that this is the
255  // cell_measure based on the new
256  // deal.II numbering. When called
257  // from inside GridReordering make
258  // sure that you reorder the
259  // vertex_indices before
260  const double x[8] = {all_vertices[vertex_indices[0]](0),
261  all_vertices[vertex_indices[1]](0),
262  all_vertices[vertex_indices[2]](0),
263  all_vertices[vertex_indices[3]](0),
264  all_vertices[vertex_indices[4]](0),
265  all_vertices[vertex_indices[5]](0),
266  all_vertices[vertex_indices[6]](0),
267  all_vertices[vertex_indices[7]](0)};
268  const double y[8] = {all_vertices[vertex_indices[0]](1),
269  all_vertices[vertex_indices[1]](1),
270  all_vertices[vertex_indices[2]](1),
271  all_vertices[vertex_indices[3]](1),
272  all_vertices[vertex_indices[4]](1),
273  all_vertices[vertex_indices[5]](1),
274  all_vertices[vertex_indices[6]](1),
275  all_vertices[vertex_indices[7]](1)};
276  const double z[8] = {all_vertices[vertex_indices[0]](2),
277  all_vertices[vertex_indices[1]](2),
278  all_vertices[vertex_indices[2]](2),
279  all_vertices[vertex_indices[3]](2),
280  all_vertices[vertex_indices[4]](2),
281  all_vertices[vertex_indices[5]](2),
282  all_vertices[vertex_indices[6]](2),
283  all_vertices[vertex_indices[7]](2)};
284 
285  /*
286  This is the same Maple script as in the barycenter method above
287  except of that here the shape functions tphi[0]-tphi[7] are ordered
288  according to the lexicographic numbering.
289 
290  x := array(0..7):
291  y := array(0..7):
292  z := array(0..7):
293  tphi[0] := (1-xi)*(1-eta)*(1-zeta):
294  tphi[1] := xi*(1-eta)*(1-zeta):
295  tphi[2] := (1-xi)* eta*(1-zeta):
296  tphi[3] := xi* eta*(1-zeta):
297  tphi[4] := (1-xi)*(1-eta)*zeta:
298  tphi[5] := xi*(1-eta)*zeta:
299  tphi[6] := (1-xi)* eta*zeta:
300  tphi[7] := xi* eta*zeta:
301  x_real := sum(x[s]*tphi[s], s=0..7):
302  y_real := sum(y[s]*tphi[s], s=0..7):
303  z_real := sum(z[s]*tphi[s], s=0..7):
304  with (linalg):
305  J := matrix(3,3, [[diff(x_real, xi), diff(x_real, eta), diff(x_real,
306  zeta)], [diff(y_real, xi), diff(y_real, eta), diff(y_real, zeta)],
307  [diff(z_real, xi), diff(z_real, eta), diff(z_real, zeta)]]):
308  detJ := det (J):
309 
310  measure := simplify ( int ( int ( int (detJ, xi=0..1), eta=0..1),
311  zeta=0..1)):
312 
313  readlib(C):
314 
315  C(measure, optimized);
316 
317  The C code produced by this maple script is further optimized by
318  hand. In particular, division by 12 is performed only once, not
319  hundred of times.
320  */
321 
322  const double t3 = y[3] * x[2];
323  const double t5 = z[1] * x[5];
324  const double t9 = z[3] * x[2];
325  const double t11 = x[1] * y[0];
326  const double t14 = x[4] * y[0];
327  const double t18 = x[5] * y[7];
328  const double t20 = y[1] * x[3];
329  const double t22 = y[5] * x[4];
330  const double t26 = z[7] * x[6];
331  const double t28 = x[0] * y[4];
332  const double t34 =
333  z[3] * x[1] * y[2] + t3 * z[1] - t5 * y[7] + y[7] * x[4] * z[6] +
334  t9 * y[6] - t11 * z[4] - t5 * y[3] - t14 * z[2] + z[1] * x[4] * y[0] -
335  t18 * z[3] + t20 * z[0] - t22 * z[0] - y[0] * x[5] * z[4] - t26 * y[3] +
336  t28 * z[2] - t9 * y[1] - y[1] * x[4] * z[0] - t11 * z[5];
337  const double t37 = y[1] * x[0];
338  const double t44 = x[1] * y[5];
339  const double t46 = z[1] * x[0];
340  const double t49 = x[0] * y[2];
341  const double t52 = y[5] * x[7];
342  const double t54 = x[3] * y[7];
343  const double t56 = x[2] * z[0];
344  const double t58 = x[3] * y[2];
345  const double t64 = -x[6] * y[4] * z[2] - t37 * z[2] + t18 * z[6] -
346  x[3] * y[6] * z[2] + t11 * z[2] + t5 * y[0] +
347  t44 * z[4] - t46 * y[4] - t20 * z[7] - t49 * z[6] -
348  t22 * z[1] + t52 * z[3] - t54 * z[2] - t56 * y[4] -
349  t58 * z[0] + y[1] * x[2] * z[0] + t9 * y[7] + t37 * z[4];
350  const double t66 = x[1] * y[7];
351  const double t68 = y[0] * x[6];
352  const double t70 = x[7] * y[6];
353  const double t73 = z[5] * x[4];
354  const double t76 = x[6] * y[7];
355  const double t90 = x[4] * z[0];
356  const double t92 = x[1] * y[3];
357  const double t95 = -t66 * z[3] - t68 * z[2] - t70 * z[2] + t26 * y[5] -
358  t73 * y[6] - t14 * z[6] + t76 * z[2] - t3 * z[6] +
359  x[6] * y[2] * z[4] - z[3] * x[6] * y[2] + t26 * y[4] -
360  t44 * z[3] - x[1] * y[2] * z[0] + x[5] * y[6] * z[4] +
361  t54 * z[5] + t90 * y[2] - t92 * z[2] + t46 * y[2];
362  const double t102 = x[2] * y[0];
363  const double t107 = y[3] * x[7];
364  const double t114 = x[0] * y[6];
365  const double t125 =
366  y[0] * x[3] * z[2] - z[7] * x[5] * y[6] - x[2] * y[6] * z[4] +
367  t102 * z[6] - t52 * z[6] + x[2] * y[4] * z[6] - t107 * z[5] - t54 * z[6] +
368  t58 * z[6] - x[7] * y[4] * z[6] + t37 * z[5] - t114 * z[4] + t102 * z[4] -
369  z[1] * x[2] * y[0] + t28 * z[6] - y[5] * x[6] * z[4] -
370  z[5] * x[1] * y[4] - t73 * y[7];
371  const double t129 = z[0] * x[6];
372  const double t133 = y[1] * x[7];
373  const double t145 = y[1] * x[5];
374  const double t156 = t90 * y[6] - t129 * y[4] + z[7] * x[2] * y[6] -
375  t133 * z[5] + x[5] * y[3] * z[7] - t26 * y[2] -
376  t70 * z[3] + t46 * y[3] + z[5] * x[7] * y[4] +
377  z[7] * x[3] * y[6] - t49 * z[4] + t145 * z[7] -
378  x[2] * y[7] * z[6] + t70 * z[5] + t66 * z[5] -
379  z[7] * x[4] * y[6] + t18 * z[4] + x[1] * y[4] * z[0];
380  const double t160 = x[5] * y[4];
381  const double t165 = z[1] * x[7];
382  const double t178 = z[1] * x[3];
383  const double t181 =
384  t107 * z[6] + t22 * z[7] + t76 * z[3] + t160 * z[1] - x[4] * y[2] * z[6] +
385  t70 * z[4] + t165 * y[5] + x[7] * y[2] * z[6] - t76 * z[5] - t76 * z[4] +
386  t133 * z[3] - t58 * z[1] + y[5] * x[0] * z[4] + t114 * z[2] - t3 * z[7] +
387  t20 * z[2] + t178 * y[7] + t129 * y[2];
388  const double t207 = t92 * z[7] + t22 * z[6] + z[3] * x[0] * y[2] -
389  x[0] * y[3] * z[2] - z[3] * x[7] * y[2] - t165 * y[3] -
390  t9 * y[0] + t58 * z[7] + y[3] * x[6] * z[2] +
391  t107 * z[2] + t73 * y[0] - x[3] * y[5] * z[7] +
392  t3 * z[0] - t56 * y[6] - z[5] * x[0] * y[4] +
393  t73 * y[1] - t160 * z[6] + t160 * z[0];
394  const double t228 = -t44 * z[7] + z[5] * x[6] * y[4] - t52 * z[4] -
395  t145 * z[4] + t68 * z[4] + t92 * z[5] - t92 * z[0] +
396  t11 * z[3] + t44 * z[0] + t178 * y[5] - t46 * y[5] -
397  t178 * y[0] - t145 * z[0] - t20 * z[5] - t37 * z[3] -
398  t160 * z[7] + t145 * z[3] + x[4] * y[6] * z[2];
399 
400  return (t34 + t64 + t95 + t125 + t156 + t181 + t207 + t228) / 12.;
401  }
402 
403 
404 
405  namespace
406  {
407  // the following class is only
408  // needed in 2d, so avoid trouble
409  // with compilers warning otherwise
410  class Rotate2d
411  {
412  public:
413  explicit Rotate2d(const double angle)
414  : angle(angle)
415  {}
416  Point<2>
417  operator()(const Point<2> &p) const
418  {
419  return {std::cos(angle) * p(0) - std::sin(angle) * p(1),
420  std::sin(angle) * p(0) + std::cos(angle) * p(1)};
421  }
422 
423  private:
424  const double angle;
425  };
426  } // namespace
427 
428 
429 
430  template <>
431  void
433  {
434  transform(Rotate2d(angle), triangulation);
435  }
436 
437 
438 
439  template <>
440  void
442  {
443  (void)angle;
444  (void)triangulation;
445 
446  AssertThrow(
447  false, ExcMessage("GridTools::rotate() is not available for dim = 3."));
448  }
449 } /* namespace GridTools */
450 
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1623
#define AssertThrow(cond, exc)
Definition: exceptions.h:1576
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:2443
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)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:394
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:393
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)