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