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