Reference documentation for deal.II version Git 3f1f337db3 2021-10-23 13:19:02 -0600
\(\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  {
1328  {
1329  // If the face is planar, the diagonal from vertex 0 to vertex 3,
1330  // v_03, should be in the plane P_012 of vertices 0, 1 and 2. Get
1331  // the normal vector of P_012 and test if v_03 is orthogonal to
1332  // that. If so, the face is planar and computing its area is simple.
1333  const Tensor<1, 3> v01 = accessor.vertex(1) - accessor.vertex(0);
1334  const Tensor<1, 3> v02 = accessor.vertex(2) - accessor.vertex(0);
1335 
1336  const Tensor<1, 3> normal = cross_product_3d(v01, v02);
1337 
1338  const Tensor<1, 3> v03 = accessor.vertex(3) - accessor.vertex(0);
1339 
1340  // check whether v03 does not lie in the plane of v01 and v02
1341  // (i.e., whether the face is not planar). we do so by checking
1342  // whether the triple product (v01 x v02) * v03 forms a positive
1343  // volume relative to |v01|*|v02|*|v03|. the test checks the
1344  // squares of these to avoid taking norms/square roots:
1345  if (std::abs((v03 * normal) * (v03 * normal) /
1346  ((v03 * v03) * (v01 * v01) * (v02 * v02))) >= 1e-24)
1347  {
1348  // If the vectors are non planar we integrate the norm of the normal
1349  // vector using a numerical Gauss scheme of order 4. In particular
1350  // we consider a bilinear quad x(u,v) = (1-v)((1-u)v_0 + u v_1) +
1351  // v((1-u)v_2 + u v_3), consequently we compute the normal vector as
1352  // n(u,v) = t_u x t_v = w_1 + u w_2 + v w_3. The integrand function
1353  // is
1354  // || n(u,v) || = sqrt(a + b u^2 + c v^2 + d u + e v + f uv).
1355  // We integrate it using a QGauss<2> (4) computed explicitly.
1356  const Tensor<1, 3> w_1 =
1357  cross_product_3d(accessor.vertex(1) - accessor.vertex(0),
1358  accessor.vertex(2) - accessor.vertex(0));
1359  const Tensor<1, 3> w_2 =
1360  cross_product_3d(accessor.vertex(1) - accessor.vertex(0),
1361  accessor.vertex(3) - accessor.vertex(2) -
1362  accessor.vertex(1) + accessor.vertex(0));
1363  const Tensor<1, 3> w_3 =
1364  cross_product_3d(accessor.vertex(3) - accessor.vertex(2) -
1365  accessor.vertex(1) + accessor.vertex(0),
1366  accessor.vertex(2) - accessor.vertex(0));
1367 
1368  double a = scalar_product(w_1, w_1);
1369  double b = scalar_product(w_2, w_2);
1370  double c = scalar_product(w_3, w_3);
1371  double d = scalar_product(w_1, w_2);
1372  double e = scalar_product(w_1, w_3);
1373  double f = scalar_product(w_2, w_3);
1374 
1375  return 0.03025074832140047 *
1376  std::sqrt(
1377  a + 0.0048207809894260144 * b +
1378  0.0048207809894260144 * c + 0.13886368840594743 * d +
1379  0.13886368840594743 * e + 0.0096415619788520288 * f) +
1380  0.056712962962962937 *
1381  std::sqrt(
1382  a + 0.0048207809894260144 * b + 0.10890625570683385 * c +
1383  0.13886368840594743 * d + 0.66001895641514374 * e +
1384  0.045826333352825557 * f) +
1385  0.056712962962962937 *
1386  std::sqrt(
1387  a + 0.0048207809894260144 * b + 0.44888729929169013 * c +
1388  0.13886368840594743 * d + 1.3399810435848563 * e +
1389  0.09303735505312187 * f) +
1390  0.03025074832140047 *
1391  std::sqrt(
1392  a + 0.0048207809894260144 * b + 0.86595709258347853 * c +
1393  0.13886368840594743 * d + 1.8611363115940525 * e +
1394  0.12922212642709538 * f) +
1395  0.056712962962962937 *
1396  std::sqrt(
1397  a + 0.10890625570683385 * b + 0.0048207809894260144 * c +
1398  0.66001895641514374 * d + 0.13886368840594743 * e +
1399  0.045826333352825557 * f) +
1400  0.10632332575267359 * std::sqrt(a + 0.10890625570683385 * b +
1401  0.10890625570683385 * c +
1402  0.66001895641514374 * d +
1403  0.66001895641514374 * e +
1404  0.2178125114136677 * f) +
1405  0.10632332575267359 * std::sqrt(a + 0.10890625570683385 * b +
1406  0.44888729929169013 * c +
1407  0.66001895641514374 * d +
1408  1.3399810435848563 * e +
1409  0.44220644500147605 * f) +
1410  0.056712962962962937 *
1411  std::sqrt(
1412  a + 0.10890625570683385 * b + 0.86595709258347853 * c +
1413  0.66001895641514374 * d + 1.8611363115940525 * e +
1414  0.61419262306231814 * f) +
1415  0.056712962962962937 *
1416  std::sqrt(
1417  a + 0.44888729929169013 * b + 0.0048207809894260144 * c +
1418  1.3399810435848563 * d + 0.13886368840594743 * e +
1419  0.09303735505312187 * f) +
1420  0.10632332575267359 * std::sqrt(a + 0.44888729929169013 * b +
1421  0.10890625570683385 * c +
1422  1.3399810435848563 * d +
1423  0.66001895641514374 * e +
1424  0.44220644500147605 * f) +
1425  0.10632332575267359 *
1426  std::sqrt(a + 0.44888729929169013 * b +
1427  0.44888729929169013 * c +
1428  1.3399810435848563 * d + 1.3399810435848563 * e +
1429  0.89777459858338027 * f) +
1430  0.056712962962962937 *
1431  std::sqrt(a + 0.44888729929169013 * b +
1432  0.86595709258347853 * c +
1433  1.3399810435848563 * d + 1.8611363115940525 * e +
1434  1.2469436885317342 * f) +
1435  0.03025074832140047 * std::sqrt(a + 0.86595709258347853 * b +
1436  0.0048207809894260144 * c +
1437  1.8611363115940525 * d +
1438  0.13886368840594743 * e +
1439  0.12922212642709538 * f) +
1440  0.056712962962962937 *
1441  std::sqrt(
1442  a + 0.86595709258347853 * b + 0.10890625570683385 * c +
1443  1.8611363115940525 * d + 0.66001895641514374 * e +
1444  0.61419262306231814 * f) +
1445  0.056712962962962937 *
1446  std::sqrt(a + 0.86595709258347853 * b +
1447  0.44888729929169013 * c +
1448  1.8611363115940525 * d + 1.3399810435848563 * e +
1449  1.2469436885317342 * f) +
1450  0.03025074832140047 *
1451  std::sqrt(a + 0.86595709258347853 * b +
1452  0.86595709258347853 * c +
1453  1.8611363115940525 * d + 1.8611363115940525 * e +
1454  1.7319141851669571 * f);
1455  }
1456 
1457  // the face is planar. then its area is 1/2 of the norm of the
1458  // cross product of the two diagonals
1459  const Tensor<1, 3> v12 = accessor.vertex(2) - accessor.vertex(1);
1460  const Tensor<1, 3> twice_area = cross_product_3d(v03, v12);
1461  return 0.5 * twice_area.norm();
1462  }
1463  else if (accessor.reference_cell() == ReferenceCells::Triangle)
1464  {
1465  // We can just use the normal triangle area formula without issue
1466  const Tensor<1, 3> v01 = accessor.vertex(1) - accessor.vertex(0);
1467  const Tensor<1, 3> v02 = accessor.vertex(2) - accessor.vertex(0);
1468  return 0.5 * cross_product_3d(v01, v02).norm();
1469  }
1470 
1471  Assert(false, ExcNotImplemented());
1472  return 0.0;
1473  }
1474 
1475 
1476 
1477  template <int structdim, int dim, int spacedim>
1478  double
1480  {
1481  // catch-all for all cases not explicitly
1482  // listed above
1483  Assert(false, ExcNotImplemented());
1484  return std::numeric_limits<double>::quiet_NaN();
1485  }
1486 
1487 
1488  template <int dim, int spacedim>
1490  get_new_point_on_object(const TriaAccessor<1, dim, spacedim> &obj)
1491  {
1493  return obj.get_manifold().get_new_point_on_line(it);
1494  }
1495 
1496  template <int dim, int spacedim>
1498  get_new_point_on_object(const TriaAccessor<2, dim, spacedim> &obj)
1499  {
1501  return obj.get_manifold().get_new_point_on_quad(it);
1502  }
1503 
1504  template <int dim, int spacedim>
1506  get_new_point_on_object(const TriaAccessor<3, dim, spacedim> &obj)
1507  {
1509  return obj.get_manifold().get_new_point_on_hex(it);
1510  }
1511 
1512  template <int structdim, int dim, int spacedim>
1514  get_new_point_on_object(const TriaAccessor<structdim, dim, spacedim> &obj,
1515  const bool use_interpolation)
1516  {
1517  if (use_interpolation)
1518  {
1520  const auto points_and_weights =
1521  Manifolds::get_default_points_and_weights(it, use_interpolation);
1522  return obj.get_manifold().get_new_point(
1523  make_array_view(points_and_weights.first.begin(),
1524  points_and_weights.first.end()),
1525  make_array_view(points_and_weights.second.begin(),
1526  points_and_weights.second.end()));
1527  }
1528  else
1529  {
1530  return get_new_point_on_object(obj);
1531  }
1532  }
1533 } // namespace
1534 
1535 
1536 
1537 /*-------------------- Static variables: TriaAccessorBase -------------------*/
1538 
1539 template <int structdim, int dim, int spacedim>
1541 
1542 template <int structdim, int dim, int spacedim>
1544 
1545 template <int structdim, int dim, int spacedim>
1546 const unsigned int
1548 
1549 
1550 /*------------------------ Functions: TriaAccessor ---------------------------*/
1551 
1552 template <int structdim, int dim, int spacedim>
1553 void
1555  const std::initializer_list<int> &new_indices) const
1556 {
1557  const ArrayView<int> bounding_object_index_ref =
1558  this->objects().get_bounding_object_indices(this->present_index);
1559 
1560  AssertIndexRange(new_indices.size(), bounding_object_index_ref.size() + 1);
1561 
1562  unsigned int i = 0;
1563  for (const auto &new_index : new_indices)
1564  {
1565  bounding_object_index_ref[i] = new_index;
1566  ++i;
1567  }
1568 }
1569 
1570 
1571 
1572 template <int structdim, int dim, int spacedim>
1573 void
1575  const std::initializer_list<unsigned int> &new_indices) const
1576 {
1577  const ArrayView<int> bounding_object_index_ref =
1578  this->objects().get_bounding_object_indices(this->present_index);
1579 
1580  AssertIndexRange(new_indices.size(), bounding_object_index_ref.size() + 1);
1581 
1582  unsigned int i = 0;
1583  for (const auto &new_index : new_indices)
1584  {
1585  bounding_object_index_ref[i] = new_index;
1586  ++i;
1587  }
1588 }
1589 
1590 
1591 
1592 template <int structdim, int dim, int spacedim>
1595 {
1596  // call the function in the anonymous
1597  // namespace above
1598  return ::barycenter(*this);
1599 }
1600 
1601 
1602 
1603 template <int structdim, int dim, int spacedim>
1604 double
1606 {
1607  // call the function in the anonymous
1608  // namespace above
1609  return ::measure(*this);
1610 }
1611 
1612 
1613 
1614 template <int structdim, int dim, int spacedim>
1617 {
1618  std::pair<Point<spacedim>, Point<spacedim>> boundary_points =
1619  std::make_pair(this->vertex(0), this->vertex(0));
1620 
1621  for (unsigned int v = 1; v < this->n_vertices(); ++v)
1622  {
1623  const Point<spacedim> &x = this->vertex(v);
1624  for (unsigned int k = 0; k < spacedim; ++k)
1625  {
1626  boundary_points.first[k] = std::min(boundary_points.first[k], x[k]);
1627  boundary_points.second[k] = std::max(boundary_points.second[k], x[k]);
1628  }
1629  }
1630 
1631  return BoundingBox<spacedim>(boundary_points);
1632 }
1633 
1634 
1635 
1636 template <int structdim, int dim, int spacedim>
1637 double
1639  const unsigned int /*axis*/) const
1640 {
1641  Assert(false, ExcNotImplemented());
1642  return std::numeric_limits<double>::signaling_NaN();
1643 }
1644 
1645 
1646 
1647 template <>
1648 double
1649 TriaAccessor<1, 1, 1>::extent_in_direction(const unsigned int axis) const
1650 {
1651  (void)axis;
1652  AssertIndexRange(axis, 1);
1653 
1654  return this->diameter();
1655 }
1656 
1657 
1658 template <>
1659 double
1660 TriaAccessor<1, 1, 2>::extent_in_direction(const unsigned int axis) const
1661 {
1662  (void)axis;
1663  AssertIndexRange(axis, 1);
1664 
1665  return this->diameter();
1666 }
1667 
1668 
1669 template <>
1670 double
1671 TriaAccessor<2, 2, 2>::extent_in_direction(const unsigned int axis) const
1672 {
1673  const unsigned int lines[2][2] = {
1674  {2, 3},
1675  {0, 1}};
1676 
1677  AssertIndexRange(axis, 2);
1678 
1679  return std::max(this->line(lines[axis][0])->diameter(),
1680  this->line(lines[axis][1])->diameter());
1681 }
1682 
1683 template <>
1684 double
1685 TriaAccessor<2, 2, 3>::extent_in_direction(const unsigned int axis) const
1686 {
1687  const unsigned int lines[2][2] = {
1688  {2, 3},
1689  {0, 1}};
1690 
1691  AssertIndexRange(axis, 2);
1692 
1693  return std::max(this->line(lines[axis][0])->diameter(),
1694  this->line(lines[axis][1])->diameter());
1695 }
1696 
1697 
1698 template <>
1699 double
1700 TriaAccessor<3, 3, 3>::extent_in_direction(const unsigned int axis) const
1701 {
1702  const unsigned int lines[3][4] = {
1703  {2, 3, 6, 7},
1704  {0, 1, 4, 5},
1705  {8, 9, 10, 11}};
1706 
1707  AssertIndexRange(axis, 3);
1708 
1709  double lengths[4] = {this->line(lines[axis][0])->diameter(),
1710  this->line(lines[axis][1])->diameter(),
1711  this->line(lines[axis][2])->diameter(),
1712  this->line(lines[axis][3])->diameter()};
1713 
1714  return std::max(std::max(lengths[0], lengths[1]),
1715  std::max(lengths[2], lengths[3]));
1716 }
1717 
1718 
1719 // Recursively set manifold ids on hex iterators.
1720 template <>
1721 void
1723  const types::manifold_id manifold_ind) const
1724 {
1725  set_manifold_id(manifold_ind);
1726 
1727  if (this->has_children())
1728  for (unsigned int c = 0; c < this->n_children(); ++c)
1729  this->child(c)->set_all_manifold_ids(manifold_ind);
1730 
1731  // for hexes also set manifold_id
1732  // of bounding quads and lines
1733 
1734  for (unsigned int i : this->face_indices())
1735  this->quad(i)->set_manifold_id(manifold_ind);
1736  for (unsigned int i : this->line_indices())
1737  this->line(i)->set_manifold_id(manifold_ind);
1738 }
1739 
1740 
1741 template <int structdim, int dim, int spacedim>
1744  const Point<structdim> &coordinates) const
1745 {
1746  // Surrounding points and weights.
1747  std::array<Point<spacedim>, GeometryInfo<structdim>::vertices_per_cell> p;
1748  std::array<double, GeometryInfo<structdim>::vertices_per_cell> w;
1749 
1750  for (const unsigned int i : this->vertex_indices())
1751  {
1752  p[i] = this->vertex(i);
1753  w[i] = GeometryInfo<structdim>::d_linear_shape_function(coordinates, i);
1754  }
1755 
1756  return this->get_manifold().get_new_point(make_array_view(p.begin(), p.end()),
1757  make_array_view(w.begin(),
1758  w.end()));
1759 }
1760 
1761 
1762 
1763 template <int structdim, int dim, int spacedim>
1766  const Point<spacedim> &point) const
1767 {
1768  std::array<Point<spacedim>, GeometryInfo<structdim>::vertices_per_cell>
1769  vertices;
1770  for (const unsigned int v : this->vertex_indices())
1771  vertices[v] = this->vertex(v);
1772 
1773  const auto A_b =
1774  GridTools::affine_cell_approximation<structdim, spacedim>(vertices);
1776  A_b.first.covariant_form().transpose();
1777  return Point<structdim>(apply_transformation(A_inv, point - A_b.second));
1778 }
1779 
1780 
1781 
1782 template <int structdim, int dim, int spacedim>
1785  const bool respect_manifold,
1786  const bool use_interpolation) const
1787 {
1788  if (respect_manifold == false)
1789  {
1790  Assert(use_interpolation == false, ExcNotImplemented());
1791  Point<spacedim> p;
1792  for (const unsigned int v : this->vertex_indices())
1793  p += vertex(v);
1794  return p / this->n_vertices();
1795  }
1796  else
1797  return get_new_point_on_object(*this, use_interpolation);
1798 }
1799 
1800 
1801 /*------------------------ Functions: CellAccessor<1> -----------------------*/
1802 
1803 
1804 
1805 template <>
1806 bool
1808 {
1809  return (this->vertex(0)[0] <= p[0]) && (p[0] <= this->vertex(1)[0]);
1810 }
1811 
1812 
1813 
1814 /*------------------------ Functions: CellAccessor<2> -----------------------*/
1815 
1816 
1817 
1818 template <>
1819 bool
1821 {
1822  // we check whether the point is
1823  // inside the cell by making sure
1824  // that it on the inner side of
1825  // each line defined by the faces,
1826  // i.e. for each of the four faces
1827  // we take the line that connects
1828  // the two vertices and subdivide
1829  // the whole domain by that in two
1830  // and check whether the point is
1831  // on the `cell-side' (rather than
1832  // the `out-side') of this line. if
1833  // the point is on the `cell-side'
1834  // for all four faces, it must be
1835  // inside the cell.
1836 
1837  // we want the faces in counter
1838  // clockwise orientation
1839  static const int direction[4] = {-1, 1, 1, -1};
1840  for (unsigned int f = 0; f < 4; ++f)
1841  {
1842  // vector from the first vertex
1843  // of the line to the point
1844  const Tensor<1, 2> to_p =
1845  p - this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 0));
1846  // vector describing the line
1847  const Tensor<1, 2> face =
1848  direction[f] *
1849  (this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 1)) -
1850  this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 0)));
1851 
1852  // if we rotate the face vector
1853  // by 90 degrees to the left
1854  // (i.e. it points to the
1855  // inside) and take the scalar
1856  // product with the vector from
1857  // the vertex to the point,
1858  // then the point is in the
1859  // `cell-side' if the scalar
1860  // product is positive. if this
1861  // is not the case, we can be
1862  // sure that the point is
1863  // outside
1864  if ((-face[1] * to_p[0] + face[0] * to_p[1]) < 0)
1865  return false;
1866  }
1867 
1868  // if we arrived here, then the
1869  // point is inside for all four
1870  // faces, and thus inside
1871  return true;
1872 }
1873 
1874 
1875 
1876 /*------------------------ Functions: CellAccessor<3> -----------------------*/
1877 
1878 
1879 
1880 template <>
1881 bool
1883 {
1884  // original implementation by Joerg
1885  // Weimar
1886 
1887  // we first eliminate points based
1888  // on the maximum and minimum of
1889  // the corner coordinates, then
1890  // transform to the unit cell, and
1891  // check there.
1892  const unsigned int dim = 3;
1893  const unsigned int spacedim = 3;
1894  Point<spacedim> maxp = this->vertex(0);
1895  Point<spacedim> minp = this->vertex(0);
1896 
1897  for (unsigned int v = 1; v < this->n_vertices(); ++v)
1898  for (unsigned int d = 0; d < dim; ++d)
1899  {
1900  maxp[d] = std::max(maxp[d], this->vertex(v)[d]);
1901  minp[d] = std::min(minp[d], this->vertex(v)[d]);
1902  }
1903 
1904  // rule out points outside the
1905  // bounding box of this cell
1906  for (unsigned int d = 0; d < dim; ++d)
1907  if ((p[d] < minp[d]) || (p[d] > maxp[d]))
1908  return false;
1909 
1910  // now we need to check more carefully: transform to the
1911  // unit cube and check there. unfortunately, this isn't
1912  // completely trivial since the transform_real_to_unit_cell
1913  // function may throw an exception that indicates that the
1914  // point given could not be inverted. we take this as a sign
1915  // that the point actually lies outside, as also documented
1916  // for that function
1917  try
1918  {
1919  const TriaRawIterator<CellAccessor<dim, spacedim>> cell_iterator(*this);
1921  reference_cell()
1922  .template get_default_linear_mapping<dim, spacedim>()
1923  .transform_real_to_unit_cell(cell_iterator, p)));
1924  }
1926  {
1927  return false;
1928  }
1929 }
1930 
1931 
1932 
1933 /*------------------- Functions: CellAccessor<dim,spacedim> -----------------*/
1934 
1935 // For codim>0 we proceed as follows:
1936 // 1) project point onto manifold and
1937 // 2) transform to the unit cell with a Q1 mapping
1938 // 3) then check if inside unit cell
1939 template <int dim, int spacedim>
1940 template <int dim_, int spacedim_>
1941 bool
1943 {
1944  const TriaRawIterator<CellAccessor<dim_, spacedim_>> cell_iterator(*this);
1945  const Point<dim_> p_unit =
1946  StaticMappingQ1<dim_, spacedim_>::mapping.transform_real_to_unit_cell(
1947  cell_iterator, p);
1948 
1950 }
1951 
1952 
1953 
1954 template <>
1955 bool
1957 {
1958  return point_inside_codim<1, 2>(p);
1959 }
1960 
1961 
1962 template <>
1963 bool
1965 {
1966  return point_inside_codim<1, 3>(p);
1967 }
1968 
1969 
1970 template <>
1971 bool
1973 {
1974  return point_inside_codim<2, 3>(p);
1975 }
1976 
1977 
1978 
1979 template <int dim, int spacedim>
1980 bool
1982 {
1983  for (const auto face : this->face_indices())
1984  if (at_boundary(face))
1985  return true;
1986 
1987  return false;
1988 }
1989 
1990 
1991 
1992 template <int dim, int spacedim>
1995 {
1997  return this->tria->levels[this->present_level]
1998  ->cells.boundary_or_material_id[this->present_index]
1999  .material_id;
2000 }
2001 
2002 
2003 
2004 template <int dim, int spacedim>
2005 void
2007  const types::material_id mat_id) const
2008 {
2011  this->tria->levels[this->present_level]
2012  ->cells.boundary_or_material_id[this->present_index]
2013  .material_id = mat_id;
2014 }
2015 
2016 
2017 
2018 template <int dim, int spacedim>
2019 void
2021  const types::material_id mat_id) const
2022 {
2023  set_material_id(mat_id);
2024 
2025  if (this->has_children())
2026  for (unsigned int c = 0; c < this->n_children(); ++c)
2027  this->child(c)->recursively_set_material_id(mat_id);
2028 }
2029 
2030 
2031 
2032 template <int dim, int spacedim>
2033 void
2035  const types::subdomain_id new_subdomain_id) const
2036 {
2038  Assert(this->is_active(),
2039  ExcMessage("set_subdomain_id() can only be called on active cells!"));
2040  this->tria->levels[this->present_level]->subdomain_ids[this->present_index] =
2041  new_subdomain_id;
2042 }
2043 
2044 
2045 
2046 template <int dim, int spacedim>
2049 {
2051  return this->tria->levels[this->present_level]
2052  ->level_subdomain_ids[this->present_index];
2053 }
2054 
2055 
2056 
2057 template <int dim, int spacedim>
2058 void
2060  const types::subdomain_id new_level_subdomain_id) const
2061 {
2063  this->tria->levels[this->present_level]
2064  ->level_subdomain_ids[this->present_index] = new_level_subdomain_id;
2065 }
2066 
2067 
2068 template <int dim, int spacedim>
2069 bool
2071 {
2073  if (dim == spacedim)
2074  return true;
2075  else
2076  return this->tria->levels[this->present_level]
2077  ->direction_flags[this->present_index];
2078 }
2079 
2080 
2081 
2082 template <int dim, int spacedim>
2083 void
2085  const bool new_direction_flag) const
2086 {
2088  if (dim < spacedim)
2089  this->tria->levels[this->present_level]
2090  ->direction_flags[this->present_index] = new_direction_flag;
2091  else
2092  Assert(new_direction_flag == true,
2093  ExcMessage("If dim==spacedim, direction flags are always true and "
2094  "can not be set to anything else."));
2095 }
2096 
2097 
2098 
2099 template <int dim, int spacedim>
2100 void
2101 CellAccessor<dim, spacedim>::set_parent(const unsigned int parent_index)
2102 {
2104  Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2105  this->tria->levels[this->present_level]->parents[this->present_index / 2] =
2106  parent_index;
2107 }
2108 
2109 
2110 
2111 template <int dim, int spacedim>
2112 int
2114 {
2115  Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2116 
2117  // the parent of two consecutive cells
2118  // is stored only once, since it is
2119  // the same
2120  return this->tria->levels[this->present_level]
2121  ->parents[this->present_index / 2];
2122 }
2123 
2124 
2125 
2126 template <int dim, int spacedim>
2127 unsigned int
2129 {
2130  Assert(this->is_active(), TriaAccessorExceptions::ExcCellNotActive());
2131  return this->tria->levels[this->present_level]
2132  ->active_cell_indices[this->present_index];
2133 }
2134 
2135 
2136 
2137 template <int dim, int spacedim>
2138 void
2140  const unsigned int active_cell_index) const
2141 {
2142  this->tria->levels[this->present_level]
2143  ->active_cell_indices[this->present_index] = active_cell_index;
2144 }
2145 
2146 
2147 
2148 template <int dim, int spacedim>
2149 void
2151  const types::global_cell_index index) const
2152 {
2153  this->tria->levels[this->present_level]
2154  ->global_active_cell_indices[this->present_index] = index;
2155 }
2156 
2157 
2158 
2159 template <int dim, int spacedim>
2162 {
2164  Assert(this->is_active(),
2165  ExcMessage(
2166  "global_active_cell_index() can only be called on active cells!"));
2167 
2168  return this->tria->levels[this->present_level]
2169  ->global_active_cell_indices[this->present_index];
2170 }
2171 
2172 
2173 
2174 template <int dim, int spacedim>
2175 void
2177  const types::global_cell_index index) const
2178 {
2179  this->tria->levels[this->present_level]
2180  ->global_level_cell_indices[this->present_index] = index;
2181 }
2182 
2183 
2184 
2185 template <int dim, int spacedim>
2188 {
2189  return this->tria->levels[this->present_level]
2190  ->global_level_cell_indices[this->present_index];
2191 }
2192 
2193 
2194 
2195 template <int dim, int spacedim>
2198 {
2200  Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2202  this->present_level - 1,
2203  parent_index());
2204 
2205  return q;
2206 }
2207 
2208 
2209 template <int dim, int spacedim>
2210 void
2212  const types::subdomain_id new_subdomain_id) const
2213 {
2214  if (this->has_children())
2215  for (unsigned int c = 0; c < this->n_children(); ++c)
2216  this->child(c)->recursively_set_subdomain_id(new_subdomain_id);
2217  else
2218  set_subdomain_id(new_subdomain_id);
2219 }
2220 
2221 
2222 
2223 template <int dim, int spacedim>
2224 void
2226  const unsigned int i,
2227  const TriaIterator<CellAccessor<dim, spacedim>> &pointer) const
2228 {
2229  AssertIndexRange(i, this->n_faces());
2230 
2231  if (pointer.state() == IteratorState::valid)
2232  {
2233  this->tria->levels[this->present_level]
2234  ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
2235  .first = pointer->present_level;
2236  this->tria->levels[this->present_level]
2237  ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
2238  .second = pointer->present_index;
2239  }
2240  else
2241  {
2242  this->tria->levels[this->present_level]
2243  ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
2244  .first = -1;
2245  this->tria->levels[this->present_level]
2246  ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
2247  .second = -1;
2248  }
2249 }
2250 
2251 
2252 
2253 template <int dim, int spacedim>
2254 CellId
2256 {
2257  std::array<unsigned char, 30> id;
2258 
2259  CellAccessor<dim, spacedim> ptr = *this;
2260  const unsigned int n_child_indices = ptr.level();
2261 
2262  while (ptr.level() > 0)
2263  {
2264  const TriaIterator<CellAccessor<dim, spacedim>> parent = ptr.parent();
2265  const unsigned int n_children = parent->n_children();
2266 
2267  // determine which child we are
2268  unsigned char v = static_cast<unsigned char>(-1);
2269  for (unsigned int c = 0; c < n_children; ++c)
2270  {
2271  if (parent->child_index(c) == ptr.index())
2272  {
2273  v = c;
2274  break;
2275  }
2276  }
2277 
2278  Assert(v != static_cast<unsigned char>(-1), ExcInternalError());
2279  id[ptr.level() - 1] = v;
2280 
2281  ptr.copy_from(*parent);
2282  }
2283 
2284  Assert(ptr.level() == 0, ExcInternalError());
2285  const unsigned int coarse_index = ptr.index();
2286 
2287  return {this->tria->coarse_cell_index_to_coarse_cell_id(coarse_index),
2288  n_child_indices,
2289  id.data()};
2290 }
2291 
2292 
2293 
2294 template <int dim, int spacedim>
2295 unsigned int
2297  const unsigned int neighbor) const
2298 {
2299  AssertIndexRange(neighbor, this->n_faces());
2300 
2301  // if we have a 1d mesh in 1d, we
2302  // can assume that the left
2303  // neighbor of the right neighbor is
2304  // the current cell. but that is an
2305  // invariant that isn't true if the
2306  // mesh is embedded in a higher
2307  // dimensional space, so we have to
2308  // fall back onto the generic code
2309  // below
2310  if ((dim == 1) && (spacedim == dim))
2311  return GeometryInfo<dim>::opposite_face[neighbor];
2312 
2313  const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2314  this->neighbor(neighbor);
2315 
2316  // usually, on regular patches of
2317  // the grid, this cell is just on
2318  // the opposite side of the
2319  // neighbor that the neighbor is of
2320  // this cell. for example in 2d, if
2321  // we want to know the
2322  // neighbor_of_neighbor if
2323  // neighbor==1 (the right
2324  // neighbor), then we will get 3
2325  // (the left neighbor) in most
2326  // cases. look up this relationship
2327  // in the table provided by
2328  // GeometryInfo and try it
2329  const unsigned int this_face_index = face_index(neighbor);
2330 
2331  const unsigned int neighbor_guess =
2333 
2334  if (neighbor_guess < neighbor_cell->n_faces() &&
2335  neighbor_cell->face_index(neighbor_guess) == this_face_index)
2336  return neighbor_guess;
2337  else
2338  // if the guess was false, then
2339  // we need to loop over all
2340  // neighbors and find the number
2341  // the hard way
2342  {
2343  for (const unsigned int face_no : neighbor_cell->face_indices())
2344  if (neighbor_cell->face_index(face_no) == this_face_index)
2345  return face_no;
2346 
2347  // running over all neighbors
2348  // faces we did not find the
2349  // present face. Thereby the
2350  // neighbor must be coarser
2351  // than the present
2352  // cell. Return an invalid
2353  // unsigned int in this case.
2355  }
2356 }
2357 
2358 
2359 
2360 template <int dim, int spacedim>
2361 unsigned int
2363  const unsigned int face_no) const
2364 {
2365  const unsigned int n2 = neighbor_of_neighbor_internal(face_no);
2368 
2369  return n2;
2370 }
2371 
2372 
2373 
2374 template <int dim, int spacedim>
2375 bool
2377  const unsigned int face_no) const
2378 {
2379  return neighbor_of_neighbor_internal(face_no) ==
2381 }
2382 
2383 
2384 
2385 template <int dim, int spacedim>
2386 std::pair<unsigned int, unsigned int>
2388  const unsigned int neighbor) const
2389 {
2390  AssertIndexRange(neighbor, this->n_faces());
2391  // make sure that the neighbor is
2392  // on a coarser level
2393  Assert(neighbor_is_coarser(neighbor),
2395 
2396  switch (dim)
2397  {
2398  case 2:
2399  {
2400  const int this_face_index = face_index(neighbor);
2401  const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2402  this->neighbor(neighbor);
2403 
2404  // usually, on regular patches of
2405  // the grid, this cell is just on
2406  // the opposite side of the
2407  // neighbor that the neighbor is of
2408  // this cell. for example in 2d, if
2409  // we want to know the
2410  // neighbor_of_neighbor if
2411  // neighbor==1 (the right
2412  // neighbor), then we will get 0
2413  // (the left neighbor) in most
2414  // cases. look up this relationship
2415  // in the table provided by
2416  // GeometryInfo and try it
2417  const unsigned int face_no_guess =
2419 
2420  const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face_guess =
2421  neighbor_cell->face(face_no_guess);
2422 
2423  if (face_guess->has_children())
2424  for (unsigned int subface_no = 0;
2425  subface_no < face_guess->n_children();
2426  ++subface_no)
2427  if (face_guess->child_index(subface_no) == this_face_index)
2428  return std::make_pair(face_no_guess, subface_no);
2429 
2430  // if the guess was false, then
2431  // we need to loop over all faces
2432  // and subfaces and find the
2433  // number the hard way
2434  for (const unsigned int face_no : neighbor_cell->face_indices())
2435  {
2436  if (face_no != face_no_guess)
2437  {
2438  const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>
2439  face = neighbor_cell->face(face_no);
2440  if (face->has_children())
2441  for (unsigned int subface_no = 0;
2442  subface_no < face->n_children();
2443  ++subface_no)
2444  if (face->child_index(subface_no) == this_face_index)
2445  return std::make_pair(face_no, subface_no);
2446  }
2447  }
2448 
2449  // we should never get here,
2450  // since then we did not find
2451  // our way back...
2452  Assert(false, ExcInternalError());
2453  return std::make_pair(numbers::invalid_unsigned_int,
2455  }
2456 
2457  case 3:
2458  {
2459  const int this_face_index = face_index(neighbor);
2460  const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2461  this->neighbor(neighbor);
2462 
2463  // usually, on regular patches of the grid, this cell is just on the
2464  // opposite side of the neighbor that the neighbor is of this cell.
2465  // for example in 2d, if we want to know the neighbor_of_neighbor if
2466  // neighbor==1 (the right neighbor), then we will get 0 (the left
2467  // neighbor) in most cases. look up this relationship in the table
2468  // provided by GeometryInfo and try it
2469  const unsigned int face_no_guess =
2471 
2472  const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face_guess =
2473  neighbor_cell->face(face_no_guess);
2474 
2475  if (face_guess->has_children())
2476  for (unsigned int subface_no = 0;
2477  subface_no < face_guess->n_children();
2478  ++subface_no)
2479  {
2480  if (face_guess->child_index(subface_no) == this_face_index)
2481  // call a helper function, that translates the current
2482  // subface number to a subface number for the current
2483  // FaceRefineCase
2484  return std::make_pair(face_no_guess,
2485  translate_subface_no(face_guess,
2486  subface_no));
2487 
2488  if (face_guess->child(subface_no)->has_children())
2489  for (unsigned int subsub_no = 0;
2490  subsub_no < face_guess->child(subface_no)->n_children();
2491  ++subsub_no)
2492  if (face_guess->child(subface_no)->child_index(subsub_no) ==
2493  this_face_index)
2494  // call a helper function, that translates the current
2495  // subface number and subsubface number to a subface
2496  // number for the current FaceRefineCase
2497  return std::make_pair(face_no_guess,
2498  translate_subface_no(face_guess,
2499  subface_no,
2500  subsub_no));
2501  }
2502 
2503  // if the guess was false, then we need to loop over all faces and
2504  // subfaces and find the number the hard way
2505  for (const unsigned int face_no : neighbor_cell->face_indices())
2506  {
2507  if (face_no == face_no_guess)
2508  continue;
2509 
2510  const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face =
2511  neighbor_cell->face(face_no);
2512 
2513  if (!face->has_children())
2514  continue;
2515 
2516  for (unsigned int subface_no = 0; subface_no < face->n_children();
2517  ++subface_no)
2518  {
2519  if (face->child_index(subface_no) == this_face_index)
2520  // call a helper function, that translates the current
2521  // subface number to a subface number for the current
2522  // FaceRefineCase
2523  return std::make_pair(face_no,
2524  translate_subface_no(face,
2525  subface_no));
2526 
2527  if (face->child(subface_no)->has_children())
2528  for (unsigned int subsub_no = 0;
2529  subsub_no < face->child(subface_no)->n_children();
2530  ++subsub_no)
2531  if (face->child(subface_no)->child_index(subsub_no) ==
2532  this_face_index)
2533  // call a helper function, that translates the current
2534  // subface number and subsubface number to a subface
2535  // number for the current FaceRefineCase
2536  return std::make_pair(face_no,
2537  translate_subface_no(face,
2538  subface_no,
2539  subsub_no));
2540  }
2541  }
2542 
2543  // we should never get here, since then we did not find our way
2544  // back...
2545  Assert(false, ExcInternalError());
2546  return std::make_pair(numbers::invalid_unsigned_int,
2548  }
2549 
2550  default:
2551  {
2552  Assert(false, ExcImpossibleInDim(1));
2553  return std::make_pair(numbers::invalid_unsigned_int,
2555  }
2556  }
2557 }
2558 
2559 
2560 
2561 template <int dim, int spacedim>
2562 bool
2564  const unsigned int i_face) const
2565 {
2566  /*
2567  * Implementation note: In all of the functions corresponding to periodic
2568  * faces we mainly use the Triangulation::periodic_face_map to find the
2569  * information about periodically connected faces. So, we actually search in
2570  * this std::map and return the cell_face on the other side of the periodic
2571  * boundary.
2572  *
2573  * We can not use operator[] as this would insert non-existing entries or
2574  * would require guarding with an extra std::map::find() or count().
2575  */
2576  AssertIndexRange(i_face, this->n_faces());
2577  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2578 
2579  cell_iterator current_cell(*this);
2580  if (this->tria->periodic_face_map.find(
2581  std::make_pair(current_cell, i_face)) !=
2582  this->tria->periodic_face_map.end())
2583  return true;
2584  return false;
2585 }
2586 
2587 
2588 
2589 template <int dim, int spacedim>
2591 CellAccessor<dim, spacedim>::periodic_neighbor(const unsigned int i_face) const
2592 {
2593  /*
2594  * To know, why we are using std::map::find() instead of [] operator, refer
2595  * to the implementation note in has_periodic_neighbor() function.
2596  *
2597  * my_it : the iterator to the current cell.
2598  * my_face_pair : the pair reported by periodic_face_map as its first pair
2599  * being the current cell_face.
2600  */
2601  AssertIndexRange(i_face, this->n_faces());
2602  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2603  cell_iterator current_cell(*this);
2604 
2605  auto my_face_pair =
2606  this->tria->periodic_face_map.find(std::make_pair(current_cell, i_face));
2607 
2608  // Make sure we are actually on a periodic boundary:
2609  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2611  return my_face_pair->second.first.first;
2612 }
2613 
2614 
2615 
2616 template <int dim, int spacedim>
2619  const unsigned int i_face) const
2620 {
2621  if (!(this->face(i_face)->at_boundary()))
2622  return this->neighbor(i_face);
2623  else if (this->has_periodic_neighbor(i_face))
2624  return this->periodic_neighbor(i_face);
2625  else
2627  // we can't come here
2628  return this->neighbor(i_face);
2629 }
2630 
2631 
2632 
2633 template <int dim, int spacedim>
2636  const unsigned int i_face,
2637  const unsigned int i_subface) const
2638 {
2639  /*
2640  * To know, why we are using std::map::find() instead of [] operator, refer
2641  * to the implementation note in has_periodic_neighbor() function.
2642  *
2643  * my_it : the iterator to the current cell.
2644  * my_face_pair : the pair reported by periodic_face_map as its first pair
2645  * being the current cell_face. nb_it : the iterator to the
2646  * neighbor of current cell at i_face. face_num_of_nb : the face number of
2647  * the periodically neighboring face in the relevant element.
2648  * nb_parent_face_it: the iterator to the parent face of the periodically
2649  * neighboring face.
2650  */
2651  AssertIndexRange(i_face, this->n_faces());
2652  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2653  cell_iterator my_it(*this);
2654 
2655  auto my_face_pair =
2656  this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2657  /*
2658  * There should be an assertion, which tells the user that this function
2659  * should not be used for a cell which is not located at a periodic boundary.
2660  */
2661  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2663  cell_iterator parent_nb_it = my_face_pair->second.first.first;
2664  unsigned int nb_face_num = my_face_pair->second.first.second;
2665  TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> nb_parent_face_it =
2666  parent_nb_it->face(nb_face_num);
2667  /*
2668  * We should check if the parent face of the neighbor has at least the same
2669  * number of children as i_subface.
2670  */
2671  AssertIndexRange(i_subface, nb_parent_face_it->n_children());
2672  unsigned int sub_neighbor_num =
2673  GeometryInfo<dim>::child_cell_on_face(parent_nb_it->refinement_case(),
2674  nb_face_num,
2675  i_subface,
2676  my_face_pair->second.second[0],
2677  my_face_pair->second.second[1],
2678  my_face_pair->second.second[2],
2679  nb_parent_face_it->refinement_case());
2680  return parent_nb_it->child(sub_neighbor_num);
2681 }
2682 
2683 
2684 
2685 template <int dim, int spacedim>
2686 std::pair<unsigned int, unsigned int>
2688  const unsigned int i_face) const
2689 {
2690  /*
2691  * To know, why we are using std::map::find() instead of [] operator, refer
2692  * to the implementation note in has_periodic_neighbor() function.
2693  *
2694  * my_it : the iterator to the current cell.
2695  * my_face_pair : the pair reported by periodic_face_map as its first pair
2696  * being the current cell_face. nb_it : the iterator to the periodic
2697  * neighbor. nb_face_pair : the pair reported by periodic_face_map as its
2698  * first pair being the periodic neighbor cell_face. p_nb_of_p_nb : the
2699  * iterator of the periodic neighbor of the periodic neighbor of the current
2700  * cell.
2701  */
2702  AssertIndexRange(i_face, this->n_faces());
2703  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2704  const int my_face_index = this->face_index(i_face);
2705  cell_iterator my_it(*this);
2706 
2707  auto my_face_pair =
2708  this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2709  /*
2710  * There should be an assertion, which tells the user that this function
2711  * should not be used for a cell which is not located at a periodic boundary.
2712  */
2713  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2715  cell_iterator nb_it = my_face_pair->second.first.first;
2716  unsigned int face_num_of_nb = my_face_pair->second.first.second;
2717 
2718  auto nb_face_pair =
2719  this->tria->periodic_face_map.find(std::make_pair(nb_it, face_num_of_nb));
2720  /*
2721  * Since, we store periodic neighbors for every cell (either active or
2722  * artificial or inactive) the nb_face_pair should also be mapped to some
2723  * cell_face pair. We assert this here.
2724  */
2725  Assert(nb_face_pair != this->tria->periodic_face_map.end(),
2727  cell_iterator p_nb_of_p_nb = nb_face_pair->second.first.first;
2728  TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> parent_face_it =
2729  p_nb_of_p_nb->face(nb_face_pair->second.first.second);
2730  for (unsigned int i_subface = 0; i_subface < parent_face_it->n_children();
2731  ++i_subface)
2732  if (parent_face_it->child_index(i_subface) == my_face_index)
2733  return std::make_pair(face_num_of_nb, i_subface);
2734  /*
2735  * Obviously, if the execution reaches to this point, some of our assumptions
2736  * should have been false. The most important one is, the user has called this
2737  * function on a face which does not have a coarser periodic neighbor.
2738  */
2740  return std::make_pair(numbers::invalid_unsigned_int,
2742 }
2743 
2744 
2745 
2746 template <int dim, int spacedim>
2747 int
2749  const unsigned int i_face) const
2750 {
2751  return periodic_neighbor(i_face)->index();
2752 }
2753 
2754 
2755 
2756 template <int dim, int spacedim>
2757 int
2759  const unsigned int i_face) const
2760 {
2761  return periodic_neighbor(i_face)->level();
2762 }
2763 
2764 
2765 
2766 template <int dim, int spacedim>
2767 unsigned int
2769  const unsigned int i_face) const
2770 {
2771  return periodic_neighbor_face_no(i_face);
2772 }
2773 
2774 
2775 
2776 template <int dim, int spacedim>
2777 unsigned int
2779  const unsigned int i_face) const
2780 {
2781  /*
2782  * To know, why we are using std::map::find() instead of [] operator, refer
2783  * to the implementation note in has_periodic_neighbor() function.
2784  *
2785  * my_it : the iterator to the current cell.
2786  * my_face_pair : the pair reported by periodic_face_map as its first pair
2787  * being the current cell_face.
2788  */
2789  AssertIndexRange(i_face, this->n_faces());
2790  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2791  cell_iterator my_it(*this);
2792 
2793  auto my_face_pair =
2794  this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2795  /*
2796  * There should be an assertion, which tells the user that this function
2797  * should not be called for a cell which is not located at a periodic boundary
2798  * !
2799  */
2800  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2802  return my_face_pair->second.first.second;
2803 }
2804 
2805 
2806 
2807 template <int dim, int spacedim>
2808 bool
2810  const unsigned int i_face) const
2811 {
2812  /*
2813  * To know, why we are using std::map::find() instead of [] operator, refer
2814  * to the implementation note in has_periodic_neighbor() function.
2815  *
2816  * Implementation note: Let p_nb_of_p_nb be the periodic neighbor of the
2817  * periodic neighbor of the current cell. Also, let p_face_of_p_nb_of_p_nb be
2818  * the periodic face of the p_nb_of_p_nb. If p_face_of_p_nb_of_p_nb has
2819  * children , then the periodic neighbor of the current cell is coarser than
2820  * itself. Although not tested, this implementation should work for
2821  * anisotropic refinement as well.
2822  *
2823  * my_it : the iterator to the current cell.
2824  * my_face_pair : the pair reported by periodic_face_map as its first pair
2825  * being the current cell_face. nb_it : the iterator to the periodic
2826  * neighbor. nb_face_pair : the pair reported by periodic_face_map as its
2827  * first pair being the periodic neighbor cell_face.
2828  */
2829  AssertIndexRange(i_face, this->n_faces());
2830  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2831  cell_iterator my_it(*this);
2832 
2833  auto my_face_pair =
2834  this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2835  /*
2836  * There should be an assertion, which tells the user that this function
2837  * should not be used for a cell which is not located at a periodic boundary.
2838  */
2839  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2841 
2842  cell_iterator nb_it = my_face_pair->second.first.first;
2843  unsigned int face_num_of_nb = my_face_pair->second.first.second;
2844 
2845  auto nb_face_pair =
2846  this->tria->periodic_face_map.find(std::make_pair(nb_it, face_num_of_nb));
2847  /*
2848  * Since, we store periodic neighbors for every cell (either active or
2849  * artificial or inactive) the nb_face_pair should also be mapped to some
2850  * cell_face pair. We assert this here.
2851  */
2852  Assert(nb_face_pair != this->tria->periodic_face_map.end(),
2854  const unsigned int my_level = this->level();
2855  const unsigned int neighbor_level = nb_face_pair->second.first.first->level();
2856  Assert(my_level >= neighbor_level, ExcInternalError());
2857  return my_level > neighbor_level;
2858 }
2859 
2860 
2861 
2862 template <int dim, int spacedim>
2863 bool
2864 CellAccessor<dim, spacedim>::at_boundary(const unsigned int i) const
2865 {
2867  AssertIndexRange(i, this->n_faces());
2868 
2869  return (neighbor_index(i) == -1);
2870 }
2871 
2872 
2873 
2874 template <int dim, int spacedim>
2875 bool
2877 {
2878  if (dim == 1)
2879  return at_boundary();
2880  else
2881  {
2882  for (unsigned int l = 0; l < this->n_lines(); ++l)
2883  if (this->line(l)->at_boundary())
2884  return true;
2885 
2886  return false;
2887  }
2888 }
2889 
2890 
2891 
2892 template <int dim, int spacedim>
2895  const unsigned int face,
2896  const unsigned int subface) const
2897 {
2898  Assert(!this->has_children(),
2899  ExcMessage("The present cell must not have children!"));
2900  Assert(!this->at_boundary(face),
2901  ExcMessage("The present cell must have a valid neighbor!"));
2902  Assert(this->neighbor(face)->has_children() == true,
2903  ExcMessage("The neighbor must have children!"));
2904 
2905  switch (dim)
2906  {
2907  case 2:
2908  {
2910  {
2911  const auto neighbor_cell = this->neighbor(face);
2912 
2913  // only for isotropic refinement at the moment
2914  Assert(neighbor_cell->refinement_case() ==
2916  ExcNotImplemented());
2917 
2918  // determine indices for this cell's subface from the perspective
2919  // of the neighboring cell
2920  const unsigned int neighbor_face =
2921  this->neighbor_of_neighbor(face);
2922  // two neighboring cells have an opposed orientation on their
2923  // shared face if both of them follow the same orientation type
2924  // (i.e., standard or non-standard).
2925  // we verify this with a XOR operation.
2926  const unsigned int neighbor_subface =
2927  (!(this->line_orientation(face)) !=
2928  !(neighbor_cell->line_orientation(neighbor_face))) ?
2929  (1 - subface) :
2930  subface;
2931 
2932  const unsigned int neighbor_child_index =
2934  neighbor_subface);
2935  const TriaIterator<CellAccessor<dim, spacedim>> sub_neighbor =
2936  neighbor_cell->child(neighbor_child_index);
2937 
2938  // neighbor's child is not allowed to be further refined for the
2939  // moment
2940  Assert(sub_neighbor->refinement_case() ==
2942  ExcNotImplemented());
2943 
2944  return sub_neighbor;
2945  }
2946  else if (this->reference_cell() == ReferenceCells::Quadrilateral)
2947  {
2948  const unsigned int neighbor_neighbor =
2949  this->neighbor_of_neighbor(face);
2950  const unsigned int neighbor_child_index =
2952  this->neighbor(face)->refinement_case(),
2953  neighbor_neighbor,
2954  subface);
2955 
2957  this->neighbor(face)->child(neighbor_child_index);
2958  // the neighbors child can have children,
2959  // which are not further refined along the
2960  // face under consideration. as we are
2961  // normally interested in one of this
2962  // child's child, search for the right one.
2963  while (sub_neighbor->has_children())
2964  {
2966  sub_neighbor->refinement_case(),
2967  neighbor_neighbor) ==
2969  ExcInternalError());
2970  sub_neighbor =
2971  sub_neighbor->child(GeometryInfo<dim>::child_cell_on_face(
2972  sub_neighbor->refinement_case(), neighbor_neighbor, 0));
2973  }
2974 
2975  return sub_neighbor;
2976  }
2977 
2978  // if no reference cell type matches
2979  Assert(false, ExcNotImplemented());
2981  }
2982 
2983 
2984  case 3:
2985  {
2987  {
2988  // this function returns the neighbor's
2989  // child on a given face and
2990  // subface.
2991 
2992  // we have to consider one other aspect here:
2993  // The face might be refined
2994  // anisotropically. In this case, the subface
2995  // number refers to the following, where we
2996  // look at the face from the current cell,
2997  // thus the subfaces are in standard
2998  // orientation concerning the cell
2999  //
3000  // for isotropic refinement
3001  //
3002  // *---*---*
3003  // | 2 | 3 |
3004  // *---*---*
3005  // | 0 | 1 |
3006  // *---*---*
3007  //
3008  // for 2*anisotropic refinement
3009  // (first cut_y, then cut_x)
3010  //
3011  // *---*---*
3012  // | 2 | 3 |
3013  // *---*---*
3014  // | 0 | 1 |
3015  // *---*---*
3016  //
3017  // for 2*anisotropic refinement
3018  // (first cut_x, then cut_y)
3019  //
3020  // *---*---*
3021  // | 1 | 3 |
3022  // *---*---*
3023  // | 0 | 2 |
3024  // *---*---*
3025  //
3026  // for purely anisotropic refinement:
3027  //
3028  // *---*---* *-------*
3029  // | | | | 1 |
3030  // | 0 | 1 | or *-------*
3031  // | | | | 0 |
3032  // *---*---* *-------*
3033  //
3034  // for "mixed" refinement:
3035  //
3036  // *---*---* *---*---* *---*---* *-------*
3037  // | | 2 | | 1 | | | 1 | 2 | | 2 |
3038  // | 0 *---* or *---* 2 | or *---*---* or *---*---*
3039  // | | 1 | | 0 | | | 0 | | 0 | 1 |
3040  // *---*---* *---*---* *-------* *---*---*
3041 
3043  mother_face = this->face(face);
3044  const unsigned int total_children =
3045  mother_face->n_active_descendants();
3046  AssertIndexRange(subface, total_children);
3048  ExcInternalError());
3049 
3050  unsigned int neighbor_neighbor;
3053  this->neighbor(face);
3054 
3055 
3056  const RefinementCase<dim - 1> mother_face_ref_case =
3057  mother_face->refinement_case();
3058  if (mother_face_ref_case ==
3059  static_cast<RefinementCase<dim - 1>>(
3060  RefinementCase<2>::cut_xy)) // total_children==4
3061  {
3062  // this case is quite easy. we are sure,
3063  // that the neighbor is not coarser.
3064 
3065  // get the neighbor's number for the given
3066  // face and the neighbor
3067  neighbor_neighbor = this->neighbor_of_neighbor(face);
3068 
3069  // now use the info provided by GeometryInfo
3070  // to extract the neighbors child number
3071  const unsigned int neighbor_child_index =
3073  neighbor->refinement_case(),
3074  neighbor_neighbor,
3075  subface,
3076  neighbor->face_orientation(neighbor_neighbor),
3077  neighbor->face_flip(neighbor_neighbor),
3078  neighbor->face_rotation(neighbor_neighbor));
3079  neighbor_child = neighbor->child(neighbor_child_index);
3080 
3081  // make sure that the neighbor child cell we
3082  // have found shares the desired subface.
3083  Assert((this->face(face)->child(subface) ==
3084  neighbor_child->face(neighbor_neighbor)),
3085  ExcInternalError());
3086  }
3087  else //-> the face is refined anisotropically
3088  {
3089  // first of all, we have to find the
3090  // neighbor at one of the anisotropic
3091  // children of the
3092  // mother_face. determine, which of
3093  // these we need.
3094  unsigned int first_child_to_find;
3095  unsigned int neighbor_child_index;
3096  if (total_children == 2)
3097  first_child_to_find = subface;
3098  else
3099  {
3100  first_child_to_find = subface / 2;
3101  if (total_children == 3 && subface == 1 &&
3102  !mother_face->child(0)->has_children())
3103  first_child_to_find = 1;
3104  }
3105  if (neighbor_is_coarser(face))
3106  {
3107  std::pair<unsigned int, unsigned int> indices =
3108  neighbor_of_coarser_neighbor(face);
3109  neighbor_neighbor = indices.first;
3110 
3111 
3112  // we have to translate our
3113  // subface_index according to the
3114  // RefineCase and subface index of
3115  // the coarser face (our face is an
3116  // anisotropic child of the coarser
3117  // face), 'a' denotes our
3118  // subface_index 0 and 'b' denotes
3119  // our subface_index 1, whereas 0...3
3120  // denote isotropic subfaces of the
3121  // coarser face
3122  //
3123  // cut_x and coarser_subface_index=0
3124  //
3125  // *---*---*
3126  // |b=2| |
3127  // | | |
3128  // |a=0| |
3129  // *---*---*
3130  //
3131  // cut_x and coarser_subface_index=1
3132  //
3133  // *---*---*
3134  // | |b=3|
3135  // | | |
3136  // | |a=1|
3137  // *---*---*
3138  //
3139  // cut_y and coarser_subface_index=0
3140  //
3141  // *-------*
3142  // | |
3143  // *-------*
3144  // |a=0 b=1|
3145  // *-------*
3146  //
3147  // cut_y and coarser_subface_index=1
3148  //
3149  // *-------*
3150  // |a=2 b=3|
3151  // *-------*
3152  // | |
3153  // *-------*
3154  unsigned int iso_subface;
3155  if (neighbor->face(neighbor_neighbor)
3156  ->refinement_case() == RefinementCase<2>::cut_x)
3157  iso_subface = 2 * first_child_to_find + indices.second;
3158  else
3159  {
3160  Assert(neighbor->face(neighbor_neighbor)
3161  ->refinement_case() ==
3163  ExcInternalError());
3164  iso_subface =
3165  first_child_to_find + 2 * indices.second;
3166  }
3167  neighbor_child_index =
3169  neighbor->refinement_case(),
3170  neighbor_neighbor,
3171  iso_subface,
3172  neighbor->face_orientation(neighbor_neighbor),
3173  neighbor->face_flip(neighbor_neighbor),
3174  neighbor->face_rotation(neighbor_neighbor));
3175  }
3176  else // neighbor is not coarser
3177  {
3178  neighbor_neighbor = neighbor_of_neighbor(face);
3179  neighbor_child_index =
3181  neighbor->refinement_case(),
3182  neighbor_neighbor,
3183  first_child_to_find,
3184  neighbor->face_orientation(neighbor_neighbor),
3185  neighbor->face_flip(neighbor_neighbor),
3186  neighbor->face_rotation(neighbor_neighbor),
3187  mother_face_ref_case);
3188  }
3189 
3190  neighbor_child = neighbor->child(neighbor_child_index);
3191  // it might be, that the neighbor_child
3192  // has children, which are not refined
3193  // along the given subface. go down that
3194  // list and deliver the last of those.
3195  while (
3196  neighbor_child->has_children() &&
3198  neighbor_child->refinement_case(), neighbor_neighbor) ==
3200  neighbor_child = neighbor_child->child(
3202  neighbor_child->refinement_case(),
3203  neighbor_neighbor,
3204  0));
3205 
3206  // if there are two total subfaces, we
3207  // are finished. if there are four we
3208  // have to get a child of our current
3209  // neighbor_child. If there are three,
3210  // we have to check which of the two
3211  // possibilities applies.
3212  if (total_children == 3)
3213  {
3214  if (mother_face->child(0)->has_children())
3215  {
3216  if (subface < 2)
3217  neighbor_child = neighbor_child->child(
3219  neighbor_child->refinement_case(),
3220  neighbor_neighbor,
3221  subface,
3222  neighbor_child->face_orientation(
3223  neighbor_neighbor),
3224  neighbor_child->face_flip(neighbor_neighbor),
3225  neighbor_child->face_rotation(
3226  neighbor_neighbor),
3227  mother_face->child(0)->refinement_case()));
3228  }
3229  else
3230  {
3231  Assert(mother_face->child(1)->has_children(),
3232  ExcInternalError());
3233  if (subface > 0)
3234  neighbor_child = neighbor_child->child(
3236  neighbor_child->refinement_case(),
3237  neighbor_neighbor,
3238  subface - 1,
3239  neighbor_child->face_orientation(
3240  neighbor_neighbor),
3241  neighbor_child->face_flip(neighbor_neighbor),
3242  neighbor_child->face_rotation(
3243  neighbor_neighbor),
3244  mother_face->child(1)->refinement_case()));
3245  }
3246  }
3247  else if (total_children == 4)
3248  {
3249  neighbor_child = neighbor_child->child(
3251  neighbor_child->refinement_case(),
3252  neighbor_neighbor,
3253  subface % 2,
3254  neighbor_child->face_orientation(neighbor_neighbor),
3255  neighbor_child->face_flip(neighbor_neighbor),
3256  neighbor_child->face_rotation(neighbor_neighbor),
3257  mother_face->child(subface / 2)->refinement_case()));
3258  }
3259  }
3260 
3261  // it might be, that the neighbor_child has
3262  // children, which are not refined along the
3263  // given subface. go down that list and
3264  // deliver the last of those.
3265  while (neighbor_child->has_children())
3266  neighbor_child =
3267  neighbor_child->child(GeometryInfo<dim>::child_cell_on_face(
3268  neighbor_child->refinement_case(), neighbor_neighbor, 0));
3269 
3270 #ifdef DEBUG
3271  // check, whether the face neighbor_child matches the requested
3272  // subface.
3274  switch (this->subface_case(face))
3275  {
3279  requested = mother_face->child(subface);
3280  break;
3283  requested =
3284  mother_face->child(subface / 2)->child(subface % 2);
3285  break;
3286 
3289  switch (subface)
3290  {
3291  case 0:
3292  case 1:
3293  requested = mother_face->child(0)->child(subface);
3294  break;
3295  case 2:
3296  requested = mother_face->child(1);
3297  break;
3298  default:
3299  Assert(false, ExcInternalError());
3300  }
3301  break;
3304  switch (subface)
3305  {
3306  case 0:
3307  requested = mother_face->child(0);
3308  break;
3309  case 1:
3310  case 2:
3311  requested = mother_face->child(1)->child(subface - 1);
3312  break;
3313  default:
3314  Assert(false, ExcInternalError());
3315  }
3316  break;
3317  default:
3318  Assert(false, ExcInternalError());
3319  break;
3320  }
3321  Assert(requested == neighbor_child->face(neighbor_neighbor),
3322  ExcInternalError());
3323 #endif
3324 
3325  return neighbor_child;
3326  }
3327 
3328  // if no reference cell type matches
3329  Assert(false, ExcNotImplemented());
3331  }
3332 
3333  default:
3334  // if 1d or more than 3d
3335  Assert(false, ExcNotImplemented());
3337  }
3338 }
3339 
3340 
3341 
3342 // explicit instantiations
3343 #include "tria_accessor.inst"
3344 
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:82
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
const ::Triangulation< dim, spacedim > & tria
#define AssertIndexRange(index, range)
Definition: exceptions.h:1720
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:2650
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:1571
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:575
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:1461
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:401
unsigned int level
Definition: grid_out.cc:4589
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:185
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:506
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:400
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:698
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