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