Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-2988-ge5f3bfba9d 2025-04-01 14:50: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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
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 const unsigned int lines[2][2] = {
1588 {2, 3}, // Lines along x-axis, see GeometryInfo
1589 {0, 1}}; // Lines along y-axis
1590
1591 AssertIndexRange(axis, 2);
1592
1593 return std::max(this->line(lines[axis][0])->diameter(),
1594 this->line(lines[axis][1])->diameter());
1595}
1596
1597template <>
1598double
1600{
1601 const unsigned int lines[2][2] = {
1602 {2, 3}, // Lines along x-axis, see GeometryInfo
1603 {0, 1}}; // Lines along y-axis
1604
1605 AssertIndexRange(axis, 2);
1606
1607 return std::max(this->line(lines[axis][0])->diameter(),
1608 this->line(lines[axis][1])->diameter());
1609}
1611
1612template <>
1613double
1615{
1616 const unsigned int lines[3][4] = {
1617 {2, 3, 6, 7}, // Lines along x-axis, see GeometryInfo
1618 {0, 1, 4, 5}, // Lines along y-axis
1619 {8, 9, 10, 11}}; // Lines along z-axis
1620
1621 AssertIndexRange(axis, 3);
1622
1623 double lengths[4] = {this->line(lines[axis][0])->diameter(),
1624 this->line(lines[axis][1])->diameter(),
1625 this->line(lines[axis][2])->diameter(),
1626 this->line(lines[axis][3])->diameter()};
1627
1628 return std::max({lengths[0], lengths[1], lengths[2], lengths[3]});
1629}
1630
1631
1632// Recursively set manifold ids on hex iterators.
1633template <>
1634void
1636 const types::manifold_id manifold_ind) const
1637{
1638 set_manifold_id(manifold_ind);
1639
1640 if (this->has_children())
1641 for (unsigned int c = 0; c < this->n_children(); ++c)
1642 this->child(c)->set_all_manifold_ids(manifold_ind);
1643
1644 // for hexes also set manifold_id
1645 // of bounding quads and lines
1646
1647 for (const unsigned int i : this->face_indices())
1648 this->quad(i)->set_manifold_id(manifold_ind);
1649 for (const unsigned int i : this->line_indices())
1650 this->line(i)->set_manifold_id(manifold_ind);
1651}
1652
1653
1654template <int structdim, int dim, int spacedim>
1657 const Point<structdim> &coordinates) const
1658{
1659 // Surrounding points and weights.
1660 std::array<Point<spacedim>, GeometryInfo<structdim>::vertices_per_cell> p;
1661 std::array<double, GeometryInfo<structdim>::vertices_per_cell> w;
1662
1663 for (const unsigned int i : this->vertex_indices())
1664 {
1665 p[i] = this->vertex(i);
1667 }
1668
1669 return this->get_manifold().get_new_point(make_array_view(p.begin(), p.end()),
1670 make_array_view(w.begin(),
1671 w.end()));
1672}
1673
1674
1675
1676template <int structdim, int dim, int spacedim>
1679 const Point<spacedim> &point) const
1680{
1681 std::array<Point<spacedim>, GeometryInfo<structdim>::vertices_per_cell>
1682 vertices;
1683 for (const unsigned int v : this->vertex_indices())
1684 vertices[v] = this->vertex(v);
1685
1686 const auto A_b =
1687 GridTools::affine_cell_approximation<structdim, spacedim>(vertices);
1689 A_b.first.covariant_form().transpose();
1690 return Point<structdim>(apply_transformation(A_inv, point - A_b.second));
1691}
1692
1693
1694
1695template <int structdim, int dim, int spacedim>
1698 const bool respect_manifold,
1699 const bool use_interpolation) const
1700{
1701 if (respect_manifold == false)
1702 {
1703 Assert(use_interpolation == false, ExcNotImplemented());
1705 for (const unsigned int v : this->vertex_indices())
1706 p += vertex(v);
1707 return p / this->n_vertices();
1708 }
1709 else
1710 return get_new_point_on_object(*this, use_interpolation);
1711}
1712
1713
1714/*---------------- Functions: TriaAccessor<0,1,spacedim> -------------------*/
1715
1716
1717template <int spacedim>
1718bool
1725
1726
1727
1728template <int spacedim>
1729void
1735
1736
1737
1738template <int spacedim>
1739void
1745
1746
1747
1748template <int spacedim>
1749void
1751{
1752 set_user_flag();
1753
1754 if (this->has_children())
1755 for (unsigned int c = 0; c < this->n_children(); ++c)
1756 this->child(c)->recursively_set_user_flag();
1757}
1758
1759
1760
1761template <int spacedim>
1762void
1764{
1765 clear_user_flag();
1766
1767 if (this->has_children())
1768 for (unsigned int c = 0; c < this->n_children(); ++c)
1769 this->child(c)->recursively_clear_user_flag();
1770}
1771
1772
1773
1774template <int spacedim>
1775void
1781
1782
1783
1784template <int spacedim>
1785void
1791
1792
1793
1794template <int spacedim>
1795void
1801
1802
1803
1804template <int spacedim>
1805void *
1807{
1810 return nullptr;
1811}
1812
1813
1814
1815template <int spacedim>
1816void
1818{
1819 set_user_pointer(p);
1820
1821 if (this->has_children())
1822 for (unsigned int c = 0; c < this->n_children(); ++c)
1823 this->child(c)->recursively_set_user_pointer(p);
1824}
1825
1826
1827
1828template <int spacedim>
1829void
1831{
1832 clear_user_pointer();
1833
1834 if (this->has_children())
1835 for (unsigned int c = 0; c < this->n_children(); ++c)
1836 this->child(c)->recursively_clear_user_pointer();
1837}
1838
1839
1840
1841template <int spacedim>
1842void
1848
1849
1850
1851template <int spacedim>
1852void
1858
1859
1860
1861template <int spacedim>
1862unsigned int
1869
1870
1871
1872template <int spacedim>
1873void
1875{
1876 set_user_index(p);
1877
1878 if (this->has_children())
1879 for (unsigned int c = 0; c < this->n_children(); ++c)
1880 this->child(c)->recursively_set_user_index(p);
1881}
1882
1883
1884
1885template <int spacedim>
1886void
1888{
1889 clear_user_index();
1890
1891 if (this->has_children())
1892 for (unsigned int c = 0; c < this->n_children(); ++c)
1893 this->child(c)->recursively_clear_user_index();
1894}
1895
1896
1897
1898/*------------------------ Functions: CellAccessor<1> -----------------------*/
1899
1900
1901
1902template <>
1903bool
1905{
1906 return (this->vertex(0)[0] <= p[0]) && (p[0] <= this->vertex(1)[0]);
1907}
1908
1909
1910
1911/*------------------------ Functions: CellAccessor<2> -----------------------*/
1912
1913
1914
1915template <>
1916bool
1918{
1919 Assert(this->reference_cell() == ReferenceCells::Quadrilateral,
1921
1922 // we check whether the point is
1923 // inside the cell by making sure
1924 // that it on the inner side of
1925 // each line defined by the faces,
1926 // i.e. for each of the four faces
1927 // we take the line that connects
1928 // the two vertices and subdivide
1929 // the whole domain by that in two
1930 // and check whether the point is
1931 // on the `cell-side' (rather than
1932 // the `out-side') of this line. if
1933 // the point is on the `cell-side'
1934 // for all four faces, it must be
1935 // inside the cell.
1936
1937 // we want the faces in counter
1938 // clockwise orientation
1939 static const int direction[4] = {-1, 1, 1, -1};
1940 for (unsigned int f = 0; f < 4; ++f)
1941 {
1942 // vector from the first vertex
1943 // of the line to the point
1944 const Tensor<1, 2> to_p =
1945 p - this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 0));
1946 // vector describing the line
1947 const Tensor<1, 2> face =
1948 direction[f] *
1949 (this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 1)) -
1950 this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 0)));
1951
1952 // if we rotate the face vector
1953 // by 90 degrees to the left
1954 // (i.e. it points to the
1955 // inside) and take the scalar
1956 // product with the vector from
1957 // the vertex to the point,
1958 // then the point is in the
1959 // `cell-side' if the scalar
1960 // product is positive. if this
1961 // is not the case, we can be
1962 // sure that the point is
1963 // outside
1964 if ((-face[1] * to_p[0] + face[0] * to_p[1]) < 0)
1965 return false;
1966 }
1967
1968 // if we arrived here, then the
1969 // point is inside for all four
1970 // faces, and thus inside
1971 return true;
1972}
1973
1974
1975
1976/*------------------------ Functions: CellAccessor<3> -----------------------*/
1977
1978
1979
1980template <>
1981bool
1983{
1984 Assert(this->reference_cell() == ReferenceCells::Hexahedron,
1986
1987 // original implementation by Joerg
1988 // Weimar
1989
1990 // we first eliminate points based
1991 // on the maximum and minimum of
1992 // the corner coordinates, then
1993 // transform to the unit cell, and
1994 // check there.
1995 const unsigned int dim = 3;
1996 const unsigned int spacedim = 3;
1997 Point<spacedim> maxp = this->vertex(0);
1998 Point<spacedim> minp = this->vertex(0);
1999
2000 for (unsigned int v = 1; v < this->n_vertices(); ++v)
2001 for (unsigned int d = 0; d < dim; ++d)
2002 {
2003 maxp[d] = std::max(maxp[d], this->vertex(v)[d]);
2004 minp[d] = std::min(minp[d], this->vertex(v)[d]);
2005 }
2006
2007 // rule out points outside the
2008 // bounding box of this cell
2009 for (unsigned int d = 0; d < dim; ++d)
2010 if ((p[d] < minp[d]) || (p[d] > maxp[d]))
2011 return false;
2012
2013 // now we need to check more carefully: transform to the
2014 // unit cube and check there. unfortunately, this isn't
2015 // completely trivial since the transform_real_to_unit_cell
2016 // function may throw an exception that indicates that the
2017 // point given could not be inverted. we take this as a sign
2018 // that the point actually lies outside, as also documented
2019 // for that function
2020 try
2021 {
2022 const TriaRawIterator<CellAccessor<dim, spacedim>> cell_iterator(*this);
2024 reference_cell()
2025 .template get_default_linear_mapping<dim, spacedim>()
2026 .transform_real_to_unit_cell(cell_iterator, p)));
2027 }
2029 {
2030 return false;
2031 }
2032}
2033
2034
2035
2036/*------------------- Functions: CellAccessor<dim,spacedim> -----------------*/
2037
2038// The return type is the same as DoFHandler<dim,spacedim>::active_cell_iterator
2039template <int dim, int spacedim>
2042 const DoFHandler<dim, spacedim> &dof_handler) const
2043{
2044 Assert(is_active(),
2045 ExcMessage("The current iterator points to an inactive cell. "
2046 "You cannot convert it to an iterator to an active cell."));
2047 Assert(&this->get_triangulation() == &dof_handler.get_triangulation(),
2048 ExcMessage("The triangulation associated with the iterator does not "
2049 "match that of the DoFHandler."));
2050
2052 &dof_handler.get_triangulation(),
2053 this->level(),
2054 this->index(),
2055 &dof_handler);
2056}
2057
2058
2059
2060template <int dim, int spacedim>
2063 const DoFHandler<dim, spacedim> &dof_handler) const
2064{
2065 Assert(&this->get_triangulation() == &dof_handler.get_triangulation(),
2066 ExcMessage("The triangulation associated with the iterator does not "
2067 "match that of the DoFHandler."));
2068
2070 &dof_handler.get_triangulation(),
2071 this->level(),
2072 this->index(),
2073 &dof_handler);
2074}
2075
2076
2077
2078// For codim>0 we proceed as follows:
2079// 1) project point onto manifold and
2080// 2) transform to the unit cell with a Q1 mapping
2081// 3) then check if inside unit cell
2082template <int dim, int spacedim>
2083template <int dim_, int spacedim_>
2084bool
2086{
2087 Assert(this->reference_cell().is_hyper_cube(), ExcNotImplemented());
2088
2089 const TriaRawIterator<CellAccessor<dim_, spacedim_>> cell_iterator(*this);
2090
2091 const Point<dim_> p_unit =
2092 this->reference_cell()
2093 .template get_default_linear_mapping<dim_, spacedim_>()
2094 .transform_real_to_unit_cell(cell_iterator, p);
2095
2097}
2098
2099
2100
2101template <>
2102bool
2104{
2105 return point_inside_codim<1, 2>(p);
2106}
2107
2108
2109template <>
2110bool
2112{
2113 return point_inside_codim<1, 3>(p);
2114}
2115
2116
2117template <>
2118bool
2120{
2121 Assert(this->reference_cell() == ReferenceCells::Quadrilateral,
2123 return point_inside_codim<2, 3>(p);
2124}
2125
2126
2127
2128template <int dim, int spacedim>
2129bool
2131{
2132 for (const auto face : this->face_indices())
2133 if (at_boundary(face))
2134 return true;
2135
2136 return false;
2137}
2138
2139
2140
2141template <int dim, int spacedim>
2144{
2146 return this->tria->levels[this->present_level]
2147 ->cells.boundary_or_material_id[this->present_index]
2148 .material_id;
2149}
2150
2151
2152
2153template <int dim, int spacedim>
2154void
2156 const types::material_id mat_id) const
2157{
2160 this->tria->levels[this->present_level]
2161 ->cells.boundary_or_material_id[this->present_index]
2162 .material_id = mat_id;
2163}
2164
2165
2166
2167template <int dim, int spacedim>
2168void
2170 const types::material_id mat_id) const
2171{
2172 set_material_id(mat_id);
2173
2174 if (this->has_children())
2175 for (unsigned int c = 0; c < this->n_children(); ++c)
2176 this->child(c)->recursively_set_material_id(mat_id);
2177}
2178
2179
2180
2181template <int dim, int spacedim>
2182void
2184 const types::subdomain_id new_subdomain_id) const
2185{
2187 Assert(this->is_active(),
2188 ExcMessage("set_subdomain_id() can only be called on active cells!"));
2189 this->tria->levels[this->present_level]->subdomain_ids[this->present_index] =
2190 new_subdomain_id;
2191}
2192
2193
2194
2195template <int dim, int spacedim>
2196void
2198 const types::subdomain_id new_level_subdomain_id) const
2199{
2201 this->tria->levels[this->present_level]
2202 ->level_subdomain_ids[this->present_index] = new_level_subdomain_id;
2203}
2204
2205
2206template <int dim, int spacedim>
2207bool
2209{
2211 if constexpr (dim == spacedim)
2212 return true;
2213 else if constexpr (dim == spacedim - 1)
2214 return this->tria->levels[this->present_level]
2215 ->direction_flags[this->present_index];
2216 else
2217 {
2218 Assert(false,
2219 ExcMessage("This function cannot be called if dim<spacedim-1."));
2220 return true;
2221 }
2222}
2223
2224
2225
2226template <int dim, int spacedim>
2227void
2229 const bool new_direction_flag) const
2230{
2231 // Some older compilers (GCC 9) print an unused variable warning about
2232 // new_direction_flag when it is only used in a subset of 'if constexpr'
2233 // statements
2235 if constexpr (dim == spacedim)
2236 Assert(new_direction_flag == true,
2237 ExcMessage("If dim==spacedim, direction flags are always true and "
2238 "can not be set to anything else."));
2239 else if constexpr (dim == spacedim - 1)
2240 this->tria->levels[this->present_level]
2241 ->direction_flags[this->present_index] = new_direction_flag;
2242 else
2243 Assert(new_direction_flag == true,
2244 ExcMessage("If dim<spacedim-1, then this function can be called "
2245 "only if the argument is 'true'."));
2246}
2247
2248
2249
2250template <int dim, int spacedim>
2251void
2252CellAccessor<dim, spacedim>::set_parent(const unsigned int parent_index)
2253{
2255 Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2256 this->tria->levels[this->present_level]->parents[this->present_index / 2] =
2257 parent_index;
2258}
2259
2260
2261
2262template <int dim, int spacedim>
2263int
2265{
2266 Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2267
2268 // the parent of two consecutive cells
2269 // is stored only once, since it is
2270 // the same
2271 return this->tria->levels[this->present_level]
2272 ->parents[this->present_index / 2];
2273}
2274
2275
2276
2277template <int dim, int spacedim>
2278void
2280 const unsigned int active_cell_index) const
2281{
2282 this->tria->levels[this->present_level]
2283 ->active_cell_indices[this->present_index] = active_cell_index;
2284}
2285
2286
2287
2288template <int dim, int spacedim>
2289void
2291 const types::global_cell_index index) const
2292{
2293 this->tria->levels[this->present_level]
2294 ->global_active_cell_indices[this->present_index] = index;
2295}
2296
2297
2298
2299template <int dim, int spacedim>
2300void
2302 const types::global_cell_index index) const
2303{
2304 this->tria->levels[this->present_level]
2305 ->global_level_cell_indices[this->present_index] = index;
2306}
2307
2308
2309
2310template <int dim, int spacedim>
2313{
2315 Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2317 this->present_level - 1,
2318 parent_index());
2319
2320 return q;
2321}
2322
2323
2324template <int dim, int spacedim>
2325void
2327 const types::subdomain_id new_subdomain_id) const
2328{
2329 if (this->has_children())
2330 for (unsigned int c = 0; c < this->n_children(); ++c)
2331 this->child(c)->recursively_set_subdomain_id(new_subdomain_id);
2332 else
2333 set_subdomain_id(new_subdomain_id);
2334}
2335
2336
2337
2338template <int dim, int spacedim>
2339void
2341 const unsigned int i,
2342 const TriaIterator<CellAccessor<dim, spacedim>> &pointer) const
2343{
2344 AssertIndexRange(i, this->n_faces());
2345
2346 auto &neighbor =
2347 this->tria->levels[this->present_level]
2348 ->neighbors[this->present_index * ReferenceCells::max_n_faces<dim>() + i];
2349 if (pointer.state() == IteratorState::valid)
2350 {
2351 neighbor.first = pointer->present_level;
2352 neighbor.second = pointer->present_index;
2353 }
2354 else
2355 {
2356 neighbor.first = -1;
2357 neighbor.second = -1;
2358 }
2359}
2360
2361
2362
2363template <int dim, int spacedim>
2364CellId
2366{
2367 std::array<unsigned char, 30> id;
2368
2369 CellAccessor<dim, spacedim> ptr = *this;
2370 const unsigned int n_child_indices = ptr.level();
2371
2372 while (ptr.level() > 0)
2373 {
2375 const unsigned int n_children = parent->n_children();
2376
2377 // determine which child we are
2378 unsigned char v = static_cast<unsigned char>(-1);
2379 for (unsigned int c = 0; c < n_children; ++c)
2380 {
2381 if (parent->child_index(c) == ptr.index())
2382 {
2383 v = c;
2384 break;
2385 }
2386 }
2387
2388 Assert(v != static_cast<unsigned char>(-1), ExcInternalError());
2389 id[ptr.level() - 1] = v;
2390
2391 ptr.copy_from(*parent);
2392 }
2393
2394 Assert(ptr.level() == 0, ExcInternalError());
2395 const unsigned int coarse_index = ptr.index();
2396
2397 return {this->tria->coarse_cell_index_to_coarse_cell_id(coarse_index),
2398 n_child_indices,
2399 id.data()};
2400}
2401
2402
2403
2404template <int dim, int spacedim>
2405unsigned int
2407 const unsigned int neighbor) const
2408{
2409 AssertIndexRange(neighbor, this->n_faces());
2410
2411 // if we have a 1d mesh in 1d, we
2412 // can assume that the left
2413 // neighbor of the right neighbor is
2414 // the current cell. but that is an
2415 // invariant that isn't true if the
2416 // mesh is embedded in a higher
2417 // dimensional space, so we have to
2418 // fall back onto the generic code
2419 // below
2420 if ((dim == 1) && (spacedim == dim))
2421 return GeometryInfo<dim>::opposite_face[neighbor];
2422
2423 const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2424 this->neighbor(neighbor);
2425
2426 // usually, on regular patches of
2427 // the grid, this cell is just on
2428 // the opposite side of the
2429 // neighbor that the neighbor is of
2430 // this cell. for example in 2d, if
2431 // we want to know the
2432 // neighbor_of_neighbor if
2433 // neighbor==1 (the right
2434 // neighbor), then we will get 3
2435 // (the left neighbor) in most
2436 // cases. look up this relationship
2437 // in the table provided by
2438 // GeometryInfo and try it
2439 const unsigned int this_face_index = face_index(neighbor);
2440
2441 const unsigned int neighbor_guess =
2443
2444 if (neighbor_guess < neighbor_cell->n_faces() &&
2445 neighbor_cell->face_index(neighbor_guess) == this_face_index)
2446 return neighbor_guess;
2447 else
2448 // if the guess was false, then
2449 // we need to loop over all
2450 // neighbors and find the number
2451 // the hard way
2452 {
2453 for (const unsigned int face_no : neighbor_cell->face_indices())
2454 if (neighbor_cell->face_index(face_no) == this_face_index)
2455 return face_no;
2456
2457 // running over all neighbors
2458 // faces we did not find the
2459 // present face. Thereby the
2460 // neighbor must be coarser
2461 // than the present
2462 // cell. Return an invalid
2463 // unsigned int in this case.
2465 }
2466}
2467
2468
2469
2470template <int dim, int spacedim>
2471unsigned int
2473 const unsigned int face_no) const
2474{
2475 const unsigned int n2 = neighbor_of_neighbor_internal(face_no);
2478
2479 return n2;
2480}
2481
2482
2483
2484template <int dim, int spacedim>
2485bool
2487 const unsigned int face_no) const
2488{
2489 return neighbor_of_neighbor_internal(face_no) ==
2491}
2492
2493
2494
2495template <int dim, int spacedim>
2496std::pair<unsigned int, unsigned int>
2498 const unsigned int neighbor) const
2499{
2500 AssertIndexRange(neighbor, this->n_faces());
2501 // make sure that the neighbor is
2502 // on a coarser level
2503 Assert(neighbor_is_coarser(neighbor),
2505
2506 switch (dim)
2507 {
2508 case 2:
2509 {
2510 const int this_face_index = face_index(neighbor);
2511 const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2512 this->neighbor(neighbor);
2513
2514 // usually, on regular patches of
2515 // the grid, this cell is just on
2516 // the opposite side of the
2517 // neighbor that the neighbor is of
2518 // this cell. for example in 2d, if
2519 // we want to know the
2520 // neighbor_of_neighbor if
2521 // neighbor==1 (the right
2522 // neighbor), then we will get 0
2523 // (the left neighbor) in most
2524 // cases. look up this relationship
2525 // in the table provided by
2526 // GeometryInfo and try it
2527 const unsigned int face_no_guess =
2529
2530 const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face_guess =
2531 neighbor_cell->face(face_no_guess);
2532
2533 if (face_guess->has_children())
2534 for (unsigned int subface_no = 0;
2535 subface_no < face_guess->n_children();
2536 ++subface_no)
2537 if (face_guess->child_index(subface_no) == this_face_index)
2538 return std::make_pair(face_no_guess, subface_no);
2539
2540 // if the guess was false, then
2541 // we need to loop over all faces
2542 // and subfaces and find the
2543 // number the hard way
2544 for (const unsigned int face_no : neighbor_cell->face_indices())
2545 {
2546 if (face_no != face_no_guess)
2547 {
2548 const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>
2549 face = neighbor_cell->face(face_no);
2550 if (face->has_children())
2551 for (unsigned int subface_no = 0;
2552 subface_no < face->n_children();
2553 ++subface_no)
2554 if (face->child_index(subface_no) == this_face_index)
2555 return std::make_pair(face_no, subface_no);
2556 }
2557 }
2558
2559 // we should never get here,
2560 // since then we did not find
2561 // our way back...
2563 return std::make_pair(numbers::invalid_unsigned_int,
2565 }
2566
2567 case 3:
2568 {
2569 const int this_face_index = face_index(neighbor);
2570 const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2571 this->neighbor(neighbor);
2572
2573 // usually, on regular patches of the grid, this cell is just on the
2574 // opposite side of the neighbor that the neighbor is of this cell.
2575 // for example in 2d, if we want to know the neighbor_of_neighbor if
2576 // neighbor==1 (the right neighbor), then we will get 0 (the left
2577 // neighbor) in most cases. look up this relationship in the table
2578 // provided by GeometryInfo and try it
2579 const unsigned int face_no_guess =
2581
2582 const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face_guess =
2583 neighbor_cell->face(face_no_guess);
2584
2585 if (face_guess->has_children())
2586 for (unsigned int subface_no = 0;
2587 subface_no < face_guess->n_children();
2588 ++subface_no)
2589 {
2590 if (face_guess->child_index(subface_no) == this_face_index)
2591 // call a helper function, that translates the current
2592 // subface number to a subface number for the current
2593 // FaceRefineCase
2594 return std::make_pair(face_no_guess,
2595 translate_subface_no(face_guess,
2596 subface_no));
2597
2598 if (face_guess->child(subface_no)->has_children())
2599 for (unsigned int subsub_no = 0;
2600 subsub_no < face_guess->child(subface_no)->n_children();
2601 ++subsub_no)
2602 if (face_guess->child(subface_no)->child_index(subsub_no) ==
2603 this_face_index)
2604 // call a helper function, that translates the current
2605 // subface number and subsubface number to a subface
2606 // number for the current FaceRefineCase
2607 return std::make_pair(face_no_guess,
2608 translate_subface_no(face_guess,
2609 subface_no,
2610 subsub_no));
2611 }
2612
2613 // if the guess was false, then we need to loop over all faces and
2614 // subfaces and find the number the hard way
2615 for (const unsigned int face_no : neighbor_cell->face_indices())
2616 {
2617 if (face_no == face_no_guess)
2618 continue;
2619
2620 const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face =
2621 neighbor_cell->face(face_no);
2622
2623 if (!face->has_children())
2624 continue;
2625
2626 for (unsigned int subface_no = 0; subface_no < face->n_children();
2627 ++subface_no)
2628 {
2629 if (face->child_index(subface_no) == this_face_index)
2630 // call a helper function, that translates the current
2631 // subface number to a subface number for the current
2632 // FaceRefineCase
2633 return std::make_pair(face_no,
2634 translate_subface_no(face,
2635 subface_no));
2636
2637 if (face->child(subface_no)->has_children())
2638 for (unsigned int subsub_no = 0;
2639 subsub_no < face->child(subface_no)->n_children();
2640 ++subsub_no)
2641 if (face->child(subface_no)->child_index(subsub_no) ==
2642 this_face_index)
2643 // call a helper function, that translates the current
2644 // subface number and subsubface number to a subface
2645 // number for the current FaceRefineCase
2646 return std::make_pair(face_no,
2647 translate_subface_no(face,
2648 subface_no,
2649 subsub_no));
2650 }
2651 }
2652
2653 // we should never get here, since then we did not find our way
2654 // back...
2656 return std::make_pair(numbers::invalid_unsigned_int,
2658 }
2659
2660 default:
2661 {
2662 Assert(false, ExcImpossibleInDim(1));
2663 return std::make_pair(numbers::invalid_unsigned_int,
2665 }
2666 }
2667}
2668
2669
2670
2671template <int dim, int spacedim>
2672bool
2674 const unsigned int i_face) const
2675{
2676 /*
2677 * Implementation note: In all of the functions corresponding to periodic
2678 * faces we mainly use the Triangulation::periodic_face_map to find the
2679 * information about periodically connected faces. So, we actually search in
2680 * this std::map and return the cell_face on the other side of the periodic
2681 * boundary.
2682 *
2683 * We can not use operator[] as this would insert non-existing entries or
2684 * would require guarding with an extra std::map::find() or count().
2685 */
2686 AssertIndexRange(i_face, this->n_faces());
2687 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2688
2689 cell_iterator current_cell(*this);
2690 if (this->tria->periodic_face_map.find(
2691 std::make_pair(current_cell, i_face)) !=
2692 this->tria->periodic_face_map.end())
2693 return true;
2694 return false;
2695}
2696
2697
2698
2699template <int dim, int spacedim>
2702{
2703 /*
2704 * To know, why we are using std::map::find() instead of [] operator, refer
2705 * to the implementation note in has_periodic_neighbor() function.
2706 *
2707 * my_it : the iterator to the current cell.
2708 * my_face_pair : the pair reported by periodic_face_map as its first pair
2709 * being the current cell_face.
2710 */
2711 AssertIndexRange(i_face, this->n_faces());
2712 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2713 cell_iterator current_cell(*this);
2714
2715 auto my_face_pair =
2716 this->tria->periodic_face_map.find(std::make_pair(current_cell, i_face));
2717
2718 // Make sure we are actually on a periodic boundary:
2719 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2721 return my_face_pair->second.first.first;
2722}
2723
2724
2725
2726template <int dim, int spacedim>
2729 const unsigned int i_face) const
2730{
2731 if (!(this->face(i_face)->at_boundary()))
2732 return this->neighbor(i_face);
2733 else if (this->has_periodic_neighbor(i_face))
2734 return this->periodic_neighbor(i_face);
2735 else
2737 // we can't come here
2738 return this->neighbor(i_face);
2739}
2740
2741
2742
2743template <int dim, int spacedim>
2746 const unsigned int i_face,
2747 const unsigned int i_subface) const
2748{
2749 /*
2750 * To know, why we are using std::map::find() instead of [] operator, refer
2751 * to the implementation note in has_periodic_neighbor() function.
2752 *
2753 * my_it : the iterator to the current cell.
2754 * my_face_pair : the pair reported by periodic_face_map as its first pair
2755 * being the current cell_face. nb_it : the iterator to the
2756 * neighbor of current cell at i_face. face_num_of_nb : the face number of
2757 * the periodically neighboring face in the relevant element.
2758 * nb_parent_face_it: the iterator to the parent face of the periodically
2759 * neighboring face.
2760 */
2761 AssertIndexRange(i_face, this->n_faces());
2762 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2763 cell_iterator my_it(*this);
2764
2765 auto my_face_pair =
2766 this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2767 /*
2768 * There should be an assertion, which tells the user that this function
2769 * should not be used for a cell which is not located at a periodic boundary.
2770 */
2771 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2773 cell_iterator parent_nb_it = my_face_pair->second.first.first;
2774 unsigned int nb_face_num = my_face_pair->second.first.second;
2775 TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> nb_parent_face_it =
2776 parent_nb_it->face(nb_face_num);
2777 /*
2778 * We should check if the parent face of the neighbor has at least the same
2779 * number of children as i_subface.
2780 */
2781 AssertIndexRange(i_subface, nb_parent_face_it->n_children());
2782
2783 const auto [orientation, rotation, flip] =
2784 internal::split_face_orientation(my_face_pair->second.second);
2785
2786 unsigned int sub_neighbor_num =
2787 GeometryInfo<dim>::child_cell_on_face(parent_nb_it->refinement_case(),
2788 nb_face_num,
2789 i_subface,
2790 orientation,
2791 flip,
2792 rotation,
2793 nb_parent_face_it->refinement_case());
2794 return parent_nb_it->child(sub_neighbor_num);
2795}
2796
2797
2798
2799template <int dim, int spacedim>
2800std::pair<unsigned int, unsigned int>
2802 const unsigned int i_face) const
2803{
2804 /*
2805 * To know, why we are using std::map::find() instead of [] operator, refer
2806 * to the implementation note in has_periodic_neighbor() function.
2807 *
2808 * my_it : the iterator to the current cell.
2809 * my_face_pair : the pair reported by periodic_face_map as its first pair
2810 * being the current cell_face. nb_it : the iterator to the periodic
2811 * neighbor. nb_face_pair : the pair reported by periodic_face_map as its
2812 * first pair being the periodic neighbor cell_face. p_nb_of_p_nb : the
2813 * iterator of the periodic neighbor of the periodic neighbor of the current
2814 * cell.
2815 */
2816 AssertIndexRange(i_face, this->n_faces());
2817 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2818 const int my_face_index = this->face_index(i_face);
2819 cell_iterator my_it(*this);
2820
2821 auto my_face_pair =
2822 this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2823 /*
2824 * There should be an assertion, which tells the user that this function
2825 * should not be used for a cell which is not located at a periodic boundary.
2826 */
2827 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2829 cell_iterator nb_it = my_face_pair->second.first.first;
2830 unsigned int face_num_of_nb = my_face_pair->second.first.second;
2831
2832 auto nb_face_pair =
2833 this->tria->periodic_face_map.find(std::make_pair(nb_it, face_num_of_nb));
2834 /*
2835 * Since, we store periodic neighbors for every cell (either active or
2836 * artificial or inactive) the nb_face_pair should also be mapped to some
2837 * cell_face pair. We assert this here.
2838 */
2839 Assert(nb_face_pair != this->tria->periodic_face_map.end(),
2841 cell_iterator p_nb_of_p_nb = nb_face_pair->second.first.first;
2842 TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> parent_face_it =
2843 p_nb_of_p_nb->face(nb_face_pair->second.first.second);
2844 for (unsigned int i_subface = 0; i_subface < parent_face_it->n_children();
2845 ++i_subface)
2846 if (parent_face_it->child_index(i_subface) == my_face_index)
2847 return std::make_pair(face_num_of_nb, i_subface);
2848 /*
2849 * Obviously, if the execution reaches to this point, some of our assumptions
2850 * should have been false. The most important one is, the user has called this
2851 * function on a face which does not have a coarser periodic neighbor.
2852 */
2854 return std::make_pair(numbers::invalid_unsigned_int,
2856}
2857
2858
2859
2860template <int dim, int spacedim>
2861int
2863 const unsigned int i_face) const
2864{
2865 return periodic_neighbor(i_face)->index();
2866}
2867
2868
2869
2870template <int dim, int spacedim>
2871int
2873 const unsigned int i_face) const
2874{
2875 return periodic_neighbor(i_face)->level();
2876}
2877
2878
2879
2880template <int dim, int spacedim>
2881unsigned int
2883 const unsigned int i_face) const
2884{
2885 return periodic_neighbor_face_no(i_face);
2886}
2887
2888
2889
2890template <int dim, int spacedim>
2891unsigned int
2893 const unsigned int i_face) const
2894{
2895 /*
2896 * To know, why we are using std::map::find() instead of [] operator, refer
2897 * to the implementation note in has_periodic_neighbor() function.
2898 *
2899 * my_it : the iterator to the current cell.
2900 * my_face_pair : the pair reported by periodic_face_map as its first pair
2901 * being the current cell_face.
2902 */
2903 AssertIndexRange(i_face, this->n_faces());
2904 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2905 cell_iterator my_it(*this);
2906
2907 auto my_face_pair =
2908 this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2909 /*
2910 * There should be an assertion, which tells the user that this function
2911 * should not be called for a cell which is not located at a periodic boundary
2912 * !
2913 */
2914 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2916 return my_face_pair->second.first.second;
2917}
2918
2919
2920
2921template <int dim, int spacedim>
2922bool
2924 const unsigned int i_face) const
2925{
2926 /*
2927 * To know, why we are using std::map::find() instead of [] operator, refer
2928 * to the implementation note in has_periodic_neighbor() function.
2929 *
2930 * Implementation note: Let p_nb_of_p_nb be the periodic neighbor of the
2931 * periodic neighbor of the current cell. Also, let p_face_of_p_nb_of_p_nb be
2932 * the periodic face of the p_nb_of_p_nb. If p_face_of_p_nb_of_p_nb has
2933 * children , then the periodic neighbor of the current cell is coarser than
2934 * itself. Although not tested, this implementation should work for
2935 * anisotropic refinement as well.
2936 *
2937 * my_it : the iterator to the current cell.
2938 * my_face_pair : the pair reported by periodic_face_map as its first pair
2939 * being the current cell_face. nb_it : the iterator to the periodic
2940 * neighbor. nb_face_pair : the pair reported by periodic_face_map as its
2941 * first pair being the periodic neighbor cell_face.
2942 */
2943 AssertIndexRange(i_face, this->n_faces());
2944 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2945 cell_iterator my_it(*this);
2946
2947 auto my_face_pair =
2948 this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2949 /*
2950 * There should be an assertion, which tells the user that this function
2951 * should not be used for a cell which is not located at a periodic boundary.
2952 */
2953 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2955
2956 cell_iterator nb_it = my_face_pair->second.first.first;
2957 unsigned int face_num_of_nb = my_face_pair->second.first.second;
2958
2959 auto nb_face_pair =
2960 this->tria->periodic_face_map.find(std::make_pair(nb_it, face_num_of_nb));
2961 /*
2962 * Since, we store periodic neighbors for every cell (either active or
2963 * artificial or inactive) the nb_face_pair should also be mapped to some
2964 * cell_face pair. We assert this here.
2965 */
2966 Assert(nb_face_pair != this->tria->periodic_face_map.end(),
2968 const unsigned int my_level = this->level();
2969 const unsigned int neighbor_level = nb_face_pair->second.first.first->level();
2970 Assert(my_level >= neighbor_level, ExcInternalError());
2971 return my_level > neighbor_level;
2972}
2973
2974
2975
2976template <int dim, int spacedim>
2977bool
2979{
2981 AssertIndexRange(i, this->n_faces());
2982
2983 return (neighbor_index(i) == -1);
2984}
2985
2986
2987
2988template <int dim, int spacedim>
2989bool
2991{
2992 if (dim == 1)
2993 return at_boundary();
2994 else
2995 {
2996 for (unsigned int l = 0; l < this->n_lines(); ++l)
2997 if (this->line(l)->at_boundary())
2998 return true;
2999
3000 return false;
3001 }
3002}
3003
3004
3005
3006template <int dim, int spacedim>
3009 const unsigned int face,
3010 const unsigned int subface) const
3011{
3012 Assert(!this->has_children(),
3013 ExcMessage("The present cell must not have children!"));
3014 Assert(!this->at_boundary(face),
3015 ExcMessage("The present cell must have a valid neighbor!"));
3016 Assert(this->neighbor(face)->has_children() == true,
3017 ExcMessage("The neighbor must have children!"));
3018
3019 switch (dim)
3020 {
3021 case 2:
3022 {
3023 if (this->reference_cell() == ReferenceCells::Triangle)
3024 {
3025 const unsigned int neighbor_neighbor =
3026 this->neighbor_of_neighbor(face);
3027 const auto neighbor = this->neighbor(face);
3028 // Triangles do not support anisotropic refinement
3029 Assert(neighbor->refinement_case() ==
3032 // Since we are looking at two faces at once, we need to check if
3033 // they have the same or opposing orientations rather than one
3034 // individual face orientation value.
3035 const auto relative_orientation =
3036 neighbor->combined_face_orientation(neighbor_neighbor) ==
3037 this->combined_face_orientation(face) ?
3040 const unsigned int neighbor_child_index =
3041 neighbor->reference_cell().child_cell_on_face(
3042 neighbor_neighbor, subface, relative_orientation);
3043 auto child = neighbor->child(neighbor_child_index);
3044 Assert(!child->has_children(), ExcInternalError());
3045 return child;
3046 }
3047 else if (this->reference_cell() == ReferenceCells::Quadrilateral)
3048 {
3049 const unsigned int neighbor_neighbor =
3050 this->neighbor_of_neighbor(face);
3051 const unsigned int neighbor_child_index =
3053 this->neighbor(face)->refinement_case(),
3054 neighbor_neighbor,
3055 subface);
3056
3058 this->neighbor(face)->child(neighbor_child_index);
3059 // the neighbors child can have children,
3060 // which are not further refined along the
3061 // face under consideration. as we are
3062 // normally interested in one of this
3063 // child's child, search for the right one.
3064 while (sub_neighbor->has_children())
3065 {
3067 sub_neighbor->refinement_case(),
3068 neighbor_neighbor) ==
3071 sub_neighbor =
3072 sub_neighbor->child(GeometryInfo<dim>::child_cell_on_face(
3073 sub_neighbor->refinement_case(), neighbor_neighbor, 0));
3074 }
3075
3076 return sub_neighbor;
3077 }
3078
3079 // if no reference cell type matches
3082 }
3083
3084
3085 case 3:
3086 {
3087 if (this->reference_cell() == ReferenceCells::Hexahedron)
3088 {
3089 // this function returns the neighbor's
3090 // child on a given face and
3091 // subface.
3092
3093 // we have to consider one other aspect here:
3094 // The face might be refined
3095 // anisotropically. In this case, the subface
3096 // number refers to the following, where we
3097 // look at the face from the current cell,
3098 // thus the subfaces are in standard
3099 // orientation concerning the cell
3100 //
3101 // for isotropic refinement
3102 //
3103 // *---*---*
3104 // | 2 | 3 |
3105 // *---*---*
3106 // | 0 | 1 |
3107 // *---*---*
3108 //
3109 // for 2*anisotropic refinement
3110 // (first cut_y, then cut_x)
3111 //
3112 // *---*---*
3113 // | 2 | 3 |
3114 // *---*---*
3115 // | 0 | 1 |
3116 // *---*---*
3117 //
3118 // for 2*anisotropic refinement
3119 // (first cut_x, then cut_y)
3120 //
3121 // *---*---*
3122 // | 1 | 3 |
3123 // *---*---*
3124 // | 0 | 2 |
3125 // *---*---*
3126 //
3127 // for purely anisotropic refinement:
3128 //
3129 // *---*---* *-------*
3130 // | | | | 1 |
3131 // | 0 | 1 | or *-------*
3132 // | | | | 0 |
3133 // *---*---* *-------*
3134 //
3135 // for "mixed" refinement:
3136 //
3137 // *---*---* *---*---* *---*---* *-------*
3138 // | | 2 | | 1 | | | 1 | 2 | | 2 |
3139 // | 0 *---* or *---* 2 | or *---*---* or *---*---*
3140 // | | 1 | | 0 | | | 0 | | 0 | 1 |
3141 // *---*---* *---*---* *-------* *---*---*
3142
3144 mother_face = this->face(face);
3145 const unsigned int total_children =
3146 mother_face->n_active_descendants();
3147 AssertIndexRange(subface, total_children);
3150
3151 unsigned int neighbor_neighbor;
3154 this->neighbor(face);
3155
3156
3157 const RefinementCase<dim - 1> mother_face_ref_case =
3158 mother_face->refinement_case();
3159 if (mother_face_ref_case ==
3160 static_cast<RefinementCase<dim - 1>>(
3161 RefinementCase<2>::cut_xy)) // total_children==4
3162 {
3163 // this case is quite easy. we are sure,
3164 // that the neighbor is not coarser.
3165
3166 // get the neighbor's number for the given
3167 // face and the neighbor
3168 neighbor_neighbor = this->neighbor_of_neighbor(face);
3169
3170 // now use the info provided by GeometryInfo
3171 // to extract the neighbors child number
3172 const unsigned int neighbor_child_index =
3174 neighbor->refinement_case(),
3175 neighbor_neighbor,
3176 subface,
3177 neighbor->face_orientation(neighbor_neighbor),
3178 neighbor->face_flip(neighbor_neighbor),
3179 neighbor->face_rotation(neighbor_neighbor));
3180 neighbor_child = neighbor->child(neighbor_child_index);
3181
3182 // make sure that the neighbor child cell we
3183 // have found shares the desired subface.
3184 Assert((this->face(face)->child(subface) ==
3185 neighbor_child->face(neighbor_neighbor)),
3187 }
3188 else //-> the face is refined anisotropically
3189 {
3190 // first of all, we have to find the
3191 // neighbor at one of the anisotropic
3192 // children of the
3193 // mother_face. determine, which of
3194 // these we need.
3195 unsigned int first_child_to_find;
3196 unsigned int neighbor_child_index;
3197 if (total_children == 2)
3198 first_child_to_find = subface;
3199 else
3200 {
3201 first_child_to_find = subface / 2;
3202 if (total_children == 3 && subface == 1 &&
3203 !mother_face->child(0)->has_children())
3204 first_child_to_find = 1;
3205 }
3206 if (neighbor_is_coarser(face))
3207 {
3208 std::pair<unsigned int, unsigned int> indices =
3209 neighbor_of_coarser_neighbor(face);
3210 neighbor_neighbor = indices.first;
3211
3212
3213 // we have to translate our
3214 // subface_index according to the
3215 // RefineCase and subface index of
3216 // the coarser face (our face is an
3217 // anisotropic child of the coarser
3218 // face), 'a' denotes our
3219 // subface_index 0 and 'b' denotes
3220 // our subface_index 1, whereas 0...3
3221 // denote isotropic subfaces of the
3222 // coarser face
3223 //
3224 // cut_x and coarser_subface_index=0
3225 //
3226 // *---*---*
3227 // |b=2| |
3228 // | | |
3229 // |a=0| |
3230 // *---*---*
3231 //
3232 // cut_x and coarser_subface_index=1
3233 //
3234 // *---*---*
3235 // | |b=3|
3236 // | | |
3237 // | |a=1|
3238 // *---*---*
3239 //
3240 // cut_y and coarser_subface_index=0
3241 //
3242 // *-------*
3243 // | |
3244 // *-------*
3245 // |a=0 b=1|
3246 // *-------*
3247 //
3248 // cut_y and coarser_subface_index=1
3249 //
3250 // *-------*
3251 // |a=2 b=3|
3252 // *-------*
3253 // | |
3254 // *-------*
3255 unsigned int iso_subface;
3256 if (neighbor->face(neighbor_neighbor)
3257 ->refinement_case() == RefinementCase<2>::cut_x)
3258 iso_subface = 2 * first_child_to_find + indices.second;
3259 else
3260 {
3261 Assert(neighbor->face(neighbor_neighbor)
3262 ->refinement_case() ==
3265 iso_subface =
3266 first_child_to_find + 2 * indices.second;
3267 }
3268 neighbor_child_index =
3270 neighbor->refinement_case(),
3271 neighbor_neighbor,
3272 iso_subface,
3273 neighbor->face_orientation(neighbor_neighbor),
3274 neighbor->face_flip(neighbor_neighbor),
3275 neighbor->face_rotation(neighbor_neighbor));
3276 }
3277 else // neighbor is not coarser
3278 {
3279 neighbor_neighbor = neighbor_of_neighbor(face);
3280 neighbor_child_index =
3282 neighbor->refinement_case(),
3283 neighbor_neighbor,
3284 first_child_to_find,
3285 neighbor->face_orientation(neighbor_neighbor),
3286 neighbor->face_flip(neighbor_neighbor),
3287 neighbor->face_rotation(neighbor_neighbor),
3288 mother_face_ref_case);
3289 }
3290
3291 neighbor_child = neighbor->child(neighbor_child_index);
3292 // it might be, that the neighbor_child
3293 // has children, which are not refined
3294 // along the given subface. go down that
3295 // list and deliver the last of those.
3296 while (
3297 neighbor_child->has_children() &&
3299 neighbor_child->refinement_case(), neighbor_neighbor) ==
3301 neighbor_child = neighbor_child->child(
3303 neighbor_child->refinement_case(),
3304 neighbor_neighbor,
3305 0));
3306
3307 // if there are two total subfaces, we
3308 // are finished. if there are four we
3309 // have to get a child of our current
3310 // neighbor_child. If there are three,
3311 // we have to check which of the two
3312 // possibilities applies.
3313 if (total_children == 3)
3314 {
3315 if (mother_face->child(0)->has_children())
3316 {
3317 if (subface < 2)
3318 neighbor_child = neighbor_child->child(
3320 neighbor_child->refinement_case(),
3321 neighbor_neighbor,
3322 subface,
3323 neighbor_child->face_orientation(
3324 neighbor_neighbor),
3325 neighbor_child->face_flip(neighbor_neighbor),
3326 neighbor_child->face_rotation(
3327 neighbor_neighbor),
3328 mother_face->child(0)->refinement_case()));
3329 }
3330 else
3331 {
3332 Assert(mother_face->child(1)->has_children(),
3334 if (subface > 0)
3335 neighbor_child = neighbor_child->child(
3337 neighbor_child->refinement_case(),
3338 neighbor_neighbor,
3339 subface - 1,
3340 neighbor_child->face_orientation(
3341 neighbor_neighbor),
3342 neighbor_child->face_flip(neighbor_neighbor),
3343 neighbor_child->face_rotation(
3344 neighbor_neighbor),
3345 mother_face->child(1)->refinement_case()));
3346 }
3347 }
3348 else if (total_children == 4)
3349 {
3350 neighbor_child = neighbor_child->child(
3352 neighbor_child->refinement_case(),
3353 neighbor_neighbor,
3354 subface % 2,
3355 neighbor_child->face_orientation(neighbor_neighbor),
3356 neighbor_child->face_flip(neighbor_neighbor),
3357 neighbor_child->face_rotation(neighbor_neighbor),
3358 mother_face->child(subface / 2)->refinement_case()));
3359 }
3360 }
3361
3362 // it might be, that the neighbor_child has
3363 // children, which are not refined along the
3364 // given subface. go down that list and
3365 // deliver the last of those.
3366 while (neighbor_child->has_children())
3367 neighbor_child =
3368 neighbor_child->child(GeometryInfo<dim>::child_cell_on_face(
3369 neighbor_child->refinement_case(), neighbor_neighbor, 0));
3370
3371 if constexpr (running_in_debug_mode())
3372 {
3373 // check, whether the face neighbor_child matches the
3374 // requested subface.
3376 requested;
3377 switch (this->subface_case(face))
3378 {
3382 requested = mother_face->child(subface);
3383 break;
3386 requested =
3387 mother_face->child(subface / 2)->child(subface % 2);
3388 break;
3389
3392 switch (subface)
3393 {
3394 case 0:
3395 case 1:
3396 requested = mother_face->child(0)->child(subface);
3397 break;
3398 case 2:
3399 requested = mother_face->child(1);
3400 break;
3401 default:
3403 }
3404 break;
3407 switch (subface)
3408 {
3409 case 0:
3410 requested = mother_face->child(0);
3411 break;
3412 case 1:
3413 case 2:
3414 requested =
3415 mother_face->child(1)->child(subface - 1);
3416 break;
3417 default:
3419 }
3420 break;
3421 default:
3423 break;
3424 }
3425 Assert(requested == neighbor_child->face(neighbor_neighbor),
3427 }
3428
3429 return neighbor_child;
3430 }
3431
3432 // if no reference cell type matches
3435 }
3436
3437 default:
3438 // if 1d or more than 3d
3441 }
3442}
3443
3444
3445
3446template <int structdim, int dim, int spacedim>
3452
3453
3454
3455template <int structdim, int dim, int spacedim>
3456int
3461
3462
3463
3464template <int structdim, int dim, int spacedim>
3465int
3470
3471
3472// explicit instantiations
3473#include "grid/tria_accessor.inst"
3474
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:232
constexpr types::geometric_orientation reverse_line_orientation
Definition types.h:359
constexpr types::material_id invalid_material_id
Definition types.h:288
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:346
::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)