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