Reference documentation for deal.II version Git 54f37b890c 2020-02-28 19:52:28 +0100
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tria_accessor.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2019 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 
16 #include <deal.II/base/geometry_info.h>
17 #include <deal.II/base/quadrature.h>
18 
19 #include <deal.II/fe/fe_q.h>
20 #include <deal.II/fe/mapping_q1.h>
21 
22 #include <deal.II/grid/grid_tools.h>
23 #include <deal.II/grid/manifold.h>
24 #include <deal.II/grid/tria.h>
25 #include <deal.II/grid/tria_accessor.h>
26 #include <deal.II/grid/tria_accessor.templates.h>
27 #include <deal.II/grid/tria_iterator.h>
28 #include <deal.II/grid/tria_iterator.templates.h>
29 #include <deal.II/grid/tria_levels.h>
30 
31 #include <array>
32 #include <cmath>
33 
34 DEAL_II_NAMESPACE_OPEN
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  {
1267  unsigned int vertex_indices[GeometryInfo<2>::vertices_per_cell];
1268  for (unsigned int i = 0; i < GeometryInfo<2>::vertices_per_cell; ++i)
1269  vertex_indices[i] = accessor.vertex_index(i);
1270 
1271  return GridTools::cell_measure<2>(
1272  accessor.get_triangulation().get_vertices(), vertex_indices);
1273  }
1274 
1275 
1276  double
1277  measure(const TriaAccessor<3, 3, 3> &accessor)
1278  {
1279  unsigned int vertex_indices[GeometryInfo<3>::vertices_per_cell];
1280  for (unsigned int i = 0; i < GeometryInfo<3>::vertices_per_cell; ++i)
1281  vertex_indices[i] = accessor.vertex_index(i);
1282 
1283  return GridTools::cell_measure<3>(
1284  accessor.get_triangulation().get_vertices(), vertex_indices);
1285  }
1286 
1287 
1288  // a 2d face in 3d space
1289  double
1290  measure(const ::TriaAccessor<2, 3, 3> &accessor)
1291  {
1292  // If the face is planar, the diagonal from vertex 0 to vertex 3,
1293  // v_03, should be in the plane P_012 of vertices 0, 1 and 2. Get
1294  // the normal vector of P_012 and test if v_03 is orthogonal to
1295  // that. If so, the face is planar and computing its area is simple.
1296  const Tensor<1, 3> v01 = accessor.vertex(1) - accessor.vertex(0);
1297  const Tensor<1, 3> v02 = accessor.vertex(2) - accessor.vertex(0);
1298 
1299  const Tensor<1, 3> normal = cross_product_3d(v01, v02);
1300 
1301  const Tensor<1, 3> v03 = accessor.vertex(3) - accessor.vertex(0);
1302 
1303  // check whether v03 does not lie in the plane of v01 and v02
1304  // (i.e., whether the face is not planar). we do so by checking
1305  // whether the triple product (v01 x v02) * v03 forms a positive
1306  // volume relative to |v01|*|v02|*|v03|. the test checks the
1307  // squares of these to avoid taking norms/square roots:
1308  if (std::abs((v03 * normal) * (v03 * normal) /
1309  ((v03 * v03) * (v01 * v01) * (v02 * v02))) >= 1e-24)
1310  {
1311  Assert(
1312  false,
1313  ExcMessage(
1314  "Computing the measure of a nonplanar face is not implemented!"));
1315  return std::numeric_limits<double>::quiet_NaN();
1316  }
1317 
1318  // the face is planar. then its area is 1/2 of the norm of the
1319  // cross product of the two diagonals
1320  const Tensor<1, 3> v12 = accessor.vertex(2) - accessor.vertex(1);
1321  const Tensor<1, 3> twice_area = cross_product_3d(v03, v12);
1322  return 0.5 * twice_area.norm();
1323  }
1324 
1325 
1326 
1327  template <int structdim, int dim, int spacedim>
1328  double
1330  {
1331  // catch-all for all cases not explicitly
1332  // listed above
1333  Assert(false, ExcNotImplemented());
1334  return std::numeric_limits<double>::quiet_NaN();
1335  }
1336 
1337 
1338  template <int dim, int spacedim>
1340  get_new_point_on_object(const TriaAccessor<1, dim, spacedim> &obj)
1341  {
1343  return obj.get_manifold().get_new_point_on_line(it);
1344  }
1345 
1346  template <int dim, int spacedim>
1348  get_new_point_on_object(const TriaAccessor<2, dim, spacedim> &obj)
1349  {
1351  return obj.get_manifold().get_new_point_on_quad(it);
1352  }
1353 
1354  template <int dim, int spacedim>
1356  get_new_point_on_object(const TriaAccessor<3, dim, spacedim> &obj)
1357  {
1359  return obj.get_manifold().get_new_point_on_hex(it);
1360  }
1361 
1362  template <int structdim, int dim, int spacedim>
1364  get_new_point_on_object(const TriaAccessor<structdim, dim, spacedim> &obj,
1365  const bool use_interpolation)
1366  {
1367  if (use_interpolation)
1368  {
1370  const auto points_and_weights =
1371  Manifolds::get_default_points_and_weights(it, use_interpolation);
1372  return obj.get_manifold().get_new_point(
1373  make_array_view(points_and_weights.first.begin(),
1374  points_and_weights.first.end()),
1375  make_array_view(points_and_weights.second.begin(),
1376  points_and_weights.second.end()));
1377  }
1378  else
1379  {
1380  return get_new_point_on_object(obj);
1381  }
1382  }
1383 } // namespace
1384 
1385 
1386 
1387 /*-------------------- Static variables: TriaAccessorBase -------------------*/
1388 
1389 template <int structdim, int dim, int spacedim>
1391 
1392 template <int structdim, int dim, int spacedim>
1394 
1395 template <int structdim, int dim, int spacedim>
1396 const unsigned int
1398 
1399 
1400 /*------------------------ Functions: TriaAccessor ---------------------------*/
1401 
1402 template <int structdim, int dim, int spacedim>
1403 void
1406  const
1407 {
1408  this->objects().cells[this->present_index] = object;
1409 }
1410 
1411 
1412 
1413 template <int structdim, int dim, int spacedim>
1416 {
1417  // call the function in the anonymous
1418  // namespace above
1419  return ::barycenter(*this);
1420 }
1421 
1422 
1423 
1424 template <int structdim, int dim, int spacedim>
1425 double
1427 {
1428  // call the function in the anonymous
1429  // namespace above
1430  return ::measure(*this);
1431 }
1432 
1433 
1434 
1435 template <int structdim, int dim, int spacedim>
1438 {
1439  std::pair<Point<spacedim>, Point<spacedim>> boundary_points =
1440  std::make_pair(this->vertex(0), this->vertex(0));
1441 
1442  for (unsigned int v = 1; v < GeometryInfo<structdim>::vertices_per_cell; ++v)
1443  {
1444  const Point<spacedim> &x = this->vertex(v);
1445  for (unsigned int k = 0; k < spacedim; ++k)
1446  {
1447  boundary_points.first[k] = std::min(boundary_points.first[k], x[k]);
1448  boundary_points.second[k] = std::max(boundary_points.second[k], x[k]);
1449  }
1450  }
1451 
1452  return BoundingBox<spacedim>(boundary_points);
1453 }
1454 
1455 
1456 
1457 template <int structdim, int dim, int spacedim>
1458 double
1460  const unsigned int /*axis*/) const
1461 {
1462  Assert(false, ExcNotImplemented());
1463  return std::numeric_limits<double>::signaling_NaN();
1464 }
1465 
1466 
1467 
1468 template <>
1469 double
1470 TriaAccessor<1, 1, 1>::extent_in_direction(const unsigned int axis) const
1471 {
1472  (void)axis;
1473  AssertIndexRange(axis, 1);
1474 
1475  return this->diameter();
1476 }
1477 
1478 
1479 template <>
1480 double
1481 TriaAccessor<1, 1, 2>::extent_in_direction(const unsigned int axis) const
1482 {
1483  (void)axis;
1484  AssertIndexRange(axis, 1);
1485 
1486  return this->diameter();
1487 }
1488 
1489 
1490 template <>
1491 double
1492 TriaAccessor<2, 2, 2>::extent_in_direction(const unsigned int axis) const
1493 {
1494  const unsigned int lines[2][2] = {
1495  {2, 3},
1496  {0, 1}};
1497 
1498  AssertIndexRange(axis, 2);
1499 
1500  return std::max(this->line(lines[axis][0])->diameter(),
1501  this->line(lines[axis][1])->diameter());
1502 }
1503 
1504 template <>
1505 double
1506 TriaAccessor<2, 2, 3>::extent_in_direction(const unsigned int axis) const
1507 {
1508  const unsigned int lines[2][2] = {
1509  {2, 3},
1510  {0, 1}};
1511 
1512  AssertIndexRange(axis, 2);
1513 
1514  return std::max(this->line(lines[axis][0])->diameter(),
1515  this->line(lines[axis][1])->diameter());
1516 }
1517 
1518 
1519 template <>
1520 double
1521 TriaAccessor<3, 3, 3>::extent_in_direction(const unsigned int axis) const
1522 {
1523  const unsigned int lines[3][4] = {
1524  {2, 3, 6, 7},
1525  {0, 1, 4, 5},
1526  {8, 9, 10, 11}};
1527 
1528  AssertIndexRange(axis, 3);
1529 
1530  double lengths[4] = {this->line(lines[axis][0])->diameter(),
1531  this->line(lines[axis][1])->diameter(),
1532  this->line(lines[axis][2])->diameter(),
1533  this->line(lines[axis][3])->diameter()};
1534 
1535  return std::max(std::max(lengths[0], lengths[1]),
1536  std::max(lengths[2], lengths[3]));
1537 }
1538 
1539 
1540 // Recursively set manifold ids on hex iterators.
1541 template <>
1542 void
1544  const types::manifold_id manifold_ind) const
1545 {
1546  set_manifold_id(manifold_ind);
1547 
1548  if (this->has_children())
1549  for (unsigned int c = 0; c < this->n_children(); ++c)
1550  this->child(c)->set_all_manifold_ids(manifold_ind);
1551 
1552  // for hexes also set manifold_id
1553  // of bounding quads and lines
1554 
1555  // Six bonding quads
1556  for (unsigned int i = 0; i < 6; ++i)
1557  this->quad(i)->set_manifold_id(manifold_ind);
1558  // Twelve bounding lines
1559  for (unsigned int i = 0; i < 12; ++i)
1560  this->line(i)->set_manifold_id(manifold_ind);
1561 }
1562 
1563 
1564 template <int structdim, int dim, int spacedim>
1567  const Point<structdim> &coordinates) const
1568 {
1569  // Surrounding points and weights.
1570  std::array<Point<spacedim>, GeometryInfo<structdim>::vertices_per_cell> p;
1571  std::array<double, GeometryInfo<structdim>::vertices_per_cell> w;
1572 
1573  for (unsigned int i = 0; i < GeometryInfo<structdim>::vertices_per_cell; ++i)
1574  {
1575  p[i] = this->vertex(i);
1576  w[i] = GeometryInfo<structdim>::d_linear_shape_function(coordinates, i);
1577  }
1578 
1579  return this->get_manifold().get_new_point(make_array_view(p.begin(), p.end()),
1580  make_array_view(w.begin(),
1581  w.end()));
1582 }
1583 
1584 
1585 namespace
1586 {
1607  template <int dim>
1608  struct TransformR2UAffine
1609  {
1610  static const double KA[GeometryInfo<dim>::vertices_per_cell][dim];
1611  static const double Kb[GeometryInfo<dim>::vertices_per_cell];
1612  };
1613 
1614 
1615  /*
1616  Octave code:
1617  M=[0 1; 1 1];
1618  K1 = transpose(M) * inverse (M*transpose(M));
1619  printf ("{%f, %f},\n", K1' );
1620  */
1621  template <>
1622  const double TransformR2UAffine<1>::KA[GeometryInfo<1>::vertices_per_cell]
1623  [1] = {{-1.000000}, {1.000000}};
1624 
1625  template <>
1626  const double TransformR2UAffine<1>::Kb[GeometryInfo<1>::vertices_per_cell] =
1627  {1.000000, 0.000000};
1628 
1629 
1630  /*
1631  Octave code:
1632  M=[0 1 0 1;0 0 1 1;1 1 1 1];
1633  K2 = transpose(M) * inverse (M*transpose(M));
1634  printf ("{%f, %f, %f},\n", K2' );
1635  */
1636  template <>
1637  const double TransformR2UAffine<2>::KA[GeometryInfo<2>::vertices_per_cell]
1638  [2] = {{-0.500000, -0.500000},
1639  {0.500000, -0.500000},
1640  {-0.500000, 0.500000},
1641  {0.500000, 0.500000}};
1642 
1643  /*
1644  Octave code:
1645  M=[0 1 0 1 0 1 0 1;0 0 1 1 0 0 1 1; 0 0 0 0 1 1 1 1; 1 1 1 1 1 1 1 1];
1646  K3 = transpose(M) * inverse (M*transpose(M))
1647  printf ("{%f, %f, %f, %f},\n", K3' );
1648  */
1649  template <>
1650  const double TransformR2UAffine<2>::Kb[GeometryInfo<2>::vertices_per_cell] =
1651  {0.750000, 0.250000, 0.250000, -0.250000};
1652 
1653 
1654  template <>
1655  const double TransformR2UAffine<3>::KA[GeometryInfo<3>::vertices_per_cell]
1656  [3] = {
1657  {-0.250000, -0.250000, -0.250000},
1658  {0.250000, -0.250000, -0.250000},
1659  {-0.250000, 0.250000, -0.250000},
1660  {0.250000, 0.250000, -0.250000},
1661  {-0.250000, -0.250000, 0.250000},
1662  {0.250000, -0.250000, 0.250000},
1663  {-0.250000, 0.250000, 0.250000},
1664  {0.250000, 0.250000, 0.250000}
1665 
1666  };
1667 
1668 
1669  template <>
1670  const double TransformR2UAffine<3>::Kb[GeometryInfo<3>::vertices_per_cell] = {
1671  0.500000,
1672  0.250000,
1673  0.250000,
1674  0.000000,
1675  0.250000,
1676  0.000000,
1677  0.000000,
1678  -0.250000};
1679 } // namespace
1680 
1681 
1682 template <int structdim, int dim, int spacedim>
1685  const Point<spacedim> &point) const
1686 {
1687  // A = vertex * KA
1689 
1690  // copy vertices to avoid expensive resolution of vertex index inside loop
1691  std::array<Point<spacedim>, GeometryInfo<structdim>::vertices_per_cell>
1692  vertices;
1693  for (unsigned int v = 0; v < GeometryInfo<structdim>::vertices_per_cell; ++v)
1694  vertices[v] = this->vertex(v);
1695  for (unsigned int d = 0; d < spacedim; ++d)
1696  for (unsigned int v = 0; v < GeometryInfo<structdim>::vertices_per_cell;
1697  ++v)
1698  for (unsigned int e = 0; e < structdim; ++e)
1699  A[d][e] += vertices[v][d] * TransformR2UAffine<structdim>::KA[v][e];
1700 
1701  // b = vertex * Kb
1702  Tensor<1, spacedim> b = point;
1703  for (unsigned int v = 0; v < GeometryInfo<structdim>::vertices_per_cell; ++v)
1704  b -= vertices[v] * TransformR2UAffine<structdim>::Kb[v];
1705 
1707  return Point<structdim>(apply_transformation(A_inv, b));
1708 }
1709 
1710 
1711 template <int structdim, int dim, int spacedim>
1714  const bool respect_manifold,
1715  const bool use_interpolation) const
1716 {
1717  if (respect_manifold == false)
1718  {
1719  Assert(use_interpolation == false, ExcNotImplemented());
1720  Point<spacedim> p;
1721  for (unsigned int v = 0; v < GeometryInfo<structdim>::vertices_per_cell;
1722  ++v)
1723  p += vertex(v);
1725  }
1726  else
1727  return get_new_point_on_object(*this, use_interpolation);
1728 }
1729 
1730 
1731 /*------------------------ Functions: CellAccessor<1> -----------------------*/
1732 
1733 
1734 
1735 template <>
1736 bool
1738 {
1739  return (this->vertex(0)[0] <= p[0]) && (p[0] <= this->vertex(1)[0]);
1740 }
1741 
1742 
1743 
1744 /*------------------------ Functions: CellAccessor<2> -----------------------*/
1745 
1746 
1747 
1748 template <>
1749 bool
1751 {
1752  // we check whether the point is
1753  // inside the cell by making sure
1754  // that it on the inner side of
1755  // each line defined by the faces,
1756  // i.e. for each of the four faces
1757  // we take the line that connects
1758  // the two vertices and subdivide
1759  // the whole domain by that in two
1760  // and check whether the point is
1761  // on the `cell-side' (rather than
1762  // the `out-side') of this line. if
1763  // the point is on the `cell-side'
1764  // for all four faces, it must be
1765  // inside the cell.
1766 
1767  // we want the faces in counter
1768  // clockwise orientation
1769  static const int direction[4] = {-1, 1, 1, -1};
1770  for (unsigned int f = 0; f < 4; ++f)
1771  {
1772  // vector from the first vertex
1773  // of the line to the point
1774  const Tensor<1, 2> to_p =
1775  p - this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 0));
1776  // vector describing the line
1777  const Tensor<1, 2> face =
1778  direction[f] *
1779  (this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 1)) -
1780  this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 0)));
1781 
1782  // if we rotate the face vector
1783  // by 90 degrees to the left
1784  // (i.e. it points to the
1785  // inside) and take the scalar
1786  // product with the vector from
1787  // the vertex to the point,
1788  // then the point is in the
1789  // `cell-side' if the scalar
1790  // product is positive. if this
1791  // is not the case, we can be
1792  // sure that the point is
1793  // outside
1794  if ((-face[1] * to_p[0] + face[0] * to_p[1]) < 0)
1795  return false;
1796  }
1797 
1798  // if we arrived here, then the
1799  // point is inside for all four
1800  // faces, and thus inside
1801  return true;
1802 }
1803 
1804 
1805 
1806 /*------------------------ Functions: CellAccessor<3> -----------------------*/
1807 
1808 
1809 
1810 template <>
1811 bool
1813 {
1814  // original implementation by Joerg
1815  // Weimar
1816 
1817  // we first eliminate points based
1818  // on the maximum and minimum of
1819  // the corner coordinates, then
1820  // transform to the unit cell, and
1821  // check there.
1822  const unsigned int dim = 3;
1823  const unsigned int spacedim = 3;
1824  Point<spacedim> maxp = this->vertex(0);
1825  Point<spacedim> minp = this->vertex(0);
1826 
1827  for (unsigned int v = 1; v < GeometryInfo<dim>::vertices_per_cell; ++v)
1828  for (unsigned int d = 0; d < dim; ++d)
1829  {
1830  maxp[d] = std::max(maxp[d], this->vertex(v)[d]);
1831  minp[d] = std::min(minp[d], this->vertex(v)[d]);
1832  }
1833 
1834  // rule out points outside the
1835  // bounding box of this cell
1836  for (unsigned int d = 0; d < dim; d++)
1837  if ((p[d] < minp[d]) || (p[d] > maxp[d]))
1838  return false;
1839 
1840  // now we need to check more carefully: transform to the
1841  // unit cube and check there. unfortunately, this isn't
1842  // completely trivial since the transform_real_to_unit_cell
1843  // function may throw an exception that indicates that the
1844  // point given could not be inverted. we take this as a sign
1845  // that the point actually lies outside, as also documented
1846  // for that function
1847  try
1848  {
1849  const TriaRawIterator<CellAccessor<dim, spacedim>> cell_iterator(*this);
1851  StaticMappingQ1<dim, spacedim>::mapping.transform_real_to_unit_cell(
1852  cell_iterator, p)));
1853  }
1855  {
1856  return false;
1857  }
1858 }
1859 
1860 
1861 
1862 /*------------------- Functions: CellAccessor<dim,spacedim> -----------------*/
1863 
1864 // For codim>0 we proceed as follows:
1865 // 1) project point onto manifold and
1866 // 2) transform to the unit cell with a Q1 mapping
1867 // 3) then check if inside unit cell
1868 template <int dim, int spacedim>
1869 template <int dim_, int spacedim_>
1870 bool
1872 {
1873  const TriaRawIterator<CellAccessor<dim_, spacedim_>> cell_iterator(*this);
1874  const Point<dim_> p_unit =
1875  StaticMappingQ1<dim_, spacedim_>::mapping.transform_real_to_unit_cell(
1876  cell_iterator, p);
1877 
1879 }
1880 
1881 
1882 
1883 template <>
1884 bool
1886 {
1887  return point_inside_codim<1, 2>(p);
1888 }
1889 
1890 
1891 template <>
1892 bool
1894 {
1895  return point_inside_codim<1, 3>(p);
1896 }
1897 
1898 
1899 template <>
1900 bool
1902 {
1903  return point_inside_codim<2, 3>(p);
1904 }
1905 
1906 
1907 
1908 template <int dim, int spacedim>
1909 bool
1911 {
1912  switch (dim)
1913  {
1914  case 1:
1915  return at_boundary(0) || at_boundary(1);
1916  case 2:
1917  return (at_boundary(0) || at_boundary(1) || at_boundary(2) ||
1918  at_boundary(3));
1919  case 3:
1920  return (at_boundary(0) || at_boundary(1) || at_boundary(2) ||
1921  at_boundary(3) || at_boundary(4) || at_boundary(5));
1922  default:
1923  Assert(false, ExcNotImplemented());
1924  return false;
1925  }
1926 }
1927 
1928 
1929 
1930 template <int dim, int spacedim>
1933 {
1935  return this->tria->levels[this->present_level]
1936  ->cells.boundary_or_material_id[this->present_index]
1937  .material_id;
1938 }
1939 
1940 
1941 
1942 template <int dim, int spacedim>
1943 void
1945  const types::material_id mat_id) const
1946 {
1949  this->tria->levels[this->present_level]
1950  ->cells.boundary_or_material_id[this->present_index]
1951  .material_id = mat_id;
1952 }
1953 
1954 
1955 
1956 template <int dim, int spacedim>
1957 void
1959  const types::material_id mat_id) const
1960 {
1961  set_material_id(mat_id);
1962 
1963  if (this->has_children())
1964  for (unsigned int c = 0; c < this->n_children(); ++c)
1965  this->child(c)->recursively_set_material_id(mat_id);
1966 }
1967 
1968 
1969 
1970 template <int dim, int spacedim>
1971 void
1973  const types::subdomain_id new_subdomain_id) const
1974 {
1976  Assert(this->is_active(),
1977  ExcMessage("set_subdomain_id() can only be called on active cells!"));
1978  this->tria->levels[this->present_level]->subdomain_ids[this->present_index] =
1979  new_subdomain_id;
1980 }
1981 
1982 
1983 
1984 template <int dim, int spacedim>
1987 {
1989  return this->tria->levels[this->present_level]
1990  ->level_subdomain_ids[this->present_index];
1991 }
1992 
1993 
1994 
1995 template <int dim, int spacedim>
1996 void
1998  const types::subdomain_id new_level_subdomain_id) const
1999 {
2001  this->tria->levels[this->present_level]
2002  ->level_subdomain_ids[this->present_index] = new_level_subdomain_id;
2003 }
2004 
2005 
2006 template <int dim, int spacedim>
2007 bool
2009 {
2011  if (dim == spacedim)
2012  return true;
2013  else
2014  return this->tria->levels[this->present_level]
2015  ->direction_flags[this->present_index];
2016 }
2017 
2018 
2019 
2020 template <int dim, int spacedim>
2021 void
2023  const bool new_direction_flag) const
2024 {
2026  if (dim < spacedim)
2027  this->tria->levels[this->present_level]
2028  ->direction_flags[this->present_index] = new_direction_flag;
2029  else
2030  Assert(new_direction_flag == true,
2031  ExcMessage("If dim==spacedim, direction flags are always true and "
2032  "can not be set to anything else."));
2033 }
2034 
2035 
2036 
2037 template <int dim, int spacedim>
2038 void
2040  const unsigned int active_cell_index)
2041 {
2042  // set the active cell index. allow setting it also for non-active (and
2043  // unused) cells to allow resetting the index after refinement
2044  this->tria->levels[this->present_level]
2045  ->active_cell_indices[this->present_index] = active_cell_index;
2046 }
2047 
2048 
2049 
2050 template <int dim, int spacedim>
2051 void
2052 CellAccessor<dim, spacedim>::set_parent(const unsigned int parent_index)
2053 {
2055  Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2056  this->tria->levels[this->present_level]->parents[this->present_index / 2] =
2057  parent_index;
2058 }
2059 
2060 
2061 
2062 template <int dim, int spacedim>
2063 int
2065 {
2066  Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2067 
2068  // the parent of two consecutive cells
2069  // is stored only once, since it is
2070  // the same
2071  return this->tria->levels[this->present_level]
2072  ->parents[this->present_index / 2];
2073 }
2074 
2075 
2076 
2077 template <int dim, int spacedim>
2078 unsigned int
2080 {
2081  Assert(this->is_active(), TriaAccessorExceptions::ExcCellNotActive());
2082  return this->tria->levels[this->present_level]
2083  ->active_cell_indices[this->present_index];
2084 }
2085 
2086 
2087 
2088 template <int dim, int spacedim>
2091 {
2093  Assert(this->present_level > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2095  this->present_level - 1,
2096  parent_index());
2097 
2098  return q;
2099 }
2100 
2101 
2102 template <int dim, int spacedim>
2103 void
2105  const types::subdomain_id new_subdomain_id) const
2106 {
2107  if (this->has_children())
2108  for (unsigned int c = 0; c < this->n_children(); ++c)
2109  this->child(c)->recursively_set_subdomain_id(new_subdomain_id);
2110  else
2111  set_subdomain_id(new_subdomain_id);
2112 }
2113 
2114 
2115 
2116 template <int dim, int spacedim>
2117 void
2119  const unsigned int i,
2120  const TriaIterator<CellAccessor<dim, spacedim>> &pointer) const
2121 {
2123 
2124  if (pointer.state() == IteratorState::valid)
2125  {
2126  this->tria->levels[this->present_level]
2127  ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
2128  .first = pointer->present_level;
2129  this->tria->levels[this->present_level]
2130  ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
2131  .second = pointer->present_index;
2132  }
2133  else
2134  {
2135  this->tria->levels[this->present_level]
2136  ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
2137  .first = -1;
2138  this->tria->levels[this->present_level]
2139  ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
2140  .second = -1;
2141  }
2142 }
2143 
2144 
2145 
2146 template <int dim, int spacedim>
2147 CellId
2149 {
2150  std::array<unsigned char, 30> id;
2151 
2152  CellAccessor<dim, spacedim> ptr = *this;
2153  const unsigned int n_child_indices = ptr.level();
2154 
2155  while (ptr.level() > 0)
2156  {
2157  const TriaIterator<CellAccessor<dim, spacedim>> parent = ptr.parent();
2158  const unsigned int n_children = parent->n_children();
2159 
2160  // determine which child we are
2161  unsigned char v = static_cast<unsigned char>(-1);
2162  for (unsigned int c = 0; c < n_children; ++c)
2163  {
2164  if (parent->child_index(c) == ptr.index())
2165  {
2166  v = c;
2167  break;
2168  }
2169  }
2170 
2171  Assert(v != static_cast<unsigned char>(-1), ExcInternalError());
2172  id[ptr.level() - 1] = v;
2173 
2174  ptr.copy_from(*parent);
2175  }
2176 
2177  Assert(ptr.level() == 0, ExcInternalError());
2178  const unsigned int coarse_index = ptr.index();
2179 
2180  return {this->tria->coarse_cell_index_to_coarse_cell_id(coarse_index),
2181  n_child_indices,
2182  id.data()};
2183 }
2184 
2185 
2186 
2187 template <int dim, int spacedim>
2188 unsigned int
2190  const unsigned int neighbor) const
2191 {
2193 
2194  // if we have a 1d mesh in 1d, we
2195  // can assume that the left
2196  // neighbor of the right neighbor is
2197  // the current cell. but that is an
2198  // invariant that isn't true if the
2199  // mesh is embedded in a higher
2200  // dimensional space, so we have to
2201  // fall back onto the generic code
2202  // below
2203  if ((dim == 1) && (spacedim == dim))
2204  return GeometryInfo<dim>::opposite_face[neighbor];
2205 
2206  const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2207  this->neighbor(neighbor);
2208 
2209  // usually, on regular patches of
2210  // the grid, this cell is just on
2211  // the opposite side of the
2212  // neighbor that the neighbor is of
2213  // this cell. for example in 2d, if
2214  // we want to know the
2215  // neighbor_of_neighbor if
2216  // neighbor==1 (the right
2217  // neighbor), then we will get 3
2218  // (the left neighbor) in most
2219  // cases. look up this relationship
2220  // in the table provided by
2221  // GeometryInfo and try it
2222  const unsigned int this_face_index = face_index(neighbor);
2223 
2224  const unsigned int neighbor_guess =
2226 
2227  if (neighbor_cell->face_index(neighbor_guess) == this_face_index)
2228  return neighbor_guess;
2229  else
2230  // if the guess was false, then
2231  // we need to loop over all
2232  // neighbors and find the number
2233  // the hard way
2234  {
2235  for (unsigned int face_no = 0;
2236  face_no < GeometryInfo<dim>::faces_per_cell;
2237  ++face_no)
2238  if (neighbor_cell->face_index(face_no) == this_face_index)
2239  return face_no;
2240 
2241  // running over all neighbors
2242  // faces we did not find the
2243  // present face. Thereby the
2244  // neighbor must be coarser
2245  // than the present
2246  // cell. Return an invalid
2247  // unsigned int in this case.
2249  }
2250 }
2251 
2252 
2253 
2254 template <int dim, int spacedim>
2255 unsigned int
2257  const unsigned int neighbor) const
2258 {
2259  const unsigned int n2 = neighbor_of_neighbor_internal(neighbor);
2262 
2263  return n2;
2264 }
2265 
2266 
2267 
2268 template <int dim, int spacedim>
2269 bool
2271  const unsigned int neighbor) const
2272 {
2273  return neighbor_of_neighbor_internal(neighbor) ==
2275 }
2276 
2277 
2278 
2279 template <int dim, int spacedim>
2280 std::pair<unsigned int, unsigned int>
2282  const unsigned int neighbor) const
2283 {
2285  // make sure that the neighbor is
2286  // on a coarser level
2287  Assert(neighbor_is_coarser(neighbor),
2289 
2290  switch (dim)
2291  {
2292  case 2:
2293  {
2294  const int this_face_index = face_index(neighbor);
2295  const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2296  this->neighbor(neighbor);
2297 
2298  // usually, on regular patches of
2299  // the grid, this cell is just on
2300  // the opposite side of the
2301  // neighbor that the neighbor is of
2302  // this cell. for example in 2d, if
2303  // we want to know the
2304  // neighbor_of_neighbor if
2305  // neighbor==1 (the right
2306  // neighbor), then we will get 0
2307  // (the left neighbor) in most
2308  // cases. look up this relationship
2309  // in the table provided by
2310  // GeometryInfo and try it
2311  const unsigned int face_no_guess =
2313 
2314  const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face_guess =
2315  neighbor_cell->face(face_no_guess);
2316 
2317  if (face_guess->has_children())
2318  for (unsigned int subface_no = 0;
2319  subface_no < face_guess->n_children();
2320  ++subface_no)
2321  if (face_guess->child_index(subface_no) == this_face_index)
2322  return std::make_pair(face_no_guess, subface_no);
2323 
2324  // if the guess was false, then
2325  // we need to loop over all faces
2326  // and subfaces and find the
2327  // number the hard way
2328  for (unsigned int face_no = 0;
2329  face_no < GeometryInfo<2>::faces_per_cell;
2330  ++face_no)
2331  {
2332  if (face_no != face_no_guess)
2333  {
2334  const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>
2335  face = neighbor_cell->face(face_no);
2336  if (face->has_children())
2337  for (unsigned int subface_no = 0;
2338  subface_no < face->n_children();
2339  ++subface_no)
2340  if (face->child_index(subface_no) == this_face_index)
2341  return std::make_pair(face_no, subface_no);
2342  }
2343  }
2344 
2345  // we should never get here,
2346  // since then we did not find
2347  // our way back...
2348  Assert(false, ExcInternalError());
2349  return std::make_pair(numbers::invalid_unsigned_int,
2351  }
2352 
2353  case 3:
2354  {
2355  const int this_face_index = face_index(neighbor);
2356  const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2357  this->neighbor(neighbor);
2358 
2359  // usually, on regular patches of the grid, this cell is just on the
2360  // opposite side of the neighbor that the neighbor is of this cell.
2361  // for example in 2d, if we want to know the neighbor_of_neighbor if
2362  // neighbor==1 (the right neighbor), then we will get 0 (the left
2363  // neighbor) in most cases. look up this relationship in the table
2364  // provided by GeometryInfo and try it
2365  const unsigned int face_no_guess =
2367 
2368  const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face_guess =
2369  neighbor_cell->face(face_no_guess);
2370 
2371  if (face_guess->has_children())
2372  for (unsigned int subface_no = 0;
2373  subface_no < face_guess->n_children();
2374  ++subface_no)
2375  {
2376  if (face_guess->child_index(subface_no) == this_face_index)
2377  // call a helper function, that translates the current
2378  // subface number to a subface number for the current
2379  // FaceRefineCase
2380  return std::make_pair(face_no_guess,
2381  translate_subface_no(face_guess,
2382  subface_no));
2383 
2384  if (face_guess->child(subface_no)->has_children())
2385  for (unsigned int subsub_no = 0;
2386  subsub_no < face_guess->child(subface_no)->n_children();
2387  ++subsub_no)
2388  if (face_guess->child(subface_no)->child_index(subsub_no) ==
2389  this_face_index)
2390  // call a helper function, that translates the current
2391  // subface number and subsubface number to a subface
2392  // number for the current FaceRefineCase
2393  return std::make_pair(face_no_guess,
2394  translate_subface_no(face_guess,
2395  subface_no,
2396  subsub_no));
2397  }
2398 
2399  // if the guess was false, then we need to loop over all faces and
2400  // subfaces and find the number the hard way
2401  for (unsigned int face_no = 0;
2402  face_no < GeometryInfo<3>::faces_per_cell;
2403  ++face_no)
2404  {
2405  if (face_no == face_no_guess)
2406  continue;
2407 
2408  const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face =
2409  neighbor_cell->face(face_no);
2410 
2411  if (!face->has_children())
2412  continue;
2413 
2414  for (unsigned int subface_no = 0; subface_no < face->n_children();
2415  ++subface_no)
2416  {
2417  if (face->child_index(subface_no) == this_face_index)
2418  // call a helper function, that translates the current
2419  // subface number to a subface number for the current
2420  // FaceRefineCase
2421  return std::make_pair(face_no,
2422  translate_subface_no(face,
2423  subface_no));
2424 
2425  if (face->child(subface_no)->has_children())
2426  for (unsigned int subsub_no = 0;
2427  subsub_no < face->child(subface_no)->n_children();
2428  ++subsub_no)
2429  if (face->child(subface_no)->child_index(subsub_no) ==
2430  this_face_index)
2431  // call a helper function, that translates the current
2432  // subface number and subsubface number to a subface
2433  // number for the current FaceRefineCase
2434  return std::make_pair(face_no,
2435  translate_subface_no(face,
2436  subface_no,
2437  subsub_no));
2438  }
2439  }
2440 
2441  // we should never get here, since then we did not find our way
2442  // back...
2443  Assert(false, ExcInternalError());
2444  return std::make_pair(numbers::invalid_unsigned_int,
2446  }
2447 
2448  default:
2449  {
2450  Assert(false, ExcImpossibleInDim(1));
2451  return std::make_pair(numbers::invalid_unsigned_int,
2453  }
2454  }
2455 }
2456 
2457 
2458 
2459 template <int dim, int spacedim>
2460 bool
2462  const unsigned int i_face) const
2463 {
2464  /*
2465  * Implementation note: In all of the functions corresponding to periodic
2466  * faces we mainly use the Triangulation::periodic_face_map to find the
2467  * information about periodically connected faces. So, we actually search in
2468  * this std::map and return the cell_face on the other side of the periodic
2469  * boundary. For this search process, we have two options:
2470  *
2471  * 1- Using the [] operator of std::map: This option results in a more
2472  * readalbe code, but requires an extra iteration in the map. Because when we
2473  * call [] on std::map, with a key which does not exist in the std::map, that
2474  * key will be created and the default value will be returned by []. This is
2475  * not desirable. So, one has to first check if the key exists in the std::map
2476  * and if it exists, then use the [] operator. The existence check is possible
2477  * using std::map::find() or std::map::count(). Using this option will result
2478  * in two iteration cycles through the map. First, existence check, then
2479  * returning the value.
2480  *
2481  * 2- Using std::map::find(): This option is less readable, but theoretically
2482  * faster, because it results in one iteration through std::map object.
2483  *
2484  * We decided to use the 2nd option.
2485  */
2487  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2488  // my_it : is the iterator to the current cell.
2489  cell_iterator my_it(*this);
2490  if (this->tria->periodic_face_map.find(
2491  std::pair<cell_iterator, unsigned int>(my_it, i_face)) !=
2492  this->tria->periodic_face_map.end())
2493  return true;
2494  return false;
2495 }
2496 
2497 
2498 
2499 template <int dim, int spacedim>
2501 CellAccessor<dim, spacedim>::periodic_neighbor(const unsigned int i_face) const
2502 {
2503  /*
2504  * To know, why we are using std::map::find() instead of [] operator, refer
2505  * to the implementation note in has_periodic_neighbor() function.
2506  *
2507  * my_it : the iterator to the current cell.
2508  * my_face_pair : the pair reported by periodic_face_map as its first pair
2509  * being the current cell_face.
2510  */
2512  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2513  cell_iterator my_it(*this);
2514 
2515  const typename std::map<std::pair<cell_iterator, unsigned int>,
2516  std::pair<std::pair<cell_iterator, unsigned int>,
2517  std::bitset<3>>>::const_iterator
2518  my_face_pair = this->tria->periodic_face_map.find(
2519  std::pair<cell_iterator, unsigned int>(my_it, i_face));
2520  // Assertion is required to check that we are actually on a periodic boundary.
2521  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2523  return my_face_pair->second.first.first;
2524 }
2525 
2526 
2527 
2528 template <int dim, int spacedim>
2531  const unsigned int i_face) const
2532 {
2533  if (!(this->face(i_face)->at_boundary()))
2534  return this->neighbor(i_face);
2535  else if (this->has_periodic_neighbor(i_face))
2536  return this->periodic_neighbor(i_face);
2537  else
2539  // we can't come here
2540  return this->neighbor(i_face);
2541 }
2542 
2543 
2544 
2545 template <int dim, int spacedim>
2548  const unsigned int i_face,
2549  const unsigned int i_subface) const
2550 {
2551  /*
2552  * To know, why we are using std::map::find() instead of [] operator, refer
2553  * to the implementation note in has_periodic_neighbor() function.
2554  *
2555  * my_it : the iterator to the current cell.
2556  * my_face_pair : the pair reported by periodic_face_map as its first pair
2557  * being the current cell_face. nb_it : the iterator to the
2558  * neighbor of current cell at i_face. face_num_of_nb : the face number of
2559  * the periodically neighboring face in the relevant element.
2560  * nb_parent_face_it: the iterator to the parent face of the periodically
2561  * neighboring face.
2562  */
2564  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2565  cell_iterator my_it(*this);
2566  const typename std::map<std::pair<cell_iterator, unsigned int>,
2567  std::pair<std::pair<cell_iterator, unsigned int>,
2568  std::bitset<3>>>::const_iterator
2569  my_face_pair = this->tria->periodic_face_map.find(
2570  std::pair<cell_iterator, unsigned int>(my_it, i_face));
2571  /*
2572  * There should be an assertion, which tells the user that this function
2573  * should not be used for a cell which is not located at a periodic boundary.
2574  */
2575  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2577  cell_iterator parent_nb_it = my_face_pair->second.first.first;
2578  unsigned int nb_face_num = my_face_pair->second.first.second;
2579  TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> nb_parent_face_it =
2580  parent_nb_it->face(nb_face_num);
2581  /*
2582  * We should check if the parent face of the neighbor has at least the same
2583  * number of children as i_subface.
2584  */
2585  AssertIndexRange(i_subface, nb_parent_face_it->n_children());
2586  unsigned int sub_neighbor_num =
2587  GeometryInfo<dim>::child_cell_on_face(parent_nb_it->refinement_case(),
2588  nb_face_num,
2589  i_subface,
2590  my_face_pair->second.second[0],
2591  my_face_pair->second.second[1],
2592  my_face_pair->second.second[2],
2593  nb_parent_face_it->refinement_case());
2594  return parent_nb_it->child(sub_neighbor_num);
2595 }
2596 
2597 
2598 
2599 template <int dim, int spacedim>
2600 std::pair<unsigned int, unsigned int>
2602  const unsigned int i_face) const
2603 {
2604  /*
2605  * To know, why we are using std::map::find() instead of [] operator, refer
2606  * to the implementation note in has_periodic_neighbor() function.
2607  *
2608  * my_it : the iterator to the current cell.
2609  * my_face_pair : the pair reported by periodic_face_map as its first pair
2610  * being the current cell_face. nb_it : the iterator to the periodic
2611  * neighbor. nb_face_pair : the pair reported by periodic_face_map as its
2612  * first pair being the periodic neighbor cell_face. p_nb_of_p_nb : the
2613  * iterator of the periodic neighbor of the periodic neighbor of the current
2614  * cell.
2615  */
2617  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2618  const int my_face_index = this->face_index(i_face);
2619  cell_iterator my_it(*this);
2620  const typename std::map<std::pair<cell_iterator, unsigned int>,
2621  std::pair<std::pair<cell_iterator, unsigned int>,
2622  std::bitset<3>>>::const_iterator
2623  my_face_pair = this->tria->periodic_face_map.find(
2624  std::pair<cell_iterator, unsigned int>(my_it, i_face));
2625  /*
2626  * There should be an assertion, which tells the user that this function
2627  * should not be used for a cell which is not located at a periodic boundary.
2628  */
2629  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2631  cell_iterator nb_it = my_face_pair->second.first.first;
2632  unsigned int face_num_of_nb = my_face_pair->second.first.second;
2633  const typename std::map<std::pair<cell_iterator, unsigned int>,
2634  std::pair<std::pair<cell_iterator, unsigned int>,
2635  std::bitset<3>>>::const_iterator
2636  nb_face_pair = this->tria->periodic_face_map.find(
2637  std::pair<cell_iterator, unsigned int>(nb_it, face_num_of_nb));
2638  /*
2639  * Since, we store periodic neighbors for every cell (either active or
2640  * artificial or inactive) the nb_face_pair should also be mapped to some
2641  * cell_face pair. We assert this here.
2642  */
2643  Assert(nb_face_pair != this->tria->periodic_face_map.end(),
2645  cell_iterator p_nb_of_p_nb = nb_face_pair->second.first.first;
2646  TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> parent_face_it =
2647  p_nb_of_p_nb->face(nb_face_pair->second.first.second);
2648  for (unsigned int i_subface = 0; i_subface < parent_face_it->n_children();
2649  ++i_subface)
2650  if (parent_face_it->child_index(i_subface) == my_face_index)
2651  return (std::pair<unsigned int, unsigned int>(face_num_of_nb, i_subface));
2652  /*
2653  * Obviously, if the execution reaches to this point, some of our assumptions
2654  * should have been false. The most important one is, the user has called this
2655  * function on a face which does not have a coarser periodic neighbor.
2656  */
2658  return std::pair<unsigned int, unsigned int>(numbers::invalid_unsigned_int,
2660 }
2661 
2662 
2663 
2664 template <int dim, int spacedim>
2665 int
2667  const unsigned int i_face) const
2668 {
2669  return periodic_neighbor(i_face)->index();
2670 }
2671 
2672 
2673 
2674 template <int dim, int spacedim>
2675 int
2677  const unsigned int i_face) const
2678 {
2679  return periodic_neighbor(i_face)->level();
2680 }
2681 
2682 
2683 
2684 template <int dim, int spacedim>
2685 unsigned int
2687  const unsigned int i_face) const
2688 {
2689  return periodic_neighbor_face_no(i_face);
2690 }
2691 
2692 
2693 
2694 template <int dim, int spacedim>
2695 unsigned int
2697  const unsigned int i_face) const
2698 {
2699  /*
2700  * To know, why we are using std::map::find() instead of [] operator, refer
2701  * to the implementation note in has_periodic_neighbor() function.
2702  *
2703  * my_it : the iterator to the current cell.
2704  * my_face_pair : the pair reported by periodic_face_map as its first pair
2705  * being the current cell_face.
2706  */
2708  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2709  cell_iterator my_it(*this);
2710  const typename std::map<std::pair<cell_iterator, unsigned int>,
2711  std::pair<std::pair<cell_iterator, unsigned int>,
2712  std::bitset<3>>>::const_iterator
2713  my_face_pair = this->tria->periodic_face_map.find(
2714  std::pair<cell_iterator, unsigned int>(my_it, i_face));
2715  /*
2716  * There should be an assertion, which tells the user that this function
2717  * should not be called for a cell which is not located at a periodic boundary
2718  * !
2719  */
2720  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2722  return my_face_pair->second.first.second;
2723 }
2724 
2725 
2726 
2727 template <int dim, int spacedim>
2728 bool
2730  const unsigned int i_face) const
2731 {
2732  /*
2733  * To know, why we are using std::map::find() instead of [] operator, refer
2734  * to the implementation note in has_periodic_neighbor() function.
2735  *
2736  * Implementation note: Let p_nb_of_p_nb be the periodic neighbor of the
2737  * periodic neighbor of the current cell. Also, let p_face_of_p_nb_of_p_nb be
2738  * the periodic face of the p_nb_of_p_nb. If p_face_of_p_nb_of_p_nb has
2739  * children , then the periodic neighbor of the current cell is coarser than
2740  * itself. Although not tested, this implementation should work for
2741  * anisotropic refinement as well.
2742  *
2743  * my_it : the iterator to the current cell.
2744  * my_face_pair : the pair reported by periodic_face_map as its first pair
2745  * being the current cell_face. nb_it : the iterator to the periodic
2746  * neighbor. nb_face_pair : the pair reported by periodic_face_map as its
2747  * first pair being the periodic neighbor cell_face.
2748  */
2750  using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2751  cell_iterator my_it(*this);
2752  const typename std::map<std::pair<cell_iterator, unsigned int>,
2753  std::pair<std::pair<cell_iterator, unsigned int>,
2754  std::bitset<3>>>::const_iterator
2755  my_face_pair = this->tria->periodic_face_map.find(
2756  std::pair<cell_iterator, unsigned int>(my_it, i_face));
2757  /*
2758  * There should be an assertion, which tells the user that this function
2759  * should not be used for a cell which is not located at a periodic boundary.
2760  */
2761  Assert(my_face_pair != this->tria->periodic_face_map.end(),
2763  cell_iterator nb_it = my_face_pair->second.first.first;
2764  unsigned int face_num_of_nb = my_face_pair->second.first.second;
2765  const typename std::map<std::pair<cell_iterator, unsigned int>,
2766  std::pair<std::pair<cell_iterator, unsigned int>,
2767  std::bitset<3>>>::const_iterator
2768  nb_face_pair = this->tria->periodic_face_map.find(
2769  std::pair<cell_iterator, unsigned int>(nb_it, face_num_of_nb));
2770  /*
2771  * Since, we store periodic neighbors for every cell (either active or
2772  * artificial or inactive) the nb_face_pair should also be mapped to some
2773  * cell_face pair. We assert this here.
2774  */
2775  Assert(nb_face_pair != this->tria->periodic_face_map.end(),
2777  const unsigned int my_level = this->level();
2778  const unsigned int neighbor_level = nb_face_pair->second.first.first->level();
2779  Assert(my_level >= neighbor_level, ExcInternalError());
2780  return my_level > neighbor_level;
2781 }
2782 
2783 
2784 
2785 template <int dim, int spacedim>
2786 bool
2787 CellAccessor<dim, spacedim>::at_boundary(const unsigned int i) const
2788 {
2791 
2792  return (neighbor_index(i) == -1);
2793 }
2794 
2795 
2796 
2797 template <int dim, int spacedim>
2798 bool
2800 {
2801  if (dim == 1)
2802  return at_boundary();
2803  else
2804  {
2805  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
2806  if (this->line(l)->at_boundary())
2807  return true;
2808 
2809  return false;
2810  }
2811 }
2812 
2813 
2814 
2815 template <int dim, int spacedim>
2818  const unsigned int face,
2819  const unsigned int subface) const
2820 {
2821  Assert(!this->has_children(),
2822  ExcMessage("The present cell must not have children!"));
2823  Assert(!this->at_boundary(face),
2824  ExcMessage("The present cell must have a valid neighbor!"));
2825  Assert(this->neighbor(face)->has_children() == true,
2826  ExcMessage("The neighbor must have children!"));
2827 
2828  switch (dim)
2829  {
2830  case 2:
2831  {
2832  const unsigned int neighbor_neighbor =
2833  this->neighbor_of_neighbor(face);
2834  const unsigned int neighbor_child_index =
2836  this->neighbor(face)->refinement_case(),
2837  neighbor_neighbor,
2838  subface);
2839 
2841  this->neighbor(face)->child(neighbor_child_index);
2842  // the neighbors child can have children,
2843  // which are not further refined along the
2844  // face under consideration. as we are
2845  // normally interested in one of this
2846  // child's child, search for the right one.
2847  while (sub_neighbor->has_children())
2848  {
2850  sub_neighbor->refinement_case(), neighbor_neighbor) ==
2852  ExcInternalError());
2853  sub_neighbor =
2854  sub_neighbor->child(GeometryInfo<dim>::child_cell_on_face(
2855  sub_neighbor->refinement_case(), neighbor_neighbor, 0));
2856  }
2857 
2858  return sub_neighbor;
2859  }
2860 
2861 
2862  case 3:
2863  {
2864  // this function returns the neighbor's
2865  // child on a given face and
2866  // subface.
2867 
2868  // we have to consider one other aspect here:
2869  // The face might be refined
2870  // anisotropically. In this case, the subface
2871  // number refers to the following, where we
2872  // look at the face from the current cell,
2873  // thus the subfaces are in standard
2874  // orientation concerning the cell
2875  //
2876  // for isotropic refinement
2877  //
2878  // *---*---*
2879  // | 2 | 3 |
2880  // *---*---*
2881  // | 0 | 1 |
2882  // *---*---*
2883  //
2884  // for 2*anisotropic refinement
2885  // (first cut_y, then cut_x)
2886  //
2887  // *---*---*
2888  // | 2 | 3 |
2889  // *---*---*
2890  // | 0 | 1 |
2891  // *---*---*
2892  //
2893  // for 2*anisotropic refinement
2894  // (first cut_x, then cut_y)
2895  //
2896  // *---*---*
2897  // | 1 | 3 |
2898  // *---*---*
2899  // | 0 | 2 |
2900  // *---*---*
2901  //
2902  // for purely anisotropic refinement:
2903  //
2904  // *---*---* *-------*
2905  // | | | | 1 |
2906  // | 0 | 1 | or *-------*
2907  // | | | | 0 |
2908  // *---*---* *-------*
2909  //
2910  // for "mixed" refinement:
2911  //
2912  // *---*---* *---*---* *---*---* *-------*
2913  // | | 2 | | 1 | | | 1 | 2 | | 2 |
2914  // | 0 *---* or *---* 2 | or *---*---* or *---*---*
2915  // | | 1 | | 0 | | | 0 | | 0 | 1 |
2916  // *---*---* *---*---* *-------* *---*---*
2917 
2919  mother_face = this->face(face);
2920  const unsigned int total_children = mother_face->number_of_children();
2921  AssertIndexRange(subface, total_children);
2923  ExcInternalError());
2924 
2925  unsigned int neighbor_neighbor;
2928  this->neighbor(face);
2929 
2930 
2931  const RefinementCase<dim - 1> mother_face_ref_case =
2932  mother_face->refinement_case();
2933  if (mother_face_ref_case ==
2934  static_cast<RefinementCase<dim - 1>>(
2935  RefinementCase<2>::cut_xy)) // total_children==4
2936  {
2937  // this case is quite easy. we are sure,
2938  // that the neighbor is not coarser.
2939 
2940  // get the neighbor's number for the given
2941  // face and the neighbor
2942  neighbor_neighbor = this->neighbor_of_neighbor(face);
2943 
2944  // now use the info provided by GeometryInfo
2945  // to extract the neighbors child number
2946  const unsigned int neighbor_child_index =
2948  neighbor->refinement_case(),
2949  neighbor_neighbor,
2950  subface,
2951  neighbor->face_orientation(neighbor_neighbor),
2952  neighbor->face_flip(neighbor_neighbor),
2953  neighbor->face_rotation(neighbor_neighbor));
2954  neighbor_child = neighbor->child(neighbor_child_index);
2955 
2956  // make sure that the neighbor child cell we
2957  // have found shares the desired subface.
2958  Assert((this->face(face)->child(subface) ==
2959  neighbor_child->face(neighbor_neighbor)),
2960  ExcInternalError());
2961  }
2962  else //-> the face is refined anisotropically
2963  {
2964  // first of all, we have to find the
2965  // neighbor at one of the anisotropic
2966  // children of the
2967  // mother_face. determine, which of
2968  // these we need.
2969  unsigned int first_child_to_find;
2970  unsigned int neighbor_child_index;
2971  if (total_children == 2)
2972  first_child_to_find = subface;
2973  else
2974  {
2975  first_child_to_find = subface / 2;
2976  if (total_children == 3 && subface == 1 &&
2977  !mother_face->child(0)->has_children())
2978  first_child_to_find = 1;
2979  }
2980  if (neighbor_is_coarser(face))
2981  {
2982  std::pair<unsigned int, unsigned int> indices =
2983  neighbor_of_coarser_neighbor(face);
2984  neighbor_neighbor = indices.first;
2985 
2986 
2987  // we have to translate our
2988  // subface_index according to the
2989  // RefineCase and subface index of
2990  // the coarser face (our face is an
2991  // anisotropic child of the coarser
2992  // face), 'a' denotes our
2993  // subface_index 0 and 'b' denotes
2994  // our subface_index 1, whereas 0...3
2995  // denote isotropic subfaces of the
2996  // coarser face
2997  //
2998  // cut_x and coarser_subface_index=0
2999  //
3000  // *---*---*
3001  // |b=2| |
3002  // | | |
3003  // |a=0| |
3004  // *---*---*
3005  //
3006  // cut_x and coarser_subface_index=1
3007  //
3008  // *---*---*
3009  // | |b=3|
3010  // | | |
3011  // | |a=1|
3012  // *---*---*
3013  //
3014  // cut_y and coarser_subface_index=0
3015  //
3016  // *-------*
3017  // | |
3018  // *-------*
3019  // |a=0 b=1|
3020  // *-------*
3021  //
3022  // cut_y and coarser_subface_index=1
3023  //
3024  // *-------*
3025  // |a=2 b=3|
3026  // *-------*
3027  // | |
3028  // *-------*
3029  unsigned int iso_subface;
3030  if (neighbor->face(neighbor_neighbor)->refinement_case() ==
3032  iso_subface = 2 * first_child_to_find + indices.second;
3033  else
3034  {
3035  Assert(
3036  neighbor->face(neighbor_neighbor)->refinement_case() ==
3038  ExcInternalError());
3039  iso_subface = first_child_to_find + 2 * indices.second;
3040  }
3041  neighbor_child_index = GeometryInfo<dim>::child_cell_on_face(
3042  neighbor->refinement_case(),
3043  neighbor_neighbor,
3044  iso_subface,
3045  neighbor->face_orientation(neighbor_neighbor),
3046  neighbor->face_flip(neighbor_neighbor),
3047  neighbor->face_rotation(neighbor_neighbor));
3048  }
3049  else // neighbor is not coarser
3050  {
3051  neighbor_neighbor = neighbor_of_neighbor(face);
3052  neighbor_child_index = GeometryInfo<dim>::child_cell_on_face(
3053  neighbor->refinement_case(),
3054  neighbor_neighbor,
3055  first_child_to_find,
3056  neighbor->face_orientation(neighbor_neighbor),
3057  neighbor->face_flip(neighbor_neighbor),
3058  neighbor->face_rotation(neighbor_neighbor),
3059  mother_face_ref_case);
3060  }
3061 
3062  neighbor_child = neighbor->child(neighbor_child_index);
3063  // it might be, that the neighbor_child
3064  // has children, which are not refined
3065  // along the given subface. go down that
3066  // list and deliver the last of those.
3067  while (neighbor_child->has_children() &&
3069  neighbor_child->refinement_case(), neighbor_neighbor) ==
3071  neighbor_child =
3072  neighbor_child->child(GeometryInfo<dim>::child_cell_on_face(
3073  neighbor_child->refinement_case(), neighbor_neighbor, 0));
3074 
3075  // if there are two total subfaces, we
3076  // are finished. if there are four we
3077  // have to get a child of our current
3078  // neighbor_child. If there are three,
3079  // we have to check which of the two
3080  // possibilities applies.
3081  if (total_children == 3)
3082  {
3083  if (mother_face->child(0)->has_children())
3084  {
3085  if (subface < 2)
3086  neighbor_child = neighbor_child->child(
3088  neighbor_child->refinement_case(),
3089  neighbor_neighbor,
3090  subface,
3091  neighbor_child->face_orientation(neighbor_neighbor),
3092  neighbor_child->face_flip(neighbor_neighbor),
3093  neighbor_child->face_rotation(neighbor_neighbor),
3094  mother_face->child(0)->refinement_case()));
3095  }
3096  else
3097  {
3098  Assert(mother_face->child(1)->has_children(),
3099  ExcInternalError());
3100  if (subface > 0)
3101  neighbor_child = neighbor_child->child(
3103  neighbor_child->refinement_case(),
3104  neighbor_neighbor,
3105  subface - 1,
3106  neighbor_child->face_orientation(neighbor_neighbor),
3107  neighbor_child->face_flip(neighbor_neighbor),
3108  neighbor_child->face_rotation(neighbor_neighbor),
3109  mother_face->child(1)->refinement_case()));
3110  }
3111  }
3112  else if (total_children == 4)
3113  {
3114  neighbor_child =
3115  neighbor_child->child(GeometryInfo<dim>::child_cell_on_face(
3116  neighbor_child->refinement_case(),
3117  neighbor_neighbor,
3118  subface % 2,
3119  neighbor_child->face_orientation(neighbor_neighbor),
3120  neighbor_child->face_flip(neighbor_neighbor),
3121  neighbor_child->face_rotation(neighbor_neighbor),
3122  mother_face->child(subface / 2)->refinement_case()));
3123  }
3124  }
3125 
3126  // it might be, that the neighbor_child has
3127  // children, which are not refined along the
3128  // given subface. go down that list and
3129  // deliver the last of those.
3130  while (neighbor_child->has_children())
3131  neighbor_child =
3132  neighbor_child->child(GeometryInfo<dim>::child_cell_on_face(
3133  neighbor_child->refinement_case(), neighbor_neighbor, 0));
3134 
3135 #ifdef DEBUG
3136  // check, whether the face neighbor_child matches the requested
3137  // subface.
3139  switch (this->subface_case(face))
3140  {
3144  requested = mother_face->child(subface);
3145  break;
3148  requested = mother_face->child(subface / 2)->child(subface % 2);
3149  break;
3150 
3153  switch (subface)
3154  {
3155  case 0:
3156  case 1:
3157  requested = mother_face->child(0)->child(subface);
3158  break;
3159  case 2:
3160  requested = mother_face->child(1);
3161  break;
3162  default:
3163  Assert(false, ExcInternalError());
3164  }
3165  break;
3168  switch (subface)
3169  {
3170  case 0:
3171  requested = mother_face->child(0);
3172  break;
3173  case 1:
3174  case 2:
3175  requested = mother_face->child(1)->child(subface - 1);
3176  break;
3177  default:
3178  Assert(false, ExcInternalError());
3179  }
3180  break;
3181  default:
3182  Assert(false, ExcInternalError());
3183  break;
3184  }
3185  Assert(requested == neighbor_child->face(neighbor_neighbor),
3186  ExcInternalError());
3187 #endif
3188 
3189  return neighbor_child;
3190  }
3191 
3192  default:
3193  // 1d or more than 3d
3194  Assert(false, ExcNotImplemented());
3196  }
3197 }
3198 
3199 
3200 
3201 // explicit instantiations
3202 #include "tria_accessor.inst"
3203 
3204 DEAL_II_NAMESPACE_CLOSE
unsigned int periodic_neighbor_of_periodic_neighbor(const unsigned int i) const
static const unsigned int invalid_unsigned_int
Definition: types.h:187
unsigned int manifold_id
Definition: types.h:137
bool point_inside_codim(const Point< spacedim_ > &p) const
std::pair< std::array< Point< MeshIteratorType::AccessorType::space_dimension >, n_default_points_per_cell< MeshIteratorType >)>, std::array< double, n_default_points_per_cell< MeshIteratorType >)> > get_default_points_and_weights(const MeshIteratorType &iterator, const bool with_interpolation=false)
TriaIterator< CellAccessor< dim, spacedim > > neighbor_or_periodic_neighbor(const unsigned int i) const
static ::ExceptionBase & ExcNeighborIsNotCoarser()
void set_active_cell_index(const unsigned int active_cell_index)
unsigned int active_cell_index() const
void set_neighbor(const unsigned int i, const TriaIterator< CellAccessor< dim, spacedim >> &pointer) const
static ::ExceptionBase & ExcNeighborIsCoarser()
unsigned int neighbor_of_neighbor_internal(const unsigned int neighbor) const
static bool is_inside_unit_cell(const Point< dim > &p)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
const Manifold< dim, spacedim > & get_manifold() const
types::material_id material_id() const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1641
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:2369
void copy_from(const TriaAccessorBase &)
double extent_in_direction(const unsigned int axis) const
unsigned int material_id
Definition: types.h:148
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
void set_all_manifold_ids(const types::manifold_id) const
#define AssertThrow(cond, exc)
Definition: exceptions.h:1523
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1443
TriaIterator< CellAccessor< dim, spacedim > > neighbor_child_on_subface(const unsigned int face_no, const unsigned int subface_no) const
bool neighbor_is_coarser(const unsigned int neighbor) const
int level() const
unsigned int vertex_index(const unsigned int i) const
static ::ExceptionBase & ExcCellNotUsed()
bool at_boundary() const
TriaIterator< CellAccessor< dim, spacedim > > parent() const
CellId id() const
static ::ExceptionBase & ExcMessage(std::string arg1)
bool has_periodic_neighbor(const unsigned int i) const
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
unsigned int subdomain_id
Definition: types.h:43
void set(const ::internal::TriangulationImplementation::TriaObject< structdim > &o) const
bool periodic_neighbor_is_coarser(const unsigned int i) const
#define Assert(cond, exc)
Definition: exceptions.h:1411
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)
Abstract base class for mapping classes.
Definition: mapping.h:302
int index() const
static ::ExceptionBase & ExcCellNotActive()
static ::ExceptionBase & ExcNoPeriodicNeighbor()
void recursively_set_material_id(const types::material_id new_material_id) const
TriaIterator< CellAccessor< dim, spacedim > > periodic_neighbor_child_on_subface(const unsigned int face_no, const unsigned int subface_no) const
Definition: cell_id.h:68
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
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:422
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
unsigned int periodic_neighbor_face_no(const unsigned int i) const
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
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:607
void set_subdomain_id(const types::subdomain_id new_subdomain_id) const
types::subdomain_id level_subdomain_id() const
double measure() const
const types::material_id invalid_material_id
Definition: types.h:217
unsigned int neighbor_of_neighbor(const unsigned int neighbor) const
bool has_boundary_lines() const
static ::ExceptionBase & ExcInternalError()