Reference documentation for deal.II version GIT 0b65fff18a 2023-09-27 19:30:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tria_accessor.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2023 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  const Point<3> x0 = accessor.vertex(0);
1333  const Point<3> x1 = accessor.vertex(1);
1334  const Point<3> x2 = accessor.vertex(2);
1335  const Point<3> x3 = accessor.vertex(3);
1336 
1337  // This is based on the approach used in libMesh (see face_quad4.C): the
1338  // primary differences are the vertex numbering and quadrature order.
1339  //
1340  // The area of a surface is the integral of the magnitude of its normal
1341  // vector, which may be computed via the cross product of two tangent
1342  // vectors. We can easily get tangent vectors from the surface
1343  // parameterization. Hence, given a bilinear surface
1344  //
1345  // X(chi, eta) = x0 + (x1 - x0) chi + (x2 - x0) eta
1346  // + (x3 + x0 - x1 - x2) chi eta
1347  //
1348  // the tangent vectors are
1349  //
1350  // t1 = (x1 - x0) + (x3 + x0 - x1 - x2) eta
1351  // t2 = (x2 - x0) + (x3 + x0 - x1 - x2) xi
1352  const Tensor<1, 3> b0 = x1 - x0;
1353  const Tensor<1, 3> b1 = x2 - x0;
1354  const Tensor<1, 3> a = x3 - x2 - b0;
1355 
1356  // The diameter is the maximum distance between any pair of vertices and
1357  // we can use it as a length scale for the cell. If all components of a
1358  // (the vector connecting x3 and the last vertex of the parallelogram
1359  // defined by the first three vertices) are zero within some tolerance,
1360  // then we have a parallelogram and can use a much simpler formula.
1361  double a_max = 0.0;
1362  for (unsigned int d = 0; d < 3; ++d)
1363  a_max = std::max(std::abs(a[d]), a_max);
1364  if (a_max < 1e-14 * accessor.diameter())
1365  return cross_product_3d(b0, b1).norm();
1366 
1367  // Otherwise, use a 4x4 quadrature to approximate the surface area.
1368  // Hard-code this in to prevent the extra overhead of always creating
1369  // the same QGauss rule.
1370  constexpr unsigned int n_qp = 4;
1371  const double c1 = 2.0 / 7.0 * std::sqrt(6.0 / 5.0);
1372  const double w0 = (18.0 - std::sqrt(30)) / 72.0;
1373  const double w1 = (18.0 + std::sqrt(30)) / 72.0;
1374 
1375  const std::array<double, n_qp> q{{
1376  0.5 - std::sqrt(3.0 / 7.0 + c1) / 2.0,
1377  0.5 - std::sqrt(3.0 / 7.0 - c1) / 2.0,
1378  0.5 + std::sqrt(3.0 / 7.0 - c1) / 2.0,
1379  0.5 + std::sqrt(3.0 / 7.0 + c1) / 2.0,
1380  }};
1381  const std::array<double, n_qp> w{{w0, w1, w1, w0}};
1382 
1383  double area = 0.;
1384  for (unsigned int i = 0; i < n_qp; ++i)
1385  for (unsigned int j = 0; j < n_qp; ++j)
1386  area += cross_product_3d(q[i] * a + b0, q[j] * a + b1).norm() *
1387  w[i] * w[j];
1388 
1389  return area;
1390  }
1391  else if (accessor.reference_cell() == ReferenceCells::Triangle)
1392  {
1393  // We can just use the normal triangle area formula without issue
1394  const Tensor<1, 3> v01 = accessor.vertex(1) - accessor.vertex(0);
1395  const Tensor<1, 3> v02 = accessor.vertex(2) - accessor.vertex(0);
1396  return 0.5 * cross_product_3d(v01, v02).norm();
1397  }
1398 
1399  Assert(false, ExcNotImplemented());
1400  return 0.0;
1401  }
1402 
1403 
1404 
1405  template <int structdim, int dim, int spacedim>
1406  double
1408  {
1409  // catch-all for all cases not explicitly
1410  // listed above
1411  Assert(false, ExcNotImplemented());
1412  return std::numeric_limits<double>::quiet_NaN();
1413  }
1414 
1415 
1416  template <int dim, int spacedim>
1418  get_new_point_on_object(const TriaAccessor<1, dim, spacedim> &obj)
1419  {
1421  return obj.get_manifold().get_new_point_on_line(it);
1422  }
1423 
1424  template <int dim, int spacedim>
1426  get_new_point_on_object(const TriaAccessor<2, dim, spacedim> &obj)
1427  {
1429  return obj.get_manifold().get_new_point_on_quad(it);
1430  }
1431 
1432  template <int dim, int spacedim>
1434  get_new_point_on_object(const TriaAccessor<3, dim, spacedim> &obj)
1435  {
1437  return obj.get_manifold().get_new_point_on_hex(it);
1438  }
1439 
1440  template <int structdim, int dim, int spacedim>
1442  get_new_point_on_object(const TriaAccessor<structdim, dim, spacedim> &obj,
1443  const bool use_interpolation)
1444  {
1445  if (use_interpolation)
1446  {
1448  const auto points_and_weights =
1449  Manifolds::get_default_points_and_weights(it, use_interpolation);
1450  return obj.get_manifold().get_new_point(
1451  make_array_view(points_and_weights.first.begin(),
1452  points_and_weights.first.end()),
1453  make_array_view(points_and_weights.second.begin(),
1454  points_and_weights.second.end()));
1455  }
1456  else
1457  {
1458  return get_new_point_on_object(obj);
1459  }
1460  }
1461 } // namespace
1462 
1463 
1464 
1465 /*-------------------- Static variables: TriaAccessorBase -------------------*/
1466 
1467 template <int structdim, int dim, int spacedim>
1469 
1470 template <int structdim, int dim, int spacedim>
1472 
1473 template <int structdim, int dim, int spacedim>
1474 const unsigned int
1476 
1477 
1478 /*------------------------ Functions: TriaAccessor ---------------------------*/
1479 
1480 template <int structdim, int dim, int spacedim>
1481 void
1483  const std::initializer_list<int> &new_indices) const
1484 {
1485  const ArrayView<int> bounding_object_index_ref =
1486  this->objects().get_bounding_object_indices(this->present_index);
1487 
1488  AssertIndexRange(new_indices.size(), bounding_object_index_ref.size() + 1);
1489 
1490  unsigned int i = 0;
1491  for (const auto &new_index : new_indices)
1492  {
1493  bounding_object_index_ref[i] = new_index;
1494  ++i;
1495  }
1496 }
1497 
1498 
1499 
1500 template <int structdim, int dim, int spacedim>
1501 void
1503  const std::initializer_list<unsigned int> &new_indices) const
1504 {
1505  const ArrayView<int> bounding_object_index_ref =
1506  this->objects().get_bounding_object_indices(this->present_index);
1507 
1508  AssertIndexRange(new_indices.size(), bounding_object_index_ref.size() + 1);
1509 
1510  unsigned int i = 0;
1511  for (const auto &new_index : new_indices)
1512  {
1513  bounding_object_index_ref[i] = new_index;
1514  ++i;
1515  }
1516 }
1517 
1518 
1519 
1520 template <int structdim, int dim, int spacedim>
1523 {
1524  // call the function in the anonymous
1525  // namespace above
1526  return ::barycenter(*this);
1527 }
1528 
1529 
1530 
1531 template <int structdim, int dim, int spacedim>
1532 double
1534 {
1535  // call the function in the anonymous
1536  // namespace above
1537  return ::measure(*this);
1538 }
1539 
1540 
1541 
1542 template <int structdim, int dim, int spacedim>
1545 {
1546  std::pair<Point<spacedim>, Point<spacedim>> boundary_points =
1547  std::make_pair(this->vertex(0), this->vertex(0));
1548 
1549  for (unsigned int v = 1; v < this->n_vertices(); ++v)
1550  {
1551  const Point<spacedim> &x = this->vertex(v);
1552  for (unsigned int k = 0; k < spacedim; ++k)
1553  {
1554  boundary_points.first[k] = std::min(boundary_points.first[k], x[k]);
1555  boundary_points.second[k] = std::max(boundary_points.second[k], x[k]);
1556  }
1557  }
1558 
1559  return BoundingBox<spacedim>(boundary_points);
1560 }
1561 
1562 
1563 
1564 template <int structdim, int dim, int spacedim>
1565 double
1567  const unsigned int /*axis*/) const
1568 {
1570  return std::numeric_limits<double>::signaling_NaN();
1571 }
1572 
1573 
1574 
1575 template <>
1576 double
1577 TriaAccessor<1, 1, 1>::extent_in_direction(const unsigned int axis) const
1578 {
1579  (void)axis;
1580  AssertIndexRange(axis, 1);
1581 
1582  return this->diameter();
1583 }
1584 
1585 
1586 template <>
1587 double
1588 TriaAccessor<1, 1, 2>::extent_in_direction(const unsigned int axis) const
1589 {
1590  (void)axis;
1592 
1593  return this->diameter();
1594 }
1595 
1596 
1597 template <>
1598 double
1599 TriaAccessor<2, 2, 2>::extent_in_direction(const unsigned int axis) const
1600 {
1601  const unsigned int lines[2][2] = {
1602  {2, 3}, // Lines along x-axis, see GeometryInfo
1603  {0, 1}}; // Lines along y-axis
1604 
1605  AssertIndexRange(axis, 2);
1606 
1607  return std::max(this->line(lines[axis][0])->diameter(),
1608  this->line(lines[axis][1])->diameter());
1609 }
1610 
1611 template <>
1612 double
1613 TriaAccessor<2, 2, 3>::extent_in_direction(const unsigned int axis) const
1614 {
1615  const unsigned int lines[2][2] = {
1616  {2, 3}, // Lines along x-axis, see GeometryInfo
1617  {0, 1}}; // Lines along y-axis
1618 
1619  AssertIndexRange(axis, 2);
1620 
1621  return std::max(this->line(lines[axis][0])->diameter(),
1622  this->line(lines[axis][1])->diameter());
1623 }
1624 
1625 
1626 template <>
1627 double
1628 TriaAccessor<3, 3, 3>::extent_in_direction(const unsigned int axis) const
1629 {
1630  const unsigned int lines[3][4] = {
1631  {2, 3, 6, 7}, // Lines along x-axis, see GeometryInfo
1632  {0, 1, 4, 5}, // Lines along y-axis
1633  {8, 9, 10, 11}}; // Lines along z-axis
1634 
1635  AssertIndexRange(axis, 3);
1636 
1637  double lengths[4] = {this->line(lines[axis][0])->diameter(),
1638  this->line(lines[axis][1])->diameter(),
1639  this->line(lines[axis][2])->diameter(),
1640  this->line(lines[axis][3])->diameter()};
1641 
1642  return std::max(std::max(lengths[0], lengths[1]),
1643  std::max(lengths[2], lengths[3]));
1644 }
1645 
1646 
1647 // Recursively set manifold ids on hex iterators.
1648 template <>
1649 void
1651  const types::manifold_id manifold_ind) const
1652 {
1653  set_manifold_id(manifold_ind);
1654 
1655  if (this->has_children())
1656  for (unsigned int c = 0; c < this->n_children(); ++c)
1657  this->child(c)->set_all_manifold_ids(manifold_ind);
1658 
1659  // for hexes also set manifold_id
1660  // of bounding quads and lines
1661 
1662  for (const unsigned int i : this->face_indices())
1663  this->quad(i)->set_manifold_id(manifold_ind);
1664  for (const unsigned int i : this->line_indices())
1665  this->line(i)->set_manifold_id(manifold_ind);
1666 }
1667 
1668 
1669 template <int structdim, int dim, int spacedim>
1672  const Point<structdim> &coordinates) const
1673 {
1674  // Surrounding points and weights.
1675  std::array<Point<spacedim>, GeometryInfo<structdim>::vertices_per_cell> p;
1676  std::array<double, GeometryInfo<structdim>::vertices_per_cell> w;
1677 
1678  for (const unsigned int i : this->vertex_indices())
1679  {
1680  p[i] = this->vertex(i);
1682  }
1683 
1684  return this->get_manifold().get_new_point(make_array_view(p.begin(), p.end()),
1685  make_array_view(w.begin(),
1686  w.end()));
1687 }
1688 
1689 
1690 
1691 template <int structdim, int dim, int spacedim>
1694  const Point<spacedim> &point) const
1695 {
1696  std::array<Point<spacedim>, GeometryInfo<structdim>::vertices_per_cell>
1697  vertices;
1698  for (const unsigned int v : this->vertex_indices())
1699  vertices[v] = this->vertex(v);
1700 
1701  const auto A_b =
1702  GridTools::affine_cell_approximation<structdim, spacedim>(vertices);
1704  A_b.first.covariant_form().transpose();
1705  return Point<structdim>(apply_transformation(A_inv, point - A_b.second));
1706 }
1708 
1709 
1710 template <int structdim, int dim, int spacedim>
1713  const bool respect_manifold,
1714  const bool use_interpolation) const
1715 {
1716  if (respect_manifold == false)
1717  {
1718  Assert(use_interpolation == false, ExcNotImplemented());
1719  Point<spacedim> p;
1720  for (const unsigned int v : this->vertex_indices())
1721  p += vertex(v);
1722  return p / this->n_vertices();
1723  }
1724  else
1725  return get_new_point_on_object(*this, use_interpolation);
1726 }
1727 
1728 
1729 /*---------------- Functions: TriaAccessor<0,1,spacedim> -------------------*/
1730 
1731 
1732 template <int spacedim>
1733 bool
1735 {
1737  Assert(false, ExcNotImplemented());
1738  return true;
1739 }
1740 
1741 
1742 
1743 template <int spacedim>
1744 void
1746 {
1748  Assert(false, ExcNotImplemented());
1749 }
1750 
1751 
1752 
1753 template <int spacedim>
1754 void
1756 {
1758  Assert(false, ExcNotImplemented());
1759 }
1760 
1761 
1762 
1763 template <int spacedim>
1764 void
1766 {
1767  set_user_flag();
1768 
1769  if (this->has_children())
1770  for (unsigned int c = 0; c < this->n_children(); ++c)
1771  this->child(c)->recursively_set_user_flag();
1772 }
1773 
1774 
1775 
1776 template <int spacedim>
1777 void
1779 {
1780  clear_user_flag();
1781 
1782  if (this->has_children())
1783  for (unsigned int c = 0; c < this->n_children(); ++c)
1784  this->child(c)->recursively_clear_user_flag();
1785 }
1786 
1787 
1788 
1789 template <int spacedim>
1790 void
1792 {
1794  Assert(false, ExcNotImplemented());
1795 }
1796 
1797 
1798 
1799 template <int spacedim>
1800 void
1802 {
1804  Assert(false, ExcNotImplemented());
1805 }
1806 
1808 
1809 template <int spacedim>
1810 void
1812 {
1814  Assert(false, ExcNotImplemented());
1815 }
1816 
1817 
1818 
1819 template <int spacedim>
1820 void *
1822 {
1824  Assert(false, ExcNotImplemented());
1825  return nullptr;
1826 }
1827 
1828 
1829 
1830 template <int spacedim>
1831 void
1833 {
1834  set_user_pointer(p);
1835 
1836  if (this->has_children())
1837  for (unsigned int c = 0; c < this->n_children(); ++c)
1838  this->child(c)->recursively_set_user_pointer(p);
1839 }
1840 
1841 
1842 
1843 template <int spacedim>
1844 void
1846 {
1848 
1849  if (this->has_children())
1850  for (unsigned int c = 0; c < this->n_children(); ++c)
1851  this->child(c)->recursively_clear_user_pointer();
1852 }
1853 
1854 
1855 
1856 template <int spacedim>
1857 void
1859 {
1861  Assert(false, ExcNotImplemented());
1862 }
1863 
1864 
1865 
1866 template <int spacedim>
1867 void
1869 {
1871  Assert(false, ExcNotImplemented());
1872 }
1873 
1874 
1875 
1876 template <int spacedim>
1877 unsigned int
1879 {
1881  Assert(false, ExcNotImplemented());
1882  return 0;
1883 }
1884 
1885 
1886 
1887 template <int spacedim>
1888 void
1890 {
1891  set_user_index(p);
1892 
1893  if (this->has_children())
1894  for (unsigned int c = 0; c < this->n_children(); ++c)
1895  this->child(c)->recursively_set_user_index(p);
1896 }
1897 
1898 
1899 
1900 template <int spacedim>
1901 void
1903 {
1904  clear_user_index();
1905 
1906  if (this->has_children())
1907  for (unsigned int c = 0; c < this->n_children(); ++c)
1908  this->child(c)->recursively_clear_user_index();
1909 }
1910 
1911 
1912 
1913 /*------------------------ Functions: CellAccessor<1> -----------------------*/
1914 
1915 
1916 
1917 template <>
1918 bool
1920 {
1921  return (this->vertex(0)[0] <= p[0]) && (p[0] <= this->vertex(1)[0]);
1922 }
1923 
1924 
1925 
1926 /*------------------------ Functions: CellAccessor<2> -----------------------*/
1927 
1928 
1929 
1930 template <>
1931 bool
1933 {
1935  ExcNotImplemented());
1936 
1937  // we check whether the point is
1938  // inside the cell by making sure
1939  // that it on the inner side of
1940  // each line defined by the faces,
1941  // i.e. for each of the four faces
1942  // we take the line that connects
1943  // the two vertices and subdivide
1944  // the whole domain by that in two
1945  // and check whether the point is
1946  // on the `cell-side' (rather than
1947  // the `out-side') of this line. if
1948  // the point is on the `cell-side'
1949  // for all four faces, it must be
1950  // inside the cell.
1951 
1952  // we want the faces in counter
1953  // clockwise orientation
1954  static const int direction[4] = {-1, 1, 1, -1};
1955  for (unsigned int f = 0; f < 4; ++f)
1956  {
1957  // vector from the first vertex
1958  // of the line to the point
1959  const Tensor<1, 2> to_p =
1960  p - this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 0));
1961  // vector describing the line
1962  const Tensor<1, 2> face =
1963  direction[f] *
1964  (this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 1)) -
1965  this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 0)));
1966 
1967  // if we rotate the face vector
1968  // by 90 degrees to the left
1969  // (i.e. it points to the
1970  // inside) and take the scalar
1971  // product with the vector from
1972  // the vertex to the point,
1973  // then the point is in the
1974  // `cell-side' if the scalar
1975  // product is positive. if this
1976  // is not the case, we can be
1977  // sure that the point is
1978  // outside
1979  if ((-face[1] * to_p[0] + face[0] * to_p[1]) < 0)
1980  return false;
1981  }
1982 
1983  // if we arrived here, then the
1984  // point is inside for all four
1985  // faces, and thus inside
1986  return true;
1987 }
1988 
1989 
1990 
1991 /*------------------------ Functions: CellAccessor<3> -----------------------*/
1992 
1993 
1994 
1995 template <>
1996 bool
1998 {
2000  ExcNotImplemented());
2001 
2002  // original implementation by Joerg
2003  // Weimar
2004 
2005  // we first eliminate points based
2006  // on the maximum and minimum of
2007  // the corner coordinates, then
2008  // transform to the unit cell, and
2009  // check there.
2010  const unsigned int dim = 3;
2011  const unsigned int spacedim = 3;
2012  Point<spacedim> maxp = this->vertex(0);
2013  Point<spacedim> minp = this->vertex(0);
2014 
2015  for (unsigned int v = 1; v < this->n_vertices(); ++v)
2016  for (unsigned int d = 0; d < dim; ++d)
2017  {
2018  maxp[d] = std::max(maxp[d], this->vertex(v)[d]);
2019  minp[d] = std::min(minp[d], this->vertex(v)[d]);
2020  }
2021 
2022  // rule out points outside the
2023  // bounding box of this cell
2024  for (unsigned int d = 0; d < dim; ++d)
2025  if ((p[d] < minp[d]) || (p[d] > maxp[d]))
2026  return false;
2027 
2028  // now we need to check more carefully: transform to the
2029  // unit cube and check there. unfortunately, this isn't
2030  // completely trivial since the transform_real_to_unit_cell
2031  // function may throw an exception that indicates that the
2032  // point given could not be inverted. we take this as a sign
2033  // that the point actually lies outside, as also documented
2034  // for that function
2035  try
2036  {
2037  const TriaRawIterator<CellAccessor<dim, spacedim>> cell_iterator(*this);
2039  reference_cell()
2040  .template get_default_linear_mapping<dim, spacedim>()
2041  .transform_real_to_unit_cell(cell_iterator, p)));
2042  }
2044  {
2045  return false;
2046  }
2047 }
2048 
2049 
2050 
2051 /*------------------- Functions: CellAccessor<dim,spacedim> -----------------*/
2052 
2053 // The return type is the same as DoFHandler<dim,spacedim>::active_cell_iterator
2054 template <int dim, int spacedim>
2057  const DoFHandler<dim, spacedim> &dof_handler) const
2058 {
2059  Assert(is_active(),
2060  ExcMessage("The current iterator points to an inactive cell. "
2061  "You cannot convert it to an iterator to an active cell."));
2062  Assert(&this->get_triangulation() == &dof_handler.get_triangulation(),
2063  ExcMessage("The triangulation associated with the iterator does not "
2064  "match that of the DoFHandler."));
2065 
2067  &dof_handler.get_triangulation(),
2068  this->level(),
2069  this->index(),
2070  &dof_handler);
2071 }
2072 
2073 
2074 
2075 template <int dim, int spacedim>
2078  const DoFHandler<dim, spacedim> &dof_handler) const
2079 {
2080  Assert(&this->get_triangulation() == &dof_handler.get_triangulation(),
2081  ExcMessage("The triangulation associated with the iterator does not "
2082  "match that of the DoFHandler."));
2083 
2085  &dof_handler.get_triangulation(),
2086  this->level(),
2087  this->index(),
2088  &dof_handler);
2089 }
2090 
2091 
2092 
2093 // For codim>0 we proceed as follows:
2094 // 1) project point onto manifold and
2095 // 2) transform to the unit cell with a Q1 mapping
2096 // 3) then check if inside unit cell
2097 template <int dim, int spacedim>
2098 template <int dim_, int spacedim_>
2099 bool
2101 {
2102  Assert(this->reference_cell().is_hyper_cube(), ExcNotImplemented());
2103 
2104  const TriaRawIterator<CellAccessor<dim_, spacedim_>> cell_iterator(*this);
2105 
2106  const Point<dim_> p_unit =
2107  this->reference_cell()
2108  .template get_default_linear_mapping<dim_, spacedim_>()
2109  .transform_real_to_unit_cell(cell_iterator, p);
2110 
2112 }
2113 
2114 
2115 
2116 template <>
2117 bool
2119 {
2120  return point_inside_codim<1, 2>(p);
2121 }
2122 
2123 
2124 template <>
2125 bool
2127 {
2128  return point_inside_codim<1, 3>(p);
2129 }
2130 
2131 
2132 template <>
2133 bool
2135 {
2137  ExcNotImplemented());
2138  return point_inside_codim<2, 3>(p);
2139 }
2140 
2141 
2142 
2143 template <int dim, int spacedim>
2144 bool
2146 {
2147  for (const auto face : this->face_indices())
2148  if (at_boundary(face))
2149  return true;
2150 
2151  return false;
2152 }
2153 
2154 
2155 
2156 template <int dim, int spacedim>
2159 {
2161  return this->tria->levels[this->present_level]
2162  ->cells.boundary_or_material_id[this->present_index]
2163  .material_id;
2164 }
2165 
2166 
2167 
2168 template <int dim, int spacedim>
2169 void
2171  const types::material_id mat_id) const
2172 {
2175  this->tria->levels[this->present_level]
2176  ->cells.boundary_or_material_id[this->present_index]
2177  .material_id = mat_id;
2178 }
2179 
2180 
2181 
2182 template <int dim, int spacedim>
2183 void
2185  const types::material_id mat_id) const
2186 {
2187  set_material_id(mat_id);
2188 
2189  if (this->has_children())
2190  for (unsigned int c = 0; c < this->n_children(); ++c)
2191  this->child(c)->recursively_set_material_id(mat_id);
2192 }
2193 
2194 
2195 
2196 template <int dim, int spacedim>
2197 void
2199  const types::subdomain_id new_subdomain_id) const
2200 {
2202  Assert(this->is_active(),
2203  ExcMessage("set_subdomain_id() can only be called on active cells!"));
2204  this->tria->levels[this->present_level]->subdomain_ids[this->present_index] =
2205  new_subdomain_id;
2206 }
2207 
2208 
2209 
2210 template <int dim, int spacedim>
2211 void
2213  const types::subdomain_id new_level_subdomain_id) const
2214 {
2216  this->tria->levels[this->present_level]
2217  ->level_subdomain_ids[this->present_index] = new_level_subdomain_id;
2218 }
2219 
2220 
2221 template <int dim, int spacedim>
2222 bool
2224 {
2226  if (dim == spacedim)
2227  return true;
2228  else
2229  return this->tria->levels[this->present_level]
2230  ->direction_flags[this->present_index];
2231 }
2232 
2233 
2234 
2235 template <int dim, int spacedim>
2236 void
2238  const bool new_direction_flag) const
2239 {
2241  if (dim < spacedim)
2242  this->tria->levels[this->present_level]
2243  ->direction_flags[this->present_index] = new_direction_flag;
2244  else
2245  Assert(new_direction_flag == true,
2246  ExcMessage("If dim==spacedim, direction flags are always true and "
2247  "can not be set to anything else."));
2248 }
2249 
2250 
2251 
2252 template <int dim, int spacedim>
2253 void
2254 CellAccessor<dim, spacedim>::set_parent(const unsigned int parent_index)
2255 {
2257  Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2258  this->tria->levels[this->present_level]->parents[this->present_index / 2] =
2259  parent_index;
2260 }
2261 
2262 
2263 
2264 template <int dim, int spacedim>
2265 int
2267 {
2268  Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2269 
2270  // the parent of two consecutive cells
2271  // is stored only once, since it is
2272  // the same
2273  return this->tria->levels[this->present_level]
2274  ->parents[this->present_index / 2];
2275 }
2276 
2277 
2278 
2279 template <int dim, int spacedim>
2280 void
2282  const unsigned int active_cell_index) const
2283 {
2284  this->tria->levels[this->present_level]
2285  ->active_cell_indices[this->present_index] = active_cell_index;
2286 }
2287 
2288 
2289 
2290 template <int dim, int spacedim>
2291 void
2293  const types::global_cell_index index) const
2294 {
2295  this->tria->levels[this->present_level]
2296  ->global_active_cell_indices[this->present_index] = index;
2297 }
2298 
2299 
2300 
2301 template <int dim, int spacedim>
2302 void
2304  const types::global_cell_index index) const
2305 {
2306  this->tria->levels[this->present_level]
2307  ->global_level_cell_indices[this->present_index] = index;
2308 }
2309 
2310 
2311 
2312 template <int dim, int spacedim>
2315 {
2317  Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2319  this->present_level - 1,
2320  parent_index());
2321 
2322  return q;
2323 }
2324 
2325 
2326 template <int dim, int spacedim>
2327 void
2329  const types::subdomain_id new_subdomain_id) const
2330 {
2331  if (this->has_children())
2332  for (unsigned int c = 0; c < this->n_children(); ++c)
2333  this->child(c)->recursively_set_subdomain_id(new_subdomain_id);
2334  else
2335  set_subdomain_id(new_subdomain_id);
2336 }
2337 
2338 
2339 
2340 template <int dim, int spacedim>
2341 void
2343  const unsigned int i,
2344  const TriaIterator<CellAccessor<dim, spacedim>> &pointer) const
2345 {
2346  AssertIndexRange(i, this->n_faces());
2347 
2348  if (pointer.state() == IteratorState::valid)
2349  {
2350  this->tria->levels[this->present_level]
2351  ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
2352  .first = pointer->present_level;
2353  this->tria->levels[this->present_level]
2354  ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
2355  .second = pointer->present_index;
2356  }
2357  else
2358  {
2359  this->tria->levels[this->present_level]
2360  ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
2361  .first = -1;
2362  this->tria->levels[this->present_level]
2363  ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
2364  .second = -1;
2365  }
2366 }
2367 
2368 
2369 
2370 template <int dim, int spacedim>
2371 CellId
2373 {
2374  std::array<unsigned char, 30> id;
2375 
2376  CellAccessor<dim, spacedim> ptr = *this;
2377  const unsigned int n_child_indices = ptr.level();
2378 
2379  while (ptr.level() > 0)
2380  {
2381  const TriaIterator<CellAccessor<dim, spacedim>> parent = ptr.parent();
2382  const unsigned int n_children = parent->n_children();
2383 
2384  // determine which child we are
2385  unsigned char v = static_cast<unsigned char>(-1);
2386  for (unsigned int c = 0; c < n_children; ++c)
2387  {
2388  if (parent->child_index(c) == ptr.index())
2389  {
2390  v = c;
2391  break;
2392  }
2393  }
2394 
2395  Assert(v != static_cast<unsigned char>(-1), ExcInternalError());
2396  id[ptr.level() - 1] = v;
2397 
2398  ptr.copy_from(*parent);
2399  }
2400 
2401  Assert(ptr.level() == 0, ExcInternalError());
2402  const unsigned int coarse_index = ptr.index();
2403 
2404  return {this->tria->coarse_cell_index_to_coarse_cell_id(coarse_index),
2405  n_child_indices,
2406  id.data()};
2407 }
2408 
2409 
2410 
2411 template <int dim, int spacedim>
2412 unsigned int
2414  const unsigned int neighbor) const
2415 {
2416  AssertIndexRange(neighbor, this->n_faces());
2417 
2418  // if we have a 1d mesh in 1d, we
2419  // can assume that the left
2420  // neighbor of the right neighbor is
2421  // the current cell. but that is an
2422  // invariant that isn't true if the
2423  // mesh is embedded in a higher
2424  // dimensional space, so we have to
2425  // fall back onto the generic code
2426  // below
2427  if ((dim == 1) && (spacedim == dim))
2428  return GeometryInfo<dim>::opposite_face[neighbor];
2429 
2430  const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2431  this->neighbor(neighbor);
2432 
2433  // usually, on regular patches of
2434  // the grid, this cell is just on
2435  // the opposite side of the
2436  // neighbor that the neighbor is of
2437  // this cell. for example in 2d, if
2438  // we want to know the
2439  // neighbor_of_neighbor if
2440  // neighbor==1 (the right
2441  // neighbor), then we will get 3
2442  // (the left neighbor) in most
2443  // cases. look up this relationship
2444  // in the table provided by
2445  // GeometryInfo and try it
2446  const unsigned int this_face_index = face_index(neighbor);
2447 
2448  const unsigned int neighbor_guess =
2450 
2451  if (neighbor_guess < neighbor_cell->n_faces() &&
2452  neighbor_cell->face_index(neighbor_guess) == this_face_index)
2453  return neighbor_guess;
2454  else
2455  // if the guess was false, then
2456  // we need to loop over all
2457  // neighbors and find the number
2458  // the hard way
2459  {
2460  for (const unsigned int face_no : neighbor_cell->face_indices())
2461  if (neighbor_cell->face_index(face_no) == this_face_index)
2462  return face_no;
2463 
2464  // running over all neighbors
2465  // faces we did not find the
2466  // present face. Thereby the
2467  // neighbor must be coarser
2468  // than the present
2469  // cell. Return an invalid
2470  // unsigned int in this case.
2472  }
2473 }
2474 
2475 
2476 
2477 template <int dim, int spacedim>
2478 unsigned int
2480  const unsigned int face_no) const
2481 {
2482  const unsigned int n2 = neighbor_of_neighbor_internal(face_no);
2485 
2486  return n2;
2487 }
2488 
2489 
2490 
2491 template <int dim, int spacedim>
2492 bool
2494  const unsigned int face_no) const
2495 {
2496  return neighbor_of_neighbor_internal(face_no) ==
2498 }
2499 
2500 
2501 
2502 template <int dim, int spacedim>
2503 std::pair<unsigned int, unsigned int>
2505  const unsigned int neighbor) const
2506 {
2507  AssertIndexRange(neighbor, this->n_faces());
2508  // make sure that the neighbor is
2509  // on a coarser level
2510  Assert(neighbor_is_coarser(neighbor),
2512 
2513  switch (dim)
2514  {
2515  case 2:
2516  {
2517  const int this_face_index = face_index(neighbor);
2518  const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2519  this->neighbor(neighbor);
2520 
2521  // usually, on regular patches of
2522  // the grid, this cell is just on
2523  // the opposite side of the
2524  // neighbor that the neighbor is of
2525  // this cell. for example in 2d, if
2526  // we want to know the
2527  // neighbor_of_neighbor if
2528  // neighbor==1 (the right
2529  // neighbor), then we will get 0
2530  // (the left neighbor) in most
2531  // cases. look up this relationship
2532  // in the table provided by
2533  // GeometryInfo and try it
2534  const unsigned int face_no_guess =
2536 
2537  const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face_guess =
2538  neighbor_cell->face(face_no_guess);
2539 
2540  if (face_guess->has_children())
2541  for (unsigned int subface_no = 0;
2542  subface_no < face_guess->n_children();
2543  ++subface_no)
2544  if (face_guess->child_index(subface_no) == this_face_index)
2545  return std::make_pair(face_no_guess, subface_no);
2546 
2547  // if the guess was false, then
2548  // we need to loop over all faces
2549  // and subfaces and find the
2550  // number the hard way
2551  for (const unsigned int face_no : neighbor_cell->face_indices())
2552  {
2553  if (face_no != face_no_guess)
2554  {
2555  const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>
2556  face = neighbor_cell->face(face_no);
2557  if (face->has_children())
2558  for (unsigned int subface_no = 0;
2559  subface_no < face->n_children();
2560  ++subface_no)
2561  if (face->child_index(subface_no) == this_face_index)
2562  return std::make_pair(face_no, subface_no);
2563  }
2564  }
2565 
2566  // we should never get here,
2567  // since then we did not find
2568  // our way back...
2569  Assert(false, ExcInternalError());
2570  return std::make_pair(numbers::invalid_unsigned_int,
2572  }
2573 
2574  case 3:
2575  {
2576  const int this_face_index = face_index(neighbor);
2577  const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2578  this->neighbor(neighbor);
2579 
2580  // usually, on regular patches of the grid, this cell is just on the
2581  // opposite side of the neighbor that the neighbor is of this cell.
2582  // for example in 2d, if we want to know the neighbor_of_neighbor if
2583  // neighbor==1 (the right neighbor), then we will get 0 (the left
2584  // neighbor) in most cases. look up this relationship in the table
2585  // provided by GeometryInfo and try it
2586  const unsigned int face_no_guess =
2588 
2589  const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face_guess =
2590  neighbor_cell->face(face_no_guess);
2591 
2592  if (face_guess->has_children())
2593  for (unsigned int subface_no = 0;
2594  subface_no < face_guess->n_children();
2595  ++subface_no)
2596  {
2597  if (face_guess->child_index(subface_no) == this_face_index)
2598  // call a helper function, that translates the current
2599  // subface number to a subface number for the current
2600  // FaceRefineCase
2601  return std::make_pair(face_no_guess,
2602  translate_subface_no(face_guess,
2603  subface_no));
2604 
2605  if (face_guess->child(subface_no)->has_children())
2606  for (unsigned int subsub_no = 0;
2607  subsub_no < face_guess->child(subface_no)->n_children();
2608  ++subsub_no)
2609  if (face_guess->child(subface_no)->child_index(subsub_no) ==
2610  this_face_index)
2611  // call a helper function, that translates the current
2612  // subface number and subsubface number to a subface
2613  // number for the current FaceRefineCase
2614  return std::make_pair(face_no_guess,
2615  translate_subface_no(face_guess,
2616  subface_no,
2617  subsub_no));
2618  }
2619 
2620  // if the guess was false, then we need to loop over all faces and
2621  // subfaces and find the number the hard way
2622  for (const unsigned int face_no : neighbor_cell->face_indices())
2623  {
2624  if (face_no == face_no_guess)
2625  continue;
2626 
2627  const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face =
2628  neighbor_cell->face(face_no);
2629 
2630  if (!face->has_children())
2631  continue;
2632 
2633  for (unsigned int subface_no = 0; subface_no < face->n_children();
2634  ++subface_no)
2635  {
2636  if (face->child_index(subface_no) == this_face_index)
2637  // call a helper function, that translates the current
2638  // subface number to a subface number for the current
2639  // FaceRefineCase
2640  return std::make_pair(face_no,
2641  translate_subface_no(face,
2642  subface_no));
2643 
2644  if (face->child(subface_no)->has_children())
2645  for (unsigned int subsub_no = 0;
2646  subsub_no < face->child(subface_no)->n_children();
2647  ++subsub_no)
2648  if (face->child(subface_no)->child_index(subsub_no) ==
2649  this_face_index)
2650  // call a helper function, that translates the current
2651  // subface number and subsubface number to a subface
2652  // number for the current FaceRefineCase
2653  return std::make_pair(face_no,
2654  translate_subface_no(face,
2655  subface_no,
2656  subsub_no));
2657  }
2658  }
2659 
2660  // we should never get here, since then we did not find our way
2661  // back...
2662  Assert(false, ExcInternalError());
2663  return std::make_pair(numbers::invalid_unsigned_int,
2665  }
2666 
2667  default:
2668  {
2669  Assert(false, ExcImpossibleInDim(1));
2670  return std::make_pair(numbers::invalid_unsigned_int,
2672  }
2673  }
2674 }
2675 
2676 
2677 
2678 template <int dim, int spacedim>
2679 bool
2681  const unsigned int i_face) const
2682 {
2683  /*
2684  * Implementation note: In all of the functions corresponding to periodic
2685  * faces we mainly use the Triangulation::periodic_face_map to find the
2686  * information about periodically connected faces. So, we actually search in
2687  * this std::map and return the cell_face on the other side of the periodic
2688  * boundary.
2689  *
2690  * We can not use operator[] as this would insert non-existing entries or
2691  * would require guarding with an extra std::map::find() or count().
2692  */
2693  AssertIndexRange(i_face, this->n_faces());
2694  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2695 
2696  cell_iterator current_cell(*this);
2697  if (this->tria->periodic_face_map.find(
2698  std::make_pair(current_cell, i_face)) !=
2699  this->tria->periodic_face_map.end())
2700  return true;
2701  return false;
2702 }
2703 
2704 
2705 
2706 template <int dim, int spacedim>
2708 CellAccessor<dim, spacedim>::periodic_neighbor(const unsigned int i_face) const
2709 {
2710  /*
2711  * To know, why we are using std::map::find() instead of [] operator, refer
2712  * to the implementation note in has_periodic_neighbor() function.
2713  *
2714  * my_it : the iterator to the current cell.
2715  * my_face_pair : the pair reported by periodic_face_map as its first pair
2716  * being the current cell_face.
2717  */
2718  AssertIndexRange(i_face, this->n_faces());
2719  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2720  cell_iterator current_cell(*this);
2721 
2722  auto my_face_pair =
2723  this->tria->periodic_face_map.find(std::make_pair(current_cell, i_face));
2724 
2725  // Make sure we are actually on a periodic boundary:
2726  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2728  return my_face_pair->second.first.first;
2729 }
2730 
2731 
2732 
2733 template <int dim, int spacedim>
2736  const unsigned int i_face) const
2737 {
2738  if (!(this->face(i_face)->at_boundary()))
2739  return this->neighbor(i_face);
2740  else if (this->has_periodic_neighbor(i_face))
2741  return this->periodic_neighbor(i_face);
2742  else
2744  // we can't come here
2745  return this->neighbor(i_face);
2746 }
2747 
2748 
2749 
2750 template <int dim, int spacedim>
2753  const unsigned int i_face,
2754  const unsigned int i_subface) const
2755 {
2756  /*
2757  * To know, why we are using std::map::find() instead of [] operator, refer
2758  * to the implementation note in has_periodic_neighbor() function.
2759  *
2760  * my_it : the iterator to the current cell.
2761  * my_face_pair : the pair reported by periodic_face_map as its first pair
2762  * being the current cell_face. nb_it : the iterator to the
2763  * neighbor of current cell at i_face. face_num_of_nb : the face number of
2764  * the periodically neighboring face in the relevant element.
2765  * nb_parent_face_it: the iterator to the parent face of the periodically
2766  * neighboring face.
2767  */
2768  AssertIndexRange(i_face, this->n_faces());
2769  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2770  cell_iterator my_it(*this);
2771 
2772  auto my_face_pair =
2773  this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2774  /*
2775  * There should be an assertion, which tells the user that this function
2776  * should not be used for a cell which is not located at a periodic boundary.
2777  */
2778  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2780  cell_iterator parent_nb_it = my_face_pair->second.first.first;
2781  unsigned int nb_face_num = my_face_pair->second.first.second;
2782  TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> nb_parent_face_it =
2783  parent_nb_it->face(nb_face_num);
2784  /*
2785  * We should check if the parent face of the neighbor has at least the same
2786  * number of children as i_subface.
2787  */
2788  AssertIndexRange(i_subface, nb_parent_face_it->n_children());
2789  unsigned int sub_neighbor_num =
2790  GeometryInfo<dim>::child_cell_on_face(parent_nb_it->refinement_case(),
2791  nb_face_num,
2792  i_subface,
2793  my_face_pair->second.second[0],
2794  my_face_pair->second.second[1],
2795  my_face_pair->second.second[2],
2796  nb_parent_face_it->refinement_case());
2797  return parent_nb_it->child(sub_neighbor_num);
2798 }
2799 
2800 
2801 
2802 template <int dim, int spacedim>
2803 std::pair<unsigned int, unsigned int>
2805  const unsigned int i_face) const
2806 {
2807  /*
2808  * To know, why we are using std::map::find() instead of [] operator, refer
2809  * to the implementation note in has_periodic_neighbor() function.
2810  *
2811  * my_it : the iterator to the current cell.
2812  * my_face_pair : the pair reported by periodic_face_map as its first pair
2813  * being the current cell_face. nb_it : the iterator to the periodic
2814  * neighbor. nb_face_pair : the pair reported by periodic_face_map as its
2815  * first pair being the periodic neighbor cell_face. p_nb_of_p_nb : the
2816  * iterator of the periodic neighbor of the periodic neighbor of the current
2817  * cell.
2818  */
2819  AssertIndexRange(i_face, this->n_faces());
2820  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2821  const int my_face_index = this->face_index(i_face);
2822  cell_iterator my_it(*this);
2823 
2824  auto my_face_pair =
2825  this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2826  /*
2827  * There should be an assertion, which tells the user that this function
2828  * should not be used for a cell which is not located at a periodic boundary.
2829  */
2830  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2832  cell_iterator nb_it = my_face_pair->second.first.first;
2833  unsigned int face_num_of_nb = my_face_pair->second.first.second;
2834 
2835  auto nb_face_pair =
2836  this->tria->periodic_face_map.find(std::make_pair(nb_it, face_num_of_nb));
2837  /*
2838  * Since, we store periodic neighbors for every cell (either active or
2839  * artificial or inactive) the nb_face_pair should also be mapped to some
2840  * cell_face pair. We assert this here.
2841  */
2842  Assert(nb_face_pair != this->tria->periodic_face_map.end(),
2844  cell_iterator p_nb_of_p_nb = nb_face_pair->second.first.first;
2845  TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> parent_face_it =
2846  p_nb_of_p_nb->face(nb_face_pair->second.first.second);
2847  for (unsigned int i_subface = 0; i_subface < parent_face_it->n_children();
2848  ++i_subface)
2849  if (parent_face_it->child_index(i_subface) == my_face_index)
2850  return std::make_pair(face_num_of_nb, i_subface);
2851  /*
2852  * Obviously, if the execution reaches to this point, some of our assumptions
2853  * should have been false. The most important one is, the user has called this
2854  * function on a face which does not have a coarser periodic neighbor.
2855  */
2857  return std::make_pair(numbers::invalid_unsigned_int,
2859 }
2860 
2861 
2862 
2863 template <int dim, int spacedim>
2864 int
2866  const unsigned int i_face) const
2867 {
2868  return periodic_neighbor(i_face)->index();
2869 }
2870 
2871 
2872 
2873 template <int dim, int spacedim>
2874 int
2876  const unsigned int i_face) const
2877 {
2878  return periodic_neighbor(i_face)->level();
2879 }
2880 
2881 
2882 
2883 template <int dim, int spacedim>
2884 unsigned int
2886  const unsigned int i_face) const
2887 {
2888  return periodic_neighbor_face_no(i_face);
2889 }
2890 
2891 
2892 
2893 template <int dim, int spacedim>
2894 unsigned int
2896  const unsigned int i_face) const
2897 {
2898  /*
2899  * To know, why we are using std::map::find() instead of [] operator, refer
2900  * to the implementation note in has_periodic_neighbor() function.
2901  *
2902  * my_it : the iterator to the current cell.
2903  * my_face_pair : the pair reported by periodic_face_map as its first pair
2904  * being the current cell_face.
2905  */
2906  AssertIndexRange(i_face, this->n_faces());
2907  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2908  cell_iterator my_it(*this);
2909 
2910  auto my_face_pair =
2911  this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2912  /*
2913  * There should be an assertion, which tells the user that this function
2914  * should not be called for a cell which is not located at a periodic boundary
2915  * !
2916  */
2917  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2919  return my_face_pair->second.first.second;
2920 }
2921 
2922 
2923 
2924 template <int dim, int spacedim>
2925 bool
2927  const unsigned int i_face) const
2928 {
2929  /*
2930  * To know, why we are using std::map::find() instead of [] operator, refer
2931  * to the implementation note in has_periodic_neighbor() function.
2932  *
2933  * Implementation note: Let p_nb_of_p_nb be the periodic neighbor of the
2934  * periodic neighbor of the current cell. Also, let p_face_of_p_nb_of_p_nb be
2935  * the periodic face of the p_nb_of_p_nb. If p_face_of_p_nb_of_p_nb has
2936  * children , then the periodic neighbor of the current cell is coarser than
2937  * itself. Although not tested, this implementation should work for
2938  * anisotropic refinement as well.
2939  *
2940  * my_it : the iterator to the current cell.
2941  * my_face_pair : the pair reported by periodic_face_map as its first pair
2942  * being the current cell_face. nb_it : the iterator to the periodic
2943  * neighbor. nb_face_pair : the pair reported by periodic_face_map as its
2944  * first pair being the periodic neighbor cell_face.
2945  */
2946  AssertIndexRange(i_face, this->n_faces());
2947  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2948  cell_iterator my_it(*this);
2949 
2950  auto my_face_pair =
2951  this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2952  /*
2953  * There should be an assertion, which tells the user that this function
2954  * should not be used for a cell which is not located at a periodic boundary.
2955  */
2956  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2958 
2959  cell_iterator nb_it = my_face_pair->second.first.first;
2960  unsigned int face_num_of_nb = my_face_pair->second.first.second;
2961 
2962  auto nb_face_pair =
2963  this->tria->periodic_face_map.find(std::make_pair(nb_it, face_num_of_nb));
2964  /*
2965  * Since, we store periodic neighbors for every cell (either active or
2966  * artificial or inactive) the nb_face_pair should also be mapped to some
2967  * cell_face pair. We assert this here.
2968  */
2969  Assert(nb_face_pair != this->tria->periodic_face_map.end(),
2971  const unsigned int my_level = this->level();
2972  const unsigned int neighbor_level = nb_face_pair->second.first.first->level();
2973  Assert(my_level >= neighbor_level, ExcInternalError());
2974  return my_level > neighbor_level;
2975 }
2976 
2977 
2978 
2979 template <int dim, int spacedim>
2980 bool
2981 CellAccessor<dim, spacedim>::at_boundary(const unsigned int i) const
2982 {
2984  AssertIndexRange(i, this->n_faces());
2985 
2986  return (neighbor_index(i) == -1);
2987 }
2988 
2989 
2990 
2991 template <int dim, int spacedim>
2992 bool
2994 {
2995  if (dim == 1)
2996  return at_boundary();
2997  else
2998  {
2999  for (unsigned int l = 0; l < this->n_lines(); ++l)
3000  if (this->line(l)->at_boundary())
3001  return true;
3002 
3003  return false;
3004  }
3005 }
3006 
3007 
3008 
3009 template <int dim, int spacedim>
3012  const unsigned int face,
3013  const unsigned int subface) const
3014 {
3015  Assert(!this->has_children(),
3016  ExcMessage("The present cell must not have children!"));
3017  Assert(!this->at_boundary(face),
3018  ExcMessage("The present cell must have a valid neighbor!"));
3019  Assert(this->neighbor(face)->has_children() == true,
3020  ExcMessage("The neighbor must have children!"));
3021 
3022  switch (dim)
3023  {
3024  case 2:
3025  {
3027  {
3028  const auto neighbor_cell = this->neighbor(face);
3029 
3030  // only for isotropic refinement at the moment
3031  Assert(neighbor_cell->refinement_case() ==
3033  ExcNotImplemented());
3034 
3035  // determine indices for this cell's subface from the perspective
3036  // of the neighboring cell
3037  const unsigned int neighbor_face =
3038  this->neighbor_of_neighbor(face);
3039  // two neighboring cells have an opposed orientation on their
3040  // shared face if both of them follow the same orientation type
3041  // (i.e., standard or non-standard).
3042  // we verify this with a XOR operation.
3043  const unsigned int neighbor_subface =
3044  (!(this->line_orientation(face)) !=
3045  !(neighbor_cell->line_orientation(neighbor_face))) ?
3046  (1 - subface) :
3047  subface;
3048 
3049  const unsigned int neighbor_child_index =
3050  neighbor_cell->reference_cell().child_cell_on_face(
3051  neighbor_face, neighbor_subface);
3052  const TriaIterator<CellAccessor<dim, spacedim>> sub_neighbor =
3053  neighbor_cell->child(neighbor_child_index);
3054 
3055  // neighbor's child is not allowed to be further refined for the
3056  // moment
3057  Assert(sub_neighbor->refinement_case() ==
3059  ExcNotImplemented());
3060 
3061  return sub_neighbor;
3062  }
3063  else if (this->reference_cell() == ReferenceCells::Quadrilateral)
3064  {
3065  const unsigned int neighbor_neighbor =
3066  this->neighbor_of_neighbor(face);
3067  const unsigned int neighbor_child_index =
3069  this->neighbor(face)->refinement_case(),
3070  neighbor_neighbor,
3071  subface);
3072 
3074  this->neighbor(face)->child(neighbor_child_index);
3075  // the neighbors child can have children,
3076  // which are not further refined along the
3077  // face under consideration. as we are
3078  // normally interested in one of this
3079  // child's child, search for the right one.
3080  while (sub_neighbor->has_children())
3081  {
3083  sub_neighbor->refinement_case(),
3084  neighbor_neighbor) ==
3086  ExcInternalError());
3087  sub_neighbor =
3088  sub_neighbor->child(GeometryInfo<dim>::child_cell_on_face(
3089  sub_neighbor->refinement_case(), neighbor_neighbor, 0));
3090  }
3091 
3092  return sub_neighbor;
3093  }
3094 
3095  // if no reference cell type matches
3096  Assert(false, ExcNotImplemented());
3098  }
3099 
3100 
3101  case 3:
3102  {
3104  {
3105  // this function returns the neighbor's
3106  // child on a given face and
3107  // subface.
3108 
3109  // we have to consider one other aspect here:
3110  // The face might be refined
3111  // anisotropically. In this case, the subface
3112  // number refers to the following, where we
3113  // look at the face from the current cell,
3114  // thus the subfaces are in standard
3115  // orientation concerning the cell
3116  //
3117  // for isotropic refinement
3118  //
3119  // *---*---*
3120  // | 2 | 3 |
3121  // *---*---*
3122  // | 0 | 1 |
3123  // *---*---*
3124  //
3125  // for 2*anisotropic refinement
3126  // (first cut_y, then cut_x)
3127  //
3128  // *---*---*
3129  // | 2 | 3 |
3130  // *---*---*
3131  // | 0 | 1 |
3132  // *---*---*
3133  //
3134  // for 2*anisotropic refinement
3135  // (first cut_x, then cut_y)
3136  //
3137  // *---*---*
3138  // | 1 | 3 |
3139  // *---*---*
3140  // | 0 | 2 |
3141  // *---*---*
3142  //
3143  // for purely anisotropic refinement:
3144  //
3145  // *---*---* *-------*
3146  // | | | | 1 |
3147  // | 0 | 1 | or *-------*
3148  // | | | | 0 |
3149  // *---*---* *-------*
3150  //
3151  // for "mixed" refinement:
3152  //
3153  // *---*---* *---*---* *---*---* *-------*
3154  // | | 2 | | 1 | | | 1 | 2 | | 2 |
3155  // | 0 *---* or *---* 2 | or *---*---* or *---*---*
3156  // | | 1 | | 0 | | | 0 | | 0 | 1 |
3157  // *---*---* *---*---* *-------* *---*---*
3158 
3160  mother_face = this->face(face);
3161  const unsigned int total_children =
3162  mother_face->n_active_descendants();
3163  AssertIndexRange(subface, total_children);
3165  ExcInternalError());
3166 
3167  unsigned int neighbor_neighbor;
3170  this->neighbor(face);
3171 
3172 
3173  const RefinementCase<dim - 1> mother_face_ref_case =
3174  mother_face->refinement_case();
3175  if (mother_face_ref_case ==
3176  static_cast<RefinementCase<dim - 1>>(
3177  RefinementCase<2>::cut_xy)) // total_children==4
3178  {
3179  // this case is quite easy. we are sure,
3180  // that the neighbor is not coarser.
3181 
3182  // get the neighbor's number for the given
3183  // face and the neighbor
3184  neighbor_neighbor = this->neighbor_of_neighbor(face);
3185 
3186  // now use the info provided by GeometryInfo
3187  // to extract the neighbors child number
3188  const unsigned int neighbor_child_index =
3190  neighbor->refinement_case(),
3191  neighbor_neighbor,
3192  subface,
3193  neighbor->face_orientation(neighbor_neighbor),
3194  neighbor->face_flip(neighbor_neighbor),
3195  neighbor->face_rotation(neighbor_neighbor));
3196  neighbor_child = neighbor->child(neighbor_child_index);
3197 
3198  // make sure that the neighbor child cell we
3199  // have found shares the desired subface.
3200  Assert((this->face(face)->child(subface) ==
3201  neighbor_child->face(neighbor_neighbor)),
3202  ExcInternalError());
3203  }
3204  else //-> the face is refined anisotropically
3205  {
3206  // first of all, we have to find the
3207  // neighbor at one of the anisotropic
3208  // children of the
3209  // mother_face. determine, which of
3210  // these we need.
3211  unsigned int first_child_to_find;
3212  unsigned int neighbor_child_index;
3213  if (total_children == 2)
3214  first_child_to_find = subface;
3215  else
3216  {
3217  first_child_to_find = subface / 2;
3218  if (total_children == 3 && subface == 1 &&
3219  !mother_face->child(0)->has_children())
3220  first_child_to_find = 1;
3221  }
3222  if (neighbor_is_coarser(face))
3223  {
3224  std::pair<unsigned int, unsigned int> indices =
3225  neighbor_of_coarser_neighbor(face);
3226  neighbor_neighbor = indices.first;
3227 
3228 
3229  // we have to translate our
3230  // subface_index according to the
3231  // RefineCase and subface index of
3232  // the coarser face (our face is an
3233  // anisotropic child of the coarser
3234  // face), 'a' denotes our
3235  // subface_index 0 and 'b' denotes
3236  // our subface_index 1, whereas 0...3
3237  // denote isotropic subfaces of the
3238  // coarser face
3239  //
3240  // cut_x and coarser_subface_index=0
3241  //
3242  // *---*---*
3243  // |b=2| |
3244  // | | |
3245  // |a=0| |
3246  // *---*---*
3247  //
3248  // cut_x and coarser_subface_index=1
3249  //
3250  // *---*---*
3251  // | |b=3|
3252  // | | |
3253  // | |a=1|
3254  // *---*---*
3255  //
3256  // cut_y and coarser_subface_index=0
3257  //
3258  // *-------*
3259  // | |
3260  // *-------*
3261  // |a=0 b=1|
3262  // *-------*
3263  //
3264  // cut_y and coarser_subface_index=1
3265  //
3266  // *-------*
3267  // |a=2 b=3|
3268  // *-------*
3269  // | |
3270  // *-------*
3271  unsigned int iso_subface;
3272  if (neighbor->face(neighbor_neighbor)
3273  ->refinement_case() == RefinementCase<2>::cut_x)
3274  iso_subface = 2 * first_child_to_find + indices.second;
3275  else
3276  {
3277  Assert(neighbor->face(neighbor_neighbor)
3278  ->refinement_case() ==
3280  ExcInternalError());
3281  iso_subface =
3282  first_child_to_find + 2 * indices.second;
3283  }
3284  neighbor_child_index =
3286  neighbor->refinement_case(),
3287  neighbor_neighbor,
3288  iso_subface,
3289  neighbor->face_orientation(neighbor_neighbor),
3290  neighbor->face_flip(neighbor_neighbor),
3291  neighbor->face_rotation(neighbor_neighbor));
3292  }
3293  else // neighbor is not coarser
3294  {
3295  neighbor_neighbor = neighbor_of_neighbor(face);
3296  neighbor_child_index =
3298  neighbor->refinement_case(),
3299  neighbor_neighbor,
3300  first_child_to_find,
3301  neighbor->face_orientation(neighbor_neighbor),
3302  neighbor->face_flip(neighbor_neighbor),
3303  neighbor->face_rotation(neighbor_neighbor),
3304  mother_face_ref_case);
3305  }
3306 
3307  neighbor_child = neighbor->child(neighbor_child_index);
3308  // it might be, that the neighbor_child
3309  // has children, which are not refined
3310  // along the given subface. go down that
3311  // list and deliver the last of those.
3312  while (
3313  neighbor_child->has_children() &&
3315  neighbor_child->refinement_case(), neighbor_neighbor) ==
3317  neighbor_child = neighbor_child->child(
3319  neighbor_child->refinement_case(),
3320  neighbor_neighbor,
3321  0));
3322 
3323  // if there are two total subfaces, we
3324  // are finished. if there are four we
3325  // have to get a child of our current
3326  // neighbor_child. If there are three,
3327  // we have to check which of the two
3328  // possibilities applies.
3329  if (total_children == 3)
3330  {
3331  if (mother_face->child(0)->has_children())
3332  {
3333  if (subface < 2)
3334  neighbor_child = neighbor_child->child(
3336  neighbor_child->refinement_case(),
3337  neighbor_neighbor,
3338  subface,
3339  neighbor_child->face_orientation(
3340  neighbor_neighbor),
3341  neighbor_child->face_flip(neighbor_neighbor),
3342  neighbor_child->face_rotation(
3343  neighbor_neighbor),
3344  mother_face->child(0)->refinement_case()));
3345  }
3346  else
3347  {
3348  Assert(mother_face->child(1)->has_children(),
3349  ExcInternalError());
3350  if (subface > 0)
3351  neighbor_child = neighbor_child->child(
3353  neighbor_child->refinement_case(),
3354  neighbor_neighbor,
3355  subface - 1,
3356  neighbor_child->face_orientation(
3357  neighbor_neighbor),
3358  neighbor_child->face_flip(neighbor_neighbor),
3359  neighbor_child->face_rotation(
3360  neighbor_neighbor),
3361  mother_face->child(1)->refinement_case()));
3362  }
3363  }
3364  else if (total_children == 4)
3365  {
3366  neighbor_child = neighbor_child->child(
3368  neighbor_child->refinement_case(),
3369  neighbor_neighbor,
3370  subface % 2,
3371  neighbor_child->face_orientation(neighbor_neighbor),
3372  neighbor_child->face_flip(neighbor_neighbor),
3373  neighbor_child->face_rotation(neighbor_neighbor),
3374  mother_face->child(subface / 2)->refinement_case()));
3375  }
3376  }
3377 
3378  // it might be, that the neighbor_child has
3379  // children, which are not refined along the
3380  // given subface. go down that list and
3381  // deliver the last of those.
3382  while (neighbor_child->has_children())
3383  neighbor_child =
3384  neighbor_child->child(GeometryInfo<dim>::child_cell_on_face(
3385  neighbor_child->refinement_case(), neighbor_neighbor, 0));
3386 
3387 #ifdef DEBUG
3388  // check, whether the face neighbor_child matches the requested
3389  // subface.
3391  switch (this->subface_case(face))
3392  {
3396  requested = mother_face->child(subface);
3397  break;
3400  requested =
3401  mother_face->child(subface / 2)->child(subface % 2);
3402  break;
3403 
3406  switch (subface)
3407  {
3408  case 0:
3409  case 1:
3410  requested = mother_face->child(0)->child(subface);
3411  break;
3412  case 2:
3413  requested = mother_face->child(1);
3414  break;
3415  default:
3416  Assert(false, ExcInternalError());
3417  }
3418  break;
3421  switch (subface)
3422  {
3423  case 0:
3424  requested = mother_face->child(0);
3425  break;
3426  case 1:
3427  case 2:
3428  requested = mother_face->child(1)->child(subface - 1);
3429  break;
3430  default:
3431  Assert(false, ExcInternalError());
3432  }
3433  break;
3434  default:
3435  Assert(false, ExcInternalError());
3436  break;
3437  }
3438  Assert(requested == neighbor_child->face(neighbor_neighbor),
3439  ExcInternalError());
3440 #endif
3441 
3442  return neighbor_child;
3443  }
3444 
3445  // if no reference cell type matches
3446  Assert(false, ExcNotImplemented());
3448  }
3449 
3450  default:
3451  // if 1d or more than 3d
3452  Assert(false, ExcNotImplemented());
3454  }
3455 }
3456 
3457 
3458 
3459 template <int structdim, int dim, int spacedim>
3462 {
3463  return IteratorState::invalid;
3464 }
3465 
3466 
3467 
3468 template <int structdim, int dim, int spacedim>
3469 int
3471 {
3472  return -1;
3473 }
3474 
3475 
3476 
3477 template <int structdim, int dim, int spacedim>
3478 int
3480 {
3481  return -1;
3482 }
3483 
3484 
3485 // explicit instantiations
3486 #include "tria_accessor.inst"
3487 
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:838
std::size_t size() const
Definition: array_view.h:573
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< DoFCellAccessor< dim, spacedim, true > > as_dof_handler_level_iterator(const DoFHandler< dim, spacedim > &dof_handler) 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
typename LevelSelector::cell_iterator level_cell_iterator
Definition: dof_handler.h:505
static int level()
static IteratorState::IteratorStates state()
static int index()
Abstract base class for mapping classes.
Definition: mapping.h:317
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
double diameter() 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:4094
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:4466
const std::vector< Point< spacedim > > & get_vertices() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
Point< 3 > vertices[4]
unsigned int level
Definition: grid_out.cc:4617
unsigned int vertex_indices[2]
Definition: grid_tools.cc:1356
static ::ExceptionBase & ExcCellHasNoParent()
static ::ExceptionBase & ExcNeighborIsNotCoarser()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcCellNotUsed()
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNoPeriodicNeighbor()
static ::ExceptionBase & ExcNeighborIsCoarser()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:441
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:111
@ 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:192
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)
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