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