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