Reference documentation for deal.II version GIT 725517848b 2022-05-23 21:10:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tria_accessor.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2022 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 + static_cast<unsigned int>(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  {
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  }
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 
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 }
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}, // Lines along x-axis, see GeometryInfo
1675  {0, 1}}; // Lines along y-axis
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}, // Lines along x-axis, see GeometryInfo
1689  {0, 1}}; // Lines along y-axis
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}, // Lines along x-axis, see GeometryInfo
1704  {0, 1, 4, 5}, // Lines along y-axis
1705  {8, 9, 10, 11}}; // Lines along z-axis
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
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);
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: TriaAccessor<0,1,spacedim> -------------------*/
1802 
1803 
1804 template <int spacedim>
1805 bool
1807 {
1809  Assert(false, ExcNotImplemented());
1810  return true;
1811 }
1812 
1813 
1814 
1815 template <int spacedim>
1816 void
1818 {
1820  Assert(false, ExcNotImplemented());
1821 }
1822 
1823 
1824 
1825 template <int spacedim>
1826 void
1828 {
1830  Assert(false, ExcNotImplemented());
1831 }
1832 
1833 
1834 
1835 template <int spacedim>
1836 void
1838 {
1839  set_user_flag();
1840 
1841  if (this->has_children())
1842  for (unsigned int c = 0; c < this->n_children(); ++c)
1843  this->child(c)->recursively_set_user_flag();
1844 }
1845 
1846 
1847 
1848 template <int spacedim>
1849 void
1851 {
1852  clear_user_flag();
1853 
1854  if (this->has_children())
1855  for (unsigned int c = 0; c < this->n_children(); ++c)
1856  this->child(c)->recursively_clear_user_flag();
1857 }
1858 
1859 
1860 
1861 template <int spacedim>
1862 void
1864 {
1866  Assert(false, ExcNotImplemented());
1867 }
1868 
1869 
1870 
1871 template <int spacedim>
1872 void
1874 {
1876  Assert(false, ExcNotImplemented());
1877 }
1878 
1879 
1880 
1881 template <int spacedim>
1882 void
1884 {
1886  Assert(false, ExcNotImplemented());
1887 }
1888 
1889 
1890 
1891 template <int spacedim>
1892 void *
1894 {
1896  Assert(false, ExcNotImplemented());
1897  return nullptr;
1898 }
1899 
1900 
1901 
1902 template <int spacedim>
1903 void
1905 {
1906  set_user_pointer(p);
1907 
1908  if (this->has_children())
1909  for (unsigned int c = 0; c < this->n_children(); ++c)
1910  this->child(c)->recursively_set_user_pointer(p);
1911 }
1912 
1913 
1914 
1915 template <int spacedim>
1916 void
1918 {
1920 
1921  if (this->has_children())
1922  for (unsigned int c = 0; c < this->n_children(); ++c)
1923  this->child(c)->recursively_clear_user_pointer();
1924 }
1925 
1926 
1927 
1928 template <int spacedim>
1929 void
1931 {
1933  Assert(false, ExcNotImplemented());
1934 }
1935 
1936 
1937 
1938 template <int spacedim>
1939 void
1941 {
1943  Assert(false, ExcNotImplemented());
1944 }
1945 
1946 
1947 
1948 template <int spacedim>
1949 unsigned int
1951 {
1953  Assert(false, ExcNotImplemented());
1954  return 0;
1955 }
1956 
1957 
1958 
1959 template <int spacedim>
1960 void
1962 {
1963  set_user_index(p);
1964 
1965  if (this->has_children())
1966  for (unsigned int c = 0; c < this->n_children(); ++c)
1967  this->child(c)->recursively_set_user_index(p);
1968 }
1969 
1970 
1971 
1972 template <int spacedim>
1973 void
1975 {
1976  clear_user_index();
1977 
1978  if (this->has_children())
1979  for (unsigned int c = 0; c < this->n_children(); ++c)
1980  this->child(c)->recursively_clear_user_index();
1981 }
1982 
1983 
1984 
1985 /*------------------------ Functions: CellAccessor<1> -----------------------*/
1986 
1987 
1988 
1989 template <>
1990 bool
1992 {
1993  return (this->vertex(0)[0] <= p[0]) && (p[0] <= this->vertex(1)[0]);
1994 }
1995 
1996 
1997 
1998 /*------------------------ Functions: CellAccessor<2> -----------------------*/
1999 
2000 
2001 
2002 template <>
2003 bool
2005 {
2006  // we check whether the point is
2007  // inside the cell by making sure
2008  // that it on the inner side of
2009  // each line defined by the faces,
2010  // i.e. for each of the four faces
2011  // we take the line that connects
2012  // the two vertices and subdivide
2013  // the whole domain by that in two
2014  // and check whether the point is
2015  // on the `cell-side' (rather than
2016  // the `out-side') of this line. if
2017  // the point is on the `cell-side'
2018  // for all four faces, it must be
2019  // inside the cell.
2020 
2021  // we want the faces in counter
2022  // clockwise orientation
2023  static const int direction[4] = {-1, 1, 1, -1};
2024  for (unsigned int f = 0; f < 4; ++f)
2025  {
2026  // vector from the first vertex
2027  // of the line to the point
2028  const Tensor<1, 2> to_p =
2029  p - this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 0));
2030  // vector describing the line
2031  const Tensor<1, 2> face =
2032  direction[f] *
2033  (this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 1)) -
2034  this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 0)));
2035 
2036  // if we rotate the face vector
2037  // by 90 degrees to the left
2038  // (i.e. it points to the
2039  // inside) and take the scalar
2040  // product with the vector from
2041  // the vertex to the point,
2042  // then the point is in the
2043  // `cell-side' if the scalar
2044  // product is positive. if this
2045  // is not the case, we can be
2046  // sure that the point is
2047  // outside
2048  if ((-face[1] * to_p[0] + face[0] * to_p[1]) < 0)
2049  return false;
2050  }
2051 
2052  // if we arrived here, then the
2053  // point is inside for all four
2054  // faces, and thus inside
2055  return true;
2056 }
2057 
2058 
2059 
2060 /*------------------------ Functions: CellAccessor<3> -----------------------*/
2061 
2062 
2063 
2064 template <>
2065 bool
2067 {
2068  // original implementation by Joerg
2069  // Weimar
2070 
2071  // we first eliminate points based
2072  // on the maximum and minimum of
2073  // the corner coordinates, then
2074  // transform to the unit cell, and
2075  // check there.
2076  const unsigned int dim = 3;
2077  const unsigned int spacedim = 3;
2078  Point<spacedim> maxp = this->vertex(0);
2079  Point<spacedim> minp = this->vertex(0);
2080 
2081  for (unsigned int v = 1; v < this->n_vertices(); ++v)
2082  for (unsigned int d = 0; d < dim; ++d)
2083  {
2084  maxp[d] = std::max(maxp[d], this->vertex(v)[d]);
2085  minp[d] = std::min(minp[d], this->vertex(v)[d]);
2086  }
2087 
2088  // rule out points outside the
2089  // bounding box of this cell
2090  for (unsigned int d = 0; d < dim; ++d)
2091  if ((p[d] < minp[d]) || (p[d] > maxp[d]))
2092  return false;
2093 
2094  // now we need to check more carefully: transform to the
2095  // unit cube and check there. unfortunately, this isn't
2096  // completely trivial since the transform_real_to_unit_cell
2097  // function may throw an exception that indicates that the
2098  // point given could not be inverted. we take this as a sign
2099  // that the point actually lies outside, as also documented
2100  // for that function
2101  try
2102  {
2103  const TriaRawIterator<CellAccessor<dim, spacedim>> cell_iterator(*this);
2105  reference_cell()
2106  .template get_default_linear_mapping<dim, spacedim>()
2107  .transform_real_to_unit_cell(cell_iterator, p)));
2108  }
2110  {
2111  return false;
2112  }
2113 }
2114 
2115 
2116 
2117 /*------------------- Functions: CellAccessor<dim,spacedim> -----------------*/
2118 
2119 // For codim>0 we proceed as follows:
2120 // 1) project point onto manifold and
2121 // 2) transform to the unit cell with a Q1 mapping
2122 // 3) then check if inside unit cell
2123 template <int dim, int spacedim>
2124 template <int dim_, int spacedim_>
2125 bool
2127 {
2128  const TriaRawIterator<CellAccessor<dim_, spacedim_>> cell_iterator(*this);
2129  const Point<dim_> p_unit =
2130  StaticMappingQ1<dim_, spacedim_>::mapping.transform_real_to_unit_cell(
2131  cell_iterator, p);
2132 
2134 }
2135 
2136 
2137 
2138 template <>
2139 bool
2141 {
2142  return point_inside_codim<1, 2>(p);
2143 }
2144 
2145 
2146 template <>
2147 bool
2149 {
2150  return point_inside_codim<1, 3>(p);
2151 }
2152 
2153 
2154 template <>
2155 bool
2157 {
2158  return point_inside_codim<2, 3>(p);
2159 }
2160 
2161 
2162 
2163 template <int dim, int spacedim>
2164 bool
2166 {
2167  for (const auto face : this->face_indices())
2168  if (at_boundary(face))
2169  return true;
2170 
2171  return false;
2172 }
2173 
2174 
2175 
2176 template <int dim, int spacedim>
2179 {
2181  return this->tria->levels[this->present_level]
2182  ->cells.boundary_or_material_id[this->present_index]
2183  .material_id;
2184 }
2185 
2186 
2187 
2188 template <int dim, int spacedim>
2189 void
2191  const types::material_id mat_id) const
2192 {
2195  this->tria->levels[this->present_level]
2196  ->cells.boundary_or_material_id[this->present_index]
2197  .material_id = mat_id;
2198 }
2199 
2200 
2201 
2202 template <int dim, int spacedim>
2203 void
2205  const types::material_id mat_id) const
2206 {
2207  set_material_id(mat_id);
2208 
2209  if (this->has_children())
2210  for (unsigned int c = 0; c < this->n_children(); ++c)
2211  this->child(c)->recursively_set_material_id(mat_id);
2212 }
2213 
2214 
2215 
2216 template <int dim, int spacedim>
2217 void
2219  const types::subdomain_id new_subdomain_id) const
2220 {
2222  Assert(this->is_active(),
2223  ExcMessage("set_subdomain_id() can only be called on active cells!"));
2224  this->tria->levels[this->present_level]->subdomain_ids[this->present_index] =
2225  new_subdomain_id;
2226 }
2227 
2228 
2229 
2230 template <int dim, int spacedim>
2233 {
2235  return this->tria->levels[this->present_level]
2236  ->level_subdomain_ids[this->present_index];
2237 }
2238 
2239 
2240 
2241 template <int dim, int spacedim>
2242 void
2244  const types::subdomain_id new_level_subdomain_id) const
2245 {
2247  this->tria->levels[this->present_level]
2248  ->level_subdomain_ids[this->present_index] = new_level_subdomain_id;
2249 }
2250 
2251 
2252 template <int dim, int spacedim>
2253 bool
2255 {
2257  if (dim == spacedim)
2258  return true;
2259  else
2260  return this->tria->levels[this->present_level]
2261  ->direction_flags[this->present_index];
2262 }
2263 
2264 
2265 
2266 template <int dim, int spacedim>
2267 void
2269  const bool new_direction_flag) const
2270 {
2272  if (dim < spacedim)
2273  this->tria->levels[this->present_level]
2274  ->direction_flags[this->present_index] = new_direction_flag;
2275  else
2276  Assert(new_direction_flag == true,
2277  ExcMessage("If dim==spacedim, direction flags are always true and "
2278  "can not be set to anything else."));
2279 }
2280 
2281 
2282 
2283 template <int dim, int spacedim>
2284 void
2285 CellAccessor<dim, spacedim>::set_parent(const unsigned int parent_index)
2286 {
2288  Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2289  this->tria->levels[this->present_level]->parents[this->present_index / 2] =
2290  parent_index;
2291 }
2292 
2293 
2294 
2295 template <int dim, int spacedim>
2296 int
2298 {
2299  Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2300 
2301  // the parent of two consecutive cells
2302  // is stored only once, since it is
2303  // the same
2304  return this->tria->levels[this->present_level]
2305  ->parents[this->present_index / 2];
2306 }
2307 
2308 
2309 
2310 template <int dim, int spacedim>
2311 void
2313  const unsigned int active_cell_index) const
2314 {
2315  this->tria->levels[this->present_level]
2316  ->active_cell_indices[this->present_index] = active_cell_index;
2317 }
2318 
2319 
2320 
2321 template <int dim, int spacedim>
2322 void
2324  const types::global_cell_index index) const
2325 {
2326  this->tria->levels[this->present_level]
2327  ->global_active_cell_indices[this->present_index] = index;
2328 }
2329 
2330 
2331 
2332 template <int dim, int spacedim>
2333 void
2335  const types::global_cell_index index) const
2336 {
2337  this->tria->levels[this->present_level]
2338  ->global_level_cell_indices[this->present_index] = index;
2339 }
2340 
2341 
2342 
2343 template <int dim, int spacedim>
2346 {
2348  Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2350  this->present_level - 1,
2351  parent_index());
2352 
2353  return q;
2354 }
2355 
2356 
2357 template <int dim, int spacedim>
2358 void
2360  const types::subdomain_id new_subdomain_id) const
2361 {
2362  if (this->has_children())
2363  for (unsigned int c = 0; c < this->n_children(); ++c)
2364  this->child(c)->recursively_set_subdomain_id(new_subdomain_id);
2365  else
2366  set_subdomain_id(new_subdomain_id);
2367 }
2368 
2369 
2370 
2371 template <int dim, int spacedim>
2372 void
2374  const unsigned int i,
2375  const TriaIterator<CellAccessor<dim, spacedim>> &pointer) const
2376 {
2377  AssertIndexRange(i, this->n_faces());
2378 
2379  if (pointer.state() == IteratorState::valid)
2380  {
2381  this->tria->levels[this->present_level]
2382  ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
2383  .first = pointer->present_level;
2384  this->tria->levels[this->present_level]
2385  ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
2386  .second = pointer->present_index;
2387  }
2388  else
2389  {
2390  this->tria->levels[this->present_level]
2391  ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
2392  .first = -1;
2393  this->tria->levels[this->present_level]
2394  ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
2395  .second = -1;
2396  }
2397 }
2398 
2399 
2400 
2401 template <int dim, int spacedim>
2402 CellId
2404 {
2405  std::array<unsigned char, 30> id;
2406 
2407  CellAccessor<dim, spacedim> ptr = *this;
2408  const unsigned int n_child_indices = ptr.level();
2409 
2410  while (ptr.level() > 0)
2411  {
2412  const TriaIterator<CellAccessor<dim, spacedim>> parent = ptr.parent();
2413  const unsigned int n_children = parent->n_children();
2414 
2415  // determine which child we are
2416  unsigned char v = static_cast<unsigned char>(-1);
2417  for (unsigned int c = 0; c < n_children; ++c)
2418  {
2419  if (parent->child_index(c) == ptr.index())
2420  {
2421  v = c;
2422  break;
2423  }
2424  }
2425 
2426  Assert(v != static_cast<unsigned char>(-1), ExcInternalError());
2427  id[ptr.level() - 1] = v;
2428 
2429  ptr.copy_from(*parent);
2430  }
2431 
2432  Assert(ptr.level() == 0, ExcInternalError());
2433  const unsigned int coarse_index = ptr.index();
2434 
2435  return {this->tria->coarse_cell_index_to_coarse_cell_id(coarse_index),
2436  n_child_indices,
2437  id.data()};
2438 }
2439 
2440 
2441 
2442 template <int dim, int spacedim>
2443 unsigned int
2445  const unsigned int neighbor) const
2446 {
2447  AssertIndexRange(neighbor, this->n_faces());
2448 
2449  // if we have a 1d mesh in 1d, we
2450  // can assume that the left
2451  // neighbor of the right neighbor is
2452  // the current cell. but that is an
2453  // invariant that isn't true if the
2454  // mesh is embedded in a higher
2455  // dimensional space, so we have to
2456  // fall back onto the generic code
2457  // below
2458  if ((dim == 1) && (spacedim == dim))
2459  return GeometryInfo<dim>::opposite_face[neighbor];
2460 
2461  const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2462  this->neighbor(neighbor);
2463 
2464  // usually, on regular patches of
2465  // the grid, this cell is just on
2466  // the opposite side of the
2467  // neighbor that the neighbor is of
2468  // this cell. for example in 2d, if
2469  // we want to know the
2470  // neighbor_of_neighbor if
2471  // neighbor==1 (the right
2472  // neighbor), then we will get 3
2473  // (the left neighbor) in most
2474  // cases. look up this relationship
2475  // in the table provided by
2476  // GeometryInfo and try it
2477  const unsigned int this_face_index = face_index(neighbor);
2478 
2479  const unsigned int neighbor_guess =
2481 
2482  if (neighbor_guess < neighbor_cell->n_faces() &&
2483  neighbor_cell->face_index(neighbor_guess) == this_face_index)
2484  return neighbor_guess;
2485  else
2486  // if the guess was false, then
2487  // we need to loop over all
2488  // neighbors and find the number
2489  // the hard way
2490  {
2491  for (const unsigned int face_no : neighbor_cell->face_indices())
2492  if (neighbor_cell->face_index(face_no) == this_face_index)
2493  return face_no;
2494 
2495  // running over all neighbors
2496  // faces we did not find the
2497  // present face. Thereby the
2498  // neighbor must be coarser
2499  // than the present
2500  // cell. Return an invalid
2501  // unsigned int in this case.
2503  }
2504 }
2505 
2506 
2507 
2508 template <int dim, int spacedim>
2509 unsigned int
2511  const unsigned int face_no) const
2512 {
2513  const unsigned int n2 = neighbor_of_neighbor_internal(face_no);
2516 
2517  return n2;
2518 }
2519 
2520 
2521 
2522 template <int dim, int spacedim>
2523 bool
2525  const unsigned int face_no) const
2526 {
2527  return neighbor_of_neighbor_internal(face_no) ==
2529 }
2530 
2531 
2532 
2533 template <int dim, int spacedim>
2534 std::pair<unsigned int, unsigned int>
2536  const unsigned int neighbor) const
2537 {
2538  AssertIndexRange(neighbor, this->n_faces());
2539  // make sure that the neighbor is
2540  // on a coarser level
2541  Assert(neighbor_is_coarser(neighbor),
2543 
2544  switch (dim)
2545  {
2546  case 2:
2547  {
2548  const int this_face_index = face_index(neighbor);
2549  const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2550  this->neighbor(neighbor);
2551 
2552  // usually, on regular patches of
2553  // the grid, this cell is just on
2554  // the opposite side of the
2555  // neighbor that the neighbor is of
2556  // this cell. for example in 2d, if
2557  // we want to know the
2558  // neighbor_of_neighbor if
2559  // neighbor==1 (the right
2560  // neighbor), then we will get 0
2561  // (the left neighbor) in most
2562  // cases. look up this relationship
2563  // in the table provided by
2564  // GeometryInfo and try it
2565  const unsigned int face_no_guess =
2567 
2568  const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face_guess =
2569  neighbor_cell->face(face_no_guess);
2570 
2571  if (face_guess->has_children())
2572  for (unsigned int subface_no = 0;
2573  subface_no < face_guess->n_children();
2574  ++subface_no)
2575  if (face_guess->child_index(subface_no) == this_face_index)
2576  return std::make_pair(face_no_guess, subface_no);
2577 
2578  // if the guess was false, then
2579  // we need to loop over all faces
2580  // and subfaces and find the
2581  // number the hard way
2582  for (const unsigned int face_no : neighbor_cell->face_indices())
2583  {
2584  if (face_no != face_no_guess)
2585  {
2586  const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>
2587  face = neighbor_cell->face(face_no);
2588  if (face->has_children())
2589  for (unsigned int subface_no = 0;
2590  subface_no < face->n_children();
2591  ++subface_no)
2592  if (face->child_index(subface_no) == this_face_index)
2593  return std::make_pair(face_no, subface_no);
2594  }
2595  }
2596 
2597  // we should never get here,
2598  // since then we did not find
2599  // our way back...
2600  Assert(false, ExcInternalError());
2601  return std::make_pair(numbers::invalid_unsigned_int,
2603  }
2604 
2605  case 3:
2606  {
2607  const int this_face_index = face_index(neighbor);
2608  const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2609  this->neighbor(neighbor);
2610 
2611  // usually, on regular patches of the grid, this cell is just on the
2612  // opposite side of the neighbor that the neighbor is of this cell.
2613  // for example in 2d, if we want to know the neighbor_of_neighbor if
2614  // neighbor==1 (the right neighbor), then we will get 0 (the left
2615  // neighbor) in most cases. look up this relationship in the table
2616  // provided by GeometryInfo and try it
2617  const unsigned int face_no_guess =
2619 
2620  const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face_guess =
2621  neighbor_cell->face(face_no_guess);
2622 
2623  if (face_guess->has_children())
2624  for (unsigned int subface_no = 0;
2625  subface_no < face_guess->n_children();
2626  ++subface_no)
2627  {
2628  if (face_guess->child_index(subface_no) == this_face_index)
2629  // call a helper function, that translates the current
2630  // subface number to a subface number for the current
2631  // FaceRefineCase
2632  return std::make_pair(face_no_guess,
2633  translate_subface_no(face_guess,
2634  subface_no));
2635 
2636  if (face_guess->child(subface_no)->has_children())
2637  for (unsigned int subsub_no = 0;
2638  subsub_no < face_guess->child(subface_no)->n_children();
2639  ++subsub_no)
2640  if (face_guess->child(subface_no)->child_index(subsub_no) ==
2641  this_face_index)
2642  // call a helper function, that translates the current
2643  // subface number and subsubface number to a subface
2644  // number for the current FaceRefineCase
2645  return std::make_pair(face_no_guess,
2646  translate_subface_no(face_guess,
2647  subface_no,
2648  subsub_no));
2649  }
2650 
2651  // if the guess was false, then we need to loop over all faces and
2652  // subfaces and find the number the hard way
2653  for (const unsigned int face_no : neighbor_cell->face_indices())
2654  {
2655  if (face_no == face_no_guess)
2656  continue;
2657 
2658  const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face =
2659  neighbor_cell->face(face_no);
2660 
2661  if (!face->has_children())
2662  continue;
2663 
2664  for (unsigned int subface_no = 0; subface_no < face->n_children();
2665  ++subface_no)
2666  {
2667  if (face->child_index(subface_no) == this_face_index)
2668  // call a helper function, that translates the current
2669  // subface number to a subface number for the current
2670  // FaceRefineCase
2671  return std::make_pair(face_no,
2672  translate_subface_no(face,
2673  subface_no));
2674 
2675  if (face->child(subface_no)->has_children())
2676  for (unsigned int subsub_no = 0;
2677  subsub_no < face->child(subface_no)->n_children();
2678  ++subsub_no)
2679  if (face->child(subface_no)->child_index(subsub_no) ==
2680  this_face_index)
2681  // call a helper function, that translates the current
2682  // subface number and subsubface number to a subface
2683  // number for the current FaceRefineCase
2684  return std::make_pair(face_no,
2685  translate_subface_no(face,
2686  subface_no,
2687  subsub_no));
2688  }
2689  }
2690 
2691  // we should never get here, since then we did not find our way
2692  // back...
2693  Assert(false, ExcInternalError());
2694  return std::make_pair(numbers::invalid_unsigned_int,
2696  }
2697 
2698  default:
2699  {
2700  Assert(false, ExcImpossibleInDim(1));
2701  return std::make_pair(numbers::invalid_unsigned_int,
2703  }
2704  }
2705 }
2706 
2707 
2708 
2709 template <int dim, int spacedim>
2710 bool
2712  const unsigned int i_face) const
2713 {
2714  /*
2715  * Implementation note: In all of the functions corresponding to periodic
2716  * faces we mainly use the Triangulation::periodic_face_map to find the
2717  * information about periodically connected faces. So, we actually search in
2718  * this std::map and return the cell_face on the other side of the periodic
2719  * boundary.
2720  *
2721  * We can not use operator[] as this would insert non-existing entries or
2722  * would require guarding with an extra std::map::find() or count().
2723  */
2724  AssertIndexRange(i_face, this->n_faces());
2725  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2726 
2727  cell_iterator current_cell(*this);
2728  if (this->tria->periodic_face_map.find(
2729  std::make_pair(current_cell, i_face)) !=
2730  this->tria->periodic_face_map.end())
2731  return true;
2732  return false;
2733 }
2734 
2735 
2736 
2737 template <int dim, int spacedim>
2739 CellAccessor<dim, spacedim>::periodic_neighbor(const unsigned int i_face) const
2740 {
2741  /*
2742  * To know, why we are using std::map::find() instead of [] operator, refer
2743  * to the implementation note in has_periodic_neighbor() function.
2744  *
2745  * my_it : the iterator to the current cell.
2746  * my_face_pair : the pair reported by periodic_face_map as its first pair
2747  * being the current cell_face.
2748  */
2749  AssertIndexRange(i_face, this->n_faces());
2750  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2751  cell_iterator current_cell(*this);
2752 
2753  auto my_face_pair =
2754  this->tria->periodic_face_map.find(std::make_pair(current_cell, i_face));
2755 
2756  // Make sure we are actually on a periodic boundary:
2757  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2759  return my_face_pair->second.first.first;
2760 }
2761 
2762 
2763 
2764 template <int dim, int spacedim>
2767  const unsigned int i_face) const
2768 {
2769  if (!(this->face(i_face)->at_boundary()))
2770  return this->neighbor(i_face);
2771  else if (this->has_periodic_neighbor(i_face))
2772  return this->periodic_neighbor(i_face);
2773  else
2775  // we can't come here
2776  return this->neighbor(i_face);
2777 }
2778 
2779 
2780 
2781 template <int dim, int spacedim>
2784  const unsigned int i_face,
2785  const unsigned int i_subface) const
2786 {
2787  /*
2788  * To know, why we are using std::map::find() instead of [] operator, refer
2789  * to the implementation note in has_periodic_neighbor() function.
2790  *
2791  * my_it : the iterator to the current cell.
2792  * my_face_pair : the pair reported by periodic_face_map as its first pair
2793  * being the current cell_face. nb_it : the iterator to the
2794  * neighbor of current cell at i_face. face_num_of_nb : the face number of
2795  * the periodically neighboring face in the relevant element.
2796  * nb_parent_face_it: the iterator to the parent face of the periodically
2797  * neighboring face.
2798  */
2799  AssertIndexRange(i_face, this->n_faces());
2800  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2801  cell_iterator my_it(*this);
2802 
2803  auto my_face_pair =
2804  this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2805  /*
2806  * There should be an assertion, which tells the user that this function
2807  * should not be used for a cell which is not located at a periodic boundary.
2808  */
2809  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2811  cell_iterator parent_nb_it = my_face_pair->second.first.first;
2812  unsigned int nb_face_num = my_face_pair->second.first.second;
2813  TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> nb_parent_face_it =
2814  parent_nb_it->face(nb_face_num);
2815  /*
2816  * We should check if the parent face of the neighbor has at least the same
2817  * number of children as i_subface.
2818  */
2819  AssertIndexRange(i_subface, nb_parent_face_it->n_children());
2820  unsigned int sub_neighbor_num =
2821  GeometryInfo<dim>::child_cell_on_face(parent_nb_it->refinement_case(),
2822  nb_face_num,
2823  i_subface,
2824  my_face_pair->second.second[0],
2825  my_face_pair->second.second[1],
2826  my_face_pair->second.second[2],
2827  nb_parent_face_it->refinement_case());
2828  return parent_nb_it->child(sub_neighbor_num);
2829 }
2830 
2831 
2832 
2833 template <int dim, int spacedim>
2834 std::pair<unsigned int, unsigned int>
2836  const unsigned int i_face) const
2837 {
2838  /*
2839  * To know, why we are using std::map::find() instead of [] operator, refer
2840  * to the implementation note in has_periodic_neighbor() function.
2841  *
2842  * my_it : the iterator to the current cell.
2843  * my_face_pair : the pair reported by periodic_face_map as its first pair
2844  * being the current cell_face. nb_it : the iterator to the periodic
2845  * neighbor. nb_face_pair : the pair reported by periodic_face_map as its
2846  * first pair being the periodic neighbor cell_face. p_nb_of_p_nb : the
2847  * iterator of the periodic neighbor of the periodic neighbor of the current
2848  * cell.
2849  */
2850  AssertIndexRange(i_face, this->n_faces());
2851  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2852  const int my_face_index = this->face_index(i_face);
2853  cell_iterator my_it(*this);
2854 
2855  auto my_face_pair =
2856  this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2857  /*
2858  * There should be an assertion, which tells the user that this function
2859  * should not be used for a cell which is not located at a periodic boundary.
2860  */
2861  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2863  cell_iterator nb_it = my_face_pair->second.first.first;
2864  unsigned int face_num_of_nb = my_face_pair->second.first.second;
2865 
2866  auto nb_face_pair =
2867  this->tria->periodic_face_map.find(std::make_pair(nb_it, face_num_of_nb));
2868  /*
2869  * Since, we store periodic neighbors for every cell (either active or
2870  * artificial or inactive) the nb_face_pair should also be mapped to some
2871  * cell_face pair. We assert this here.
2872  */
2873  Assert(nb_face_pair != this->tria->periodic_face_map.end(),
2875  cell_iterator p_nb_of_p_nb = nb_face_pair->second.first.first;
2876  TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> parent_face_it =
2877  p_nb_of_p_nb->face(nb_face_pair->second.first.second);
2878  for (unsigned int i_subface = 0; i_subface < parent_face_it->n_children();
2879  ++i_subface)
2880  if (parent_face_it->child_index(i_subface) == my_face_index)
2881  return std::make_pair(face_num_of_nb, i_subface);
2882  /*
2883  * Obviously, if the execution reaches to this point, some of our assumptions
2884  * should have been false. The most important one is, the user has called this
2885  * function on a face which does not have a coarser periodic neighbor.
2886  */
2888  return std::make_pair(numbers::invalid_unsigned_int,
2890 }
2891 
2892 
2893 
2894 template <int dim, int spacedim>
2895 int
2897  const unsigned int i_face) const
2898 {
2899  return periodic_neighbor(i_face)->index();
2900 }
2901 
2902 
2903 
2904 template <int dim, int spacedim>
2905 int
2907  const unsigned int i_face) const
2908 {
2909  return periodic_neighbor(i_face)->level();
2910 }
2911 
2912 
2913 
2914 template <int dim, int spacedim>
2915 unsigned int
2917  const unsigned int i_face) const
2918 {
2919  return periodic_neighbor_face_no(i_face);
2920 }
2921 
2922 
2923 
2924 template <int dim, int spacedim>
2925 unsigned int
2927  const unsigned int i_face) const
2928 {
2929  /*
2930  * To know, why we are using std::map::find() instead of [] operator, refer
2931  * to the implementation note in has_periodic_neighbor() function.
2932  *
2933  * my_it : the iterator to the current cell.
2934  * my_face_pair : the pair reported by periodic_face_map as its first pair
2935  * being the current cell_face.
2936  */
2937  AssertIndexRange(i_face, this->n_faces());
2938  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2939  cell_iterator my_it(*this);
2940 
2941  auto my_face_pair =
2942  this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2943  /*
2944  * There should be an assertion, which tells the user that this function
2945  * should not be called for a cell which is not located at a periodic boundary
2946  * !
2947  */
2948  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2950  return my_face_pair->second.first.second;
2951 }
2952 
2953 
2954 
2955 template <int dim, int spacedim>
2956 bool
2958  const unsigned int i_face) const
2959 {
2960  /*
2961  * To know, why we are using std::map::find() instead of [] operator, refer
2962  * to the implementation note in has_periodic_neighbor() function.
2963  *
2964  * Implementation note: Let p_nb_of_p_nb be the periodic neighbor of the
2965  * periodic neighbor of the current cell. Also, let p_face_of_p_nb_of_p_nb be
2966  * the periodic face of the p_nb_of_p_nb. If p_face_of_p_nb_of_p_nb has
2967  * children , then the periodic neighbor of the current cell is coarser than
2968  * itself. Although not tested, this implementation should work for
2969  * anisotropic refinement as well.
2970  *
2971  * my_it : the iterator to the current cell.
2972  * my_face_pair : the pair reported by periodic_face_map as its first pair
2973  * being the current cell_face. nb_it : the iterator to the periodic
2974  * neighbor. nb_face_pair : the pair reported by periodic_face_map as its
2975  * first pair being the periodic neighbor cell_face.
2976  */
2977  AssertIndexRange(i_face, this->n_faces());
2978  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2979  cell_iterator my_it(*this);
2980 
2981  auto my_face_pair =
2982  this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2983  /*
2984  * There should be an assertion, which tells the user that this function
2985  * should not be used for a cell which is not located at a periodic boundary.
2986  */
2987  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2989 
2990  cell_iterator nb_it = my_face_pair->second.first.first;
2991  unsigned int face_num_of_nb = my_face_pair->second.first.second;
2992 
2993  auto nb_face_pair =
2994  this->tria->periodic_face_map.find(std::make_pair(nb_it, face_num_of_nb));
2995  /*
2996  * Since, we store periodic neighbors for every cell (either active or
2997  * artificial or inactive) the nb_face_pair should also be mapped to some
2998  * cell_face pair. We assert this here.
2999  */
3000  Assert(nb_face_pair != this->tria->periodic_face_map.end(),
3002  const unsigned int my_level = this->level();
3003  const unsigned int neighbor_level = nb_face_pair->second.first.first->level();
3004  Assert(my_level >= neighbor_level, ExcInternalError());
3005  return my_level > neighbor_level;
3006 }
3007 
3008 
3009 
3010 template <int dim, int spacedim>
3011 bool
3012 CellAccessor<dim, spacedim>::at_boundary(const unsigned int i) const
3013 {
3015  AssertIndexRange(i, this->n_faces());
3016 
3017  return (neighbor_index(i) == -1);
3018 }
3019 
3020 
3021 
3022 template <int dim, int spacedim>
3023 bool
3025 {
3026  if (dim == 1)
3027  return at_boundary();
3028  else
3029  {
3030  for (unsigned int l = 0; l < this->n_lines(); ++l)
3031  if (this->line(l)->at_boundary())
3032  return true;
3033 
3034  return false;
3035  }
3036 }
3037 
3038 
3039 
3040 template <int dim, int spacedim>
3043  const unsigned int face,
3044  const unsigned int subface) const
3045 {
3046  Assert(!this->has_children(),
3047  ExcMessage("The present cell must not have children!"));
3048  Assert(!this->at_boundary(face),
3049  ExcMessage("The present cell must have a valid neighbor!"));
3050  Assert(this->neighbor(face)->has_children() == true,
3051  ExcMessage("The neighbor must have children!"));
3052 
3053  switch (dim)
3054  {
3055  case 2:
3056  {
3058  {
3059  const auto neighbor_cell = this->neighbor(face);
3060 
3061  // only for isotropic refinement at the moment
3062  Assert(neighbor_cell->refinement_case() ==
3064  ExcNotImplemented());
3065 
3066  // determine indices for this cell's subface from the perspective
3067  // of the neighboring cell
3068  const unsigned int neighbor_face =
3069  this->neighbor_of_neighbor(face);
3070  // two neighboring cells have an opposed orientation on their
3071  // shared face if both of them follow the same orientation type
3072  // (i.e., standard or non-standard).
3073  // we verify this with a XOR operation.
3074  const unsigned int neighbor_subface =
3075  (!(this->line_orientation(face)) !=
3076  !(neighbor_cell->line_orientation(neighbor_face))) ?
3077  (1 - subface) :
3078  subface;
3079 
3080  const unsigned int neighbor_child_index =
3081  neighbor_cell->reference_cell().child_cell_on_face(
3082  neighbor_face, neighbor_subface);
3083  const TriaIterator<CellAccessor<dim, spacedim>> sub_neighbor =
3084  neighbor_cell->child(neighbor_child_index);
3085 
3086  // neighbor's child is not allowed to be further refined for the
3087  // moment
3088  Assert(sub_neighbor->refinement_case() ==
3090  ExcNotImplemented());
3091 
3092  return sub_neighbor;
3093  }
3094  else if (this->reference_cell() == ReferenceCells::Quadrilateral)
3095  {
3096  const unsigned int neighbor_neighbor =
3097  this->neighbor_of_neighbor(face);
3098  const unsigned int neighbor_child_index =
3100  this->neighbor(face)->refinement_case(),
3101  neighbor_neighbor,
3102  subface);
3103 
3105  this->neighbor(face)->child(neighbor_child_index);
3106  // the neighbors child can have children,
3107  // which are not further refined along the
3108  // face under consideration. as we are
3109  // normally interested in one of this
3110  // child's child, search for the right one.
3111  while (sub_neighbor->has_children())
3112  {
3114  sub_neighbor->refinement_case(),
3115  neighbor_neighbor) ==
3117  ExcInternalError());
3118  sub_neighbor =
3119  sub_neighbor->child(GeometryInfo<dim>::child_cell_on_face(
3120  sub_neighbor->refinement_case(), neighbor_neighbor, 0));
3121  }
3122 
3123  return sub_neighbor;
3124  }
3125 
3126  // if no reference cell type matches
3127  Assert(false, ExcNotImplemented());
3129  }
3130 
3131 
3132  case 3:
3133  {
3135  {
3136  // this function returns the neighbor's
3137  // child on a given face and
3138  // subface.
3139 
3140  // we have to consider one other aspect here:
3141  // The face might be refined
3142  // anisotropically. In this case, the subface
3143  // number refers to the following, where we
3144  // look at the face from the current cell,
3145  // thus the subfaces are in standard
3146  // orientation concerning the cell
3147  //
3148  // for isotropic refinement
3149  //
3150  // *---*---*
3151  // | 2 | 3 |
3152  // *---*---*
3153  // | 0 | 1 |
3154  // *---*---*
3155  //
3156  // for 2*anisotropic refinement
3157  // (first cut_y, then cut_x)
3158  //
3159  // *---*---*
3160  // | 2 | 3 |
3161  // *---*---*
3162  // | 0 | 1 |
3163  // *---*---*
3164  //
3165  // for 2*anisotropic refinement
3166  // (first cut_x, then cut_y)
3167  //
3168  // *---*---*
3169  // | 1 | 3 |
3170  // *---*---*
3171  // | 0 | 2 |
3172  // *---*---*
3173  //
3174  // for purely anisotropic refinement:
3175  //
3176  // *---*---* *-------*
3177  // | | | | 1 |
3178  // | 0 | 1 | or *-------*
3179  // | | | | 0 |
3180  // *---*---* *-------*
3181  //
3182  // for "mixed" refinement:
3183  //
3184  // *---*---* *---*---* *---*---* *-------*
3185  // | | 2 | | 1 | | | 1 | 2 | | 2 |
3186  // | 0 *---* or *---* 2 | or *---*---* or *---*---*
3187  // | | 1 | | 0 | | | 0 | | 0 | 1 |
3188  // *---*---* *---*---* *-------* *---*---*
3189 
3191  mother_face = this->face(face);
3192  const unsigned int total_children =
3193  mother_face->n_active_descendants();
3194  AssertIndexRange(subface, total_children);
3196  ExcInternalError());
3197 
3198  unsigned int neighbor_neighbor;
3201  this->neighbor(face);
3202 
3203 
3204  const RefinementCase<dim - 1> mother_face_ref_case =
3205  mother_face->refinement_case();
3206  if (mother_face_ref_case ==
3207  static_cast<RefinementCase<dim - 1>>(
3208  RefinementCase<2>::cut_xy)) // total_children==4
3209  {
3210  // this case is quite easy. we are sure,
3211  // that the neighbor is not coarser.
3212 
3213  // get the neighbor's number for the given
3214  // face and the neighbor
3215  neighbor_neighbor = this->neighbor_of_neighbor(face);
3216 
3217  // now use the info provided by GeometryInfo
3218  // to extract the neighbors child number
3219  const unsigned int neighbor_child_index =
3221  neighbor->refinement_case(),
3222  neighbor_neighbor,
3223  subface,
3224  neighbor->face_orientation(neighbor_neighbor),
3225  neighbor->face_flip(neighbor_neighbor),
3226  neighbor->face_rotation(neighbor_neighbor));
3227  neighbor_child = neighbor->child(neighbor_child_index);
3228 
3229  // make sure that the neighbor child cell we
3230  // have found shares the desired subface.
3231  Assert((this->face(face)->child(subface) ==
3232  neighbor_child->face(neighbor_neighbor)),
3233  ExcInternalError());
3234  }
3235  else //-> the face is refined anisotropically
3236  {
3237  // first of all, we have to find the
3238  // neighbor at one of the anisotropic
3239  // children of the
3240  // mother_face. determine, which of
3241  // these we need.
3242  unsigned int first_child_to_find;
3243  unsigned int neighbor_child_index;
3244  if (total_children == 2)
3245  first_child_to_find = subface;
3246  else
3247  {
3248  first_child_to_find = subface / 2;
3249  if (total_children == 3 && subface == 1 &&
3250  !mother_face->child(0)->has_children())
3251  first_child_to_find = 1;
3252  }
3253  if (neighbor_is_coarser(face))
3254  {
3255  std::pair<unsigned int, unsigned int> indices =
3256  neighbor_of_coarser_neighbor(face);
3257  neighbor_neighbor = indices.first;
3258 
3259 
3260  // we have to translate our
3261  // subface_index according to the
3262  // RefineCase and subface index of
3263  // the coarser face (our face is an
3264  // anisotropic child of the coarser
3265  // face), 'a' denotes our
3266  // subface_index 0 and 'b' denotes
3267  // our subface_index 1, whereas 0...3
3268  // denote isotropic subfaces of the
3269  // coarser face
3270  //
3271  // cut_x and coarser_subface_index=0
3272  //
3273  // *---*---*
3274  // |b=2| |
3275  // | | |
3276  // |a=0| |
3277  // *---*---*
3278  //
3279  // cut_x and coarser_subface_index=1
3280  //
3281  // *---*---*
3282  // | |b=3|
3283  // | | |
3284  // | |a=1|
3285  // *---*---*
3286  //
3287  // cut_y and coarser_subface_index=0
3288  //
3289  // *-------*
3290  // | |
3291  // *-------*
3292  // |a=0 b=1|
3293  // *-------*
3294  //
3295  // cut_y and coarser_subface_index=1
3296  //
3297  // *-------*
3298  // |a=2 b=3|
3299  // *-------*
3300  // | |
3301  // *-------*
3302  unsigned int iso_subface;
3303  if (neighbor->face(neighbor_neighbor)
3304  ->refinement_case() == RefinementCase<2>::cut_x)
3305  iso_subface = 2 * first_child_to_find + indices.second;
3306  else
3307  {
3308  Assert(neighbor->face(neighbor_neighbor)
3309  ->refinement_case() ==
3311  ExcInternalError());
3312  iso_subface =
3313  first_child_to_find + 2 * indices.second;
3314  }
3315  neighbor_child_index =
3317  neighbor->refinement_case(),
3318  neighbor_neighbor,
3319  iso_subface,
3320  neighbor->face_orientation(neighbor_neighbor),
3321  neighbor->face_flip(neighbor_neighbor),
3322  neighbor->face_rotation(neighbor_neighbor));
3323  }
3324  else // neighbor is not coarser
3325  {
3326  neighbor_neighbor = neighbor_of_neighbor(face);
3327  neighbor_child_index =
3329  neighbor->refinement_case(),
3330  neighbor_neighbor,
3331  first_child_to_find,
3332  neighbor->face_orientation(neighbor_neighbor),
3333  neighbor->face_flip(neighbor_neighbor),
3334  neighbor->face_rotation(neighbor_neighbor),
3335  mother_face_ref_case);
3336  }
3337 
3338  neighbor_child = neighbor->child(neighbor_child_index);
3339  // it might be, that the neighbor_child
3340  // has children, which are not refined
3341  // along the given subface. go down that
3342  // list and deliver the last of those.
3343  while (
3344  neighbor_child->has_children() &&
3346  neighbor_child->refinement_case(), neighbor_neighbor) ==
3348  neighbor_child = neighbor_child->child(
3350  neighbor_child->refinement_case(),
3351  neighbor_neighbor,
3352  0));
3353 
3354  // if there are two total subfaces, we
3355  // are finished. if there are four we
3356  // have to get a child of our current
3357  // neighbor_child. If there are three,
3358  // we have to check which of the two
3359  // possibilities applies.
3360  if (total_children == 3)
3361  {
3362  if (mother_face->child(0)->has_children())
3363  {
3364  if (subface < 2)
3365  neighbor_child = neighbor_child->child(
3367  neighbor_child->refinement_case(),
3368  neighbor_neighbor,
3369  subface,
3370  neighbor_child->face_orientation(
3371  neighbor_neighbor),
3372  neighbor_child->face_flip(neighbor_neighbor),
3373  neighbor_child->face_rotation(
3374  neighbor_neighbor),
3375  mother_face->child(0)->refinement_case()));
3376  }
3377  else
3378  {
3379  Assert(mother_face->child(1)->has_children(),
3380  ExcInternalError());
3381  if (subface > 0)
3382  neighbor_child = neighbor_child->child(
3384  neighbor_child->refinement_case(),
3385  neighbor_neighbor,
3386  subface - 1,
3387  neighbor_child->face_orientation(
3388  neighbor_neighbor),
3389  neighbor_child->face_flip(neighbor_neighbor),
3390  neighbor_child->face_rotation(
3391  neighbor_neighbor),
3392  mother_face->child(1)->refinement_case()));
3393  }
3394  }
3395  else if (total_children == 4)
3396  {
3397  neighbor_child = neighbor_child->child(
3399  neighbor_child->refinement_case(),
3400  neighbor_neighbor,
3401  subface % 2,
3402  neighbor_child->face_orientation(neighbor_neighbor),
3403  neighbor_child->face_flip(neighbor_neighbor),
3404  neighbor_child->face_rotation(neighbor_neighbor),
3405  mother_face->child(subface / 2)->refinement_case()));
3406  }
3407  }
3408 
3409  // it might be, that the neighbor_child has
3410  // children, which are not refined along the
3411  // given subface. go down that list and
3412  // deliver the last of those.
3413  while (neighbor_child->has_children())
3414  neighbor_child =
3415  neighbor_child->child(GeometryInfo<dim>::child_cell_on_face(
3416  neighbor_child->refinement_case(), neighbor_neighbor, 0));
3417 
3418 #ifdef DEBUG
3419  // check, whether the face neighbor_child matches the requested
3420  // subface.
3422  switch (this->subface_case(face))
3423  {
3427  requested = mother_face->child(subface);
3428  break;
3431  requested =
3432  mother_face->child(subface / 2)->child(subface % 2);
3433  break;
3434 
3437  switch (subface)
3438  {
3439  case 0:
3440  case 1:
3441  requested = mother_face->child(0)->child(subface);
3442  break;
3443  case 2:
3444  requested = mother_face->child(1);
3445  break;
3446  default:
3447  Assert(false, ExcInternalError());
3448  }
3449  break;
3452  switch (subface)
3453  {
3454  case 0:
3455  requested = mother_face->child(0);
3456  break;
3457  case 1:
3458  case 2:
3459  requested = mother_face->child(1)->child(subface - 1);
3460  break;
3461  default:
3462  Assert(false, ExcInternalError());
3463  }
3464  break;
3465  default:
3466  Assert(false, ExcInternalError());
3467  break;
3468  }
3469  Assert(requested == neighbor_child->face(neighbor_neighbor),
3470  ExcInternalError());
3471 #endif
3472 
3473  return neighbor_child;
3474  }
3475 
3476  // if no reference cell type matches
3477  Assert(false, ExcNotImplemented());
3479  }
3480 
3481  default:
3482  // if 1d or more than 3d
3483  Assert(false, ExcNotImplemented());
3485  }
3486 }
3487 
3488 
3489 
3490 // explicit instantiations
3491 #include "tria_accessor.inst"
3492 
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:699
std::size_t size() const
Definition: array_view.h:576
void recursively_set_subdomain_id(const types::subdomain_id new_subdomain_id) const
TriaIterator< CellAccessor< dim, spacedim > > parent() const
void set_active_cell_index(const unsigned int active_cell_index) const
unsigned int neighbor_of_neighbor_internal(const unsigned int neighbor) const
TriaIterator< CellAccessor< dim, spacedim > > periodic_neighbor(const unsigned int i) const
void set_direction_flag(const bool new_direction_flag) const
void recursively_set_material_id(const types::material_id new_material_id) const
void set_level_subdomain_id(const types::subdomain_id new_level_subdomain_id) const
types::subdomain_id level_subdomain_id() const
TriaIterator< CellAccessor< dim, spacedim > > neighbor_child_on_subface(const unsigned int face_no, const unsigned int subface_no) const
void set_subdomain_id(const types::subdomain_id new_subdomain_id) const
bool neighbor_is_coarser(const unsigned int face_no) const
void set_global_level_cell_index(const types::global_cell_index index) const
bool has_periodic_neighbor(const unsigned int i) const
int periodic_neighbor_level(const unsigned int i) const
std::pair< unsigned int, unsigned int > neighbor_of_coarser_neighbor(const unsigned int neighbor) const
unsigned int neighbor_of_neighbor(const unsigned int face_no) const
void set_material_id(const types::material_id new_material_id) const
bool point_inside_codim(const Point< spacedim_ > &p) const
bool has_boundary_lines() const
TriaIterator< CellAccessor< dim, spacedim > > periodic_neighbor_child_on_subface(const unsigned int face_no, const unsigned int subface_no) const
int periodic_neighbor_index(const unsigned int i) const
bool periodic_neighbor_is_coarser(const unsigned int i) const
void set_global_active_cell_index(const types::global_cell_index index) const
void set_neighbor(const unsigned int i, const TriaIterator< CellAccessor< dim, spacedim >> &pointer) const
void set_parent(const unsigned int parent_index)
std::pair< unsigned int, unsigned int > periodic_neighbor_of_coarser_periodic_neighbor(const unsigned face_no) const
bool at_boundary() const
bool point_inside(const Point< spacedim > &p) const
bool direction_flag() const
types::material_id material_id() const
CellId id() const
TriaIterator< CellAccessor< dim, spacedim > > neighbor_or_periodic_neighbor(const unsigned int i) const
int parent_index() const
unsigned int periodic_neighbor_of_periodic_neighbor(const unsigned int i) const
unsigned int periodic_neighbor_face_no(const unsigned int i) const
Definition: cell_id.h:71
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)
DerivativeForm< 1, dim, spacedim, Number > covariant_form() const
Abstract base class for mapping classes.
Definition: mapping.h:311
constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
Definition: tensor.h:503
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:2660
numbers::NumberTraits< Number >::real_type norm() const
void copy_from(const TriaAccessorBase &)
int index() const
int level() const
const Triangulation< dim, spacedim > & get_triangulation() const
void set_user_index(const unsigned int p) const
void clear_user_pointer() const
void recursively_set_user_index(const unsigned int p) const
void clear_user_data() const
Point< spacedim > & vertex(const unsigned int i) const
Point< structdim > real_to_unit_cell_affine_approximation(const Point< spacedim > &point) const
void recursively_clear_user_index() const
void recursively_set_user_pointer(void *p) const
double extent_in_direction(const unsigned int axis) const
Point< spacedim > intermediate_point(const Point< structdim > &coordinates) const
unsigned int n_vertices() const
bool has_children() const
void recursively_clear_user_flag() const
Point< spacedim > barycenter() const
BoundingBox< spacedim > bounding_box() const
void clear_user_flag() const
unsigned int n_children() const
void recursively_set_user_flag() const
bool user_flag_set() const
void set_user_flag() const
unsigned int vertex_index(const unsigned int i) const
Point< spacedim > center(const bool respect_manifold=false, const bool interpolate_from_surrounding=false) const
void clear_user_index() const
double measure() const
void set_bounding_object_indices(const std::initializer_list< int > &new_indices) const
unsigned int user_index() const
void set_user_pointer(void *p) const
void recursively_clear_user_pointer() const
ReferenceCell reference_cell() const
void * user_pointer() const
TriaIterator< TriaAccessor< structdim, dim, spacedim > > child(const unsigned int i) const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices() const
const Manifold< dim, spacedim > & get_manifold() const
bool used() const
std::map< std::pair< cell_iterator, unsigned int >, std::pair< std::pair< cell_iterator, unsigned int >, std::bitset< 3 > > > periodic_face_map
Definition: tria.h:3775
virtual types::coarse_cell_id coarse_cell_index_to_coarse_cell_id(const unsigned int coarse_cell_index) const
std::vector< std::unique_ptr<::internal::TriangulationImplementation::TriaLevel > > levels
Definition: tria.h:4153
const std::vector< Point< spacedim > > & get_vertices() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:416
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:417
Point< 3 > vertices[4]
unsigned int level
Definition: grid_out.cc:4606
unsigned int vertex_indices[2]
Definition: grid_tools.cc:1294
static ::ExceptionBase & ExcCellHasNoParent()
static ::ExceptionBase & ExcNeighborIsNotCoarser()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcCellNotUsed()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNoPeriodicNeighbor()
static ::ExceptionBase & ExcNeighborIsCoarser()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
void set_all_manifold_ids(const types::manifold_id) const
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
double cell_measure< 3 >(const std::vector< Point< 3 >> &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
double cell_measure< 2 >(const std::vector< Point< 2 >> &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
double diameter(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:84
@ valid
Iterator points to a valid object.
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)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:185
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Quadrilateral
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Pyramid
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Hexahedron
Point< 1 > transform_real_to_unit_cell(const std::array< Point< spacedim >, GeometryInfo< 1 >::vertices_per_cell > &vertices, const Point< spacedim > &p)
const types::material_id invalid_material_id
Definition: types.h:233
static const unsigned int invalid_unsigned_int
Definition: types.h:201
unsigned int manifold_id
Definition: types.h:141
unsigned int subdomain_id
Definition: types.h:43
unsigned int global_cell_index
Definition: types.h:105
unsigned int material_id
Definition: types.h:152
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
static unsigned int child_cell_on_face(const RefinementCase< dim > &ref_case, const unsigned int face, const unsigned int subface, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const RefinementCase< dim - 1 > &face_refinement_case=RefinementCase< dim - 1 >::isotropic_refinement)
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
static bool is_inside_unit_cell(const Point< dim > &p)
const ::Triangulation< dim, spacedim > & tria