deal.II version GIT relicensing-3683-g2da9d4bdba 2025-07-07 18:10:00+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
tria_accessor.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1998 - 2025 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
17
20
21#include <deal.II/fe/fe_q.h>
22#include <deal.II/fe/mapping.h>
23
26#include <deal.II/grid/tria.h>
30
31#include <algorithm>
32#include <array>
33#include <cmath>
34#include <limits>
35
37
38// anonymous namespace for helper functions
39namespace
40{
41 // given the number of face's child
42 // (subface_no), return the number of the
43 // subface concerning the FaceRefineCase of
44 // the face
45 inline unsigned int
46 translate_subface_no(const TriaIterator<TriaAccessor<2, 3, 3>> &face,
47 const unsigned int subface_no)
48 {
49 Assert(face->has_children(), ExcInternalError());
50 Assert(subface_no < face->n_children(), ExcInternalError());
51
52 if (face->child(subface_no)->has_children())
53 // although the subface is refine, it
54 // still matches the face of the cell
55 // invoking the
56 // neighbor_of_coarser_neighbor
57 // function. this means that we are
58 // looking from one cell (anisotropic
59 // child) to a coarser neighbor which is
60 // refined stronger than we are
61 // (isotropically). So we won't be able
62 // to use the neighbor_child_on_subface
63 // function anyway, as the neighbor is
64 // not active. In this case, simply
65 // return the subface_no.
66 return subface_no;
67
68 const bool first_child_has_children = face->child(0)->has_children();
69 // if the first child has children
70 // (FaceRefineCase case_x1y or case_y1x),
71 // then the current subface_no needs to be
72 // 1 and the result of this function is 2,
73 // else simply return the given number,
74 // which is 0 or 1 in an anisotropic case
75 // (case_x, case_y, casex2y or casey2x) or
76 // 0...3 in an isotropic case (case_xy)
77 return subface_no + static_cast<unsigned int>(first_child_has_children);
78 }
79
80
81
82 // given the number of face's child
83 // (subface_no) and grandchild
84 // (subsubface_no), return the number of the
85 // subface concerning the FaceRefineCase of
86 // the face
87 inline unsigned int
88 translate_subface_no(const TriaIterator<TriaAccessor<2, 3, 3>> &face,
89 const unsigned int subface_no,
90 const unsigned int subsubface_no)
91 {
92 Assert(face->has_children(), ExcInternalError());
93 // the subface must be refined, otherwise
94 // we would have ended up in the second
95 // function of this name...
96 Assert(face->child(subface_no)->has_children(), ExcInternalError());
97 Assert(subsubface_no < face->child(subface_no)->n_children(),
99 // This can only be an anisotropic refinement case
100 Assert(face->refinement_case() < RefinementCase<2>::isotropic_refinement,
102
103 const bool first_child_has_children = face->child(0)->has_children();
104
105 static const unsigned int e = numbers::invalid_unsigned_int;
106
107 // array containing the translation of the
108 // numbers,
109 //
110 // first index: subface_no
111 // second index: subsubface_no
112 // third index: does the first subface have children? -> no and yes
113 static const unsigned int translated_subface_no[2][2][2] = {
114 {{e, 0}, // first subface, first subsubface,
115 // first_child_has_children==no and yes
116 {e, 1}}, // first subface, second subsubface,
117 // first_child_has_children==no and yes
118 {{1, 2}, // second subface, first subsubface,
119 // first_child_has_children==no and yes
120 {2, 3}}}; // second subface, second subsubface,
121 // first_child_has_children==no and yes
122
123 Assert(translated_subface_no[subface_no][subsubface_no]
124 [first_child_has_children] != e,
126
127 return translated_subface_no[subface_no][subsubface_no]
128 [first_child_has_children];
129 }
130
131
132 template <int dim, int spacedim>
134 barycenter(const TriaAccessor<1, dim, spacedim> &accessor)
135 {
136 return (accessor.vertex(1) + accessor.vertex(0)) / 2.;
137 }
138
139
141 barycenter(const TriaAccessor<2, 2, 2> &accessor)
142 {
144 {
145 // We define the center in the same way as a simplex barycenter
146 return accessor.center();
147 }
148 else if (accessor.reference_cell() == ReferenceCells::Quadrilateral)
149 {
150 // the evaluation of the formulae
151 // is a bit tricky when done dimension
152 // independently, so we write this function
153 // for 2d and 3d separately
154 /*
155 Get the computation of the barycenter by this little Maple script. We
156 use the bilinear mapping of the unit quad to the real quad. However,
157 every transformation mapping the unit faces to straight lines should
158 do.
159
160 Remember that the area of the quad is given by
161 |K| = \int_K 1 dx dy = \int_{\hat K} |det J| d(xi) d(eta)
162 and that the barycenter is given by
163 \vec x_s = 1/|K| \int_K \vec x dx dy
164 = 1/|K| \int_{\hat K} \vec x(xi,eta) |det J| d(xi) d(eta)
165
166 # x and y are arrays holding the x- and y-values of the four vertices
167 # of this cell in real space.
168 x := array(0..3);
169 y := array(0..3);
170 tphi[0] := (1-xi)*(1-eta):
171 tphi[1] := xi*(1-eta):
172 tphi[2] := (1-xi)*eta:
173 tphi[3] := xi*eta:
174 x_real := sum(x[s]*tphi[s], s=0..3):
175 y_real := sum(y[s]*tphi[s], s=0..3):
176 detJ := diff(x_real,xi)*diff(y_real,eta) -
177 diff(x_real,eta)*diff(y_real,xi):
178
179 measure := simplify ( int ( int (detJ, xi=0..1), eta=0..1)):
180
181 xs := simplify (1/measure * int ( int (x_real * detJ, xi=0..1),
182 eta=0..1)): ys := simplify (1/measure * int ( int (y_real * detJ,
183 xi=0..1), eta=0..1)): readlib(C):
184
185 C(array(1..2, [xs, ys]), optimized);
186 */
187
188 const double x[4] = {accessor.vertex(0)[0],
189 accessor.vertex(1)[0],
190 accessor.vertex(2)[0],
191 accessor.vertex(3)[0]};
192 const double y[4] = {accessor.vertex(0)[1],
193 accessor.vertex(1)[1],
194 accessor.vertex(2)[1],
195 accessor.vertex(3)[1]};
196 const double t1 = x[0] * x[1];
197 const double t3 = x[0] * x[0];
198 const double t5 = x[1] * x[1];
199 const double t9 = y[0] * x[0];
200 const double t11 = y[1] * x[1];
201 const double t14 = x[2] * x[2];
202 const double t16 = x[3] * x[3];
203 const double t20 = x[2] * x[3];
204 const double t27 = t1 * y[1] + t3 * y[1] - t5 * y[0] - t3 * y[2] +
205 t5 * y[3] + t9 * x[2] - t11 * x[3] - t1 * y[0] -
206 t14 * y[3] + t16 * y[2] - t16 * y[1] + t14 * y[0] -
207 t20 * y[3] - x[0] * x[2] * y[2] +
208 x[1] * x[3] * y[3] + t20 * y[2];
209 const double t37 =
210 1 / (-x[1] * y[0] + x[1] * y[3] + y[0] * x[2] + x[0] * y[1] -
211 x[0] * y[2] - y[1] * x[3] - x[2] * y[3] + x[3] * y[2]);
212 const double t39 = y[2] * y[2];
213 const double t51 = y[0] * y[0];
214 const double t53 = y[1] * y[1];
215 const double t59 = y[3] * y[3];
216 const double t63 =
217 t39 * x[3] + y[2] * y[0] * x[2] + y[3] * x[3] * y[2] -
218 y[2] * x[2] * y[3] - y[3] * y[1] * x[3] - t9 * y[2] + t11 * y[3] +
219 t51 * x[2] - t53 * x[3] - x[1] * t51 + t9 * y[1] - t11 * y[0] +
220 x[0] * t53 - t59 * x[2] + t59 * x[1] - t39 * x[0];
221
222 return {t27 * t37 / 3, t63 * t37 / 3};
223 }
224 else
225 {
227 return {};
228 }
229 }
230
231
232
234 barycenter(const TriaAccessor<3, 3, 3> &accessor)
235 {
237 {
238 // We define the center in the same way as a simplex barycenter
239 return accessor.center();
240 }
241 else if (accessor.reference_cell() == ReferenceCells::Hexahedron)
242 {
243 /*
244 Get the computation of the barycenter by this little Maple script. We
245 use the trilinear mapping of the unit hex to the real hex.
246
247 Remember that the area of the hex is given by
248 |K| = \int_K 1 dx dy dz = \int_{\hat K} |det J| d(xi) d(eta) d(zeta)
249 and that the barycenter is given by
250 \vec x_s = 1/|K| \int_K \vec x dx dy dz
251 = 1/|K| \int_{\hat K} \vec x(xi,eta,zeta) |det J| d(xi) d(eta) d(zeta)
252
253 Note, that in the ordering of the shape functions tphi[0]-tphi[7]
254 below, eta and zeta have been exchanged (zeta belongs to the y, and
255 eta to the z direction). However, the resulting Jacobian determinant
256 detJ should be the same, as a matrix and the matrix created from it
257 by exchanging two consecutive lines and two neighboring columns have
258 the same determinant.
259
260 # x, y and z are arrays holding the x-, y- and z-values of the four
261 vertices # of this cell in real space. x := array(0..7): y :=
262 array(0..7): z := array(0..7): tphi[0] := (1-xi)*(1-eta)*(1-zeta):
263 tphi[1] := xi*(1-eta)*(1-zeta):
264 tphi[2] := xi*eta*(1-zeta):
265 tphi[3] := (1-xi)*eta*(1-zeta):
266 tphi[4] := (1-xi)*(1-eta)*zeta:
267 tphi[5] := xi*(1-eta)*zeta:
268 tphi[6] := xi*eta*zeta:
269 tphi[7] := (1-xi)*eta*zeta:
270 x_real := sum(x[s]*tphi[s], s=0..7):
271 y_real := sum(y[s]*tphi[s], s=0..7):
272 z_real := sum(z[s]*tphi[s], s=0..7):
273 with (linalg):
274 J := matrix(3,3, [[diff(x_real, xi), diff(x_real, eta), diff(x_real,
275 zeta)], [diff(y_real, xi), diff(y_real, eta), diff(y_real, zeta)],
276 [diff(z_real, xi), diff(z_real, eta), diff(z_real, zeta)]]):
277 detJ := det (J):
278
279 measure := simplify ( int ( int ( int (detJ, xi=0..1), eta=0..1),
280 zeta=0..1)):
281
282 xs := simplify (1/measure * int ( int ( int (x_real * detJ, xi=0..1),
283 eta=0..1), zeta=0..1)): ys := simplify (1/measure * int ( int ( int
284 (y_real * detJ, xi=0..1), eta=0..1), zeta=0..1)): zs := simplify
285 (1/measure * int ( int ( int (z_real * detJ, xi=0..1), eta=0..1),
286 zeta=0..1)):
287
288 readlib(C):
289
290 C(array(1..3, [xs, ys, zs]));
291
292
293 This script takes more than several hours when using an old version
294 of maple on an old and slow computer. Therefore, when changing to
295 the new deal.II numbering scheme (lexicographic numbering) the code
296 lines below have not been reproduced with maple but only the
297 ordering of points in the definitions of x[], y[] and z[] have been
298 changed.
299
300 For the case, someone is willing to rerun the maple script, he/she
301 should use following ordering of shape functions:
302
303 tphi[0] := (1-xi)*(1-eta)*(1-zeta):
304 tphi[1] := xi*(1-eta)*(1-zeta):
305 tphi[2] := (1-xi)* eta*(1-zeta):
306 tphi[3] := xi* eta*(1-zeta):
307 tphi[4] := (1-xi)*(1-eta)*zeta:
308 tphi[5] := xi*(1-eta)*zeta:
309 tphi[6] := (1-xi)* eta*zeta:
310 tphi[7] := xi* eta*zeta:
311
312 and change the ordering of points in the definitions of x[], y[] and
313 z[] back to the standard ordering.
314 */
315
316 const double x[8] = {accessor.vertex(0)[0],
317 accessor.vertex(1)[0],
318 accessor.vertex(5)[0],
319 accessor.vertex(4)[0],
320 accessor.vertex(2)[0],
321 accessor.vertex(3)[0],
322 accessor.vertex(7)[0],
323 accessor.vertex(6)[0]};
324 const double y[8] = {accessor.vertex(0)[1],
325 accessor.vertex(1)[1],
326 accessor.vertex(5)[1],
327 accessor.vertex(4)[1],
328 accessor.vertex(2)[1],
329 accessor.vertex(3)[1],
330 accessor.vertex(7)[1],
331 accessor.vertex(6)[1]};
332 const double z[8] = {accessor.vertex(0)[2],
333 accessor.vertex(1)[2],
334 accessor.vertex(5)[2],
335 accessor.vertex(4)[2],
336 accessor.vertex(2)[2],
337 accessor.vertex(3)[2],
338 accessor.vertex(7)[2],
339 accessor.vertex(6)[2]};
340
341 double s1, s2, s3, s4, s5, s6, s7, s8;
342
343 s1 = 1.0 / 6.0;
344 s8 = -x[2] * x[2] * y[0] * z[3] - 2.0 * z[6] * x[7] * x[7] * y[4] -
345 z[5] * x[7] * x[7] * y[4] - z[6] * x[7] * x[7] * y[5] +
346 2.0 * y[6] * x[7] * x[7] * z[4] - z[5] * x[6] * x[6] * y[4] +
347 x[6] * x[6] * y[4] * z[7] - z[1] * x[0] * x[0] * y[2] -
348 x[6] * x[6] * y[7] * z[4] + 2.0 * x[6] * x[6] * y[5] * z[7] -
349 2.0 * x[6] * x[6] * y[7] * z[5] + y[5] * x[6] * x[6] * z[4] +
350 2.0 * x[5] * x[5] * y[4] * z[6] + x[0] * x[0] * y[7] * z[4] -
351 2.0 * x[5] * x[5] * y[6] * z[4];
352 s7 = s8 - y[6] * x[5] * x[5] * z[7] + z[6] * x[5] * x[5] * y[7] -
353 y[1] * x[0] * x[0] * z[5] + x[7] * z[5] * x[4] * y[7] -
354 x[7] * y[6] * x[5] * z[7] - 2.0 * x[7] * x[6] * y[7] * z[4] +
355 2.0 * x[7] * x[6] * y[4] * z[7] - x[7] * x[5] * y[7] * z[4] -
356 2.0 * x[7] * y[6] * x[4] * z[7] - x[7] * y[5] * x[4] * z[7] +
357 x[2] * x[2] * y[3] * z[0] - x[7] * x[6] * y[7] * z[5] +
358 x[7] * x[6] * y[5] * z[7] + 2.0 * x[1] * x[1] * y[0] * z[5] +
359 x[7] * z[6] * x[5] * y[7];
360 s8 = -2.0 * x[1] * x[1] * y[5] * z[0] + z[1] * x[0] * x[0] * y[5] +
361 2.0 * x[2] * x[2] * y[3] * z[1] - z[5] * x[4] * x[4] * y[1] +
362 y[5] * x[4] * x[4] * z[1] - 2.0 * x[5] * x[5] * y[4] * z[1] +
363 2.0 * x[5] * x[5] * y[1] * z[4] - 2.0 * x[2] * x[2] * y[1] * z[3] -
364 y[1] * x[2] * x[2] * z[0] + x[7] * y[2] * x[3] * z[7] +
365 x[7] * z[2] * x[6] * y[3] + 2.0 * x[7] * z[6] * x[4] * y[7] +
366 z[5] * x[1] * x[1] * y[4] + z[1] * x[2] * x[2] * y[0] -
367 2.0 * y[0] * x[3] * x[3] * z[7];
368 s6 = s8 + 2.0 * z[0] * x[3] * x[3] * y[7] - x[7] * x[2] * y[3] * z[7] -
369 x[7] * z[2] * x[3] * y[7] + x[7] * x[2] * y[7] * z[3] -
370 x[7] * y[2] * x[6] * z[3] + x[4] * x[5] * y[1] * z[4] -
371 x[4] * x[5] * y[4] * z[1] + x[4] * z[5] * x[1] * y[4] -
372 x[4] * y[5] * x[1] * z[4] - 2.0 * x[5] * z[5] * x[4] * y[1] -
373 2.0 * x[5] * y[5] * x[1] * z[4] + 2.0 * x[5] * z[5] * x[1] * y[4] +
374 2.0 * x[5] * y[5] * x[4] * z[1] - x[6] * z[5] * x[7] * y[4] -
375 z[2] * x[3] * x[3] * y[6] + s7;
376 s8 = -2.0 * x[6] * z[6] * x[7] * y[5] - x[6] * y[6] * x[4] * z[7] +
377 y[2] * x[3] * x[3] * z[6] + x[6] * y[6] * x[7] * z[4] +
378 2.0 * y[2] * x[3] * x[3] * z[7] + x[0] * x[1] * y[0] * z[5] +
379 x[0] * y[1] * x[5] * z[0] - x[0] * z[1] * x[5] * y[0] -
380 2.0 * z[2] * x[3] * x[3] * y[7] + 2.0 * x[6] * z[6] * x[5] * y[7] -
381 x[0] * x[1] * y[5] * z[0] - x[6] * y[5] * x[4] * z[6] -
382 2.0 * x[3] * z[0] * x[7] * y[3] - x[6] * z[6] * x[7] * y[4] -
383 2.0 * x[1] * z[1] * x[5] * y[0];
384 s7 = s8 + 2.0 * x[1] * y[1] * x[5] * z[0] +
385 2.0 * x[1] * z[1] * x[0] * y[5] + 2.0 * x[3] * y[0] * x[7] * z[3] +
386 2.0 * x[3] * x[0] * y[3] * z[7] - 2.0 * x[3] * x[0] * y[7] * z[3] -
387 2.0 * x[1] * y[1] * x[0] * z[5] - 2.0 * x[6] * y[6] * x[5] * z[7] +
388 s6 - y[5] * x[1] * x[1] * z[4] + x[6] * z[6] * x[4] * y[7] -
389 2.0 * x[2] * y[2] * x[3] * z[1] + x[6] * z[5] * x[4] * y[6] +
390 x[6] * x[5] * y[4] * z[6] - y[6] * x[7] * x[7] * z[2] -
391 x[6] * x[5] * y[6] * z[4];
392 s8 = x[3] * x[3] * y[7] * z[4] - 2.0 * y[6] * x[7] * x[7] * z[3] +
393 z[6] * x[7] * x[7] * y[2] + 2.0 * z[6] * x[7] * x[7] * y[3] +
394 2.0 * y[1] * x[0] * x[0] * z[3] + 2.0 * x[0] * x[1] * y[3] * z[0] -
395 2.0 * x[0] * y[0] * x[3] * z[4] - 2.0 * x[0] * z[1] * x[4] * y[0] -
396 2.0 * x[0] * y[1] * x[3] * z[0] + 2.0 * x[0] * y[0] * x[4] * z[3] -
397 2.0 * x[0] * z[0] * x[4] * y[3] + 2.0 * x[0] * x[1] * y[0] * z[4] +
398 2.0 * x[0] * z[1] * x[3] * y[0] - 2.0 * x[0] * x[1] * y[0] * z[3] -
399 2.0 * x[0] * x[1] * y[4] * z[0] + 2.0 * x[0] * y[1] * x[4] * z[0];
400 s5 = s8 + 2.0 * x[0] * z[0] * x[3] * y[4] + x[1] * y[1] * x[0] * z[3] -
401 x[1] * z[1] * x[4] * y[0] - x[1] * y[1] * x[0] * z[4] +
402 x[1] * z[1] * x[0] * y[4] - x[1] * y[1] * x[3] * z[0] -
403 x[1] * z[1] * x[0] * y[3] - x[0] * z[5] * x[4] * y[1] +
404 x[0] * y[5] * x[4] * z[1] - 2.0 * x[4] * x[0] * y[4] * z[7] -
405 2.0 * x[4] * y[5] * x[0] * z[4] + 2.0 * x[4] * z[5] * x[0] * y[4] -
406 2.0 * x[4] * x[5] * y[4] * z[0] - 2.0 * x[4] * y[0] * x[7] * z[4] -
407 x[5] * y[5] * x[0] * z[4] + s7;
408 s8 = x[5] * z[5] * x[0] * y[4] - x[5] * z[5] * x[4] * y[0] +
409 x[1] * z[5] * x[0] * y[4] + x[5] * y[5] * x[4] * z[0] -
410 x[0] * y[0] * x[7] * z[4] - x[0] * z[5] * x[4] * y[0] -
411 x[1] * y[5] * x[0] * z[4] + x[0] * z[0] * x[7] * y[4] +
412 x[0] * y[5] * x[4] * z[0] - x[0] * z[0] * x[4] * y[7] +
413 x[0] * x[5] * y[0] * z[4] + x[0] * y[0] * x[4] * z[7] -
414 x[0] * x[5] * y[4] * z[0] - x[3] * x[3] * y[4] * z[7] +
415 2.0 * x[2] * z[2] * x[3] * y[1];
416 s7 = s8 - x[5] * x[5] * y[4] * z[0] + 2.0 * y[5] * x[4] * x[4] * z[0] -
417 2.0 * z[0] * x[4] * x[4] * y[7] + 2.0 * y[0] * x[4] * x[4] * z[7] -
418 2.0 * z[5] * x[4] * x[4] * y[0] + x[5] * x[5] * y[4] * z[7] -
419 x[5] * x[5] * y[7] * z[4] - 2.0 * y[5] * x[4] * x[4] * z[7] +
420 2.0 * z[5] * x[4] * x[4] * y[7] - x[0] * x[0] * y[7] * z[3] +
421 y[2] * x[0] * x[0] * z[3] + x[0] * x[0] * y[3] * z[7] -
422 x[5] * x[1] * y[4] * z[0] + x[5] * y[1] * x[4] * z[0] -
423 x[4] * y[0] * x[3] * z[4];
424 s8 = -x[4] * y[1] * x[0] * z[4] + x[4] * z[1] * x[0] * y[4] +
425 x[4] * x[0] * y[3] * z[4] - x[4] * x[0] * y[4] * z[3] +
426 x[4] * x[1] * y[0] * z[4] - x[4] * x[1] * y[4] * z[0] +
427 x[4] * z[0] * x[3] * y[4] + x[5] * x[1] * y[0] * z[4] +
428 x[1] * z[1] * x[3] * y[0] + x[1] * y[1] * x[4] * z[0] -
429 x[5] * z[1] * x[4] * y[0] - 2.0 * y[1] * x[0] * x[0] * z[4] +
430 2.0 * z[1] * x[0] * x[0] * y[4] + 2.0 * x[0] * x[0] * y[3] * z[4] -
431 2.0 * z[1] * x[0] * x[0] * y[3];
432 s6 = s8 - 2.0 * x[0] * x[0] * y[4] * z[3] + x[1] * x[1] * y[3] * z[0] +
433 x[1] * x[1] * y[0] * z[4] - x[1] * x[1] * y[0] * z[3] -
434 x[1] * x[1] * y[4] * z[0] - z[1] * x[4] * x[4] * y[0] +
435 y[0] * x[4] * x[4] * z[3] - z[0] * x[4] * x[4] * y[3] +
436 y[1] * x[4] * x[4] * z[0] - x[0] * x[0] * y[4] * z[7] -
437 y[5] * x[0] * x[0] * z[4] + z[5] * x[0] * x[0] * y[4] +
438 x[5] * x[5] * y[0] * z[4] - x[0] * y[0] * x[3] * z[7] +
439 x[0] * z[0] * x[3] * y[7] + s7;
440 s8 = s6 + x[0] * x[2] * y[3] * z[0] - x[0] * x[2] * y[0] * z[3] +
441 x[0] * y[0] * x[7] * z[3] - x[0] * y[2] * x[3] * z[0] +
442 x[0] * z[2] * x[3] * y[0] - x[0] * z[0] * x[7] * y[3] +
443 x[1] * x[2] * y[3] * z[0] - z[2] * x[0] * x[0] * y[3] +
444 x[3] * z[2] * x[6] * y[3] - x[3] * x[2] * y[3] * z[6] +
445 x[3] * x[2] * y[6] * z[3] - x[3] * y[2] * x[6] * z[3] -
446 2.0 * x[3] * y[2] * x[7] * z[3] + 2.0 * x[3] * z[2] * x[7] * y[3];
447 s7 = s8 + 2.0 * x[4] * y[5] * x[7] * z[4] +
448 2.0 * x[4] * x[5] * y[4] * z[7] - 2.0 * x[4] * z[5] * x[7] * y[4] -
449 2.0 * x[4] * x[5] * y[7] * z[4] + x[5] * y[5] * x[7] * z[4] -
450 x[5] * z[5] * x[7] * y[4] - x[5] * y[5] * x[4] * z[7] +
451 x[5] * z[5] * x[4] * y[7] + 2.0 * x[3] * x[2] * y[7] * z[3] -
452 2.0 * x[2] * z[2] * x[1] * y[3] + 2.0 * x[4] * z[0] * x[7] * y[4] +
453 2.0 * x[4] * x[0] * y[7] * z[4] + 2.0 * x[4] * x[5] * y[0] * z[4] -
454 x[7] * x[6] * y[2] * z[7] - 2.0 * x[3] * x[2] * y[3] * z[7] -
455 x[0] * x[4] * y[7] * z[3];
456 s8 = x[0] * x[3] * y[7] * z[4] - x[0] * x[3] * y[4] * z[7] +
457 x[0] * x[4] * y[3] * z[7] - 2.0 * x[7] * z[6] * x[3] * y[7] +
458 x[3] * x[7] * y[4] * z[3] - x[3] * x[4] * y[7] * z[3] -
459 x[3] * x[7] * y[3] * z[4] + x[3] * x[4] * y[3] * z[7] +
460 2.0 * x[2] * y[2] * x[1] * z[3] + y[6] * x[3] * x[3] * z[7] -
461 z[6] * x[3] * x[3] * y[7] - x[1] * z[5] * x[4] * y[1] -
462 x[1] * x[5] * y[4] * z[1] - x[1] * z[2] * x[0] * y[3] -
463 x[1] * x[2] * y[0] * z[3] + x[1] * y[2] * x[0] * z[3];
464 s4 = s8 + x[1] * x[5] * y[1] * z[4] + x[1] * y[5] * x[4] * z[1] +
465 x[4] * y[0] * x[7] * z[3] - x[4] * z[0] * x[7] * y[3] -
466 x[4] * x[4] * y[7] * z[3] + x[4] * x[4] * y[3] * z[7] +
467 x[3] * z[6] * x[7] * y[3] - x[3] * x[6] * y[3] * z[7] +
468 x[3] * x[6] * y[7] * z[3] - x[3] * z[6] * x[2] * y[7] -
469 x[3] * y[6] * x[7] * z[3] + x[3] * z[6] * x[7] * y[2] +
470 x[3] * y[6] * x[2] * z[7] + 2.0 * x[5] * z[5] * x[4] * y[6] + s5 +
471 s7;
472 s8 = s4 - 2.0 * x[5] * z[5] * x[6] * y[4] - x[5] * z[6] * x[7] * y[5] +
473 x[5] * x[6] * y[5] * z[7] - x[5] * x[6] * y[7] * z[5] -
474 2.0 * x[5] * y[5] * x[4] * z[6] + 2.0 * x[5] * y[5] * x[6] * z[4] -
475 x[3] * y[6] * x[7] * z[2] + x[4] * x[7] * y[4] * z[3] +
476 x[4] * x[3] * y[7] * z[4] - x[4] * x[7] * y[3] * z[4] -
477 x[4] * x[3] * y[4] * z[7] - z[1] * x[5] * x[5] * y[0] +
478 y[1] * x[5] * x[5] * z[0] + x[4] * y[6] * x[7] * z[4];
479 s7 = s8 - x[4] * x[6] * y[7] * z[4] + x[4] * x[6] * y[4] * z[7] -
480 x[4] * z[6] * x[7] * y[4] - x[5] * y[6] * x[4] * z[7] -
481 x[5] * x[6] * y[7] * z[4] + x[5] * x[6] * y[4] * z[7] +
482 x[5] * z[6] * x[4] * y[7] - y[6] * x[4] * x[4] * z[7] +
483 z[6] * x[4] * x[4] * y[7] + x[7] * x[5] * y[4] * z[7] -
484 y[2] * x[7] * x[7] * z[3] + z[2] * x[7] * x[7] * y[3] -
485 y[0] * x[3] * x[3] * z[4] - y[1] * x[3] * x[3] * z[0] +
486 z[1] * x[3] * x[3] * y[0];
487 s8 = z[0] * x[3] * x[3] * y[4] - x[2] * y[1] * x[3] * z[0] +
488 x[2] * z[1] * x[3] * y[0] + x[3] * y[1] * x[0] * z[3] +
489 x[3] * x[1] * y[3] * z[0] + x[3] * x[0] * y[3] * z[4] -
490 x[3] * z[1] * x[0] * y[3] - x[3] * x[0] * y[4] * z[3] +
491 x[3] * y[0] * x[4] * z[3] - x[3] * z[0] * x[4] * y[3] -
492 x[3] * x[1] * y[0] * z[3] + x[3] * z[0] * x[7] * y[4] -
493 x[3] * y[0] * x[7] * z[4] + z[0] * x[7] * x[7] * y[4] -
494 y[0] * x[7] * x[7] * z[4];
495 s6 = s8 + y[1] * x[0] * x[0] * z[2] - 2.0 * y[2] * x[3] * x[3] * z[0] +
496 2.0 * z[2] * x[3] * x[3] * y[0] - 2.0 * x[1] * x[1] * y[0] * z[2] +
497 2.0 * x[1] * x[1] * y[2] * z[0] - y[2] * x[3] * x[3] * z[1] +
498 z[2] * x[3] * x[3] * y[1] - y[5] * x[4] * x[4] * z[6] +
499 z[5] * x[4] * x[4] * y[6] + x[7] * x[0] * y[7] * z[4] -
500 x[7] * z[0] * x[4] * y[7] - x[7] * x[0] * y[4] * z[7] +
501 x[7] * y[0] * x[4] * z[7] - x[0] * x[1] * y[0] * z[2] +
502 x[0] * z[1] * x[2] * y[0] + s7;
503 s8 = s6 + x[0] * x[1] * y[2] * z[0] - x[0] * y[1] * x[2] * z[0] -
504 x[3] * z[1] * x[0] * y[2] + 2.0 * x[3] * x[2] * y[3] * z[0] +
505 y[0] * x[7] * x[7] * z[3] - z[0] * x[7] * x[7] * y[3] -
506 2.0 * x[3] * z[2] * x[0] * y[3] - 2.0 * x[3] * x[2] * y[0] * z[3] +
507 2.0 * x[3] * y[2] * x[0] * z[3] + x[3] * x[2] * y[3] * z[1] -
508 x[3] * x[2] * y[1] * z[3] - x[5] * y[1] * x[0] * z[5] +
509 x[3] * y[1] * x[0] * z[2] + x[4] * y[6] * x[7] * z[5];
510 s7 = s8 - x[5] * x[1] * y[5] * z[0] + 2.0 * x[1] * z[1] * x[2] * y[0] -
511 2.0 * x[1] * z[1] * x[0] * y[2] + x[1] * x[2] * y[3] * z[1] -
512 x[1] * x[2] * y[1] * z[3] + 2.0 * x[1] * y[1] * x[0] * z[2] -
513 2.0 * x[1] * y[1] * x[2] * z[0] - z[2] * x[1] * x[1] * y[3] +
514 y[2] * x[1] * x[1] * z[3] + y[5] * x[7] * x[7] * z[4] +
515 y[6] * x[7] * x[7] * z[5] + x[7] * x[6] * y[7] * z[2] +
516 x[7] * y[6] * x[2] * z[7] - x[7] * z[6] * x[2] * y[7] -
517 2.0 * x[7] * x[6] * y[3] * z[7];
518 s8 = s7 + 2.0 * x[7] * x[6] * y[7] * z[3] +
519 2.0 * x[7] * y[6] * x[3] * z[7] - x[3] * z[2] * x[1] * y[3] +
520 x[3] * y[2] * x[1] * z[3] + x[5] * x[1] * y[0] * z[5] +
521 x[4] * y[5] * x[6] * z[4] + x[5] * z[1] * x[0] * y[5] -
522 x[4] * z[6] * x[7] * y[5] - x[4] * x[5] * y[6] * z[4] +
523 x[4] * x[5] * y[4] * z[6] - x[4] * z[5] * x[6] * y[4] -
524 x[1] * y[2] * x[3] * z[1] + x[1] * z[2] * x[3] * y[1] -
525 x[2] * x[1] * y[0] * z[2] - x[2] * z[1] * x[0] * y[2];
526 s5 = s8 + x[2] * x[1] * y[2] * z[0] - x[2] * z[2] * x[0] * y[3] +
527 x[2] * y[2] * x[0] * z[3] - x[2] * y[2] * x[3] * z[0] +
528 x[2] * z[2] * x[3] * y[0] + x[2] * y[1] * x[0] * z[2] +
529 x[5] * y[6] * x[7] * z[5] + x[6] * y[5] * x[7] * z[4] +
530 2.0 * x[6] * y[6] * x[7] * z[5] - x[7] * y[0] * x[3] * z[7] +
531 x[7] * z[0] * x[3] * y[7] - x[7] * x[0] * y[7] * z[3] +
532 x[7] * x[0] * y[3] * z[7] + 2.0 * x[7] * x[7] * y[4] * z[3] -
533 2.0 * x[7] * x[7] * y[3] * z[4] - 2.0 * x[1] * x[1] * y[2] * z[5];
534 s8 = s5 - 2.0 * x[7] * x[4] * y[7] * z[3] +
535 2.0 * x[7] * x[3] * y[7] * z[4] - 2.0 * x[7] * x[3] * y[4] * z[7] +
536 2.0 * x[7] * x[4] * y[3] * z[7] + 2.0 * x[1] * x[1] * y[5] * z[2] -
537 x[1] * x[1] * y[2] * z[6] + x[1] * x[1] * y[6] * z[2] +
538 z[1] * x[5] * x[5] * y[2] - y[1] * x[5] * x[5] * z[2] -
539 x[1] * x[1] * y[6] * z[5] + x[1] * x[1] * y[5] * z[6] +
540 x[5] * x[5] * y[6] * z[2] - x[5] * x[5] * y[2] * z[6] -
541 2.0 * y[1] * x[5] * x[5] * z[6];
542 s7 = s8 + 2.0 * z[1] * x[5] * x[5] * y[6] +
543 2.0 * x[1] * z[1] * x[5] * y[2] + 2.0 * x[1] * y[1] * x[2] * z[5] -
544 2.0 * x[1] * z[1] * x[2] * y[5] - 2.0 * x[1] * y[1] * x[5] * z[2] -
545 x[1] * y[1] * x[6] * z[2] - x[1] * z[1] * x[2] * y[6] +
546 x[1] * z[1] * x[6] * y[2] + x[1] * y[1] * x[2] * z[6] -
547 x[5] * x[1] * y[2] * z[5] + x[5] * y[1] * x[2] * z[5] -
548 x[5] * z[1] * x[2] * y[5] + x[5] * x[1] * y[5] * z[2] -
549 x[5] * y[1] * x[6] * z[2] - x[5] * x[1] * y[2] * z[6];
550 s8 = s7 + x[5] * x[1] * y[6] * z[2] + x[5] * z[1] * x[6] * y[2] +
551 x[1] * x[2] * y[5] * z[6] - x[1] * x[2] * y[6] * z[5] -
552 x[1] * z[1] * x[6] * y[5] - x[1] * y[1] * x[5] * z[6] +
553 x[1] * z[1] * x[5] * y[6] + x[1] * y[1] * x[6] * z[5] -
554 x[5] * x[6] * y[5] * z[2] + x[5] * x[2] * y[5] * z[6] -
555 x[5] * x[2] * y[6] * z[5] + x[5] * x[6] * y[2] * z[5] -
556 2.0 * x[5] * z[1] * x[6] * y[5] - 2.0 * x[5] * x[1] * y[6] * z[5] +
557 2.0 * x[5] * x[1] * y[5] * z[6];
558 s6 = s8 + 2.0 * x[5] * y[1] * x[6] * z[5] +
559 2.0 * x[2] * x[1] * y[6] * z[2] + 2.0 * x[2] * z[1] * x[6] * y[2] -
560 2.0 * x[2] * x[1] * y[2] * z[6] + x[2] * x[5] * y[6] * z[2] +
561 x[2] * x[6] * y[2] * z[5] - x[2] * x[5] * y[2] * z[6] +
562 y[1] * x[2] * x[2] * z[5] - z[1] * x[2] * x[2] * y[5] -
563 2.0 * x[2] * y[1] * x[6] * z[2] - x[2] * x[6] * y[5] * z[2] -
564 2.0 * z[1] * x[2] * x[2] * y[6] + x[2] * x[2] * y[5] * z[6] -
565 x[2] * x[2] * y[6] * z[5] + 2.0 * y[1] * x[2] * x[2] * z[6] +
566 x[2] * z[1] * x[5] * y[2];
567 s8 = s6 - x[2] * x[1] * y[2] * z[5] + x[2] * x[1] * y[5] * z[2] -
568 x[2] * y[1] * x[5] * z[2] + x[6] * y[1] * x[2] * z[5] -
569 x[6] * z[1] * x[2] * y[5] - z[1] * x[6] * x[6] * y[5] +
570 y[1] * x[6] * x[6] * z[5] - y[1] * x[6] * x[6] * z[2] -
571 2.0 * x[6] * x[6] * y[5] * z[2] + 2.0 * x[6] * x[6] * y[2] * z[5] +
572 z[1] * x[6] * x[6] * y[2] - x[6] * x[1] * y[6] * z[5] -
573 x[6] * y[1] * x[5] * z[6] + x[6] * x[1] * y[5] * z[6];
574 s7 = s8 + x[6] * z[1] * x[5] * y[6] - x[6] * z[1] * x[2] * y[6] -
575 x[6] * x[1] * y[2] * z[6] + 2.0 * x[6] * x[5] * y[6] * z[2] +
576 2.0 * x[6] * x[2] * y[5] * z[6] - 2.0 * x[6] * x[2] * y[6] * z[5] -
577 2.0 * x[6] * x[5] * y[2] * z[6] + x[6] * x[1] * y[6] * z[2] +
578 x[6] * y[1] * x[2] * z[6] - x[2] * x[2] * y[3] * z[7] +
579 x[2] * x[2] * y[7] * z[3] - x[2] * z[2] * x[3] * y[7] -
580 x[2] * y[2] * x[7] * z[3] + x[2] * z[2] * x[7] * y[3] +
581 x[2] * y[2] * x[3] * z[7] - x[6] * x[6] * y[3] * z[7];
582 s8 = s7 + x[6] * x[6] * y[7] * z[3] - x[6] * x[2] * y[3] * z[7] +
583 x[6] * x[2] * y[7] * z[3] - x[6] * y[6] * x[7] * z[3] +
584 x[6] * y[6] * x[3] * z[7] - x[6] * z[6] * x[3] * y[7] +
585 x[6] * z[6] * x[7] * y[3] + y[6] * x[2] * x[2] * z[7] -
586 z[6] * x[2] * x[2] * y[7] + 2.0 * x[2] * x[2] * y[6] * z[3] -
587 x[2] * y[6] * x[7] * z[2] - 2.0 * x[2] * y[2] * x[6] * z[3] -
588 2.0 * x[2] * x[2] * y[3] * z[6] + 2.0 * x[2] * y[2] * x[3] * z[6] -
589 x[2] * x[6] * y[2] * z[7];
590 s3 = s8 + x[2] * x[6] * y[7] * z[2] + x[2] * z[6] * x[7] * y[2] +
591 2.0 * x[2] * z[2] * x[6] * y[3] - 2.0 * x[2] * z[2] * x[3] * y[6] -
592 y[2] * x[6] * x[6] * z[3] - 2.0 * x[6] * x[6] * y[2] * z[7] +
593 2.0 * x[6] * x[6] * y[7] * z[2] + z[2] * x[6] * x[6] * y[3] -
594 2.0 * x[6] * y[6] * x[7] * z[2] + x[6] * y[2] * x[3] * z[6] -
595 x[6] * x[2] * y[3] * z[6] + 2.0 * x[6] * z[6] * x[7] * y[2] +
596 2.0 * x[6] * y[6] * x[2] * z[7] - 2.0 * x[6] * z[6] * x[2] * y[7] +
597 x[6] * x[2] * y[6] * z[3] - x[6] * z[2] * x[3] * y[6];
598 s8 = y[1] * x[0] * z[3] + x[1] * y[3] * z[0] - y[0] * x[3] * z[7] -
599 x[1] * y[5] * z[0] - y[0] * x[3] * z[4] - x[1] * y[0] * z[2] +
600 z[1] * x[2] * y[0] - y[1] * x[0] * z[5] - z[1] * x[0] * y[2] -
601 y[1] * x[0] * z[4] + z[1] * x[5] * y[2] + z[0] * x[7] * y[4] +
602 z[0] * x[3] * y[7] + z[1] * x[0] * y[4] - x[1] * y[2] * z[5] +
603 x[2] * y[3] * z[0] + y[1] * x[2] * z[5] - x[2] * y[3] * z[7];
604 s7 = s8 - z[1] * x[2] * y[5] - y[1] * x[3] * z[0] - x[0] * y[7] * z[3] -
605 z[1] * x[0] * y[3] + y[5] * x[4] * z[0] - x[0] * y[4] * z[3] +
606 y[5] * x[7] * z[4] - z[0] * x[4] * y[3] + x[1] * y[0] * z[4] -
607 z[2] * x[3] * y[7] - y[6] * x[7] * z[2] + x[1] * y[5] * z[2] +
608 y[6] * x[7] * z[5] + x[0] * y[7] * z[4] + x[1] * y[2] * z[0] -
609 z[1] * x[4] * y[0] - z[0] * x[4] * y[7] - z[2] * x[0] * y[3];
610 s8 = x[5] * y[0] * z[4] + z[1] * x[0] * y[5] - x[2] * y[0] * z[3] -
611 z[1] * x[5] * y[0] + y[1] * x[5] * z[0] - x[1] * y[0] * z[3] -
612 x[1] * y[4] * z[0] - y[1] * x[5] * z[2] + x[2] * y[7] * z[3] +
613 y[0] * x[4] * z[3] - x[0] * y[4] * z[7] + x[1] * y[0] * z[5] -
614 y[1] * x[6] * z[2] - y[2] * x[6] * z[3] + y[0] * x[7] * z[3] -
615 y[2] * x[7] * z[3] + z[2] * x[7] * y[3] + y[2] * x[0] * z[3];
616 s6 = s8 + y[2] * x[3] * z[7] - y[2] * x[3] * z[0] - x[6] * y[5] * z[2] -
617 y[5] * x[0] * z[4] + z[2] * x[3] * y[0] + x[2] * y[3] * z[1] +
618 x[0] * y[3] * z[7] - x[2] * y[1] * z[3] + y[1] * x[4] * z[0] +
619 y[1] * x[0] * z[2] - z[1] * x[2] * y[6] + y[2] * x[3] * z[6] -
620 y[1] * x[2] * z[0] + z[1] * x[3] * y[0] - x[1] * y[2] * z[6] -
621 x[2] * y[3] * z[6] + x[0] * y[3] * z[4] + z[0] * x[3] * y[4] + s7;
622 s8 = x[5] * y[4] * z[7] + s6 + y[5] * x[6] * z[4] - y[5] * x[4] * z[6] +
623 z[6] * x[5] * y[7] - x[6] * y[2] * z[7] - x[6] * y[7] * z[5] +
624 x[5] * y[6] * z[2] + x[6] * y[5] * z[7] + x[6] * y[7] * z[2] +
625 y[6] * x[7] * z[4] - y[6] * x[4] * z[7] - y[6] * x[7] * z[3] +
626 z[6] * x[7] * y[2] + x[2] * y[5] * z[6] - x[2] * y[6] * z[5] +
627 y[6] * x[2] * z[7] + x[6] * y[2] * z[5];
628 s7 = s8 - x[5] * y[2] * z[6] - z[6] * x[7] * y[5] - z[5] * x[7] * y[4] +
629 z[5] * x[0] * y[4] - y[5] * x[4] * z[7] + y[0] * x[4] * z[7] -
630 z[6] * x[2] * y[7] - x[5] * y[4] * z[0] - x[5] * y[7] * z[4] -
631 y[0] * x[7] * z[4] + y[5] * x[4] * z[1] - x[6] * y[7] * z[4] +
632 x[7] * y[4] * z[3] - x[4] * y[7] * z[3] + x[3] * y[7] * z[4] -
633 x[7] * y[3] * z[4] - x[6] * y[3] * z[7] + x[6] * y[4] * z[7];
634 s8 = -x[3] * y[4] * z[7] + x[4] * y[3] * z[7] - z[6] * x[7] * y[4] -
635 z[1] * x[6] * y[5] + x[6] * y[7] * z[3] - x[1] * y[6] * z[5] -
636 y[1] * x[5] * z[6] + z[5] * x[4] * y[7] - z[5] * x[4] * y[0] +
637 x[1] * y[5] * z[6] - y[6] * x[5] * z[7] - y[2] * x[3] * z[1] +
638 z[1] * x[5] * y[6] - y[5] * x[1] * z[4] + z[6] * x[4] * y[7] +
639 x[5] * y[1] * z[4] - x[5] * y[6] * z[4] + y[6] * x[3] * z[7] -
640 x[5] * y[4] * z[1];
641 s5 = s8 + x[5] * y[4] * z[6] + z[5] * x[1] * y[4] + y[1] * x[6] * z[5] -
642 z[6] * x[3] * y[7] + z[6] * x[7] * y[3] - z[5] * x[6] * y[4] -
643 z[5] * x[4] * y[1] + z[5] * x[4] * y[6] + x[1] * y[6] * z[2] +
644 x[2] * y[6] * z[3] + z[2] * x[6] * y[3] + z[1] * x[6] * y[2] +
645 z[2] * x[3] * y[1] - z[2] * x[1] * y[3] - z[2] * x[3] * y[6] +
646 y[2] * x[1] * z[3] + y[1] * x[2] * z[6] - z[0] * x[7] * y[3] + s7;
647 s4 = 1 / s5;
648 s2 = s3 * s4;
649 const double unknown0 = s1 * s2;
650 s1 = 1.0 / 6.0;
651 s8 = 2.0 * x[1] * y[0] * y[0] * z[4] + x[5] * y[0] * y[0] * z[4] -
652 x[1] * y[4] * y[4] * z[0] + z[1] * x[0] * y[4] * y[4] +
653 x[1] * y[0] * y[0] * z[5] - z[1] * x[5] * y[0] * y[0] -
654 2.0 * z[1] * x[4] * y[0] * y[0] + 2.0 * z[1] * x[3] * y[0] * y[0] +
655 z[2] * x[3] * y[0] * y[0] + y[0] * y[0] * x[7] * z[3] +
656 2.0 * y[0] * y[0] * x[4] * z[3] - 2.0 * x[1] * y[0] * y[0] * z[3] -
657 2.0 * x[5] * y[4] * y[4] * z[0] + 2.0 * z[5] * x[0] * y[4] * y[4] +
658 2.0 * y[4] * y[5] * x[7] * z[4];
659 s7 = s8 - x[3] * y[4] * y[4] * z[7] + x[7] * y[4] * y[4] * z[3] +
660 z[0] * x[3] * y[4] * y[4] - 2.0 * x[0] * y[4] * y[4] * z[7] -
661 y[1] * x[1] * y[4] * z[0] - x[0] * y[4] * y[4] * z[3] +
662 2.0 * z[0] * x[7] * y[4] * y[4] + y[4] * z[6] * x[4] * y[7] -
663 y[0] * y[0] * x[7] * z[4] + y[0] * y[0] * x[4] * z[7] +
664 2.0 * y[4] * z[5] * x[4] * y[7] - 2.0 * y[4] * x[5] * y[7] * z[4] -
665 y[4] * x[6] * y[7] * z[4] - y[4] * y[6] * x[4] * z[7] -
666 2.0 * y[4] * y[5] * x[4] * z[7];
667 s8 = y[4] * y[6] * x[7] * z[4] - y[7] * y[2] * x[7] * z[3] +
668 y[7] * z[2] * x[7] * y[3] + y[7] * y[2] * x[3] * z[7] +
669 2.0 * x[5] * y[4] * y[4] * z[7] - y[7] * x[2] * y[3] * z[7] -
670 y[0] * z[0] * x[4] * y[7] + z[6] * x[7] * y[3] * y[3] -
671 y[0] * x[0] * y[4] * z[7] + y[0] * x[0] * y[7] * z[4] -
672 2.0 * x[2] * y[3] * y[3] * z[7] - z[5] * x[4] * y[0] * y[0] +
673 y[0] * z[0] * x[7] * y[4] - 2.0 * z[6] * x[3] * y[7] * y[7] +
674 z[1] * x[2] * y[0] * y[0];
675 s6 = s8 + y[4] * y[0] * x[4] * z[3] - 2.0 * y[4] * z[0] * x[4] * y[7] +
676 2.0 * y[4] * x[0] * y[7] * z[4] - y[4] * z[0] * x[4] * y[3] -
677 y[4] * x[0] * y[7] * z[3] + y[4] * z[0] * x[3] * y[7] -
678 y[4] * y[0] * x[3] * z[4] + y[0] * x[4] * y[3] * z[7] -
679 y[0] * x[7] * y[3] * z[4] - y[0] * x[3] * y[4] * z[7] +
680 y[0] * x[7] * y[4] * z[3] + x[2] * y[7] * y[7] * z[3] -
681 z[2] * x[3] * y[7] * y[7] - 2.0 * z[2] * x[0] * y[3] * y[3] +
682 2.0 * y[0] * z[1] * x[0] * y[4] + s7;
683 s8 = -2.0 * y[0] * y[1] * x[0] * z[4] - y[0] * y[1] * x[0] * z[5] -
684 y[0] * y[0] * x[3] * z[7] - z[1] * x[0] * y[3] * y[3] -
685 y[0] * x[1] * y[5] * z[0] - 2.0 * z[0] * x[7] * y[3] * y[3] +
686 x[0] * y[3] * y[3] * z[4] + 2.0 * x[0] * y[3] * y[3] * z[7] -
687 z[0] * x[4] * y[3] * y[3] + 2.0 * x[2] * y[3] * y[3] * z[0] +
688 x[1] * y[3] * y[3] * z[0] + 2.0 * y[7] * z[6] * x[7] * y[3] +
689 2.0 * y[7] * y[6] * x[3] * z[7] - 2.0 * y[7] * y[6] * x[7] * z[3] -
690 2.0 * y[7] * x[6] * y[3] * z[7];
691 s7 = s8 + y[4] * x[4] * y[3] * z[7] - y[4] * x[4] * y[7] * z[3] +
692 y[4] * x[3] * y[7] * z[4] - y[4] * x[7] * y[3] * z[4] +
693 2.0 * y[4] * y[0] * x[4] * z[7] - 2.0 * y[4] * y[0] * x[7] * z[4] +
694 2.0 * x[6] * y[7] * y[7] * z[3] + y[4] * x[0] * y[3] * z[4] +
695 y[0] * y[1] * x[5] * z[0] + y[0] * z[1] * x[0] * y[5] -
696 x[2] * y[0] * y[0] * z[3] + x[4] * y[3] * y[3] * z[7] -
697 x[7] * y[3] * y[3] * z[4] - x[5] * y[4] * y[4] * z[1] +
698 y[3] * z[0] * x[3] * y[4];
699 s8 = y[3] * y[0] * x[4] * z[3] + 2.0 * y[3] * y[0] * x[7] * z[3] +
700 2.0 * y[3] * y[2] * x[0] * z[3] - 2.0 * y[3] * y[2] * x[3] * z[0] +
701 2.0 * y[3] * z[2] * x[3] * y[0] + y[3] * z[1] * x[3] * y[0] -
702 2.0 * y[3] * x[2] * y[0] * z[3] - y[3] * x[1] * y[0] * z[3] -
703 y[3] * y[1] * x[3] * z[0] - 2.0 * y[3] * x[0] * y[7] * z[3] -
704 y[3] * x[0] * y[4] * z[3] - 2.0 * y[3] * y[0] * x[3] * z[7] -
705 y[3] * y[0] * x[3] * z[4] + 2.0 * y[3] * z[0] * x[3] * y[7] +
706 y[3] * y[1] * x[0] * z[3] + z[5] * x[1] * y[4] * y[4];
707 s5 = s8 - 2.0 * y[0] * y[0] * x[3] * z[4] -
708 2.0 * y[0] * x[1] * y[4] * z[0] + y[3] * x[7] * y[4] * z[3] -
709 y[3] * x[4] * y[7] * z[3] + y[3] * x[3] * y[7] * z[4] -
710 y[3] * x[3] * y[4] * z[7] + y[3] * x[0] * y[7] * z[4] -
711 y[3] * z[0] * x[4] * y[7] - 2.0 * y[4] * y[5] * x[0] * z[4] + s6 +
712 y[7] * x[0] * y[3] * z[7] - y[7] * z[0] * x[7] * y[3] +
713 y[7] * y[0] * x[7] * z[3] - y[7] * y[0] * x[3] * z[7] +
714 2.0 * y[0] * y[1] * x[4] * z[0] + s7;
715 s8 = -2.0 * y[7] * x[7] * y[3] * z[4] -
716 2.0 * y[7] * x[3] * y[4] * z[7] + 2.0 * y[7] * x[4] * y[3] * z[7] +
717 y[7] * y[0] * x[4] * z[7] - y[7] * y[0] * x[7] * z[4] +
718 2.0 * y[7] * x[7] * y[4] * z[3] - y[7] * x[0] * y[4] * z[7] +
719 y[7] * z[0] * x[7] * y[4] + z[5] * x[4] * y[7] * y[7] +
720 2.0 * z[6] * x[4] * y[7] * y[7] - x[5] * y[7] * y[7] * z[4] -
721 2.0 * x[6] * y[7] * y[7] * z[4] + 2.0 * y[7] * x[6] * y[4] * z[7] -
722 2.0 * y[7] * z[6] * x[7] * y[4] + 2.0 * y[7] * y[6] * x[7] * z[4];
723 s7 = s8 - 2.0 * y[7] * y[6] * x[4] * z[7] - y[7] * z[5] * x[7] * y[4] -
724 y[7] * y[5] * x[4] * z[7] - x[0] * y[7] * y[7] * z[3] +
725 z[0] * x[3] * y[7] * y[7] + y[7] * x[5] * y[4] * z[7] +
726 y[7] * y[5] * x[7] * z[4] - y[4] * x[1] * y[5] * z[0] -
727 x[1] * y[0] * y[0] * z[2] - y[4] * y[5] * x[1] * z[4] -
728 2.0 * y[4] * z[5] * x[4] * y[0] - y[4] * y[1] * x[0] * z[4] +
729 y[4] * y[5] * x[4] * z[1] + y[0] * z[0] * x[3] * y[7] -
730 y[0] * z[1] * x[0] * y[2];
731 s8 = 2.0 * y[0] * x[1] * y[3] * z[0] + y[4] * y[1] * x[4] * z[0] +
732 2.0 * y[0] * y[1] * x[0] * z[3] + y[4] * x[1] * y[0] * z[5] -
733 y[4] * z[1] * x[5] * y[0] + y[4] * z[1] * x[0] * y[5] -
734 y[4] * z[1] * x[4] * y[0] + y[4] * x[1] * y[0] * z[4] -
735 y[4] * z[5] * x[4] * y[1] + x[5] * y[4] * y[4] * z[6] -
736 z[5] * x[6] * y[4] * y[4] + y[4] * x[5] * y[1] * z[4] -
737 y[0] * z[2] * x[0] * y[3] + y[0] * y[5] * x[4] * z[0] +
738 y[0] * x[1] * y[2] * z[0];
739 s6 = s8 - 2.0 * y[0] * z[0] * x[4] * y[3] -
740 2.0 * y[0] * x[0] * y[4] * z[3] - 2.0 * y[0] * z[1] * x[0] * y[3] -
741 y[0] * x[0] * y[7] * z[3] - 2.0 * y[0] * y[1] * x[3] * z[0] +
742 y[0] * x[2] * y[3] * z[0] - y[0] * y[1] * x[2] * z[0] +
743 y[0] * y[1] * x[0] * z[2] - y[0] * x[2] * y[1] * z[3] +
744 y[0] * x[0] * y[3] * z[7] + y[0] * x[2] * y[3] * z[1] -
745 y[0] * y[2] * x[3] * z[0] + y[0] * y[2] * x[0] * z[3] -
746 y[0] * y[5] * x[0] * z[4] - y[4] * y[5] * x[4] * z[6] + s7;
747 s8 = s6 + y[4] * z[6] * x[5] * y[7] - y[4] * x[6] * y[7] * z[5] +
748 y[4] * x[6] * y[5] * z[7] - y[4] * z[6] * x[7] * y[5] -
749 y[4] * x[5] * y[6] * z[4] + y[4] * z[5] * x[4] * y[6] +
750 y[4] * y[5] * x[6] * z[4] - 2.0 * y[1] * y[1] * x[0] * z[5] +
751 2.0 * y[1] * y[1] * x[5] * z[0] - 2.0 * y[2] * y[2] * x[6] * z[3] +
752 x[5] * y[1] * y[1] * z[4] - z[5] * x[4] * y[1] * y[1] -
753 x[6] * y[2] * y[2] * z[7] + z[6] * x[7] * y[2] * y[2];
754 s7 = s8 - x[1] * y[5] * y[5] * z[0] + z[1] * x[0] * y[5] * y[5] +
755 y[1] * y[5] * x[4] * z[1] - y[1] * y[5] * x[1] * z[4] -
756 2.0 * y[2] * z[2] * x[3] * y[6] + 2.0 * y[1] * z[1] * x[0] * y[5] -
757 2.0 * y[1] * z[1] * x[5] * y[0] + 2.0 * y[1] * x[1] * y[0] * z[5] -
758 y[2] * x[2] * y[3] * z[7] - y[2] * z[2] * x[3] * y[7] +
759 y[2] * x[2] * y[7] * z[3] + y[2] * z[2] * x[7] * y[3] -
760 2.0 * y[2] * x[2] * y[3] * z[6] + 2.0 * y[2] * x[2] * y[6] * z[3] +
761 2.0 * y[2] * z[2] * x[6] * y[3] - y[3] * y[2] * x[6] * z[3];
762 s8 = y[3] * y[2] * x[3] * z[6] + y[3] * x[2] * y[6] * z[3] -
763 y[3] * z[2] * x[3] * y[6] - y[2] * y[2] * x[7] * z[3] +
764 2.0 * y[2] * y[2] * x[3] * z[6] + y[2] * y[2] * x[3] * z[7] -
765 2.0 * y[1] * x[1] * y[5] * z[0] - x[2] * y[3] * y[3] * z[6] +
766 z[2] * x[6] * y[3] * y[3] + 2.0 * y[6] * x[2] * y[5] * z[6] +
767 2.0 * y[6] * x[6] * y[2] * z[5] - 2.0 * y[6] * x[5] * y[2] * z[6] +
768 2.0 * y[3] * x[2] * y[7] * z[3] - 2.0 * y[3] * z[2] * x[3] * y[7] -
769 y[0] * z[0] * x[7] * y[3] - y[0] * z[2] * x[1] * y[3];
770 s4 = s8 - y[2] * y[6] * x[7] * z[2] + y[0] * z[2] * x[3] * y[1] +
771 y[1] * z[5] * x[1] * y[4] - y[1] * x[5] * y[4] * z[1] +
772 2.0 * y[0] * z[0] * x[3] * y[4] + 2.0 * y[0] * x[0] * y[3] * z[4] +
773 2.0 * z[2] * x[7] * y[3] * y[3] - 2.0 * z[5] * x[7] * y[4] * y[4] +
774 x[6] * y[4] * y[4] * z[7] - z[6] * x[7] * y[4] * y[4] +
775 y[1] * y[1] * x[0] * z[3] + y[3] * x[6] * y[7] * z[2] -
776 y[3] * z[6] * x[2] * y[7] + 2.0 * y[3] * y[2] * x[3] * z[7] + s5 +
777 s7;
778 s8 = s4 + y[2] * x[6] * y[7] * z[2] - y[2] * y[6] * x[7] * z[3] +
779 y[2] * y[6] * x[2] * z[7] - y[2] * z[6] * x[2] * y[7] -
780 y[2] * x[6] * y[3] * z[7] + y[2] * y[6] * x[3] * z[7] +
781 y[2] * z[6] * x[7] * y[3] - 2.0 * y[3] * y[2] * x[7] * z[3] -
782 x[6] * y[3] * y[3] * z[7] + y[1] * y[1] * x[4] * z[0] -
783 y[1] * y[1] * x[3] * z[0] + x[2] * y[6] * y[6] * z[3] -
784 z[2] * x[3] * y[6] * y[6] - y[1] * y[1] * x[0] * z[4];
785 s7 = s8 + y[5] * x[1] * y[0] * z[5] + y[6] * x[2] * y[7] * z[3] -
786 y[6] * y[2] * x[6] * z[3] + y[6] * y[2] * x[3] * z[6] -
787 y[6] * x[2] * y[3] * z[6] + y[6] * z[2] * x[6] * y[3] -
788 y[5] * y[1] * x[0] * z[5] - y[5] * z[1] * x[5] * y[0] +
789 y[5] * y[1] * x[5] * z[0] - y[6] * z[2] * x[3] * y[7] -
790 y[7] * y[6] * x[7] * z[2] + 2.0 * y[6] * y[6] * x[2] * z[7] +
791 y[6] * y[6] * x[3] * z[7] + x[6] * y[7] * y[7] * z[2] -
792 z[6] * x[2] * y[7] * y[7];
793 s8 = -x[2] * y[1] * y[1] * z[3] + 2.0 * y[1] * y[1] * x[0] * z[2] -
794 2.0 * y[1] * y[1] * x[2] * z[0] + z[2] * x[3] * y[1] * y[1] -
795 z[1] * x[0] * y[2] * y[2] + x[1] * y[2] * y[2] * z[0] +
796 y[2] * y[2] * x[0] * z[3] - y[2] * y[2] * x[3] * z[0] -
797 2.0 * y[2] * y[2] * x[3] * z[1] + y[1] * x[1] * y[3] * z[0] -
798 2.0 * y[6] * y[6] * x[7] * z[2] + 2.0 * y[5] * y[5] * x[4] * z[1] -
799 2.0 * y[5] * y[5] * x[1] * z[4] - y[6] * y[6] * x[7] * z[3] -
800 2.0 * y[1] * x[1] * y[0] * z[2];
801 s6 = s8 + 2.0 * y[1] * z[1] * x[2] * y[0] -
802 2.0 * y[1] * z[1] * x[0] * y[2] + 2.0 * y[1] * x[1] * y[2] * z[0] +
803 y[1] * x[2] * y[3] * z[1] - y[1] * y[2] * x[3] * z[1] -
804 y[1] * z[2] * x[1] * y[3] + y[1] * y[2] * x[1] * z[3] -
805 y[2] * x[1] * y[0] * z[2] + y[2] * z[1] * x[2] * y[0] +
806 y[2] * x[2] * y[3] * z[0] - y[7] * x[6] * y[2] * z[7] +
807 y[7] * z[6] * x[7] * y[2] + y[7] * y[6] * x[2] * z[7] -
808 y[6] * x[6] * y[3] * z[7] + y[6] * x[6] * y[7] * z[3] + s7;
809 s8 = s6 - y[6] * z[6] * x[3] * y[7] + y[6] * z[6] * x[7] * y[3] +
810 2.0 * y[2] * y[2] * x[1] * z[3] + x[2] * y[3] * y[3] * z[1] -
811 z[2] * x[1] * y[3] * y[3] + y[1] * x[1] * y[0] * z[4] +
812 y[1] * z[1] * x[3] * y[0] - y[1] * x[1] * y[0] * z[3] +
813 2.0 * y[5] * x[5] * y[1] * z[4] - 2.0 * y[5] * x[5] * y[4] * z[1] +
814 2.0 * y[5] * z[5] * x[1] * y[4] - 2.0 * y[5] * z[5] * x[4] * y[1] -
815 2.0 * y[6] * x[6] * y[2] * z[7] + 2.0 * y[6] * x[6] * y[7] * z[2];
816 s7 = s8 + 2.0 * y[6] * z[6] * x[7] * y[2] -
817 2.0 * y[6] * z[6] * x[2] * y[7] - y[1] * z[1] * x[4] * y[0] +
818 y[1] * z[1] * x[0] * y[4] - y[1] * z[1] * x[0] * y[3] +
819 2.0 * y[6] * y[6] * x[7] * z[5] + 2.0 * y[5] * y[5] * x[6] * z[4] -
820 2.0 * y[5] * y[5] * x[4] * z[6] + x[6] * y[5] * y[5] * z[7] -
821 y[3] * x[2] * y[1] * z[3] - y[3] * y[2] * x[3] * z[1] +
822 y[3] * z[2] * x[3] * y[1] + y[3] * y[2] * x[1] * z[3] -
823 y[2] * x[2] * y[0] * z[3] + y[2] * z[2] * x[3] * y[0];
824 s8 = s7 + 2.0 * y[2] * x[2] * y[3] * z[1] -
825 2.0 * y[2] * x[2] * y[1] * z[3] + y[2] * y[1] * x[0] * z[2] -
826 y[2] * y[1] * x[2] * z[0] + 2.0 * y[2] * z[2] * x[3] * y[1] -
827 2.0 * y[2] * z[2] * x[1] * y[3] - y[2] * z[2] * x[0] * y[3] +
828 y[5] * z[6] * x[5] * y[7] - y[5] * x[6] * y[7] * z[5] -
829 y[5] * y[6] * x[4] * z[7] - y[5] * y[6] * x[5] * z[7] -
830 2.0 * y[5] * x[5] * y[6] * z[4] + 2.0 * y[5] * x[5] * y[4] * z[6] -
831 2.0 * y[5] * z[5] * x[6] * y[4] + 2.0 * y[5] * z[5] * x[4] * y[6];
832 s5 = s8 - y[1] * y[5] * x[0] * z[4] - z[6] * x[7] * y[5] * y[5] +
833 y[6] * y[6] * x[7] * z[4] - y[6] * y[6] * x[4] * z[7] -
834 2.0 * y[6] * y[6] * x[5] * z[7] - x[5] * y[6] * y[6] * z[4] +
835 z[5] * x[4] * y[6] * y[6] + z[6] * x[5] * y[7] * y[7] -
836 x[6] * y[7] * y[7] * z[5] + y[1] * y[5] * x[4] * z[0] +
837 y[7] * y[6] * x[7] * z[5] + y[6] * y[5] * x[7] * z[4] +
838 y[5] * y[6] * x[7] * z[5] + y[6] * y[5] * x[6] * z[4] -
839 y[6] * y[5] * x[4] * z[6] + 2.0 * y[6] * z[6] * x[5] * y[7];
840 s8 = s5 - 2.0 * y[6] * x[6] * y[7] * z[5] +
841 2.0 * y[6] * x[6] * y[5] * z[7] - 2.0 * y[6] * z[6] * x[7] * y[5] -
842 y[6] * x[5] * y[7] * z[4] - y[6] * x[6] * y[7] * z[4] +
843 y[6] * x[6] * y[4] * z[7] - y[6] * z[6] * x[7] * y[4] +
844 y[6] * z[5] * x[4] * y[7] + y[6] * z[6] * x[4] * y[7] +
845 y[6] * x[5] * y[4] * z[6] - y[6] * z[5] * x[6] * y[4] +
846 y[7] * x[6] * y[5] * z[7] - y[7] * z[6] * x[7] * y[5] -
847 2.0 * y[6] * x[6] * y[5] * z[2];
848 s7 = s8 - y[7] * y[6] * x[5] * z[7] + 2.0 * y[4] * y[5] * x[4] * z[0] +
849 2.0 * x[3] * y[7] * y[7] * z[4] - 2.0 * x[4] * y[7] * y[7] * z[3] -
850 z[0] * x[4] * y[7] * y[7] + x[0] * y[7] * y[7] * z[4] -
851 y[0] * z[5] * x[4] * y[1] + y[0] * x[5] * y[1] * z[4] -
852 y[0] * x[5] * y[4] * z[0] + y[0] * z[5] * x[0] * y[4] -
853 y[5] * y[5] * x[0] * z[4] + y[5] * y[5] * x[4] * z[0] +
854 2.0 * y[1] * y[1] * x[2] * z[5] - 2.0 * y[1] * y[1] * x[5] * z[2] +
855 z[1] * x[5] * y[2] * y[2];
856 s8 = s7 - x[1] * y[2] * y[2] * z[5] - y[5] * z[5] * x[4] * y[0] +
857 y[5] * z[5] * x[0] * y[4] - y[5] * x[5] * y[4] * z[0] -
858 y[2] * x[1] * y[6] * z[5] - y[2] * y[1] * x[5] * z[6] +
859 y[2] * z[1] * x[5] * y[6] + y[2] * y[1] * x[6] * z[5] -
860 y[1] * z[1] * x[6] * y[5] - y[1] * x[1] * y[6] * z[5] +
861 y[1] * x[1] * y[5] * z[6] + y[1] * z[1] * x[5] * y[6] +
862 y[5] * x[5] * y[0] * z[4] + y[2] * y[1] * x[2] * z[5] -
863 y[2] * z[1] * x[2] * y[5];
864 s6 = s8 + y[2] * x[1] * y[5] * z[2] - y[2] * y[1] * x[5] * z[2] -
865 y[1] * y[1] * x[5] * z[6] + y[1] * y[1] * x[6] * z[5] -
866 z[1] * x[2] * y[5] * y[5] + x[1] * y[5] * y[5] * z[2] +
867 2.0 * y[1] * z[1] * x[5] * y[2] - 2.0 * y[1] * x[1] * y[2] * z[5] -
868 2.0 * y[1] * z[1] * x[2] * y[5] + 2.0 * y[1] * x[1] * y[5] * z[2] -
869 y[1] * y[1] * x[6] * z[2] + y[1] * y[1] * x[2] * z[6] -
870 2.0 * y[5] * x[1] * y[6] * z[5] - 2.0 * y[5] * y[1] * x[5] * z[6] +
871 2.0 * y[5] * z[1] * x[5] * y[6] + 2.0 * y[5] * y[1] * x[6] * z[5];
872 s8 = s6 - y[6] * z[1] * x[6] * y[5] - y[6] * y[1] * x[5] * z[6] +
873 y[6] * x[1] * y[5] * z[6] + y[6] * y[1] * x[6] * z[5] -
874 2.0 * z[1] * x[6] * y[5] * y[5] + 2.0 * x[1] * y[5] * y[5] * z[6] -
875 x[1] * y[6] * y[6] * z[5] + z[1] * x[5] * y[6] * y[6] +
876 y[5] * z[1] * x[5] * y[2] - y[5] * x[1] * y[2] * z[5] +
877 y[5] * y[1] * x[2] * z[5] - y[5] * y[1] * x[5] * z[2] -
878 y[6] * z[1] * x[2] * y[5] + y[6] * x[1] * y[5] * z[2];
879 s7 = s8 - y[1] * z[1] * x[2] * y[6] - y[1] * x[1] * y[2] * z[6] +
880 y[1] * x[1] * y[6] * z[2] + y[1] * z[1] * x[6] * y[2] +
881 y[5] * x[5] * y[6] * z[2] - y[5] * x[2] * y[6] * z[5] +
882 y[5] * x[6] * y[2] * z[5] - y[5] * x[5] * y[2] * z[6] -
883 x[6] * y[5] * y[5] * z[2] + x[2] * y[5] * y[5] * z[6] -
884 y[5] * y[5] * x[4] * z[7] + y[5] * y[5] * x[7] * z[4] -
885 y[1] * x[6] * y[5] * z[2] + y[1] * x[2] * y[5] * z[6] -
886 y[2] * x[6] * y[5] * z[2] - 2.0 * y[2] * y[1] * x[6] * z[2];
887 s8 = s7 - 2.0 * y[2] * z[1] * x[2] * y[6] +
888 2.0 * y[2] * x[1] * y[6] * z[2] + 2.0 * y[2] * y[1] * x[2] * z[6] -
889 2.0 * x[1] * y[2] * y[2] * z[6] + 2.0 * z[1] * x[6] * y[2] * y[2] +
890 x[6] * y[2] * y[2] * z[5] - x[5] * y[2] * y[2] * z[6] +
891 2.0 * x[5] * y[6] * y[6] * z[2] - 2.0 * x[2] * y[6] * y[6] * z[5] -
892 z[1] * x[2] * y[6] * y[6] - y[6] * y[1] * x[6] * z[2] -
893 y[6] * x[1] * y[2] * z[6] + y[6] * z[1] * x[6] * y[2] +
894 y[6] * y[1] * x[2] * z[6] + x[1] * y[6] * y[6] * z[2];
895 s3 = s8 + y[2] * x[5] * y[6] * z[2] + y[2] * x[2] * y[5] * z[6] -
896 y[2] * x[2] * y[6] * z[5] + y[5] * z[5] * x[4] * y[7] +
897 y[5] * x[5] * y[4] * z[7] - y[5] * z[5] * x[7] * y[4] -
898 y[5] * x[5] * y[7] * z[4] + 2.0 * y[4] * x[5] * y[0] * z[4] -
899 y[3] * z[6] * x[3] * y[7] + y[3] * y[6] * x[3] * z[7] +
900 y[3] * x[6] * y[7] * z[3] - y[3] * y[6] * x[7] * z[3] -
901 y[2] * y[1] * x[3] * z[0] - y[2] * z[1] * x[0] * y[3] +
902 y[2] * y[1] * x[0] * z[3] + y[2] * x[1] * y[3] * z[0];
903 s8 = y[1] * x[0] * z[3] + x[1] * y[3] * z[0] - y[0] * x[3] * z[7] -
904 x[1] * y[5] * z[0] - y[0] * x[3] * z[4] - x[1] * y[0] * z[2] +
905 z[1] * x[2] * y[0] - y[1] * x[0] * z[5] - z[1] * x[0] * y[2] -
906 y[1] * x[0] * z[4] + z[1] * x[5] * y[2] + z[0] * x[7] * y[4] +
907 z[0] * x[3] * y[7] + z[1] * x[0] * y[4] - x[1] * y[2] * z[5] +
908 x[2] * y[3] * z[0] + y[1] * x[2] * z[5] - x[2] * y[3] * z[7];
909 s7 = s8 - z[1] * x[2] * y[5] - y[1] * x[3] * z[0] - x[0] * y[7] * z[3] -
910 z[1] * x[0] * y[3] + y[5] * x[4] * z[0] - x[0] * y[4] * z[3] +
911 y[5] * x[7] * z[4] - z[0] * x[4] * y[3] + x[1] * y[0] * z[4] -
912 z[2] * x[3] * y[7] - y[6] * x[7] * z[2] + x[1] * y[5] * z[2] +
913 y[6] * x[7] * z[5] + x[0] * y[7] * z[4] + x[1] * y[2] * z[0] -
914 z[1] * x[4] * y[0] - z[0] * x[4] * y[7] - z[2] * x[0] * y[3];
915 s8 = x[5] * y[0] * z[4] + z[1] * x[0] * y[5] - x[2] * y[0] * z[3] -
916 z[1] * x[5] * y[0] + y[1] * x[5] * z[0] - x[1] * y[0] * z[3] -
917 x[1] * y[4] * z[0] - y[1] * x[5] * z[2] + x[2] * y[7] * z[3] +
918 y[0] * x[4] * z[3] - x[0] * y[4] * z[7] + x[1] * y[0] * z[5] -
919 y[1] * x[6] * z[2] - y[2] * x[6] * z[3] + y[0] * x[7] * z[3] -
920 y[2] * x[7] * z[3] + z[2] * x[7] * y[3] + y[2] * x[0] * z[3];
921 s6 = s8 + y[2] * x[3] * z[7] - y[2] * x[3] * z[0] - x[6] * y[5] * z[2] -
922 y[5] * x[0] * z[4] + z[2] * x[3] * y[0] + x[2] * y[3] * z[1] +
923 x[0] * y[3] * z[7] - x[2] * y[1] * z[3] + y[1] * x[4] * z[0] +
924 y[1] * x[0] * z[2] - z[1] * x[2] * y[6] + y[2] * x[3] * z[6] -
925 y[1] * x[2] * z[0] + z[1] * x[3] * y[0] - x[1] * y[2] * z[6] -
926 x[2] * y[3] * z[6] + x[0] * y[3] * z[4] + z[0] * x[3] * y[4] + s7;
927 s8 = x[5] * y[4] * z[7] + s6 + y[5] * x[6] * z[4] - y[5] * x[4] * z[6] +
928 z[6] * x[5] * y[7] - x[6] * y[2] * z[7] - x[6] * y[7] * z[5] +
929 x[5] * y[6] * z[2] + x[6] * y[5] * z[7] + x[6] * y[7] * z[2] +
930 y[6] * x[7] * z[4] - y[6] * x[4] * z[7] - y[6] * x[7] * z[3] +
931 z[6] * x[7] * y[2] + x[2] * y[5] * z[6] - x[2] * y[6] * z[5] +
932 y[6] * x[2] * z[7] + x[6] * y[2] * z[5];
933 s7 = s8 - x[5] * y[2] * z[6] - z[6] * x[7] * y[5] - z[5] * x[7] * y[4] +
934 z[5] * x[0] * y[4] - y[5] * x[4] * z[7] + y[0] * x[4] * z[7] -
935 z[6] * x[2] * y[7] - x[5] * y[4] * z[0] - x[5] * y[7] * z[4] -
936 y[0] * x[7] * z[4] + y[5] * x[4] * z[1] - x[6] * y[7] * z[4] +
937 x[7] * y[4] * z[3] - x[4] * y[7] * z[3] + x[3] * y[7] * z[4] -
938 x[7] * y[3] * z[4] - x[6] * y[3] * z[7] + x[6] * y[4] * z[7];
939 s8 = -x[3] * y[4] * z[7] + x[4] * y[3] * z[7] - z[6] * x[7] * y[4] -
940 z[1] * x[6] * y[5] + x[6] * y[7] * z[3] - x[1] * y[6] * z[5] -
941 y[1] * x[5] * z[6] + z[5] * x[4] * y[7] - z[5] * x[4] * y[0] +
942 x[1] * y[5] * z[6] - y[6] * x[5] * z[7] - y[2] * x[3] * z[1] +
943 z[1] * x[5] * y[6] - y[5] * x[1] * z[4] + z[6] * x[4] * y[7] +
944 x[5] * y[1] * z[4] - x[5] * y[6] * z[4] + y[6] * x[3] * z[7] -
945 x[5] * y[4] * z[1];
946 s5 = s8 + x[5] * y[4] * z[6] + z[5] * x[1] * y[4] + y[1] * x[6] * z[5] -
947 z[6] * x[3] * y[7] + z[6] * x[7] * y[3] - z[5] * x[6] * y[4] -
948 z[5] * x[4] * y[1] + z[5] * x[4] * y[6] + x[1] * y[6] * z[2] +
949 x[2] * y[6] * z[3] + z[2] * x[6] * y[3] + z[1] * x[6] * y[2] +
950 z[2] * x[3] * y[1] - z[2] * x[1] * y[3] - z[2] * x[3] * y[6] +
951 y[2] * x[1] * z[3] + y[1] * x[2] * z[6] - z[0] * x[7] * y[3] + s7;
952 s4 = 1 / s5;
953 s2 = s3 * s4;
954 const double unknown1 = s1 * s2;
955 s1 = 1.0 / 6.0;
956 s8 = -z[2] * x[1] * y[2] * z[5] + z[2] * y[1] * x[2] * z[5] -
957 z[2] * z[1] * x[2] * y[5] + z[2] * z[1] * x[5] * y[2] +
958 2.0 * y[5] * x[7] * z[4] * z[4] - y[1] * x[2] * z[0] * z[0] +
959 x[0] * y[3] * z[7] * z[7] - 2.0 * z[5] * z[5] * x[4] * y[1] +
960 2.0 * z[5] * z[5] * x[1] * y[4] + z[5] * z[5] * x[0] * y[4] -
961 2.0 * z[2] * z[2] * x[1] * y[3] + 2.0 * z[2] * z[2] * x[3] * y[1] -
962 x[0] * y[4] * z[7] * z[7] - y[0] * x[3] * z[7] * z[7] +
963 x[1] * y[0] * z[5] * z[5];
964 s7 = s8 - y[1] * x[0] * z[5] * z[5] + z[1] * y[1] * x[2] * z[6] +
965 y[1] * x[0] * z[2] * z[2] + z[2] * z[2] * x[3] * y[0] -
966 z[2] * z[2] * x[0] * y[3] - x[1] * y[0] * z[2] * z[2] +
967 2.0 * z[5] * z[5] * x[4] * y[6] - 2.0 * z[5] * z[5] * x[6] * y[4] -
968 z[5] * z[5] * x[7] * y[4] - x[6] * y[7] * z[5] * z[5] +
969 2.0 * z[2] * y[1] * x[2] * z[6] - 2.0 * z[2] * x[1] * y[2] * z[6] +
970 2.0 * z[2] * z[1] * x[6] * y[2] - y[6] * x[5] * z[7] * z[7] +
971 2.0 * x[6] * y[4] * z[7] * z[7];
972 s8 = -2.0 * y[6] * x[4] * z[7] * z[7] + x[6] * y[5] * z[7] * z[7] -
973 2.0 * z[2] * z[1] * x[2] * y[6] + z[4] * y[6] * x[7] * z[5] +
974 x[5] * y[4] * z[6] * z[6] + z[6] * z[6] * x[4] * y[7] -
975 z[6] * z[6] * x[7] * y[4] - 2.0 * z[6] * z[6] * x[7] * y[5] +
976 2.0 * z[6] * z[6] * x[5] * y[7] - y[5] * x[4] * z[6] * z[6] +
977 2.0 * z[0] * z[0] * x[3] * y[4] - x[6] * y[5] * z[2] * z[2] +
978 z[1] * z[1] * x[5] * y[6] - z[1] * z[1] * x[6] * y[5] -
979 z[5] * z[5] * x[4] * y[0];
980 s6 = s8 + 2.0 * x[1] * y[3] * z[0] * z[0] +
981 2.0 * x[1] * y[6] * z[2] * z[2] - 2.0 * y[1] * x[6] * z[2] * z[2] -
982 y[1] * x[5] * z[2] * z[2] - z[1] * z[1] * x[2] * y[6] -
983 2.0 * z[1] * z[1] * x[2] * y[5] + 2.0 * z[1] * z[1] * x[5] * y[2] +
984 z[1] * y[1] * x[6] * z[5] + y[1] * x[2] * z[5] * z[5] +
985 z[2] * z[1] * x[2] * y[0] + z[1] * x[1] * y[5] * z[6] -
986 z[1] * x[1] * y[6] * z[5] - z[1] * y[1] * x[5] * z[6] -
987 z[1] * x[2] * y[6] * z[5] + z[1] * x[6] * y[2] * z[5] + s7;
988 s8 = -x[1] * y[2] * z[5] * z[5] + z[1] * x[5] * y[6] * z[2] -
989 2.0 * z[2] * z[2] * x[3] * y[6] + 2.0 * z[2] * z[2] * x[6] * y[3] +
990 z[2] * z[2] * x[7] * y[3] - z[2] * z[2] * x[3] * y[7] -
991 z[1] * x[6] * y[5] * z[2] + 2.0 * z[1] * x[1] * y[5] * z[2] -
992 2.0 * x[3] * y[4] * z[7] * z[7] + 2.0 * x[4] * y[3] * z[7] * z[7] +
993 x[5] * y[6] * z[2] * z[2] + y[1] * x[2] * z[6] * z[6] +
994 y[0] * x[4] * z[7] * z[7] + z[2] * x[2] * y[3] * z[0] -
995 x[1] * y[2] * z[6] * z[6];
996 s7 = s8 - z[7] * z[2] * x[3] * y[7] + x[2] * y[6] * z[3] * z[3] -
997 y[2] * x[6] * z[3] * z[3] - z[6] * x[2] * y[3] * z[7] -
998 z[2] * z[1] * x[0] * y[2] + z[6] * z[2] * x[6] * y[3] -
999 z[6] * z[2] * x[3] * y[6] + z[6] * x[2] * y[6] * z[3] +
1000 z[2] * x[1] * y[2] * z[0] + z[6] * y[2] * x[3] * z[7] -
1001 z[4] * z[5] * x[6] * y[4] + z[4] * z[5] * x[4] * y[6] -
1002 z[4] * y[6] * x[5] * z[7] + z[4] * z[6] * x[4] * y[7] +
1003 z[4] * x[5] * y[4] * z[6];
1004 s8 = -z[6] * y[2] * x[6] * z[3] - z[4] * y[5] * x[4] * z[6] -
1005 z[2] * y[1] * x[5] * z[6] + z[2] * x[1] * y[5] * z[6] +
1006 z[4] * x[6] * y[4] * z[7] + 2.0 * z[4] * z[5] * x[4] * y[7] -
1007 z[4] * z[6] * x[7] * y[4] + x[6] * y[7] * z[3] * z[3] -
1008 2.0 * z[4] * z[5] * x[7] * y[4] - 2.0 * z[4] * y[5] * x[4] * z[7] -
1009 z[4] * y[6] * x[4] * z[7] + z[4] * x[6] * y[5] * z[7] -
1010 z[4] * x[6] * y[7] * z[5] + 2.0 * z[4] * x[5] * y[4] * z[7] +
1011 z[2] * x[2] * y[5] * z[6] - z[2] * x[2] * y[6] * z[5];
1012 s5 = s8 + z[2] * x[6] * y[2] * z[5] - z[2] * x[5] * y[2] * z[6] -
1013 z[2] * x[2] * y[3] * z[7] - x[2] * y[3] * z[7] * z[7] +
1014 2.0 * z[2] * x[2] * y[3] * z[1] - z[2] * y[2] * x[3] * z[0] +
1015 z[2] * y[2] * x[0] * z[3] - z[2] * x[2] * y[0] * z[3] -
1016 z[7] * y[2] * x[7] * z[3] + z[7] * z[2] * x[7] * y[3] +
1017 z[7] * x[2] * y[7] * z[3] + z[6] * y[1] * x[2] * z[5] -
1018 z[6] * x[1] * y[2] * z[5] + z[5] * x[1] * y[5] * z[2] + s6 + s7;
1019 s8 = z[5] * z[1] * x[5] * y[2] - z[5] * z[1] * x[2] * y[5] -
1020 y[6] * x[7] * z[2] * z[2] + 2.0 * z[2] * x[2] * y[6] * z[3] -
1021 2.0 * z[2] * x[2] * y[3] * z[6] + 2.0 * z[2] * y[2] * x[3] * z[6] +
1022 y[2] * x[3] * z[6] * z[6] + y[6] * x[7] * z[5] * z[5] +
1023 z[2] * y[2] * x[3] * z[7] - z[2] * y[2] * x[7] * z[3] -
1024 2.0 * z[2] * y[2] * x[6] * z[3] + z[2] * x[2] * y[7] * z[3] +
1025 x[6] * y[2] * z[5] * z[5] - 2.0 * z[2] * x[2] * y[1] * z[3] -
1026 x[2] * y[6] * z[5] * z[5];
1027 s7 = s8 - y[1] * x[5] * z[6] * z[6] + z[6] * x[1] * y[6] * z[2] -
1028 z[3] * z[2] * x[3] * y[6] + z[6] * z[1] * x[6] * y[2] -
1029 z[6] * z[1] * x[2] * y[6] - z[6] * y[1] * x[6] * z[2] -
1030 2.0 * x[5] * y[2] * z[6] * z[6] + z[4] * z[1] * x[0] * y[4] -
1031 z[3] * x[2] * y[3] * z[6] - z[5] * y[1] * x[5] * z[2] +
1032 z[3] * y[2] * x[3] * z[6] + 2.0 * x[2] * y[5] * z[6] * z[6] -
1033 z[5] * x[1] * y[5] * z[0] + y[2] * x[3] * z[7] * z[7] -
1034 x[2] * y[3] * z[6] * z[6];
1035 s8 = z[5] * y[5] * x[4] * z[0] + z[3] * z[2] * x[6] * y[3] +
1036 x[1] * y[5] * z[6] * z[6] + z[5] * y[5] * x[7] * z[4] -
1037 z[1] * x[1] * y[2] * z[6] + z[1] * x[1] * y[6] * z[2] +
1038 2.0 * z[6] * y[6] * x[7] * z[5] - z[7] * y[6] * x[7] * z[2] -
1039 z[3] * y[6] * x[7] * z[2] + x[6] * y[7] * z[2] * z[2] -
1040 2.0 * z[6] * y[6] * x[7] * z[2] - 2.0 * x[6] * y[3] * z[7] * z[7] -
1041 x[6] * y[2] * z[7] * z[7] - z[5] * x[6] * y[5] * z[2] +
1042 y[6] * x[2] * z[7] * z[7];
1043 s6 = s8 + 2.0 * y[6] * x[3] * z[7] * z[7] + z[6] * z[6] * x[7] * y[3] -
1044 y[6] * x[7] * z[3] * z[3] + z[5] * x[5] * y[0] * z[4] +
1045 2.0 * z[6] * z[6] * x[7] * y[2] - 2.0 * z[6] * z[6] * x[2] * y[7] -
1046 z[6] * z[6] * x[3] * y[7] + z[7] * y[6] * x[7] * z[5] +
1047 z[7] * y[5] * x[7] * z[4] - 2.0 * z[7] * x[7] * y[3] * z[4] +
1048 2.0 * z[7] * x[3] * y[7] * z[4] - 2.0 * z[7] * x[4] * y[7] * z[3] +
1049 2.0 * z[7] * x[7] * y[4] * z[3] - z[7] * y[0] * x[7] * z[4] -
1050 2.0 * z[7] * z[6] * x[3] * y[7] + s7;
1051 s8 = s6 + 2.0 * z[7] * z[6] * x[7] * y[3] +
1052 2.0 * z[7] * x[6] * y[7] * z[3] + z[7] * x[6] * y[7] * z[2] -
1053 2.0 * z[7] * y[6] * x[7] * z[3] + z[7] * z[6] * x[7] * y[2] -
1054 z[7] * z[6] * x[2] * y[7] + z[5] * y[1] * x[5] * z[0] -
1055 z[5] * z[1] * x[5] * y[0] + 2.0 * y[1] * x[6] * z[5] * z[5] -
1056 2.0 * x[1] * y[6] * z[5] * z[5] + z[5] * z[1] * x[0] * y[5] +
1057 z[6] * y[6] * x[3] * z[7] + 2.0 * z[6] * x[6] * y[7] * z[2] -
1058 z[6] * y[6] * x[7] * z[3];
1059 s7 = s8 + 2.0 * z[6] * y[6] * x[2] * z[7] - z[6] * x[6] * y[3] * z[7] +
1060 z[6] * x[6] * y[7] * z[3] - 2.0 * z[6] * x[6] * y[2] * z[7] -
1061 2.0 * z[1] * y[1] * x[5] * z[2] - z[1] * y[1] * x[6] * z[2] -
1062 z[7] * z[0] * x[7] * y[3] - 2.0 * z[6] * x[6] * y[5] * z[2] -
1063 z[2] * z[6] * x[3] * y[7] + z[2] * x[6] * y[7] * z[3] -
1064 z[2] * z[6] * x[2] * y[7] + y[5] * x[6] * z[4] * z[4] +
1065 z[2] * y[6] * x[2] * z[7] + y[6] * x[7] * z[4] * z[4] +
1066 z[2] * z[6] * x[7] * y[2] - 2.0 * x[5] * y[7] * z[4] * z[4];
1067 s8 = -x[6] * y[7] * z[4] * z[4] - z[5] * y[5] * x[0] * z[4] -
1068 z[2] * x[6] * y[2] * z[7] - x[5] * y[6] * z[4] * z[4] -
1069 2.0 * z[5] * y[1] * x[5] * z[6] + 2.0 * z[5] * z[1] * x[5] * y[6] +
1070 2.0 * z[5] * x[1] * y[5] * z[6] - 2.0 * z[5] * z[1] * x[6] * y[5] -
1071 z[5] * x[5] * y[2] * z[6] + z[5] * x[5] * y[6] * z[2] +
1072 z[5] * x[2] * y[5] * z[6] + z[5] * z[5] * x[4] * y[7] -
1073 y[5] * x[4] * z[7] * z[7] + x[5] * y[4] * z[7] * z[7] +
1074 z[6] * z[1] * x[5] * y[6] + z[6] * y[1] * x[6] * z[5];
1075 s4 = s8 - z[6] * z[1] * x[6] * y[5] - z[6] * x[1] * y[6] * z[5] +
1076 z[2] * z[6] * x[7] * y[3] + 2.0 * z[6] * x[6] * y[2] * z[5] +
1077 2.0 * z[6] * x[5] * y[6] * z[2] - 2.0 * z[6] * x[2] * y[6] * z[5] +
1078 z[7] * z[0] * x[3] * y[7] + z[7] * z[0] * x[7] * y[4] +
1079 z[3] * z[6] * x[7] * y[3] - z[3] * z[6] * x[3] * y[7] -
1080 z[3] * x[6] * y[3] * z[7] + z[3] * y[6] * x[2] * z[7] -
1081 z[3] * x[6] * y[2] * z[7] + z[5] * x[5] * y[4] * z[7] + s5 + s7;
1082 s8 = s4 + z[3] * y[6] * x[3] * z[7] - z[7] * x[0] * y[7] * z[3] +
1083 z[6] * x[5] * y[4] * z[7] + z[7] * y[0] * x[7] * z[3] +
1084 z[5] * z[6] * x[4] * y[7] - 2.0 * z[5] * x[5] * y[6] * z[4] +
1085 2.0 * z[5] * x[5] * y[4] * z[6] - z[5] * x[5] * y[7] * z[4] -
1086 z[5] * y[6] * x[5] * z[7] - z[5] * z[6] * x[7] * y[4] -
1087 z[7] * z[0] * x[4] * y[7] - z[5] * z[6] * x[7] * y[5] -
1088 z[5] * y[5] * x[4] * z[7] + z[7] * x[0] * y[7] * z[4];
1089 s7 = s8 - 2.0 * z[5] * y[5] * x[4] * z[6] + z[5] * z[6] * x[5] * y[7] +
1090 z[5] * x[6] * y[5] * z[7] + 2.0 * z[5] * y[5] * x[6] * z[4] +
1091 z[6] * z[5] * x[4] * y[6] - z[6] * x[5] * y[6] * z[4] -
1092 z[6] * z[5] * x[6] * y[4] - z[6] * x[6] * y[7] * z[4] -
1093 2.0 * z[6] * y[6] * x[5] * z[7] + z[6] * x[6] * y[4] * z[7] -
1094 z[6] * y[5] * x[4] * z[7] - z[6] * y[6] * x[4] * z[7] +
1095 z[6] * y[6] * x[7] * z[4] + z[6] * y[5] * x[6] * z[4] +
1096 2.0 * z[6] * x[6] * y[5] * z[7];
1097 s8 = -2.0 * z[6] * x[6] * y[7] * z[5] - z[2] * y[1] * x[2] * z[0] +
1098 2.0 * z[7] * z[6] * x[4] * y[7] - 2.0 * z[7] * x[6] * y[7] * z[4] -
1099 2.0 * z[7] * z[6] * x[7] * y[4] + z[7] * z[5] * x[4] * y[7] -
1100 z[7] * z[5] * x[7] * y[4] - z[7] * x[5] * y[7] * z[4] +
1101 2.0 * z[7] * y[6] * x[7] * z[4] - z[7] * z[6] * x[7] * y[5] +
1102 z[7] * z[6] * x[5] * y[7] - z[7] * x[6] * y[7] * z[5] +
1103 z[1] * z[1] * x[6] * y[2] + s7 + x[1] * y[5] * z[2] * z[2];
1104 s6 = s8 + 2.0 * z[2] * y[2] * x[1] * z[3] -
1105 2.0 * z[2] * y[2] * x[3] * z[1] - 2.0 * x[1] * y[4] * z[0] * z[0] +
1106 2.0 * y[1] * x[4] * z[0] * z[0] + 2.0 * x[2] * y[7] * z[3] * z[3] -
1107 2.0 * y[2] * x[7] * z[3] * z[3] - x[1] * y[5] * z[0] * z[0] +
1108 z[0] * z[0] * x[7] * y[4] + z[0] * z[0] * x[3] * y[7] +
1109 x[2] * y[3] * z[0] * z[0] - 2.0 * y[1] * x[3] * z[0] * z[0] +
1110 y[5] * x[4] * z[0] * z[0] - 2.0 * z[0] * z[0] * x[4] * y[3] +
1111 x[1] * y[2] * z[0] * z[0] - z[0] * z[0] * x[4] * y[7] +
1112 y[1] * x[5] * z[0] * z[0];
1113 s8 = s6 - y[2] * x[3] * z[0] * z[0] + y[1] * x[0] * z[3] * z[3] -
1114 2.0 * x[0] * y[7] * z[3] * z[3] - x[0] * y[4] * z[3] * z[3] -
1115 2.0 * x[2] * y[0] * z[3] * z[3] - x[1] * y[0] * z[3] * z[3] +
1116 y[0] * x[4] * z[3] * z[3] - 2.0 * z[0] * y[1] * x[0] * z[4] +
1117 2.0 * z[0] * z[1] * x[0] * y[4] + 2.0 * z[0] * x[1] * y[0] * z[4] -
1118 2.0 * z[0] * z[1] * x[4] * y[0] - 2.0 * z[3] * x[2] * y[3] * z[7] -
1119 2.0 * z[3] * z[2] * x[3] * y[7] + 2.0 * z[3] * z[2] * x[7] * y[3];
1120 s7 = s8 + 2.0 * z[3] * y[2] * x[3] * z[7] +
1121 2.0 * z[5] * y[5] * x[4] * z[1] + 2.0 * z[0] * y[1] * x[0] * z[3] -
1122 z[0] * y[0] * x[3] * z[7] - 2.0 * z[0] * y[0] * x[3] * z[4] -
1123 z[0] * x[1] * y[0] * z[2] + z[0] * z[1] * x[2] * y[0] -
1124 z[0] * y[1] * x[0] * z[5] - z[0] * z[1] * x[0] * y[2] -
1125 z[0] * x[0] * y[7] * z[3] - 2.0 * z[0] * z[1] * x[0] * y[3] -
1126 z[5] * x[5] * y[4] * z[0] - 2.0 * z[0] * x[0] * y[4] * z[3] +
1127 z[0] * x[0] * y[7] * z[4] - z[0] * z[2] * x[0] * y[3];
1128 s8 = s7 + z[0] * x[5] * y[0] * z[4] + z[0] * z[1] * x[0] * y[5] -
1129 z[0] * x[2] * y[0] * z[3] - z[0] * z[1] * x[5] * y[0] -
1130 2.0 * z[0] * x[1] * y[0] * z[3] + 2.0 * z[0] * y[0] * x[4] * z[3] -
1131 z[0] * x[0] * y[4] * z[7] + z[0] * x[1] * y[0] * z[5] +
1132 z[0] * y[0] * x[7] * z[3] + z[0] * y[2] * x[0] * z[3] -
1133 z[0] * y[5] * x[0] * z[4] + z[0] * z[2] * x[3] * y[0] +
1134 z[0] * x[2] * y[3] * z[1] + z[0] * x[0] * y[3] * z[7] -
1135 z[0] * x[2] * y[1] * z[3];
1136 s5 = s8 + z[0] * y[1] * x[0] * z[2] + z[3] * x[1] * y[3] * z[0] -
1137 2.0 * z[3] * y[0] * x[3] * z[7] - z[3] * y[0] * x[3] * z[4] -
1138 z[3] * x[1] * y[0] * z[2] + z[3] * z[0] * x[7] * y[4] +
1139 2.0 * z[3] * z[0] * x[3] * y[7] + 2.0 * z[3] * x[2] * y[3] * z[0] -
1140 z[3] * y[1] * x[3] * z[0] - z[3] * z[1] * x[0] * y[3] -
1141 z[3] * z[0] * x[4] * y[3] + z[3] * x[1] * y[2] * z[0] -
1142 z[3] * z[0] * x[4] * y[7] - 2.0 * z[3] * z[2] * x[0] * y[3] -
1143 z[3] * x[0] * y[4] * z[7] - 2.0 * z[3] * y[2] * x[3] * z[0];
1144 s8 = s5 + 2.0 * z[3] * z[2] * x[3] * y[0] + z[3] * x[2] * y[3] * z[1] +
1145 2.0 * z[3] * x[0] * y[3] * z[7] + z[3] * y[1] * x[0] * z[2] -
1146 z[4] * y[0] * x[3] * z[7] - z[4] * x[1] * y[5] * z[0] -
1147 z[4] * y[1] * x[0] * z[5] + 2.0 * z[4] * z[0] * x[7] * y[4] +
1148 z[4] * z[0] * x[3] * y[7] + 2.0 * z[4] * y[5] * x[4] * z[0] +
1149 2.0 * y[0] * x[7] * z[3] * z[3] + 2.0 * y[2] * x[0] * z[3] * z[3] -
1150 x[2] * y[1] * z[3] * z[3] - y[0] * x[3] * z[4] * z[4];
1151 s7 = s8 - y[1] * x[0] * z[4] * z[4] + x[1] * y[0] * z[4] * z[4] +
1152 2.0 * x[0] * y[7] * z[4] * z[4] + 2.0 * x[5] * y[0] * z[4] * z[4] -
1153 2.0 * y[5] * x[0] * z[4] * z[4] + 2.0 * z[1] * z[1] * x[2] * y[0] -
1154 2.0 * z[1] * z[1] * x[0] * y[2] + z[1] * z[1] * x[0] * y[4] -
1155 z[1] * z[1] * x[0] * y[3] - z[1] * z[1] * x[4] * y[0] +
1156 2.0 * z[1] * z[1] * x[0] * y[5] - 2.0 * z[1] * z[1] * x[5] * y[0] +
1157 x[2] * y[3] * z[1] * z[1] - x[5] * y[4] * z[0] * z[0] -
1158 z[0] * z[0] * x[7] * y[3];
1159 s8 = s7 + x[7] * y[4] * z[3] * z[3] - x[4] * y[7] * z[3] * z[3] +
1160 y[2] * x[1] * z[3] * z[3] + x[0] * y[3] * z[4] * z[4] -
1161 2.0 * y[0] * x[7] * z[4] * z[4] + x[3] * y[7] * z[4] * z[4] -
1162 x[7] * y[3] * z[4] * z[4] - y[5] * x[1] * z[4] * z[4] +
1163 x[5] * y[1] * z[4] * z[4] + z[1] * z[1] * x[3] * y[0] +
1164 y[5] * x[4] * z[1] * z[1] - y[2] * x[3] * z[1] * z[1] -
1165 x[5] * y[4] * z[1] * z[1] - z[4] * x[0] * y[4] * z[3] -
1166 z[4] * z[0] * x[4] * y[3];
1167 s6 = s8 - z[4] * z[1] * x[4] * y[0] - 2.0 * z[4] * z[0] * x[4] * y[7] +
1168 z[4] * y[1] * x[5] * z[0] - 2.0 * z[5] * x[5] * y[4] * z[1] -
1169 z[4] * x[1] * y[4] * z[0] + z[4] * y[0] * x[4] * z[3] -
1170 2.0 * z[4] * x[0] * y[4] * z[7] + z[4] * x[1] * y[0] * z[5] -
1171 2.0 * z[1] * x[1] * y[2] * z[5] + z[4] * x[0] * y[3] * z[7] +
1172 2.0 * z[5] * x[5] * y[1] * z[4] + z[4] * y[1] * x[4] * z[0] +
1173 z[1] * y[1] * x[0] * z[3] + z[1] * x[1] * y[3] * z[0] -
1174 2.0 * z[1] * x[1] * y[5] * z[0] - 2.0 * z[1] * x[1] * y[0] * z[2];
1175 s8 = s6 - 2.0 * z[1] * y[1] * x[0] * z[5] - z[1] * y[1] * x[0] * z[4] +
1176 2.0 * z[1] * y[1] * x[2] * z[5] - z[1] * y[1] * x[3] * z[0] -
1177 2.0 * z[5] * y[5] * x[1] * z[4] + z[1] * y[5] * x[4] * z[0] +
1178 z[1] * x[1] * y[0] * z[4] + 2.0 * z[1] * x[1] * y[2] * z[0] -
1179 z[1] * z[2] * x[0] * y[3] + 2.0 * z[1] * y[1] * x[5] * z[0] -
1180 z[1] * x[1] * y[0] * z[3] - z[1] * x[1] * y[4] * z[0] +
1181 2.0 * z[1] * x[1] * y[0] * z[5] - z[1] * y[2] * x[3] * z[0];
1182 s7 = s8 + z[1] * z[2] * x[3] * y[0] - z[1] * x[2] * y[1] * z[3] +
1183 z[1] * y[1] * x[4] * z[0] + 2.0 * z[1] * y[1] * x[0] * z[2] +
1184 2.0 * z[0] * z[1] * x[3] * y[0] + 2.0 * z[0] * x[0] * y[3] * z[4] +
1185 z[0] * z[5] * x[0] * y[4] + z[0] * y[0] * x[4] * z[7] -
1186 z[0] * y[0] * x[7] * z[4] - z[0] * x[7] * y[3] * z[4] -
1187 z[0] * z[5] * x[4] * y[0] - z[0] * x[5] * y[4] * z[1] +
1188 z[3] * z[1] * x[3] * y[0] + z[3] * x[0] * y[3] * z[4] +
1189 z[3] * z[0] * x[3] * y[4] + z[3] * y[0] * x[4] * z[7];
1190 s8 = s7 + z[3] * x[3] * y[7] * z[4] - z[3] * x[7] * y[3] * z[4] -
1191 z[3] * x[3] * y[4] * z[7] + z[3] * x[4] * y[3] * z[7] -
1192 z[3] * y[2] * x[3] * z[1] + z[3] * z[2] * x[3] * y[1] -
1193 z[3] * z[2] * x[1] * y[3] - 2.0 * z[3] * z[0] * x[7] * y[3] +
1194 z[4] * z[0] * x[3] * y[4] + 2.0 * z[4] * z[5] * x[0] * y[4] +
1195 2.0 * z[4] * y[0] * x[4] * z[7] - 2.0 * z[4] * x[5] * y[4] * z[0] +
1196 z[4] * y[5] * x[4] * z[1] + z[4] * x[7] * y[4] * z[3] -
1197 z[4] * x[4] * y[7] * z[3];
1198 s3 = s8 - z[4] * x[3] * y[4] * z[7] + z[4] * x[4] * y[3] * z[7] -
1199 2.0 * z[4] * z[5] * x[4] * y[0] - z[4] * x[5] * y[4] * z[1] +
1200 z[4] * z[5] * x[1] * y[4] - z[4] * z[5] * x[4] * y[1] -
1201 2.0 * z[1] * y[1] * x[2] * z[0] + z[1] * z[5] * x[0] * y[4] -
1202 z[1] * z[5] * x[4] * y[0] - z[1] * y[5] * x[1] * z[4] +
1203 z[1] * x[5] * y[1] * z[4] + z[1] * z[5] * x[1] * y[4] -
1204 z[1] * z[5] * x[4] * y[1] + z[1] * z[2] * x[3] * y[1] -
1205 z[1] * z[2] * x[1] * y[3] + z[1] * y[2] * x[1] * z[3];
1206 s8 = y[1] * x[0] * z[3] + x[1] * y[3] * z[0] - y[0] * x[3] * z[7] -
1207 x[1] * y[5] * z[0] - y[0] * x[3] * z[4] - x[1] * y[0] * z[2] +
1208 z[1] * x[2] * y[0] - y[1] * x[0] * z[5] - z[1] * x[0] * y[2] -
1209 y[1] * x[0] * z[4] + z[1] * x[5] * y[2] + z[0] * x[7] * y[4] +
1210 z[0] * x[3] * y[7] + z[1] * x[0] * y[4] - x[1] * y[2] * z[5] +
1211 x[2] * y[3] * z[0] + y[1] * x[2] * z[5] - x[2] * y[3] * z[7];
1212 s7 = s8 - z[1] * x[2] * y[5] - y[1] * x[3] * z[0] - x[0] * y[7] * z[3] -
1213 z[1] * x[0] * y[3] + y[5] * x[4] * z[0] - x[0] * y[4] * z[3] +
1214 y[5] * x[7] * z[4] - z[0] * x[4] * y[3] + x[1] * y[0] * z[4] -
1215 z[2] * x[3] * y[7] - y[6] * x[7] * z[2] + x[1] * y[5] * z[2] +
1216 y[6] * x[7] * z[5] + x[0] * y[7] * z[4] + x[1] * y[2] * z[0] -
1217 z[1] * x[4] * y[0] - z[0] * x[4] * y[7] - z[2] * x[0] * y[3];
1218 s8 = x[5] * y[0] * z[4] + z[1] * x[0] * y[5] - x[2] * y[0] * z[3] -
1219 z[1] * x[5] * y[0] + y[1] * x[5] * z[0] - x[1] * y[0] * z[3] -
1220 x[1] * y[4] * z[0] - y[1] * x[5] * z[2] + x[2] * y[7] * z[3] +
1221 y[0] * x[4] * z[3] - x[0] * y[4] * z[7] + x[1] * y[0] * z[5] -
1222 y[1] * x[6] * z[2] - y[2] * x[6] * z[3] + y[0] * x[7] * z[3] -
1223 y[2] * x[7] * z[3] + z[2] * x[7] * y[3] + y[2] * x[0] * z[3];
1224 s6 = s8 + y[2] * x[3] * z[7] - y[2] * x[3] * z[0] - x[6] * y[5] * z[2] -
1225 y[5] * x[0] * z[4] + z[2] * x[3] * y[0] + x[2] * y[3] * z[1] +
1226 x[0] * y[3] * z[7] - x[2] * y[1] * z[3] + y[1] * x[4] * z[0] +
1227 y[1] * x[0] * z[2] - z[1] * x[2] * y[6] + y[2] * x[3] * z[6] -
1228 y[1] * x[2] * z[0] + z[1] * x[3] * y[0] - x[1] * y[2] * z[6] -
1229 x[2] * y[3] * z[6] + x[0] * y[3] * z[4] + z[0] * x[3] * y[4] + s7;
1230 s8 = x[5] * y[4] * z[7] + s6 + y[5] * x[6] * z[4] - y[5] * x[4] * z[6] +
1231 z[6] * x[5] * y[7] - x[6] * y[2] * z[7] - x[6] * y[7] * z[5] +
1232 x[5] * y[6] * z[2] + x[6] * y[5] * z[7] + x[6] * y[7] * z[2] +
1233 y[6] * x[7] * z[4] - y[6] * x[4] * z[7] - y[6] * x[7] * z[3] +
1234 z[6] * x[7] * y[2] + x[2] * y[5] * z[6] - x[2] * y[6] * z[5] +
1235 y[6] * x[2] * z[7] + x[6] * y[2] * z[5];
1236 s7 = s8 - x[5] * y[2] * z[6] - z[6] * x[7] * y[5] - z[5] * x[7] * y[4] +
1237 z[5] * x[0] * y[4] - y[5] * x[4] * z[7] + y[0] * x[4] * z[7] -
1238 z[6] * x[2] * y[7] - x[5] * y[4] * z[0] - x[5] * y[7] * z[4] -
1239 y[0] * x[7] * z[4] + y[5] * x[4] * z[1] - x[6] * y[7] * z[4] +
1240 x[7] * y[4] * z[3] - x[4] * y[7] * z[3] + x[3] * y[7] * z[4] -
1241 x[7] * y[3] * z[4] - x[6] * y[3] * z[7] + x[6] * y[4] * z[7];
1242 s8 = -x[3] * y[4] * z[7] + x[4] * y[3] * z[7] - z[6] * x[7] * y[4] -
1243 z[1] * x[6] * y[5] + x[6] * y[7] * z[3] - x[1] * y[6] * z[5] -
1244 y[1] * x[5] * z[6] + z[5] * x[4] * y[7] - z[5] * x[4] * y[0] +
1245 x[1] * y[5] * z[6] - y[6] * x[5] * z[7] - y[2] * x[3] * z[1] +
1246 z[1] * x[5] * y[6] - y[5] * x[1] * z[4] + z[6] * x[4] * y[7] +
1247 x[5] * y[1] * z[4] - x[5] * y[6] * z[4] + y[6] * x[3] * z[7] -
1248 x[5] * y[4] * z[1];
1249 s5 = s8 + x[5] * y[4] * z[6] + z[5] * x[1] * y[4] + y[1] * x[6] * z[5] -
1250 z[6] * x[3] * y[7] + z[6] * x[7] * y[3] - z[5] * x[6] * y[4] -
1251 z[5] * x[4] * y[1] + z[5] * x[4] * y[6] + x[1] * y[6] * z[2] +
1252 x[2] * y[6] * z[3] + z[2] * x[6] * y[3] + z[1] * x[6] * y[2] +
1253 z[2] * x[3] * y[1] - z[2] * x[1] * y[3] - z[2] * x[3] * y[6] +
1254 y[2] * x[1] * z[3] + y[1] * x[2] * z[6] - z[0] * x[7] * y[3] + s7;
1255 s4 = 1 / s5;
1256 s2 = s3 * s4;
1257 const double unknown2 = s1 * s2;
1258
1259 return {unknown0, unknown1, unknown2};
1260 }
1261 else
1262 {
1263 // Be somewhat particular in which exception we throw
1268
1269 return {};
1270 }
1271 }
1272
1273
1274
1275 template <int structdim, int dim, int spacedim>
1277 barycenter(const TriaAccessor<structdim, dim, spacedim> &)
1278 {
1279 // this function catches all the cases not
1280 // explicitly handled above
1282 return {};
1283 }
1284
1285
1286
1287 template <int dim, int spacedim>
1288 double
1289 measure(const TriaAccessor<1, dim, spacedim> &accessor)
1290 {
1291 // remember that we use (dim-)linear
1292 // mappings
1293 return (accessor.vertex(1) - accessor.vertex(0)).norm();
1294 }
1295
1296
1297
1298 double
1299 measure(const TriaAccessor<2, 2, 2> &accessor)
1300 {
1302 for (const unsigned int i : accessor.vertex_indices())
1303 vertex_indices[i] = accessor.vertex_index(i);
1304
1306 accessor.get_triangulation().get_vertices(),
1308 }
1309
1310
1311 double
1312 measure(const TriaAccessor<3, 3, 3> &accessor)
1313 {
1315 for (const unsigned int i : accessor.vertex_indices())
1316 vertex_indices[i] = accessor.vertex_index(i);
1317
1319 accessor.get_triangulation().get_vertices(),
1321 }
1322
1323
1324 // a 2d face in 3d space
1325 template <int dim>
1326 double
1327 measure(const TriaAccessor<2, dim, 3> &accessor)
1328 {
1330 {
1331 const Point<3> x0 = accessor.vertex(0);
1332 const Point<3> x1 = accessor.vertex(1);
1333 const Point<3> x2 = accessor.vertex(2);
1334 const Point<3> x3 = accessor.vertex(3);
1335
1336 // This is based on the approach used in libMesh (see face_quad4.C): the
1337 // primary differences are the vertex numbering and quadrature order.
1338 //
1339 // The area of a surface is the integral of the magnitude of its normal
1340 // vector, which may be computed via the cross product of two tangent
1341 // vectors. We can easily get tangent vectors from the surface
1342 // parameterization. Hence, given a bilinear surface
1343 //
1344 // X(chi, eta) = x0 + (x1 - x0) chi + (x2 - x0) eta
1345 // + (x3 + x0 - x1 - x2) chi eta
1346 //
1347 // the tangent vectors are
1348 //
1349 // t1 = (x1 - x0) + (x3 + x0 - x1 - x2) eta
1350 // t2 = (x2 - x0) + (x3 + x0 - x1 - x2) xi
1351 const Tensor<1, 3> b0 = x1 - x0;
1352 const Tensor<1, 3> b1 = x2 - x0;
1353 const Tensor<1, 3> a = x3 - x2 - b0;
1354
1355 // The diameter is the maximum distance between any pair of vertices and
1356 // we can use it as a length scale for the cell. If all components of a
1357 // (the vector connecting x3 and the last vertex of the parallelogram
1358 // defined by the first three vertices) are zero within some tolerance,
1359 // then we have a parallelogram and can use a much simpler formula.
1360 double a_max = 0.0;
1361 for (unsigned int d = 0; d < 3; ++d)
1362 a_max = std::max(std::abs(a[d]), a_max);
1363 if (a_max < 1e-14 * accessor.diameter())
1364 return cross_product_3d(b0, b1).norm();
1365
1366 // Otherwise, use a 4x4 quadrature to approximate the surface area.
1367 // Hard-code this in to prevent the extra overhead of always creating
1368 // the same QGauss rule.
1369 constexpr unsigned int n_qp = 4;
1370 const double c1 = 2.0 / 7.0 * std::sqrt(6.0 / 5.0);
1371 const double w0 = (18.0 - std::sqrt(30)) / 72.0;
1372 const double w1 = (18.0 + std::sqrt(30)) / 72.0;
1373
1374 const std::array<double, n_qp> q{{
1375 0.5 - std::sqrt(3.0 / 7.0 + c1) / 2.0,
1376 0.5 - std::sqrt(3.0 / 7.0 - c1) / 2.0,
1377 0.5 + std::sqrt(3.0 / 7.0 - c1) / 2.0,
1378 0.5 + std::sqrt(3.0 / 7.0 + c1) / 2.0,
1379 }};
1380 const std::array<double, n_qp> w{{w0, w1, w1, w0}};
1381
1382 double area = 0.;
1383 for (unsigned int i = 0; i < n_qp; ++i)
1384 for (unsigned int j = 0; j < n_qp; ++j)
1385 area += cross_product_3d(q[i] * a + b0, q[j] * a + b1).norm() *
1386 w[i] * w[j];
1387
1388 return area;
1389 }
1390 else if (accessor.reference_cell() == ReferenceCells::Triangle)
1391 {
1392 // We can just use the normal triangle area formula without issue
1393 const Tensor<1, 3> v01 = accessor.vertex(1) - accessor.vertex(0);
1394 const Tensor<1, 3> v02 = accessor.vertex(2) - accessor.vertex(0);
1395 return 0.5 * cross_product_3d(v01, v02).norm();
1396 }
1397
1399 return 0.0;
1400 }
1401
1402
1403
1404 template <int structdim, int dim, int spacedim>
1405 double
1407 {
1408 // catch-all for all cases not explicitly
1409 // listed above
1411 return std::numeric_limits<double>::quiet_NaN();
1412 }
1413
1414
1415 template <int structdim, int dim, int spacedim>
1417 get_new_point_on_object(const TriaAccessor<structdim, dim, spacedim> &obj,
1418 const bool use_interpolation)
1419 {
1420 if (use_interpolation)
1421 {
1423 const auto points_and_weights =
1424 Manifolds::get_default_points_and_weights(it, use_interpolation);
1425 return obj.get_manifold().get_new_point(
1426 make_array_view(points_and_weights.first.begin(),
1427 points_and_weights.first.end()),
1428 make_array_view(points_and_weights.second.begin(),
1429 points_and_weights.second.end()));
1430 }
1431 else
1432 {
1434 if constexpr (structdim == 1)
1435 return obj.get_manifold().get_new_point_on_line(it);
1436 else if constexpr (structdim == 2)
1437 return obj.get_manifold().get_new_point_on_quad(it);
1438 else if constexpr (structdim == 3)
1439 return obj.get_manifold().get_new_point_on_hex(it);
1440 else
1442
1443 return {};
1444 }
1445 }
1446} // namespace
1447
1448
1449
1450/*-------------------- Static variables: TriaAccessorBase -------------------*/
1451#ifndef DOXYGEN
1452
1453template <int structdim, int dim, int spacedim>
1455
1456template <int structdim, int dim, int spacedim>
1458
1459template <int structdim, int dim, int spacedim>
1460const unsigned int
1462
1463#endif
1464/*------------------------ Functions: TriaAccessor ---------------------------*/
1465#ifndef DOXYGEN
1466
1467template <int structdim, int dim, int spacedim>
1468void
1470 const std::initializer_list<int> &new_indices) const
1471{
1472 const ArrayView<int> bounding_object_index_ref =
1473 this->objects().get_bounding_object_indices(this->present_index);
1474
1475 AssertIndexRange(new_indices.size(), bounding_object_index_ref.size() + 1);
1476
1477 unsigned int i = 0;
1478 for (const auto &new_index : new_indices)
1479 {
1480 bounding_object_index_ref[i] = new_index;
1481 ++i;
1482 }
1483}
1484
1485
1486
1487template <int structdim, int dim, int spacedim>
1488void
1490 const std::initializer_list<unsigned int> &new_indices) const
1491{
1492 const ArrayView<int> bounding_object_index_ref =
1493 this->objects().get_bounding_object_indices(this->present_index);
1494
1495 AssertIndexRange(new_indices.size(), bounding_object_index_ref.size() + 1);
1496
1497 unsigned int i = 0;
1498 for (const auto &new_index : new_indices)
1499 {
1500 bounding_object_index_ref[i] = new_index;
1501 ++i;
1502 }
1503}
1504
1505
1506
1507template <int structdim, int dim, int spacedim>
1510{
1511 // call the function in the anonymous
1512 // namespace above
1513 return ::barycenter(*this);
1514}
1515
1516
1517
1518template <int structdim, int dim, int spacedim>
1519double
1521{
1522 // call the function in the anonymous
1523 // namespace above
1524 return ::measure(*this);
1525}
1526
1527
1528
1529template <int structdim, int dim, int spacedim>
1532{
1533 std::pair<Point<spacedim>, Point<spacedim>> boundary_points =
1534 std::make_pair(this->vertex(0), this->vertex(0));
1535
1536 const unsigned int n_vertices = this->n_vertices();
1537 for (unsigned int v = 1; v < n_vertices; ++v)
1538 {
1539 const Point<spacedim> x = this->vertex(v);
1540 for (unsigned int k = 0; k < spacedim; ++k)
1541 {
1542 boundary_points.first[k] = std::min(boundary_points.first[k], x[k]);
1543 boundary_points.second[k] = std::max(boundary_points.second[k], x[k]);
1544 }
1545 }
1546
1547 return BoundingBox<spacedim>(boundary_points);
1548}
1549
1550
1551
1552template <int structdim, int dim, int spacedim>
1553double
1555 const unsigned int /*axis*/) const
1556{
1558 return std::numeric_limits<double>::signaling_NaN();
1559}
1560
1561#endif
1562
1563template <>
1564double
1566{
1567 AssertIndexRange(axis, 1);
1568
1569 return this->diameter();
1570}
1571
1572
1573template <>
1574double
1576{
1577 AssertIndexRange(axis, 1);
1578
1579 return this->diameter();
1580}
1581
1582
1583template <>
1584double
1586{
1587 Assert(this->reference_cell() == ReferenceCells::Quadrilateral,
1589
1590 constexpr unsigned int lines[2][2] = {
1591 {2, 3}, // Lines along x-axis, see GeometryInfo
1592 {0, 1}}; // Lines along y-axis
1593
1594 AssertIndexRange(axis, 2);
1595
1596 return std::max(this->line(lines[axis][0])->diameter(),
1597 this->line(lines[axis][1])->diameter());
1598}
1599
1600template <>
1601double
1603{
1604 Assert(this->reference_cell() == ReferenceCells::Quadrilateral,
1606
1607 constexpr unsigned int lines[2][2] = {
1608 {2, 3}, // Lines along x-axis, see GeometryInfo
1609 {0, 1}}; // Lines along y-axis
1611 AssertIndexRange(axis, 2);
1612
1613 return std::max(this->line(lines[axis][0])->diameter(),
1614 this->line(lines[axis][1])->diameter());
1615}
1616
1617
1618template <>
1619double
1621{
1622 Assert(this->reference_cell() == ReferenceCells::Hexahedron,
1624
1625 constexpr unsigned int lines[3][4] = {
1626 {2, 3, 6, 7}, // Lines along x-axis, see GeometryInfo
1627 {0, 1, 4, 5}, // Lines along y-axis
1628 {8, 9, 10, 11}}; // Lines along z-axis
1629
1630 AssertIndexRange(axis, 3);
1631
1632 const double lengths[4] = {this->line(lines[axis][0])->diameter(),
1633 this->line(lines[axis][1])->diameter(),
1634 this->line(lines[axis][2])->diameter(),
1635 this->line(lines[axis][3])->diameter()};
1636
1637 return std::max({lengths[0], lengths[1], lengths[2], lengths[3]});
1638}
1639
1640
1641// Recursively set manifold ids on hex iterators.
1642template <>
1643void
1645 const types::manifold_id manifold_ind) const
1646{
1647 set_manifold_id(manifold_ind);
1648
1649 if (this->has_children())
1650 for (unsigned int c = 0; c < this->n_children(); ++c)
1651 this->child(c)->set_all_manifold_ids(manifold_ind);
1652
1653 // for hexes also set manifold_id
1654 // of bounding quads and lines
1655
1656 for (const unsigned int i : this->face_indices())
1657 this->quad(i)->set_manifold_id(manifold_ind);
1658 for (const unsigned int i : this->line_indices())
1659 this->line(i)->set_manifold_id(manifold_ind);
1660}
1661
1662
1663template <int structdim, int dim, int spacedim>
1666 const Point<structdim> &coordinates) const
1667{
1668 // Surrounding points and weights.
1669 std::array<Point<spacedim>, GeometryInfo<structdim>::vertices_per_cell> p;
1670 std::array<double, GeometryInfo<structdim>::vertices_per_cell> w;
1671
1672 for (const unsigned int i : this->vertex_indices())
1673 {
1674 p[i] = this->vertex(i);
1676 }
1677
1678 return this->get_manifold().get_new_point(make_array_view(p.begin(), p.end()),
1679 make_array_view(w.begin(),
1680 w.end()));
1681}
1682
1683
1684
1685template <int structdim, int dim, int spacedim>
1688 const Point<spacedim> &point) const
1689{
1690 std::array<Point<spacedim>, GeometryInfo<structdim>::vertices_per_cell>
1691 vertices;
1692 for (const unsigned int v : this->vertex_indices())
1693 vertices[v] = this->vertex(v);
1694
1695 const auto A_b =
1696 GridTools::affine_cell_approximation<structdim, spacedim>(vertices);
1698 A_b.first.covariant_form().transpose();
1699 return Point<structdim>(apply_transformation(A_inv, point - A_b.second));
1700}
1701
1702
1703
1704template <int structdim, int dim, int spacedim>
1707 const bool respect_manifold,
1708 const bool use_interpolation) const
1709{
1710 if (respect_manifold == false)
1711 {
1712 Assert(use_interpolation == false, ExcNotImplemented());
1714 for (const unsigned int v : this->vertex_indices())
1715 p += vertex(v);
1716 return p / this->n_vertices();
1717 }
1718 else
1719 return get_new_point_on_object(*this, use_interpolation);
1720}
1721
1722
1723/*---------------- Functions: TriaAccessor<0,1,spacedim> -------------------*/
1724
1725
1726template <int spacedim>
1727bool
1734
1735
1736
1737template <int spacedim>
1738void
1744
1745
1746
1747template <int spacedim>
1748void
1754
1755
1756
1757template <int spacedim>
1758void
1760{
1761 set_user_flag();
1762
1763 if (this->has_children())
1764 for (unsigned int c = 0; c < this->n_children(); ++c)
1765 this->child(c)->recursively_set_user_flag();
1766}
1767
1768
1769
1770template <int spacedim>
1771void
1773{
1774 clear_user_flag();
1775
1776 if (this->has_children())
1777 for (unsigned int c = 0; c < this->n_children(); ++c)
1778 this->child(c)->recursively_clear_user_flag();
1779}
1780
1781
1782
1783template <int spacedim>
1784void
1790
1791
1792
1793template <int spacedim>
1794void
1800
1801
1802
1803template <int spacedim>
1804void
1810
1811
1812
1813template <int spacedim>
1814void *
1816{
1819 return nullptr;
1820}
1821
1822
1823
1824template <int spacedim>
1825void
1827{
1828 set_user_pointer(p);
1829
1830 if (this->has_children())
1831 for (unsigned int c = 0; c < this->n_children(); ++c)
1832 this->child(c)->recursively_set_user_pointer(p);
1833}
1834
1835
1836
1837template <int spacedim>
1838void
1840{
1841 clear_user_pointer();
1842
1843 if (this->has_children())
1844 for (unsigned int c = 0; c < this->n_children(); ++c)
1845 this->child(c)->recursively_clear_user_pointer();
1846}
1847
1848
1849
1850template <int spacedim>
1851void
1857
1858
1859
1860template <int spacedim>
1861void
1867
1868
1869
1870template <int spacedim>
1871unsigned int
1878
1879
1880
1881template <int spacedim>
1882void
1884{
1885 set_user_index(p);
1886
1887 if (this->has_children())
1888 for (unsigned int c = 0; c < this->n_children(); ++c)
1889 this->child(c)->recursively_set_user_index(p);
1890}
1891
1892
1893
1894template <int spacedim>
1895void
1897{
1898 clear_user_index();
1899
1900 if (this->has_children())
1901 for (unsigned int c = 0; c < this->n_children(); ++c)
1902 this->child(c)->recursively_clear_user_index();
1903}
1904
1905
1906
1907/*------------------------ Functions: CellAccessor<1> -----------------------*/
1908
1909
1910
1911template <>
1912bool
1914{
1915 return (this->vertex(0)[0] <= p[0]) && (p[0] <= this->vertex(1)[0]);
1916}
1917
1918
1919
1920/*------------------------ Functions: CellAccessor<2> -----------------------*/
1921
1922
1923
1924template <>
1925bool
1927{
1928 Assert(this->reference_cell() == ReferenceCells::Quadrilateral,
1930
1931 // we check whether the point is
1932 // inside the cell by making sure
1933 // that it on the inner side of
1934 // each line defined by the faces,
1935 // i.e. for each of the four faces
1936 // we take the line that connects
1937 // the two vertices and subdivide
1938 // the whole domain by that in two
1939 // and check whether the point is
1940 // on the `cell-side' (rather than
1941 // the `out-side') of this line. if
1942 // the point is on the `cell-side'
1943 // for all four faces, it must be
1944 // inside the cell.
1945
1946 // we want the faces in counter
1947 // clockwise orientation
1948 static const int direction[4] = {-1, 1, 1, -1};
1949 for (unsigned int f = 0; f < 4; ++f)
1950 {
1951 // vector from the first vertex
1952 // of the line to the point
1953 const Tensor<1, 2> to_p =
1954 p - this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 0));
1955 // vector describing the line
1956 const Tensor<1, 2> face =
1957 direction[f] *
1958 (this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 1)) -
1959 this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 0)));
1960
1961 // if we rotate the face vector
1962 // by 90 degrees to the left
1963 // (i.e. it points to the
1964 // inside) and take the scalar
1965 // product with the vector from
1966 // the vertex to the point,
1967 // then the point is in the
1968 // `cell-side' if the scalar
1969 // product is positive. if this
1970 // is not the case, we can be
1971 // sure that the point is
1972 // outside
1973 if ((-face[1] * to_p[0] + face[0] * to_p[1]) < 0)
1974 return false;
1975 }
1976
1977 // if we arrived here, then the
1978 // point is inside for all four
1979 // faces, and thus inside
1980 return true;
1981}
1982
1983
1984
1985/*------------------------ Functions: CellAccessor<3> -----------------------*/
1986
1987
1988
1989template <>
1990bool
1992{
1993 Assert(this->reference_cell() == ReferenceCells::Hexahedron,
1995
1996 // original implementation by Joerg
1997 // Weimar
1998
1999 // we first eliminate points based
2000 // on the maximum and minimum of
2001 // the corner coordinates, then
2002 // transform to the unit cell, and
2003 // check there.
2004 const unsigned int dim = 3;
2005 const unsigned int spacedim = 3;
2006 Point<spacedim> maxp = this->vertex(0);
2007 Point<spacedim> minp = this->vertex(0);
2008
2009 for (unsigned int v = 1; v < this->n_vertices(); ++v)
2010 for (unsigned int d = 0; d < dim; ++d)
2011 {
2012 maxp[d] = std::max(maxp[d], this->vertex(v)[d]);
2013 minp[d] = std::min(minp[d], this->vertex(v)[d]);
2014 }
2015
2016 // rule out points outside the
2017 // bounding box of this cell
2018 for (unsigned int d = 0; d < dim; ++d)
2019 if ((p[d] < minp[d]) || (p[d] > maxp[d]))
2020 return false;
2021
2022 // now we need to check more carefully: transform to the
2023 // unit cube and check there. unfortunately, this isn't
2024 // completely trivial since the transform_real_to_unit_cell
2025 // function may throw an exception that indicates that the
2026 // point given could not be inverted. we take this as a sign
2027 // that the point actually lies outside, as also documented
2028 // for that function
2029 try
2030 {
2031 const TriaRawIterator<CellAccessor<dim, spacedim>> cell_iterator(*this);
2033 reference_cell()
2034 .template get_default_linear_mapping<dim, spacedim>()
2035 .transform_real_to_unit_cell(cell_iterator, p)));
2036 }
2038 {
2039 return false;
2040 }
2041}
2042
2043
2044
2045/*------------------- Functions: CellAccessor<dim,spacedim> -----------------*/
2046
2047// The return type is the same as DoFHandler<dim,spacedim>::active_cell_iterator
2048template <int dim, int spacedim>
2051 const DoFHandler<dim, spacedim> &dof_handler) const
2052{
2053 Assert(is_active(),
2054 ExcMessage("The current iterator points to an inactive cell. "
2055 "You cannot convert it to an iterator to an active cell."));
2056 Assert(&this->get_triangulation() == &dof_handler.get_triangulation(),
2057 ExcMessage("The triangulation associated with the iterator does not "
2058 "match that of the DoFHandler."));
2059
2061 &dof_handler.get_triangulation(),
2062 this->level(),
2063 this->index(),
2064 &dof_handler);
2065}
2066
2067
2068
2069template <int dim, int spacedim>
2072 const DoFHandler<dim, spacedim> &dof_handler) const
2073{
2074 Assert(&this->get_triangulation() == &dof_handler.get_triangulation(),
2075 ExcMessage("The triangulation associated with the iterator does not "
2076 "match that of the DoFHandler."));
2077
2079 &dof_handler.get_triangulation(),
2080 this->level(),
2081 this->index(),
2082 &dof_handler);
2083}
2084
2085
2086
2087// For codim>0 we proceed as follows:
2088// 1) project point onto manifold and
2089// 2) transform to the unit cell with a Q1 mapping
2090// 3) then check if inside unit cell
2091template <int dim, int spacedim>
2092template <int dim_, int spacedim_>
2093bool
2095{
2096 Assert(this->reference_cell().is_hyper_cube(), ExcNotImplemented());
2097
2098 const TriaRawIterator<CellAccessor<dim_, spacedim_>> cell_iterator(*this);
2099
2100 const Point<dim_> p_unit =
2101 this->reference_cell()
2102 .template get_default_linear_mapping<dim_, spacedim_>()
2103 .transform_real_to_unit_cell(cell_iterator, p);
2104
2106}
2107
2108
2109
2110template <>
2111bool
2113{
2114 return point_inside_codim<1, 2>(p);
2115}
2116
2117
2118template <>
2119bool
2121{
2122 return point_inside_codim<1, 3>(p);
2123}
2124
2125
2126template <>
2127bool
2129{
2130 Assert(this->reference_cell() == ReferenceCells::Quadrilateral,
2132 return point_inside_codim<2, 3>(p);
2133}
2134
2135
2136
2137template <int dim, int spacedim>
2138bool
2140{
2141 for (const auto face : this->face_indices())
2142 if (at_boundary(face))
2143 return true;
2144
2145 return false;
2146}
2147
2148
2149
2150template <int dim, int spacedim>
2153{
2155 return this->tria->levels[this->present_level]
2156 ->cells.boundary_or_material_id[this->present_index]
2157 .material_id;
2158}
2159
2160
2161
2162template <int dim, int spacedim>
2163void
2165 const types::material_id mat_id) const
2166{
2169 this->tria->levels[this->present_level]
2170 ->cells.boundary_or_material_id[this->present_index]
2171 .material_id = mat_id;
2172}
2173
2174
2175
2176template <int dim, int spacedim>
2177void
2179 const types::material_id mat_id) const
2180{
2181 set_material_id(mat_id);
2182
2183 if (this->has_children())
2184 for (unsigned int c = 0; c < this->n_children(); ++c)
2185 this->child(c)->recursively_set_material_id(mat_id);
2186}
2187
2188
2189
2190template <int dim, int spacedim>
2191void
2193 const types::subdomain_id new_subdomain_id) const
2194{
2196 Assert(this->is_active(),
2197 ExcMessage("set_subdomain_id() can only be called on active cells!"));
2198 this->tria->levels[this->present_level]->subdomain_ids[this->present_index] =
2199 new_subdomain_id;
2200}
2201
2202
2203
2204template <int dim, int spacedim>
2205void
2207 const types::subdomain_id new_level_subdomain_id) const
2208{
2210 this->tria->levels[this->present_level]
2211 ->level_subdomain_ids[this->present_index] = new_level_subdomain_id;
2212}
2213
2214
2215template <int dim, int spacedim>
2216bool
2218{
2220 if constexpr (dim == spacedim)
2221 return true;
2222 else if constexpr (dim == spacedim - 1)
2223 return this->tria->levels[this->present_level]
2224 ->direction_flags[this->present_index];
2225 else
2226 {
2227 Assert(false,
2228 ExcMessage("This function cannot be called if dim<spacedim-1."));
2229 return true;
2230 }
2231}
2232
2233
2234
2235template <int dim, int spacedim>
2236void
2238 const bool new_direction_flag) const
2239{
2240 // Some older compilers (GCC 9) print an unused variable warning about
2241 // new_direction_flag when it is only used in a subset of 'if constexpr'
2242 // statements
2244 if constexpr (dim == spacedim)
2245 Assert(new_direction_flag == true,
2246 ExcMessage("If dim==spacedim, direction flags are always true and "
2247 "can not be set to anything else."));
2248 else if constexpr (dim == spacedim - 1)
2249 this->tria->levels[this->present_level]
2250 ->direction_flags[this->present_index] = new_direction_flag;
2251 else
2252 Assert(new_direction_flag == true,
2253 ExcMessage("If dim<spacedim-1, then this function can be called "
2254 "only if the argument is 'true'."));
2255}
2256
2257
2258
2259template <int dim, int spacedim>
2260void
2261CellAccessor<dim, spacedim>::set_parent(const unsigned int parent_index)
2262{
2264 Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2265 this->tria->levels[this->present_level]->parents[this->present_index / 2] =
2266 parent_index;
2267}
2268
2269
2270
2271template <int dim, int spacedim>
2272int
2274{
2275 Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2276
2277 // the parent of two consecutive cells
2278 // is stored only once, since it is
2279 // the same
2280 return this->tria->levels[this->present_level]
2281 ->parents[this->present_index / 2];
2282}
2283
2284
2285
2286template <int dim, int spacedim>
2287void
2289 const unsigned int active_cell_index) const
2290{
2291 this->tria->levels[this->present_level]
2292 ->active_cell_indices[this->present_index] = active_cell_index;
2293}
2294
2295
2296
2297template <int dim, int spacedim>
2298void
2300 const types::global_cell_index index) const
2301{
2302 this->tria->levels[this->present_level]
2303 ->global_active_cell_indices[this->present_index] = index;
2304}
2305
2306
2307
2308template <int dim, int spacedim>
2309void
2311 const types::global_cell_index index) const
2312{
2313 this->tria->levels[this->present_level]
2314 ->global_level_cell_indices[this->present_index] = index;
2315}
2316
2317
2318
2319template <int dim, int spacedim>
2322{
2324 Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2326 this->present_level - 1,
2327 parent_index());
2328
2329 return q;
2330}
2331
2332
2333template <int dim, int spacedim>
2334void
2336 const types::subdomain_id new_subdomain_id) const
2337{
2338 if (this->has_children())
2339 for (unsigned int c = 0; c < this->n_children(); ++c)
2340 this->child(c)->recursively_set_subdomain_id(new_subdomain_id);
2341 else
2342 set_subdomain_id(new_subdomain_id);
2343}
2344
2345
2346
2347template <int dim, int spacedim>
2348void
2350 const unsigned int i,
2351 const TriaIterator<CellAccessor<dim, spacedim>> &pointer) const
2352{
2353 AssertIndexRange(i, this->n_faces());
2354
2355 auto &neighbor =
2356 this->tria->levels[this->present_level]
2357 ->neighbors[this->present_index * ReferenceCells::max_n_faces<dim>() + i];
2358 if (pointer.state() == IteratorState::valid)
2359 {
2360 neighbor.first = pointer->present_level;
2361 neighbor.second = pointer->present_index;
2362 }
2363 else
2364 {
2365 neighbor.first = -1;
2366 neighbor.second = -1;
2367 }
2368}
2369
2370
2371
2372template <int dim, int spacedim>
2373CellId
2375{
2376 std::array<unsigned char, 30> id;
2377
2378 CellAccessor<dim, spacedim> ptr = *this;
2379 const unsigned int n_child_indices = ptr.level();
2380
2381 while (ptr.level() > 0)
2382 {
2384 const unsigned int n_children = parent->n_children();
2385
2386 // determine which child we are
2387 unsigned char v = static_cast<unsigned char>(-1);
2388 for (unsigned int c = 0; c < n_children; ++c)
2389 {
2390 if (parent->child_index(c) == ptr.index())
2391 {
2392 v = c;
2393 break;
2394 }
2395 }
2396
2397 Assert(v != static_cast<unsigned char>(-1), ExcInternalError());
2398 id[ptr.level() - 1] = v;
2399
2400 ptr.copy_from(*parent);
2401 }
2402
2403 Assert(ptr.level() == 0, ExcInternalError());
2404 const unsigned int coarse_index = ptr.index();
2405
2406 return {this->tria->coarse_cell_index_to_coarse_cell_id(coarse_index),
2407 n_child_indices,
2408 id.data()};
2409}
2410
2411
2412
2413template <int dim, int spacedim>
2414unsigned int
2416 const unsigned int neighbor) const
2417{
2418 AssertIndexRange(neighbor, this->n_faces());
2419
2420 // if we have a 1d mesh in 1d, we
2421 // can assume that the left
2422 // neighbor of the right neighbor is
2423 // the current cell. but that is an
2424 // invariant that isn't true if the
2425 // mesh is embedded in a higher
2426 // dimensional space, so we have to
2427 // fall back onto the generic code
2428 // below
2429 if ((dim == 1) && (spacedim == dim))
2430 return GeometryInfo<dim>::opposite_face[neighbor];
2431
2432 const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2433 this->neighbor(neighbor);
2434
2435 // usually, on regular patches of
2436 // the grid, this cell is just on
2437 // the opposite side of the
2438 // neighbor that the neighbor is of
2439 // this cell. for example in 2d, if
2440 // we want to know the
2441 // neighbor_of_neighbor if
2442 // neighbor==1 (the right
2443 // neighbor), then we will get 3
2444 // (the left neighbor) in most
2445 // cases. look up this relationship
2446 // in the table provided by
2447 // GeometryInfo and try it
2448 const unsigned int this_face_index = face_index(neighbor);
2449
2450 const unsigned int neighbor_guess =
2452
2453 if (neighbor_guess < neighbor_cell->n_faces() &&
2454 neighbor_cell->face_index(neighbor_guess) == this_face_index)
2455 return neighbor_guess;
2456 else
2457 // if the guess was false, then
2458 // we need to loop over all
2459 // neighbors and find the number
2460 // the hard way
2461 {
2462 for (const unsigned int face_no : neighbor_cell->face_indices())
2463 if (neighbor_cell->face_index(face_no) == this_face_index)
2464 return face_no;
2465
2466 // running over all neighbors
2467 // faces we did not find the
2468 // present face. Thereby the
2469 // neighbor must be coarser
2470 // than the present
2471 // cell. Return an invalid
2472 // unsigned int in this case.
2474 }
2475}
2476
2477
2478
2479template <int dim, int spacedim>
2480unsigned int
2482 const unsigned int face_no) const
2483{
2484 const unsigned int n2 = neighbor_of_neighbor_internal(face_no);
2487
2488 return n2;
2489}
2490
2491
2492
2493template <int dim, int spacedim>
2494bool
2496 const unsigned int face_no) const
2497{
2498 return neighbor_of_neighbor_internal(face_no) ==
2500}
2501
2502
2503
2504template <int dim, int spacedim>
2505std::pair<unsigned int, unsigned int>
2507 const unsigned int neighbor) const
2508{
2509 AssertIndexRange(neighbor, this->n_faces());
2510 // make sure that the neighbor is
2511 // on a coarser level
2512 Assert(neighbor_is_coarser(neighbor),
2514
2515 switch (dim)
2516 {
2517 case 2:
2518 {
2519 const int this_face_index = face_index(neighbor);
2520 const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2521 this->neighbor(neighbor);
2522
2523 // usually, on regular patches of
2524 // the grid, this cell is just on
2525 // the opposite side of the
2526 // neighbor that the neighbor is of
2527 // this cell. for example in 2d, if
2528 // we want to know the
2529 // neighbor_of_neighbor if
2530 // neighbor==1 (the right
2531 // neighbor), then we will get 0
2532 // (the left neighbor) in most
2533 // cases. look up this relationship
2534 // in the table provided by
2535 // GeometryInfo and try it
2536 const unsigned int face_no_guess =
2538
2539 const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face_guess =
2540 neighbor_cell->face(face_no_guess);
2541
2542 if (face_guess->has_children())
2543 for (unsigned int subface_no = 0;
2544 subface_no < face_guess->n_children();
2545 ++subface_no)
2546 if (face_guess->child_index(subface_no) == this_face_index)
2547 return std::make_pair(face_no_guess, subface_no);
2548
2549 // if the guess was false, then
2550 // we need to loop over all faces
2551 // and subfaces and find the
2552 // number the hard way
2553 for (const unsigned int face_no : neighbor_cell->face_indices())
2554 {
2555 if (face_no != face_no_guess)
2556 {
2557 const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>
2558 face = neighbor_cell->face(face_no);
2559 if (face->has_children())
2560 for (unsigned int subface_no = 0;
2561 subface_no < face->n_children();
2562 ++subface_no)
2563 if (face->child_index(subface_no) == this_face_index)
2564 return std::make_pair(face_no, subface_no);
2565 }
2566 }
2567
2568 // we should never get here,
2569 // since then we did not find
2570 // our way back...
2572 return std::make_pair(numbers::invalid_unsigned_int,
2574 }
2575
2576 case 3:
2577 {
2578 const int this_face_index = face_index(neighbor);
2579 const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2580 this->neighbor(neighbor);
2581
2582 // usually, on regular patches of the grid, this cell is just on the
2583 // opposite side of the neighbor that the neighbor is of this cell.
2584 // for example in 2d, if we want to know the neighbor_of_neighbor if
2585 // neighbor==1 (the right neighbor), then we will get 0 (the left
2586 // neighbor) in most cases. look up this relationship in the table
2587 // provided by GeometryInfo and try it
2588 const unsigned int face_no_guess =
2590
2591 const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face_guess =
2592 neighbor_cell->face(face_no_guess);
2593
2594 if (face_guess->has_children())
2595 for (unsigned int subface_no = 0;
2596 subface_no < face_guess->n_children();
2597 ++subface_no)
2598 {
2599 if (face_guess->child_index(subface_no) == this_face_index)
2600 // call a helper function, that translates the current
2601 // subface number to a subface number for the current
2602 // FaceRefineCase
2603 return std::make_pair(face_no_guess,
2604 translate_subface_no(face_guess,
2605 subface_no));
2606
2607 if (face_guess->child(subface_no)->has_children())
2608 for (unsigned int subsub_no = 0;
2609 subsub_no < face_guess->child(subface_no)->n_children();
2610 ++subsub_no)
2611 if (face_guess->child(subface_no)->child_index(subsub_no) ==
2612 this_face_index)
2613 // call a helper function, that translates the current
2614 // subface number and subsubface number to a subface
2615 // number for the current FaceRefineCase
2616 return std::make_pair(face_no_guess,
2617 translate_subface_no(face_guess,
2618 subface_no,
2619 subsub_no));
2620 }
2621
2622 // if the guess was false, then we need to loop over all faces and
2623 // subfaces and find the number the hard way
2624 for (const unsigned int face_no : neighbor_cell->face_indices())
2625 {
2626 if (face_no == face_no_guess)
2627 continue;
2628
2629 const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face =
2630 neighbor_cell->face(face_no);
2631
2632 if (!face->has_children())
2633 continue;
2634
2635 for (unsigned int subface_no = 0; subface_no < face->n_children();
2636 ++subface_no)
2637 {
2638 if (face->child_index(subface_no) == this_face_index)
2639 // call a helper function, that translates the current
2640 // subface number to a subface number for the current
2641 // FaceRefineCase
2642 return std::make_pair(face_no,
2643 translate_subface_no(face,
2644 subface_no));
2645
2646 if (face->child(subface_no)->has_children())
2647 for (unsigned int subsub_no = 0;
2648 subsub_no < face->child(subface_no)->n_children();
2649 ++subsub_no)
2650 if (face->child(subface_no)->child_index(subsub_no) ==
2651 this_face_index)
2652 // call a helper function, that translates the current
2653 // subface number and subsubface number to a subface
2654 // number for the current FaceRefineCase
2655 return std::make_pair(face_no,
2656 translate_subface_no(face,
2657 subface_no,
2658 subsub_no));
2659 }
2660 }
2661
2662 // we should never get here, since then we did not find our way
2663 // back...
2665 return std::make_pair(numbers::invalid_unsigned_int,
2667 }
2668
2669 default:
2670 {
2671 Assert(false, ExcImpossibleInDim(1));
2672 return std::make_pair(numbers::invalid_unsigned_int,
2674 }
2675 }
2676}
2677
2678
2679
2680template <int dim, int spacedim>
2681bool
2683 const unsigned int i_face) const
2684{
2685 /*
2686 * Implementation note: In all of the functions corresponding to periodic
2687 * faces we mainly use the Triangulation::periodic_face_map to find the
2688 * information about periodically connected faces. So, we actually search in
2689 * this std::map and return the cell_face on the other side of the periodic
2690 * boundary.
2691 *
2692 * We can not use operator[] as this would insert non-existing entries or
2693 * would require guarding with an extra std::map::find() or count().
2694 */
2695 AssertIndexRange(i_face, this->n_faces());
2696 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2697
2698 if (at_boundary(i_face) && this->tria->periodic_face_map.find(
2699 std::make_pair(cell_iterator(*this), i_face)) !=
2700 this->tria->periodic_face_map.end())
2701 return true;
2702 return false;
2703}
2704
2705
2706
2707template <int dim, int spacedim>
2710{
2711 /*
2712 * To know, why we are using std::map::find() instead of [] operator, refer
2713 * to the implementation note in has_periodic_neighbor() function.
2714 *
2715 * my_it : the iterator to the current cell.
2716 * my_face_pair : the pair reported by periodic_face_map as its first pair
2717 * being the current cell_face.
2718 */
2719 AssertIndexRange(i_face, this->n_faces());
2720 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2721 cell_iterator current_cell(*this);
2722
2723 auto my_face_pair =
2724 this->tria->periodic_face_map.find(std::make_pair(current_cell, i_face));
2725
2726 // Make sure we are actually on a periodic boundary:
2727 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2729 return my_face_pair->second.first.first;
2730}
2731
2732
2733
2734template <int dim, int spacedim>
2737 const unsigned int i_face) const
2738{
2739 if (!(this->face(i_face)->at_boundary()))
2740 return this->neighbor(i_face);
2741 else if (this->has_periodic_neighbor(i_face))
2742 return this->periodic_neighbor(i_face);
2743 else
2745 // we can't come here
2746 return this->neighbor(i_face);
2747}
2748
2749
2750
2751template <int dim, int spacedim>
2754 const unsigned int i_face,
2755 const unsigned int i_subface) const
2756{
2757 /*
2758 * To know, why we are using std::map::find() instead of [] operator, refer
2759 * to the implementation note in has_periodic_neighbor() function.
2760 *
2761 * my_it : the iterator to the current cell.
2762 * my_face_pair : the pair reported by periodic_face_map as its first pair
2763 * being the current cell_face. nb_it : the iterator to the
2764 * neighbor of current cell at i_face. face_num_of_nb : the face number of
2765 * the periodically neighboring face in the relevant element.
2766 * nb_parent_face_it: the iterator to the parent face of the periodically
2767 * neighboring face.
2768 */
2769 AssertIndexRange(i_face, this->n_faces());
2770 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2771 cell_iterator my_it(*this);
2772
2773 auto my_face_pair =
2774 this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2775 /*
2776 * There should be an assertion, which tells the user that this function
2777 * should not be used for a cell which is not located at a periodic boundary.
2778 */
2779 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2781 cell_iterator parent_nb_it = my_face_pair->second.first.first;
2782 unsigned int nb_face_num = my_face_pair->second.first.second;
2783 TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> nb_parent_face_it =
2784 parent_nb_it->face(nb_face_num);
2785 /*
2786 * We should check if the parent face of the neighbor has at least the same
2787 * number of children as i_subface.
2788 */
2789 AssertIndexRange(i_subface, nb_parent_face_it->n_children());
2790
2791 const auto [orientation, rotation, flip] =
2792 internal::split_face_orientation(my_face_pair->second.second);
2793
2794 unsigned int sub_neighbor_num =
2795 GeometryInfo<dim>::child_cell_on_face(parent_nb_it->refinement_case(),
2796 nb_face_num,
2797 i_subface,
2798 orientation,
2799 flip,
2800 rotation,
2801 nb_parent_face_it->refinement_case());
2802 return parent_nb_it->child(sub_neighbor_num);
2803}
2804
2805
2806
2807template <int dim, int spacedim>
2808std::pair<unsigned int, unsigned int>
2810 const unsigned int i_face) const
2811{
2812 /*
2813 * To know, why we are using std::map::find() instead of [] operator, refer
2814 * to the implementation note in has_periodic_neighbor() function.
2815 *
2816 * my_it : the iterator to the current cell.
2817 * my_face_pair : the pair reported by periodic_face_map as its first pair
2818 * being the current cell_face. nb_it : the iterator to the periodic
2819 * neighbor. nb_face_pair : the pair reported by periodic_face_map as its
2820 * first pair being the periodic neighbor cell_face. p_nb_of_p_nb : the
2821 * iterator of the periodic neighbor of the periodic neighbor of the current
2822 * cell.
2823 */
2824 AssertIndexRange(i_face, this->n_faces());
2825 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2826 const int my_face_index = this->face_index(i_face);
2827 cell_iterator my_it(*this);
2828
2829 auto my_face_pair =
2830 this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2831 /*
2832 * There should be an assertion, which tells the user that this function
2833 * should not be used for a cell which is not located at a periodic boundary.
2834 */
2835 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2837 cell_iterator nb_it = my_face_pair->second.first.first;
2838 unsigned int face_num_of_nb = my_face_pair->second.first.second;
2839
2840 auto nb_face_pair =
2841 this->tria->periodic_face_map.find(std::make_pair(nb_it, face_num_of_nb));
2842 /*
2843 * Since, we store periodic neighbors for every cell (either active or
2844 * artificial or inactive) the nb_face_pair should also be mapped to some
2845 * cell_face pair. We assert this here.
2846 */
2847 Assert(nb_face_pair != this->tria->periodic_face_map.end(),
2849 cell_iterator p_nb_of_p_nb = nb_face_pair->second.first.first;
2850 TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> parent_face_it =
2851 p_nb_of_p_nb->face(nb_face_pair->second.first.second);
2852 for (unsigned int i_subface = 0; i_subface < parent_face_it->n_children();
2853 ++i_subface)
2854 if (parent_face_it->child_index(i_subface) == my_face_index)
2855 return std::make_pair(face_num_of_nb, i_subface);
2856 /*
2857 * Obviously, if the execution reaches to this point, some of our assumptions
2858 * should have been false. The most important one is, the user has called this
2859 * function on a face which does not have a coarser periodic neighbor.
2860 */
2862 return std::make_pair(numbers::invalid_unsigned_int,
2864}
2865
2866
2867
2868template <int dim, int spacedim>
2869int
2871 const unsigned int i_face) const
2872{
2873 return periodic_neighbor(i_face)->index();
2874}
2875
2876
2877
2878template <int dim, int spacedim>
2879int
2881 const unsigned int i_face) const
2882{
2883 return periodic_neighbor(i_face)->level();
2884}
2885
2886
2887
2888template <int dim, int spacedim>
2889unsigned int
2891 const unsigned int i_face) const
2892{
2893 return periodic_neighbor_face_no(i_face);
2894}
2895
2896
2897
2898template <int dim, int spacedim>
2899unsigned int
2901 const unsigned int i_face) const
2902{
2903 /*
2904 * To know, why we are using std::map::find() instead of [] operator, refer
2905 * to the implementation note in has_periodic_neighbor() function.
2906 *
2907 * my_it : the iterator to the current cell.
2908 * my_face_pair : the pair reported by periodic_face_map as its first pair
2909 * being the current cell_face.
2910 */
2911 AssertIndexRange(i_face, this->n_faces());
2912 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2913 cell_iterator my_it(*this);
2914
2915 auto my_face_pair =
2916 this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2917 /*
2918 * There should be an assertion, which tells the user that this function
2919 * should not be called for a cell which is not located at a periodic boundary
2920 * !
2921 */
2922 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2924 return my_face_pair->second.first.second;
2925}
2926
2927
2928
2929template <int dim, int spacedim>
2930bool
2932 const unsigned int i_face) const
2933{
2934 /*
2935 * To know, why we are using std::map::find() instead of [] operator, refer
2936 * to the implementation note in has_periodic_neighbor() function.
2937 *
2938 * Implementation note: Let p_nb_of_p_nb be the periodic neighbor of the
2939 * periodic neighbor of the current cell. Also, let p_face_of_p_nb_of_p_nb be
2940 * the periodic face of the p_nb_of_p_nb. If p_face_of_p_nb_of_p_nb has
2941 * children , then the periodic neighbor of the current cell is coarser than
2942 * itself. Although not tested, this implementation should work for
2943 * anisotropic refinement as well.
2944 *
2945 * my_it : the iterator to the current cell.
2946 * my_face_pair : the pair reported by periodic_face_map as its first pair
2947 * being the current cell_face. nb_it : the iterator to the periodic
2948 * neighbor. nb_face_pair : the pair reported by periodic_face_map as its
2949 * first pair being the periodic neighbor cell_face.
2950 */
2951 AssertIndexRange(i_face, this->n_faces());
2952 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2953 cell_iterator my_it(*this);
2954
2955 auto my_face_pair =
2956 this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2957 /*
2958 * There should be an assertion, which tells the user that this function
2959 * should not be used for a cell which is not located at a periodic boundary.
2960 */
2961 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2963
2964 cell_iterator nb_it = my_face_pair->second.first.first;
2965 unsigned int face_num_of_nb = my_face_pair->second.first.second;
2966
2967 auto nb_face_pair =
2968 this->tria->periodic_face_map.find(std::make_pair(nb_it, face_num_of_nb));
2969 /*
2970 * Since, we store periodic neighbors for every cell (either active or
2971 * artificial or inactive) the nb_face_pair should also be mapped to some
2972 * cell_face pair. We assert this here.
2973 */
2974 Assert(nb_face_pair != this->tria->periodic_face_map.end(),
2976 const unsigned int my_level = this->level();
2977 const unsigned int neighbor_level = nb_face_pair->second.first.first->level();
2978 Assert(my_level >= neighbor_level, ExcInternalError());
2979 return my_level > neighbor_level;
2980}
2981
2982
2983
2984template <int dim, int spacedim>
2985bool
2987{
2989 AssertIndexRange(i, this->n_faces());
2990
2991 return (neighbor_index(i) == -1);
2992}
2993
2994
2995
2996template <int dim, int spacedim>
2997bool
2999{
3000 if (dim == 1)
3001 return at_boundary();
3002 else
3003 {
3004 for (unsigned int l = 0; l < this->n_lines(); ++l)
3005 if (this->line(l)->at_boundary())
3006 return true;
3007
3008 return false;
3009 }
3010}
3011
3012
3013
3014template <int dim, int spacedim>
3017 const unsigned int face,
3018 const unsigned int subface) const
3019{
3020 Assert(!this->has_children(),
3021 ExcMessage("The present cell must not have children!"));
3022 Assert(!this->at_boundary(face),
3023 ExcMessage("The present cell must have a valid neighbor!"));
3024 Assert(this->neighbor(face)->has_children() == true,
3025 ExcMessage("The neighbor must have children!"));
3026
3027 switch (dim)
3028 {
3029 case 2:
3030 {
3031 if (this->reference_cell() == ReferenceCells::Triangle)
3032 {
3033 const unsigned int neighbor_neighbor =
3034 this->neighbor_of_neighbor(face);
3035 const auto neighbor = this->neighbor(face);
3036 // Triangles do not support anisotropic refinement
3037 Assert(neighbor->refinement_case() ==
3040 // Since we are looking at two faces at once, we need to check if
3041 // they have the same or opposing orientations rather than one
3042 // individual face orientation value.
3043 const auto relative_orientation =
3044 neighbor->combined_face_orientation(neighbor_neighbor) ==
3045 this->combined_face_orientation(face) ?
3048 const unsigned int neighbor_child_index =
3049 neighbor->reference_cell().child_cell_on_face(
3050 neighbor_neighbor, subface, relative_orientation);
3051 auto child = neighbor->child(neighbor_child_index);
3052 Assert(!child->has_children(), ExcInternalError());
3053 return child;
3054 }
3055 else if (this->reference_cell() == ReferenceCells::Quadrilateral)
3056 {
3057 const unsigned int neighbor_neighbor =
3058 this->neighbor_of_neighbor(face);
3059 const unsigned int neighbor_child_index =
3061 this->neighbor(face)->refinement_case(),
3062 neighbor_neighbor,
3063 subface);
3064
3066 this->neighbor(face)->child(neighbor_child_index);
3067 // the neighbors child can have children,
3068 // which are not further refined along the
3069 // face under consideration. as we are
3070 // normally interested in one of this
3071 // child's child, search for the right one.
3072 while (sub_neighbor->has_children())
3073 {
3075 sub_neighbor->refinement_case(),
3076 neighbor_neighbor) ==
3079 sub_neighbor =
3080 sub_neighbor->child(GeometryInfo<dim>::child_cell_on_face(
3081 sub_neighbor->refinement_case(), neighbor_neighbor, 0));
3082 }
3083
3084 return sub_neighbor;
3085 }
3086
3087 // if no reference cell type matches
3090 }
3091
3092
3093 case 3:
3094 {
3095 if (this->reference_cell() == ReferenceCells::Hexahedron)
3096 {
3097 // this function returns the neighbor's
3098 // child on a given face and
3099 // subface.
3100
3101 // we have to consider one other aspect here:
3102 // The face might be refined
3103 // anisotropically. In this case, the subface
3104 // number refers to the following, where we
3105 // look at the face from the current cell,
3106 // thus the subfaces are in standard
3107 // orientation concerning the cell
3108 //
3109 // for isotropic refinement
3110 //
3111 // *---*---*
3112 // | 2 | 3 |
3113 // *---*---*
3114 // | 0 | 1 |
3115 // *---*---*
3116 //
3117 // for 2*anisotropic refinement
3118 // (first cut_y, then cut_x)
3119 //
3120 // *---*---*
3121 // | 2 | 3 |
3122 // *---*---*
3123 // | 0 | 1 |
3124 // *---*---*
3125 //
3126 // for 2*anisotropic refinement
3127 // (first cut_x, then cut_y)
3128 //
3129 // *---*---*
3130 // | 1 | 3 |
3131 // *---*---*
3132 // | 0 | 2 |
3133 // *---*---*
3134 //
3135 // for purely anisotropic refinement:
3136 //
3137 // *---*---* *-------*
3138 // | | | | 1 |
3139 // | 0 | 1 | or *-------*
3140 // | | | | 0 |
3141 // *---*---* *-------*
3142 //
3143 // for "mixed" refinement:
3144 //
3145 // *---*---* *---*---* *---*---* *-------*
3146 // | | 2 | | 1 | | | 1 | 2 | | 2 |
3147 // | 0 *---* or *---* 2 | or *---*---* or *---*---*
3148 // | | 1 | | 0 | | | 0 | | 0 | 1 |
3149 // *---*---* *---*---* *-------* *---*---*
3150
3152 mother_face = this->face(face);
3153 const unsigned int total_children =
3154 mother_face->n_active_descendants();
3155 AssertIndexRange(subface, total_children);
3158
3159 unsigned int neighbor_neighbor;
3162 this->neighbor(face);
3163
3164
3165 const RefinementCase<dim - 1> mother_face_ref_case =
3166 mother_face->refinement_case();
3167 if (mother_face_ref_case ==
3168 static_cast<RefinementCase<dim - 1>>(
3169 RefinementCase<2>::cut_xy)) // total_children==4
3170 {
3171 // this case is quite easy. we are sure,
3172 // that the neighbor is not coarser.
3173
3174 // get the neighbor's number for the given
3175 // face and the neighbor
3176 neighbor_neighbor = this->neighbor_of_neighbor(face);
3177
3178 // now use the info provided by GeometryInfo
3179 // to extract the neighbors child number
3180 const unsigned int neighbor_child_index =
3182 neighbor->refinement_case(),
3183 neighbor_neighbor,
3184 subface,
3185 neighbor->face_orientation(neighbor_neighbor),
3186 neighbor->face_flip(neighbor_neighbor),
3187 neighbor->face_rotation(neighbor_neighbor));
3188 neighbor_child = neighbor->child(neighbor_child_index);
3189
3190 // make sure that the neighbor child cell we
3191 // have found shares the desired subface.
3192 Assert((this->face(face)->child(subface) ==
3193 neighbor_child->face(neighbor_neighbor)),
3195 }
3196 else //-> the face is refined anisotropically
3197 {
3198 // first of all, we have to find the
3199 // neighbor at one of the anisotropic
3200 // children of the
3201 // mother_face. determine, which of
3202 // these we need.
3203 unsigned int first_child_to_find;
3204 unsigned int neighbor_child_index;
3205 if (total_children == 2)
3206 first_child_to_find = subface;
3207 else
3208 {
3209 first_child_to_find = subface / 2;
3210 if (total_children == 3 && subface == 1 &&
3211 !mother_face->child(0)->has_children())
3212 first_child_to_find = 1;
3213 }
3214 if (neighbor_is_coarser(face))
3215 {
3216 std::pair<unsigned int, unsigned int> indices =
3217 neighbor_of_coarser_neighbor(face);
3218 neighbor_neighbor = indices.first;
3219
3220
3221 // we have to translate our
3222 // subface_index according to the
3223 // RefineCase and subface index of
3224 // the coarser face (our face is an
3225 // anisotropic child of the coarser
3226 // face), 'a' denotes our
3227 // subface_index 0 and 'b' denotes
3228 // our subface_index 1, whereas 0...3
3229 // denote isotropic subfaces of the
3230 // coarser face
3231 //
3232 // cut_x and coarser_subface_index=0
3233 //
3234 // *---*---*
3235 // |b=2| |
3236 // | | |
3237 // |a=0| |
3238 // *---*---*
3239 //
3240 // cut_x and coarser_subface_index=1
3241 //
3242 // *---*---*
3243 // | |b=3|
3244 // | | |
3245 // | |a=1|
3246 // *---*---*
3247 //
3248 // cut_y and coarser_subface_index=0
3249 //
3250 // *-------*
3251 // | |
3252 // *-------*
3253 // |a=0 b=1|
3254 // *-------*
3255 //
3256 // cut_y and coarser_subface_index=1
3257 //
3258 // *-------*
3259 // |a=2 b=3|
3260 // *-------*
3261 // | |
3262 // *-------*
3263 unsigned int iso_subface;
3264 if (neighbor->face(neighbor_neighbor)
3265 ->refinement_case() == RefinementCase<2>::cut_x)
3266 iso_subface = 2 * first_child_to_find + indices.second;
3267 else
3268 {
3269 Assert(neighbor->face(neighbor_neighbor)
3270 ->refinement_case() ==
3273 iso_subface =
3274 first_child_to_find + 2 * indices.second;
3275 }
3276 neighbor_child_index =
3278 neighbor->refinement_case(),
3279 neighbor_neighbor,
3280 iso_subface,
3281 neighbor->face_orientation(neighbor_neighbor),
3282 neighbor->face_flip(neighbor_neighbor),
3283 neighbor->face_rotation(neighbor_neighbor));
3284 }
3285 else // neighbor is not coarser
3286 {
3287 neighbor_neighbor = neighbor_of_neighbor(face);
3288 neighbor_child_index =
3290 neighbor->refinement_case(),
3291 neighbor_neighbor,
3292 first_child_to_find,
3293 neighbor->face_orientation(neighbor_neighbor),
3294 neighbor->face_flip(neighbor_neighbor),
3295 neighbor->face_rotation(neighbor_neighbor),
3296 mother_face_ref_case);
3297 }
3298
3299 neighbor_child = neighbor->child(neighbor_child_index);
3300 // it might be, that the neighbor_child
3301 // has children, which are not refined
3302 // along the given subface. go down that
3303 // list and deliver the last of those.
3304 while (
3305 neighbor_child->has_children() &&
3307 neighbor_child->refinement_case(), neighbor_neighbor) ==
3309 neighbor_child = neighbor_child->child(
3311 neighbor_child->refinement_case(),
3312 neighbor_neighbor,
3313 0));
3314
3315 // if there are two total subfaces, we
3316 // are finished. if there are four we
3317 // have to get a child of our current
3318 // neighbor_child. If there are three,
3319 // we have to check which of the two
3320 // possibilities applies.
3321 if (total_children == 3)
3322 {
3323 if (mother_face->child(0)->has_children())
3324 {
3325 if (subface < 2)
3326 neighbor_child = neighbor_child->child(
3328 neighbor_child->refinement_case(),
3329 neighbor_neighbor,
3330 subface,
3331 neighbor_child->face_orientation(
3332 neighbor_neighbor),
3333 neighbor_child->face_flip(neighbor_neighbor),
3334 neighbor_child->face_rotation(
3335 neighbor_neighbor),
3336 mother_face->child(0)->refinement_case()));
3337 }
3338 else
3339 {
3340 Assert(mother_face->child(1)->has_children(),
3342 if (subface > 0)
3343 neighbor_child = neighbor_child->child(
3345 neighbor_child->refinement_case(),
3346 neighbor_neighbor,
3347 subface - 1,
3348 neighbor_child->face_orientation(
3349 neighbor_neighbor),
3350 neighbor_child->face_flip(neighbor_neighbor),
3351 neighbor_child->face_rotation(
3352 neighbor_neighbor),
3353 mother_face->child(1)->refinement_case()));
3354 }
3355 }
3356 else if (total_children == 4)
3357 {
3358 neighbor_child = neighbor_child->child(
3360 neighbor_child->refinement_case(),
3361 neighbor_neighbor,
3362 subface % 2,
3363 neighbor_child->face_orientation(neighbor_neighbor),
3364 neighbor_child->face_flip(neighbor_neighbor),
3365 neighbor_child->face_rotation(neighbor_neighbor),
3366 mother_face->child(subface / 2)->refinement_case()));
3367 }
3368 }
3369
3370 // it might be, that the neighbor_child has
3371 // children, which are not refined along the
3372 // given subface. go down that list and
3373 // deliver the last of those.
3374 while (neighbor_child->has_children())
3375 neighbor_child =
3376 neighbor_child->child(GeometryInfo<dim>::child_cell_on_face(
3377 neighbor_child->refinement_case(), neighbor_neighbor, 0));
3378
3379 if constexpr (running_in_debug_mode())
3380 {
3381 // check, whether the face neighbor_child matches the
3382 // requested subface.
3384 requested;
3385 switch (this->subface_case(face))
3386 {
3390 requested = mother_face->child(subface);
3391 break;
3394 requested =
3395 mother_face->child(subface / 2)->child(subface % 2);
3396 break;
3397
3400 switch (subface)
3401 {
3402 case 0:
3403 case 1:
3404 requested = mother_face->child(0)->child(subface);
3405 break;
3406 case 2:
3407 requested = mother_face->child(1);
3408 break;
3409 default:
3411 }
3412 break;
3415 switch (subface)
3416 {
3417 case 0:
3418 requested = mother_face->child(0);
3419 break;
3420 case 1:
3421 case 2:
3422 requested =
3423 mother_face->child(1)->child(subface - 1);
3424 break;
3425 default:
3427 }
3428 break;
3429 default:
3431 break;
3432 }
3433 Assert(requested == neighbor_child->face(neighbor_neighbor),
3435 }
3436
3437 return neighbor_child;
3438 }
3439
3440 // if no reference cell type matches
3443 }
3444
3445 default:
3446 // if 1d or more than 3d
3449 }
3450}
3451
3452
3453
3454template <int structdim, int dim, int spacedim>
3460
3461
3462
3463template <int structdim, int dim, int spacedim>
3464int
3469
3470
3471
3472template <int structdim, int dim, int spacedim>
3473int
3478
3479
3480// explicit instantiations
3481#include "grid/tria_accessor.inst"
3482
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:954
std::size_t size() const
Definition array_view.h:689
void recursively_set_subdomain_id(const types::subdomain_id new_subdomain_id) const
TriaIterator< CellAccessor< dim, spacedim > > parent() const
void set_active_cell_index(const unsigned int active_cell_index) const
unsigned int neighbor_of_neighbor_internal(const unsigned int neighbor) const
TriaIterator< CellAccessor< dim, spacedim > > periodic_neighbor(const unsigned int i) const
void set_neighbor(const unsigned int i, const TriaIterator< CellAccessor< dim, spacedim > > &pointer) const
void set_direction_flag(const bool new_direction_flag) const
void recursively_set_material_id(const types::material_id new_material_id) const
void set_level_subdomain_id(const types::subdomain_id new_level_subdomain_id) const
TriaActiveIterator< DoFCellAccessor< dim, spacedim, false > > as_dof_handler_iterator(const DoFHandler< dim, spacedim > &dof_handler) const
TriaIterator< CellAccessor< dim, spacedim > > neighbor_child_on_subface(const unsigned int face_no, const unsigned int subface_no) const
void set_subdomain_id(const types::subdomain_id new_subdomain_id) const
bool neighbor_is_coarser(const unsigned int face_no) const
void set_global_level_cell_index(const types::global_cell_index index) const
bool has_periodic_neighbor(const unsigned int i) const
int periodic_neighbor_level(const unsigned int i) const
std::pair< unsigned int, unsigned int > neighbor_of_coarser_neighbor(const unsigned int neighbor) const
unsigned int neighbor_of_neighbor(const unsigned int face_no) const
void set_material_id(const types::material_id new_material_id) const
bool point_inside_codim(const Point< spacedim_ > &p) const
bool has_boundary_lines() const
TriaIterator< CellAccessor< dim, spacedim > > periodic_neighbor_child_on_subface(const unsigned int face_no, const unsigned int subface_no) const
int periodic_neighbor_index(const unsigned int i) const
bool periodic_neighbor_is_coarser(const unsigned int i) const
void set_global_active_cell_index(const types::global_cell_index index) const
void set_parent(const unsigned int parent_index)
std::pair< unsigned int, unsigned int > periodic_neighbor_of_coarser_periodic_neighbor(const unsigned face_no) const
bool at_boundary() const
bool point_inside(const Point< spacedim > &p) const
bool direction_flag() const
types::material_id material_id() const
CellId id() const
TriaIterator< DoFCellAccessor< dim, spacedim, true > > as_dof_handler_level_iterator(const DoFHandler< dim, spacedim > &dof_handler) const
TriaIterator< CellAccessor< dim, spacedim > > neighbor_or_periodic_neighbor(const unsigned int i) const
int parent_index() const
unsigned int periodic_neighbor_of_periodic_neighbor(const unsigned int i) const
unsigned int periodic_neighbor_face_no(const unsigned int i) const
DerivativeForm< 1, dim, spacedim, Number > covariant_form() const
DerivativeForm< 1, spacedim, dim, Number > transpose() const
const Triangulation< dim, spacedim > & get_triangulation() const
typename LevelSelector::cell_iterator level_cell_iterator
static int level()
static IteratorState::IteratorStates state()
static int index()
virtual Point< spacedim > get_new_point_on_hex(const typename Triangulation< dim, spacedim >::hex_iterator &hex) const
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
virtual Point< spacedim > get_new_point_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad) const
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim > > &surrounding_points, const ArrayView< const double > &weights) const
Abstract base class for mapping classes.
Definition mapping.h:320
Definition point.h:113
numbers::NumberTraits< Number >::real_type norm() const
void copy_from(const TriaAccessorBase &)
const Triangulation< dim, spacedim > & get_triangulation() const
int index() const
int level() const
void set_user_index(const unsigned int p) const
void clear_user_pointer() const
void recursively_set_user_index(const unsigned int p) const
void clear_user_data() const
Point< structdim > real_to_unit_cell_affine_approximation(const Point< spacedim > &point) const
void recursively_clear_user_index() const
const Manifold< dim, spacedim > & get_manifold() const
void recursively_set_user_pointer(void *p) const
double extent_in_direction(const unsigned int axis) const
Point< spacedim > intermediate_point(const Point< structdim > &coordinates) const
unsigned int n_vertices() const
void recursively_clear_user_flag() const
Point< spacedim > barycenter() const
BoundingBox< spacedim > bounding_box() const
void clear_user_flag() const
void recursively_set_user_flag() const
bool user_flag_set() const
void set_user_flag() const
void * user_pointer() const
Point< spacedim > center(const bool respect_manifold=false, const bool interpolate_from_surrounding=false) const
void clear_user_index() const
double measure() const
void set_bounding_object_indices(const std::initializer_list< int > &new_indices) const
Point< spacedim > & vertex(const unsigned int i) const
unsigned int user_index() const
void set_user_pointer(void *p) const
void recursively_clear_user_pointer() const
ReferenceCell reference_cell() const
double diameter() const
const std::vector< Point< spacedim > > & get_vertices() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
constexpr bool running_in_debug_mode()
Definition config.h:73
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
unsigned int level
Definition grid_out.cc:4635
unsigned int vertex_indices[2]
static ::ExceptionBase & ExcCellHasNoParent()
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcCellNotUsed()
static ::ExceptionBase & ExcNoPeriodicNeighbor()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
static ::ExceptionBase & ExcNeighborIsNotCoarser()
static ::ExceptionBase & ExcNeighborIsCoarser()
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void set_all_manifold_ids(const types::manifold_id) const
double cell_measure< 2 >(const std::vector< Point< 2 > > &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)
@ valid
Iterator points to a valid object.
@ invalid
Iterator is invalid, probably due to an error.
std::pair< std::array< Point< MeshIteratorType::AccessorType::space_dimension >, n_default_points_per_cell< MeshIteratorType >()>, std::array< double, n_default_points_per_cell< MeshIteratorType >()> > get_default_points_and_weights(const MeshIteratorType &iterator, const bool with_interpolation=false)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr ReferenceCell Triangle
constexpr ReferenceCell Hexahedron
constexpr ReferenceCell Pyramid
constexpr ReferenceCell Wedge
constexpr ReferenceCell Quadrilateral
constexpr ReferenceCell Tetrahedron
std::tuple< bool, bool, bool > split_face_orientation(const types::geometric_orientation combined_face_orientation)
constexpr unsigned int invalid_unsigned_int
Definition types.h:238
constexpr types::geometric_orientation reverse_line_orientation
Definition types.h:365
constexpr types::material_id invalid_material_id
Definition types.h:294
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:352
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement)
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
static bool is_inside_unit_cell(const Point< dim > &p)