Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
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
30namespace 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} /* namespace GridTools */
420
Definition point.h:112
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
unsigned int vertex_indices[2]
#define AssertDimension(dim1, dim2)
double cell_measure< 2 >(const std::vector< Point< 2 > > &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
double cell_measure< 1 >(const std::vector< Point< 1 > > &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
double cell_measure< 3 >(const std::vector< Point< 3 > > &all_vertices, const ArrayView< const unsigned int > &vertex_indices)