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