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