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