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