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