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