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