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