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