deal.II version GIT relicensing-6834-g5b78e6bcdf 2026-10-01 11:20:01+00:00
\(\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\}}\)
Loading...
Searching...
No Matches
tria_accessor.cc
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 1998 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13#include <deal.II/base/geometry_info.h>
14#include <deal.II/base/quadrature.h>
15
16#include <deal.II/dofs/dof_accessor.h>
17#include <deal.II/dofs/dof_handler.h>
18
19#include <deal.II/fe/fe_q.h>
20#include <deal.II/fe/mapping.h>
21
22#include <deal.II/grid/grid_tools_geometry.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_iterator.h>
27#include <deal.II/grid/tria_levels.h>
28
29#include <algorithm>
30#include <array>
31#include <cmath>
32#include <limits>
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 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 + static_cast<unsigned int>(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 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 DEAL_II_ASSERT_UNREACHABLE();
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 DEAL_II_ASSERT_UNREACHABLE();
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 DEAL_II_NOT_IMPLEMENTED();
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 (accessor.reference_cell() == ReferenceCells::Quadrilateral)
1328 {
1329 const Point<3> x0 = accessor.vertex(0);
1330 const Point<3> x1 = accessor.vertex(1);
1331 const Point<3> x2 = accessor.vertex(2);
1332 const Point<3> x3 = accessor.vertex(3);
1333
1334 // This is based on the approach used in libMesh (see face_quad4.C): the
1335 // primary differences are the vertex numbering and quadrature order.
1336 //
1337 // The area of a surface is the integral of the magnitude of its normal
1338 // vector, which may be computed via the cross product of two tangent
1339 // vectors. We can easily get tangent vectors from the surface
1340 // parameterization. Hence, given a bilinear surface
1341 //
1342 // X(chi, eta) = x0 + (x1 - x0) chi + (x2 - x0) eta
1343 // + (x3 + x0 - x1 - x2) chi eta
1344 //
1345 // the tangent vectors are
1346 //
1347 // t1 = (x1 - x0) + (x3 + x0 - x1 - x2) eta
1348 // t2 = (x2 - x0) + (x3 + x0 - x1 - x2) xi
1349 const Tensor<1, 3> b0 = x1 - x0;
1350 const Tensor<1, 3> b1 = x2 - x0;
1351 const Tensor<1, 3> a = x3 - x2 - b0;
1352
1353 // The diameter is the maximum distance between any pair of vertices and
1354 // we can use it as a length scale for the cell. If all components of a
1355 // (the vector connecting x3 and the last vertex of the parallelogram
1356 // defined by the first three vertices) are zero within some tolerance,
1357 // then we have a parallelogram and can use a much simpler formula.
1358 double a_max = 0.0;
1359 for (unsigned int d = 0; d < 3; ++d)
1360 a_max = std::max(std::abs(a[d]), a_max);
1361 if (a_max < 1e-14 * accessor.diameter())
1362 return cross_product_3d(b0, b1).norm();
1363
1364 // Otherwise, use a 4x4 quadrature to approximate the surface area.
1365 // Hard-code this in to prevent the extra overhead of always creating
1366 // the same QGauss rule.
1367 constexpr unsigned int n_qp = 4;
1368 const double c1 = 2.0 / 7.0 * std::sqrt(6.0 / 5.0);
1369 const double w0 = (18.0 - std::sqrt(30)) / 72.0;
1370 const double w1 = (18.0 + std::sqrt(30)) / 72.0;
1371
1372 const std::array<double, n_qp> q{{
1373 0.5 - std::sqrt(3.0 / 7.0 + c1) / 2.0,
1374 0.5 - std::sqrt(3.0 / 7.0 - c1) / 2.0,
1375 0.5 + std::sqrt(3.0 / 7.0 - c1) / 2.0,
1376 0.5 + std::sqrt(3.0 / 7.0 + c1) / 2.0,
1377 }};
1378 const std::array<double, n_qp> w{{w0, w1, w1, w0}};
1379
1380 double area = 0.;
1381 for (unsigned int i = 0; i < n_qp; ++i)
1382 for (unsigned int j = 0; j < n_qp; ++j)
1383 area += cross_product_3d(q[i] * a + b0, q[j] * a + b1).norm() *
1384 w[i] * w[j];
1385
1386 return area;
1387 }
1388 else if (accessor.reference_cell() == ReferenceCells::Triangle)
1389 {
1390 // We can just use the normal triangle area formula without issue
1391 const Tensor<1, 3> v01 = accessor.vertex(1) - accessor.vertex(0);
1392 const Tensor<1, 3> v02 = accessor.vertex(2) - accessor.vertex(0);
1393 return 0.5 * cross_product_3d(v01, v02).norm();
1394 }
1395
1396 DEAL_II_NOT_IMPLEMENTED();
1397 return 0.0;
1398 }
1399
1400
1401
1402 template <int structdim, int dim, int spacedim>
1403 Point<spacedim>
1404 get_new_point_on_object(const TriaAccessor<structdim, dim, spacedim> &obj,
1405 const bool use_interpolation)
1406 {
1407 if (use_interpolation)
1408 {
1409 const TriaRawIterator<TriaAccessor<structdim, dim, spacedim>> it(obj);
1410 const auto points_and_weights =
1411 Manifolds::get_default_points_and_weights(it, use_interpolation);
1412 return obj.get_manifold().get_new_point(
1413 make_array_view(points_and_weights.first.begin(),
1414 points_and_weights.first.end()),
1415 make_array_view(points_and_weights.second.begin(),
1416 points_and_weights.second.end()));
1417 }
1418 else
1419 {
1420 const TriaIterator<TriaAccessor<structdim, dim, spacedim>> it(obj);
1421 if constexpr (structdim == 1)
1422 return obj.get_manifold().get_new_point_on_line(it);
1423 else if constexpr (structdim == 2)
1424 return obj.get_manifold().get_new_point_on_quad(it);
1425 else if constexpr (structdim == 3)
1426 return obj.get_manifold().get_new_point_on_hex(it);
1427 else
1428 DEAL_II_ASSERT_UNREACHABLE();
1429
1430 return {};
1431 }
1432 }
1433} // namespace
1434
1435
1436
1437/*-------------------- Static variables: TriaAccessorBase -------------------*/
1438#ifndef DOXYGEN
1439
1440template <int structdim, int dim, int spacedim>
1441const unsigned int TriaAccessorBase<structdim, dim, spacedim>::dimension;
1442
1443template <int structdim, int dim, int spacedim>
1444const unsigned int TriaAccessorBase<structdim, dim, spacedim>::space_dimension;
1445
1446template <int structdim, int dim, int spacedim>
1447const unsigned int
1448 TriaAccessorBase<structdim, dim, spacedim>::structure_dimension;
1449
1450#endif
1451/*------------------------ Functions: TriaAccessor ---------------------------*/
1452#ifndef DOXYGEN
1453
1454template <int structdim, int dim, int spacedim>
1455void
1456TriaAccessor<structdim, dim, spacedim>::set_bounding_object_indices(
1457 const std::initializer_list<int> &new_indices) const
1458{
1459 ArrayView<int> bounding_object_index_ref =
1460 this->objects().get_bounding_object_indices(this->present_index);
1461
1462 if constexpr (running_in_debug_mode())
1463 for (const auto &v : new_indices)
1464 Assert(v >= 0, ExcInternalError());
1465
1466 AssertIndexRange(new_indices.size(), bounding_object_index_ref.size() + 1);
1467 std::copy(new_indices.begin(),
1468 new_indices.end(),
1469 bounding_object_index_ref.begin());
1470}
1471
1472
1473
1474template <int structdim, int dim, int spacedim>
1475void
1476TriaAccessor<structdim, dim, spacedim>::set_bounding_object_indices(
1477 const std::initializer_list<unsigned int> &new_indices) const
1478{
1479 const ArrayView<int> bounding_object_index_ref =
1480 this->objects().get_bounding_object_indices(this->present_index);
1481
1482 AssertIndexRange(new_indices.size(), bounding_object_index_ref.size() + 1);
1483 std::copy(new_indices.begin(),
1484 new_indices.end(),
1485 bounding_object_index_ref.begin());
1486}
1487
1488
1489
1490template <int structdim, int dim, int spacedim>
1491Point<spacedim>
1492TriaAccessor<structdim, dim, spacedim>::barycenter() const
1493{
1494 // call the function in the anonymous
1495 // namespace above
1496 return ::barycenter(*this);
1497}
1498
1499
1500
1501template <int structdim, int dim, int spacedim>
1502double
1503TriaAccessor<structdim, dim, spacedim>::measure() const
1504{
1505 // call the function in the anonymous
1506 // namespace above
1507 return ::measure(*this);
1508}
1509
1510
1511
1512template <int structdim, int dim, int spacedim>
1513BoundingBox<spacedim>
1514TriaAccessor<structdim, dim, spacedim>::bounding_box() const
1515{
1516 std::pair<Point<spacedim>, Point<spacedim>> boundary_points =
1517 std::make_pair(this->vertex(0), this->vertex(0));
1518
1519 const unsigned int n_vertices = this->n_vertices();
1520 for (unsigned int v = 1; v < n_vertices; ++v)
1521 {
1522 const Point<spacedim> x = this->vertex(v);
1523 for (unsigned int k = 0; k < spacedim; ++k)
1524 {
1525 boundary_points.first[k] = std::min(boundary_points.first[k], x[k]);
1526 boundary_points.second[k] = std::max(boundary_points.second[k], x[k]);
1527 }
1528 }
1529
1530 return BoundingBox<spacedim>(boundary_points);
1531}
1532
1533
1534
1535template <int structdim, int dim, int spacedim>
1536double
1537TriaAccessor<structdim, dim, spacedim>::extent_in_direction(
1538 const unsigned int /*axis*/) const
1539{
1540 DEAL_II_NOT_IMPLEMENTED();
1541 return std::numeric_limits<double>::signaling_NaN();
1542}
1543
1544#endif
1545
1546template <>
1547double
1548TriaAccessor<1, 1, 1>::extent_in_direction(const unsigned int axis) const
1549{
1550 AssertIndexRange(axis, 1);
1551
1552 return this->diameter();
1553}
1554
1555
1556template <>
1557double
1558TriaAccessor<1, 1, 2>::extent_in_direction(const unsigned int axis) const
1559{
1560 AssertIndexRange(axis, 1);
1561
1562 return this->diameter();
1563}
1564
1565
1566template <>
1567double
1568TriaAccessor<2, 2, 2>::extent_in_direction(const unsigned int axis) const
1569{
1570 Assert(this->reference_cell() == ReferenceCells::Quadrilateral,
1571 ExcNotImplemented());
1572
1573 constexpr unsigned int lines[2][2] = {
1574 {2, 3}, // Lines along x-axis, see GeometryInfo
1575 {0, 1}}; // Lines along y-axis
1576
1577 AssertIndexRange(axis, 2);
1578
1579 return std::max(this->line(lines[axis][0])->diameter(),
1580 this->line(lines[axis][1])->diameter());
1581}
1582
1583template <>
1584double
1585TriaAccessor<2, 2, 3>::extent_in_direction(const unsigned int axis) const
1586{
1587 Assert(this->reference_cell() == ReferenceCells::Quadrilateral,
1588 ExcNotImplemented());
1589
1590 constexpr unsigned int lines[2][2] = {
1591 {2, 3}, // Lines along x-axis, see GeometryInfo
1592 {0, 1}}; // Lines along y-axis
1593
1594 AssertIndexRange(axis, 2);
1595
1596 return std::max(this->line(lines[axis][0])->diameter(),
1597 this->line(lines[axis][1])->diameter());
1598}
1599
1600
1601template <>
1602double
1603TriaAccessor<3, 3, 3>::extent_in_direction(const unsigned int axis) const
1604{
1605 Assert(this->reference_cell() == ReferenceCells::Hexahedron,
1606 ExcNotImplemented());
1607
1608 constexpr unsigned int lines[3][4] = {
1609 {2, 3, 6, 7}, // Lines along x-axis, see GeometryInfo
1610 {0, 1, 4, 5}, // Lines along y-axis
1611 {8, 9, 10, 11}}; // Lines along z-axis
1612
1613 AssertIndexRange(axis, 3);
1614
1615 const double lengths[4] = {this->line(lines[axis][0])->diameter(),
1616 this->line(lines[axis][1])->diameter(),
1617 this->line(lines[axis][2])->diameter(),
1618 this->line(lines[axis][3])->diameter()};
1619
1620 return std::max({lengths[0], lengths[1], lengths[2], lengths[3]});
1621}
1622
1623
1624// Recursively set manifold ids on hex iterators.
1625template <>
1626void
1627TriaAccessor<3, 3, 3>::set_all_manifold_ids(
1628 const types::manifold_id manifold_ind) const
1629{
1630 set_manifold_id(manifold_ind);
1631
1632 if (this->has_children())
1633 for (unsigned int c = 0; c < this->n_children(); ++c)
1634 this->child(c)->set_all_manifold_ids(manifold_ind);
1635
1636 // for hexes also set manifold_id
1637 // of bounding quads and lines
1638
1639 for (const unsigned int i : this->face_indices())
1640 this->quad(i)->set_manifold_id(manifold_ind);
1641 for (const unsigned int i : this->line_indices())
1642 this->line(i)->set_manifold_id(manifold_ind);
1643}
1644
1645
1646template <int structdim, int dim, int spacedim>
1647Point<spacedim>
1648TriaAccessor<structdim, dim, spacedim>::intermediate_point(
1649 const Point<structdim> &coordinates) const
1650{
1651 // Surrounding points and weights.
1652 std::array<Point<spacedim>, GeometryInfo<structdim>::vertices_per_cell> p;
1653 std::array<double, GeometryInfo<structdim>::vertices_per_cell> w;
1654
1655 for (const unsigned int i : this->vertex_indices())
1656 {
1657 p[i] = this->vertex(i);
1658 w[i] = GeometryInfo<structdim>::d_linear_shape_function(coordinates, i);
1659 }
1660
1661 return this->get_manifold().get_new_point(make_array_view(p.begin(), p.end()),
1662 make_array_view(w.begin(),
1663 w.end()));
1664}
1665
1666
1667
1668template <int structdim, int dim, int spacedim>
1669Point<structdim>
1670TriaAccessor<structdim, dim, spacedim>::real_to_unit_cell_affine_approximation(
1671 const Point<spacedim> &point) const
1672{
1673 Assert(this->reference_cell().is_hyper_cube(), ExcNotImplemented());
1674
1675 std::array<Point<spacedim>, GeometryInfo<structdim>::vertices_per_cell>
1676 vertices;
1677 for (const unsigned int v : this->vertex_indices())
1678 vertices[v] = this->vertex(v);
1679
1680 const auto A_b =
1681 GridTools::affine_cell_approximation<structdim, spacedim>(vertices);
1682 DerivativeForm<1, spacedim, structdim> A_inv =
1683 A_b.first.covariant_form().transpose();
1684 return Point<structdim>(apply_transformation(A_inv, point - A_b.second));
1685}
1686
1687
1688
1689template <int structdim, int dim, int spacedim>
1690Point<spacedim>
1691TriaAccessor<structdim, dim, spacedim>::center(
1692 const bool respect_manifold,
1693 const bool use_interpolation) const
1694{
1695 if (respect_manifold == false)
1696 {
1697 Assert(use_interpolation == false, ExcNotImplemented());
1698 Point<spacedim> p;
1699 for (const unsigned int v : this->vertex_indices())
1700 p += vertex(v);
1701 return p / this->n_vertices();
1702 }
1703 else
1704 return get_new_point_on_object(*this, use_interpolation);
1705}
1706
1707
1708/*---------------- Functions: TriaAccessor<0,1,spacedim> -------------------*/
1709
1710
1711template <int spacedim>
1712bool
1719
1720
1721
1722template <int spacedim>
1723void
1729
1730
1731
1732template <int spacedim>
1733void
1739
1740
1741
1742template <int spacedim>
1743void
1744TriaAccessor<0, 1, spacedim>::recursively_set_user_flag() const
1745{
1746 set_user_flag();
1747
1748 if (this->has_children())
1749 for (unsigned int c = 0; c < this->n_children(); ++c)
1750 this->child(c)->recursively_set_user_flag();
1751}
1752
1753
1754
1755template <int spacedim>
1756void
1757TriaAccessor<0, 1, spacedim>::recursively_clear_user_flag() const
1758{
1759 clear_user_flag();
1760
1761 if (this->has_children())
1762 for (unsigned int c = 0; c < this->n_children(); ++c)
1763 this->child(c)->recursively_clear_user_flag();
1764}
1765
1766
1767
1768template <int spacedim>
1769void
1775
1776
1777
1778template <int spacedim>
1779void
1785
1786
1787
1788template <int spacedim>
1789void
1795
1796
1797
1798template <int spacedim>
1799void *
1800TriaAccessor<0, 1, spacedim>::user_pointer() const
1801{
1802 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
1803 DEAL_II_NOT_IMPLEMENTED();
1804 return nullptr;
1805}
1806
1807
1808
1809template <int spacedim>
1810void
1811TriaAccessor<0, 1, spacedim>::recursively_set_user_pointer(void *p) const
1812{
1813 set_user_pointer(p);
1814
1815 if (this->has_children())
1816 for (unsigned int c = 0; c < this->n_children(); ++c)
1817 this->child(c)->recursively_set_user_pointer(p);
1818}
1819
1820
1821
1822template <int spacedim>
1823void
1824TriaAccessor<0, 1, spacedim>::recursively_clear_user_pointer() const
1825{
1826 clear_user_pointer();
1827
1828 if (this->has_children())
1829 for (unsigned int c = 0; c < this->n_children(); ++c)
1830 this->child(c)->recursively_clear_user_pointer();
1831}
1832
1833
1834
1835template <int spacedim>
1836void
1842
1843
1844
1845template <int spacedim>
1846void
1852
1853
1854
1855template <int spacedim>
1856unsigned int
1863
1864
1865
1866template <int spacedim>
1867void
1868TriaAccessor<0, 1, spacedim>::recursively_set_user_index(unsigned int p) const
1869{
1870 set_user_index(p);
1871
1872 if (this->has_children())
1873 for (unsigned int c = 0; c < this->n_children(); ++c)
1874 this->child(c)->recursively_set_user_index(p);
1875}
1876
1877
1878
1879template <int spacedim>
1880void
1881TriaAccessor<0, 1, spacedim>::recursively_clear_user_index() const
1882{
1883 clear_user_index();
1884
1885 if (this->has_children())
1886 for (unsigned int c = 0; c < this->n_children(); ++c)
1887 this->child(c)->recursively_clear_user_index();
1888}
1889
1890
1891
1892/*------------------------ Functions: CellAccessor<1> -----------------------*/
1893
1894
1895
1896template <>
1897bool
1898CellAccessor<1>::point_inside(const Point<1> &p) const
1899{
1900 return (this->vertex(0)[0] <= p[0]) && (p[0] <= this->vertex(1)[0]);
1901}
1902
1903
1904
1905/*------------------------ Functions: CellAccessor<2> -----------------------*/
1906
1907
1908
1909template <>
1910bool
1911CellAccessor<2>::point_inside(const Point<2> &p) const
1912{
1913 Assert(this->reference_cell() == ReferenceCells::Quadrilateral,
1914 ExcNotImplemented());
1915
1916 // we check whether the point is
1917 // inside the cell by making sure
1918 // that it on the inner side of
1919 // each line defined by the faces,
1920 // i.e. for each of the four faces
1921 // we take the line that connects
1922 // the two vertices and subdivide
1923 // the whole domain by that in two
1924 // and check whether the point is
1925 // on the `cell-side' (rather than
1926 // the `out-side') of this line. if
1927 // the point is on the `cell-side'
1928 // for all four faces, it must be
1929 // inside the cell.
1930
1931 // we want the faces in counter
1932 // clockwise orientation
1933 static const int direction[4] = {-1, 1, 1, -1};
1934 for (unsigned int f = 0; f < 4; ++f)
1935 {
1936 // vector from the first vertex
1937 // of the line to the point
1938 const Tensor<1, 2> to_p =
1939 p - this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 0));
1940 // vector describing the line
1941 const Tensor<1, 2> face =
1942 direction[f] *
1943 (this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 1)) -
1944 this->vertex(GeometryInfo<2>::face_to_cell_vertices(f, 0)));
1945
1946 // if we rotate the face vector
1947 // by 90 degrees to the left
1948 // (i.e. it points to the
1949 // inside) and take the scalar
1950 // product with the vector from
1951 // the vertex to the point,
1952 // then the point is in the
1953 // `cell-side' if the scalar
1954 // product is positive. if this
1955 // is not the case, we can be
1956 // sure that the point is
1957 // outside
1958 if ((-face[1] * to_p[0] + face[0] * to_p[1]) < 0)
1959 return false;
1960 }
1961
1962 // if we arrived here, then the
1963 // point is inside for all four
1964 // faces, and thus inside
1965 return true;
1966}
1967
1968
1969
1970/*------------------------ Functions: CellAccessor<3> -----------------------*/
1971
1972
1973
1974template <>
1975bool
1976CellAccessor<3>::point_inside(const Point<3> &p) const
1977{
1978 Assert(this->reference_cell() == ReferenceCells::Hexahedron,
1979 ExcNotImplemented());
1980
1981 // original implementation by Joerg
1982 // Weimar
1983
1984 // we first eliminate points based
1985 // on the maximum and minimum of
1986 // the corner coordinates, then
1987 // transform to the unit cell, and
1988 // check there.
1989 const unsigned int dim = 3;
1990 const unsigned int spacedim = 3;
1991 Point<spacedim> maxp = this->vertex(0);
1992 Point<spacedim> minp = this->vertex(0);
1993
1994 for (unsigned int v = 1; v < this->n_vertices(); ++v)
1995 for (unsigned int d = 0; d < dim; ++d)
1996 {
1997 maxp[d] = std::max(maxp[d], this->vertex(v)[d]);
1998 minp[d] = std::min(minp[d], this->vertex(v)[d]);
1999 }
2000
2001 // rule out points outside the
2002 // bounding box of this cell
2003 for (unsigned int d = 0; d < dim; ++d)
2004 if ((p[d] < minp[d]) || (p[d] > maxp[d]))
2005 return false;
2006
2007 // now we need to check more carefully: transform to the
2008 // unit cube and check there. unfortunately, this isn't
2009 // completely trivial since the transform_real_to_unit_cell
2010 // function may throw an exception that indicates that the
2011 // point given could not be inverted. we take this as a sign
2012 // that the point actually lies outside, as also documented
2013 // for that function
2014 try
2015 {
2016 const TriaRawIterator<CellAccessor<dim, spacedim>> cell_iterator(*this);
2017 return (GeometryInfo<dim>::is_inside_unit_cell(
2018 reference_cell()
2019 .template get_default_linear_mapping<spacedim>()
2020 .transform_real_to_unit_cell(cell_iterator, p)));
2021 }
2022 catch (const Mapping<dim, spacedim>::ExcTransformationFailed &)
2023 {
2024 return false;
2025 }
2026}
2027
2028
2029
2030/*------------------- Functions: CellAccessor<dim,spacedim> -----------------*/
2031
2032// The return type is the same as DoFHandler<dim,spacedim>::active_cell_iterator
2033template <int dim, int spacedim>
2034TriaActiveIterator<DoFCellAccessor<dim, spacedim, false>>
2035CellAccessor<dim, spacedim>::as_dof_handler_iterator(
2036 const DoFHandler<dim, spacedim> &dof_handler) const
2037{
2038 Assert(is_active(),
2039 ExcMessage("The current iterator points to an inactive cell. "
2040 "You cannot convert it to an iterator to an active cell."));
2041 Assert(&this->get_triangulation() == &dof_handler.get_triangulation(),
2042 ExcMessage("The triangulation associated with the iterator does not "
2043 "match that of the DoFHandler."));
2044
2045 return typename DoFHandler<dim, spacedim>::active_cell_iterator(
2046 &dof_handler.get_triangulation(),
2047 this->level(),
2048 this->index(),
2049 &dof_handler);
2050}
2051
2052
2053
2054template <int dim, int spacedim>
2055TriaIterator<DoFCellAccessor<dim, spacedim, true>>
2056CellAccessor<dim, spacedim>::as_dof_handler_level_iterator(
2057 const DoFHandler<dim, spacedim> &dof_handler) const
2058{
2059 Assert(&this->get_triangulation() == &dof_handler.get_triangulation(),
2060 ExcMessage("The triangulation associated with the iterator does not "
2061 "match that of the DoFHandler."));
2062
2063 return typename DoFHandler<dim, spacedim>::level_cell_iterator(
2064 &dof_handler.get_triangulation(),
2065 this->level(),
2066 this->index(),
2067 &dof_handler);
2068}
2069
2070
2071
2072// For codim>0 we proceed as follows:
2073// 1) project point onto manifold and
2074// 2) transform to the unit cell with a Q1 mapping
2075// 3) then check if inside unit cell
2076template <int dim, int spacedim>
2077template <int dim_, int spacedim_>
2078bool
2079CellAccessor<dim, spacedim>::point_inside_codim(const Point<spacedim_> &p) const
2080{
2081 Assert(this->reference_cell().is_hyper_cube(), ExcNotImplemented());
2082
2083 const TriaRawIterator<CellAccessor<dim_, spacedim_>> cell_iterator(*this);
2084
2085 const Point<dim_> p_unit = this->reference_cell()
2086 .template get_default_linear_mapping<spacedim_>()
2087 .transform_real_to_unit_cell(cell_iterator, p);
2088
2089 return GeometryInfo<dim_>::is_inside_unit_cell(p_unit);
2090}
2091
2092
2093
2094template <>
2095bool
2096CellAccessor<1, 2>::point_inside(const Point<2> &p) const
2097{
2098 return point_inside_codim<1, 2>(p);
2099}
2100
2101
2102template <>
2103bool
2104CellAccessor<1, 3>::point_inside(const Point<3> &p) const
2105{
2106 return point_inside_codim<1, 3>(p);
2107}
2108
2109
2110template <>
2111bool
2112CellAccessor<2, 3>::point_inside(const Point<3> &p) const
2113{
2114 Assert(this->reference_cell() == ReferenceCells::Quadrilateral,
2115 ExcNotImplemented());
2116 return point_inside_codim<2, 3>(p);
2117}
2118
2119
2120
2121template <int dim, int spacedim>
2122bool
2123CellAccessor<dim, spacedim>::at_boundary() const
2124{
2125 for (const auto face : this->face_indices())
2126 if (at_boundary(face))
2127 return true;
2128
2129 return false;
2130}
2131
2132
2133
2134template <int dim, int spacedim>
2135void
2136CellAccessor<dim, spacedim>::set_material_id(
2137 const types::material_id mat_id) const
2138{
2139 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
2140 AssertIndexRange(mat_id, numbers::invalid_material_id);
2141 this->tria->levels[this->level()]
2142 ->cells.boundary_or_material_id[this->present_index]
2143 .material_id = mat_id;
2144}
2145
2146
2147
2148template <int dim, int spacedim>
2149void
2150CellAccessor<dim, spacedim>::recursively_set_material_id(
2151 const types::material_id mat_id) const
2152{
2153 set_material_id(mat_id);
2154
2155 if (this->has_children())
2156 for (unsigned int c = 0; c < this->n_children(); ++c)
2157 this->child(c)->recursively_set_material_id(mat_id);
2158}
2159
2160
2161
2162template <int dim, int spacedim>
2163void
2164CellAccessor<dim, spacedim>::set_subdomain_id(
2165 const types::subdomain_id new_subdomain_id) const
2166{
2167 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
2168 Assert(this->is_active(),
2169 ExcMessage("set_subdomain_id() can only be called on active cells!"));
2170 this->tria->levels[this->level()]->subdomain_ids[this->present_index] =
2171 new_subdomain_id;
2172}
2173
2174
2175
2176template <int dim, int spacedim>
2177void
2178CellAccessor<dim, spacedim>::set_level_subdomain_id(
2179 const types::subdomain_id new_level_subdomain_id) const
2180{
2181 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
2182 this->tria->levels[this->level()]->level_subdomain_ids[this->present_index] =
2183 new_level_subdomain_id;
2184}
2185
2186
2187
2188template <int dim, int spacedim>
2189void
2190CellAccessor<dim, spacedim>::set_direction_flag(
2191 const bool new_direction_flag) const
2192{
2193 // Some older compilers (GCC 9) print an unused variable warning about
2194 // new_direction_flag when it is only used in a subset of 'if constexpr'
2195 // statements
2196 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
2197 if constexpr (dim == spacedim)
2198 Assert(new_direction_flag == true,
2199 ExcMessage("If dim==spacedim, direction flags are always true and "
2200 "can not be set to anything else."));
2201 else if constexpr (dim == spacedim - 1)
2202 this->tria->levels[this->level()]->direction_flags[this->present_index] =
2203 new_direction_flag;
2204 else
2205 Assert(new_direction_flag == true,
2206 ExcMessage("If dim<spacedim-1, then this function can be called "
2207 "only if the argument is 'true'."));
2208}
2209
2210
2211
2212template <int dim, int spacedim>
2213void
2214CellAccessor<dim, spacedim>::set_parent(const unsigned int parent_index)
2215{
2216 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
2217 Assert(this->level() > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2218
2219 // We only store the parent for every second cell. That's because cells are
2220 // created during refinement in multiples of two, and so two successive
2221 // cells always share the same parent.
2222 this->tria->levels[this->level()]->parents[this->present_index / 2] =
2223 parent_index;
2224}
2225
2226
2227
2228template <int dim, int spacedim>
2229int
2230CellAccessor<dim, spacedim>::parent_index() const
2231{
2232 Assert(this->level() > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2233
2234 // We only store the parent for every second cell. That's because cells are
2235 // created during refinement in multiples of two, and so two successive
2236 // cells always share the same parent.
2237 return this->tria->levels[this->level()]->parents[this->present_index / 2];
2238}
2239
2240
2241
2242template <int dim, int spacedim>
2243void
2244CellAccessor<dim, spacedim>::set_active_cell_index(
2245 const unsigned int active_cell_index) const
2246{
2247 this->tria->levels[this->level()]->active_cell_indices[this->present_index] =
2248 active_cell_index;
2249}
2250
2251
2252
2253template <int dim, int spacedim>
2254void
2255CellAccessor<dim, spacedim>::set_global_active_cell_index(
2256 const types::global_cell_index index) const
2257{
2258 this->tria->levels[this->level()]
2259 ->global_active_cell_indices[this->present_index] = index;
2260}
2261
2262
2263
2264template <int dim, int spacedim>
2265void
2266CellAccessor<dim, spacedim>::set_global_level_cell_index(
2267 const types::global_cell_index index) const
2268{
2269 this->tria->levels[this->level()]
2270 ->global_level_cell_indices[this->present_index] = index;
2271}
2272
2273
2274
2275template <int dim, int spacedim>
2276TriaIterator<CellAccessor<dim, spacedim>>
2277CellAccessor<dim, spacedim>::parent() const
2278{
2279 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
2280 Assert(this->level() > 0, TriaAccessorExceptions::ExcCellHasNoParent());
2281 TriaIterator<CellAccessor<dim, spacedim>> q(this->tria,
2282 this->level() - 1,
2283 parent_index());
2284
2285 return q;
2286}
2287
2288
2289template <int dim, int spacedim>
2290void
2291CellAccessor<dim, spacedim>::recursively_set_subdomain_id(
2292 const types::subdomain_id new_subdomain_id) const
2293{
2294 if (this->has_children())
2295 for (unsigned int c = 0; c < this->n_children(); ++c)
2296 this->child(c)->recursively_set_subdomain_id(new_subdomain_id);
2297 else
2298 set_subdomain_id(new_subdomain_id);
2299}
2300
2301
2302template <int dim, int spacedim>
2303std::set<TriaActiveIterator<CellAccessor<dim, spacedim>>>
2304CellAccessor<dim, spacedim>::get_cells_adjacent_to_line(
2305 const unsigned int i) const
2306{
2307 AssertIndexRange(i, this->n_lines());
2308
2309 // We explicitly allow the call of the get_cells_adjacent_to_line() before
2310 // compute_line_to_adjacent_cells_map() was called, as this happens during
2311 // the initialization of the FE element (e.g., in FE_NedelecSZ).
2312 if (!this->tria->line_to_adjacent_cells_map.has_value())
2313 return std::set<TriaActiveIterator<CellAccessor<dim, spacedim>>>();
2314
2315 const auto &map = this->tria->line_to_adjacent_cells_map;
2316 return map.value()[this->active_cell_index()][i];
2317}
2318
2319
2320template <int dim, int spacedim>
2321CellId
2322CellAccessor<dim, spacedim>::id() const
2323{
2324 std::array<unsigned char, 30> id;
2325
2326 CellAccessor<dim, spacedim> ptr = *this;
2327 const unsigned int n_child_indices = ptr.level();
2328
2329 while (ptr.level() > 0)
2330 {
2331 const TriaIterator<CellAccessor<dim, spacedim>> parent = ptr.parent();
2332 const unsigned int n_children = parent->n_children();
2333
2334 // determine which child we are
2335 unsigned char v = static_cast<unsigned char>(-1);
2336 for (unsigned int c = 0; c < n_children; ++c)
2337 {
2338 if (parent->child_index(c) == ptr.index())
2339 {
2340 v = c;
2341 break;
2342 }
2343 }
2344
2345 Assert(v != static_cast<unsigned char>(-1), ExcInternalError());
2346 id[ptr.level() - 1] = v;
2347
2348 ptr.copy_from(*parent);
2349 }
2350
2351 Assert(ptr.level() == 0, ExcInternalError());
2352 const unsigned int coarse_index = ptr.index();
2353
2354 return {this->tria->coarse_cell_index_to_coarse_cell_id(coarse_index),
2355 n_child_indices,
2356 id.data()};
2357}
2358
2359
2360
2361template <int dim, int spacedim>
2362unsigned int
2363CellAccessor<dim, spacedim>::neighbor_of_neighbor_internal(
2364 const unsigned int neighbor) const
2365{
2366 AssertIndexRange(neighbor, this->n_faces());
2367
2368 // if we have a 1d mesh in 1d, we
2369 // can assume that the left
2370 // neighbor of the right neighbor is
2371 // the current cell. but that is an
2372 // invariant that isn't true if the
2373 // mesh is embedded in a higher
2374 // dimensional space, so we have to
2375 // fall back onto the generic code
2376 // below
2377 if ((dim == 1) && (spacedim == dim))
2378 return GeometryInfo<dim>::opposite_face[neighbor];
2379
2380 const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2381 this->neighbor(neighbor);
2382
2383 // usually, on regular patches of
2384 // the grid, this cell is just on
2385 // the opposite side of the
2386 // neighbor that the neighbor is of
2387 // this cell. on structured grids
2388 // for example in 2d, if
2389 // we want to know the
2390 // neighbor_of_neighbor if
2391 // neighbor==1 (the right
2392 // neighbor), then we will get 3
2393 // (the left neighbor) in most
2394 // cases. look up this relationship
2395 // in the table provided by
2396 // reference cell and try it.
2397 // on mixed meshes we cannot get any
2398 // meaningful information about the opposite
2399 // face, so guess simply zero
2400 const unsigned int this_face_index = face_index(neighbor);
2401
2402 const unsigned int neighbor_guess =
2403 this->reference_cell() == neighbor_cell->reference_cell() ?
2404 this->reference_cell().opposite_face_index(neighbor) :
2405 0;
2406
2407 if (neighbor_cell->face_index(neighbor_guess) == this_face_index)
2408 return neighbor_guess;
2409
2410 // if the guess was false, then
2411 // we need to loop over all
2412 // neighbors and find the number
2413 // the hard way
2414 for (const unsigned int face_no : neighbor_cell->face_indices())
2415 if (face_no != neighbor_guess)
2416 if (neighbor_cell->face_index(face_no) == this_face_index)
2417 return face_no;
2418
2419 // running over all neighbors
2420 // faces we did not find the
2421 // present face. Thereby the
2422 // neighbor must be coarser
2423 // than the present
2424 // cell. Return an invalid
2425 // unsigned int in this case.
2426 return numbers::invalid_unsigned_int;
2427}
2428
2429
2430
2431template <int dim, int spacedim>
2432unsigned int
2433CellAccessor<dim, spacedim>::neighbor_of_neighbor(
2434 const unsigned int face_no) const
2435{
2436 const unsigned int n2 = neighbor_of_neighbor_internal(face_no);
2437 Assert(n2 != numbers::invalid_unsigned_int,
2438 TriaAccessorExceptions::ExcNeighborIsCoarser());
2439
2440 return n2;
2441}
2442
2443
2444
2445template <int dim, int spacedim>
2446bool
2447CellAccessor<dim, spacedim>::neighbor_is_coarser(
2448 const unsigned int face_no) const
2449{
2450 return neighbor_of_neighbor_internal(face_no) ==
2451 numbers::invalid_unsigned_int;
2452}
2453
2454
2455
2456template <int dim, int spacedim>
2457std::pair<unsigned int, unsigned int>
2458CellAccessor<dim, spacedim>::neighbor_of_coarser_neighbor(
2459 const unsigned int neighbor) const
2460{
2461 AssertIndexRange(neighbor, this->n_faces());
2462 // make sure that the neighbor is
2463 // on a coarser level
2464 Assert(neighbor_is_coarser(neighbor),
2465 TriaAccessorExceptions::ExcNeighborIsNotCoarser());
2466
2467 switch (dim)
2468 {
2469 case 2:
2470 {
2471 const int this_face_index = face_index(neighbor);
2472 const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2473 this->neighbor(neighbor);
2474
2475 // usually, on regular patches of
2476 // the grid, this cell is just on
2477 // the opposite side of the
2478 // neighbor that the neighbor is of
2479 // this cell. for example in 2d, if
2480 // we want to know the
2481 // neighbor_of_neighbor if
2482 // neighbor==1 (the right
2483 // neighbor), then we will get 0
2484 // (the left neighbor) in most
2485 // cases. look up this relationship
2486 // in the table provided by
2487 // the reference cell and try it
2488 // on mixed meshes we cannot get any
2489 // meaningful information about the
2490 // opposite face, so guess simply zero
2491 const unsigned int face_no_guess =
2492 this->reference_cell() == neighbor_cell->reference_cell() ?
2493 this->reference_cell().opposite_face_index(neighbor) :
2494 0;
2495
2496 const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face_guess =
2497 neighbor_cell->face(face_no_guess);
2498
2499 if (face_guess->has_children())
2500 for (unsigned int subface_no = 0;
2501 subface_no < face_guess->n_children();
2502 ++subface_no)
2503 if (face_guess->child_index(subface_no) == this_face_index)
2504 return std::make_pair(face_no_guess, subface_no);
2505
2506 // if the guess was false, then
2507 // we need to loop over all faces
2508 // and subfaces and find the
2509 // number the hard way
2510 for (const unsigned int face_no : neighbor_cell->face_indices())
2511 {
2512 if (face_no != face_no_guess)
2513 {
2514 const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>
2515 face = neighbor_cell->face(face_no);
2516 if (face->has_children())
2517 for (unsigned int subface_no = 0;
2518 subface_no < face->n_children();
2519 ++subface_no)
2520 if (face->child_index(subface_no) == this_face_index)
2521 return std::make_pair(face_no, subface_no);
2522 }
2523 }
2524
2525 // we should never get here,
2526 // since then we did not find
2527 // our way back...
2528 DEAL_II_ASSERT_UNREACHABLE();
2529 return std::make_pair(numbers::invalid_unsigned_int,
2530 numbers::invalid_unsigned_int);
2531 }
2532
2533 case 3:
2534 {
2535 const int this_face_index = face_index(neighbor);
2536 const TriaIterator<CellAccessor<dim, spacedim>> neighbor_cell =
2537 this->neighbor(neighbor);
2538
2539 // usually, on regular patches of the grid, this cell is just on the
2540 // opposite side of the neighbor that the neighbor is of this cell.
2541 // for example in 2d, if we want to know the neighbor_of_neighbor if
2542 // neighbor==1 (the right neighbor), then we will get 0 (the left
2543 // neighbor) in most cases. look up this relationship in the table
2544 // provided by the reference cell and try it
2545 // on mixed meshes we cannot get any meaningful information about the
2546 // opposite face, so guess simply zero
2547 const unsigned int face_no_guess =
2548 this->reference_cell() == neighbor_cell->reference_cell() ?
2549 this->reference_cell().opposite_face_index(neighbor) :
2550 0;
2551
2552 const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face_guess =
2553 neighbor_cell->face(face_no_guess);
2554
2555 if (face_guess->has_children())
2556 for (unsigned int subface_no = 0;
2557 subface_no < face_guess->n_children();
2558 ++subface_no)
2559 {
2560 if (face_guess->child_index(subface_no) == this_face_index)
2561 // call a helper function, that translates the current
2562 // subface number to a subface number for the current
2563 // FaceRefineCase
2564 return std::make_pair(face_no_guess,
2565 translate_subface_no(face_guess,
2566 subface_no));
2567
2568 if (face_guess->child(subface_no)->has_children())
2569 for (unsigned int subsub_no = 0;
2570 subsub_no < face_guess->child(subface_no)->n_children();
2571 ++subsub_no)
2572 if (face_guess->child(subface_no)->child_index(subsub_no) ==
2573 this_face_index)
2574 // call a helper function, that translates the current
2575 // subface number and subsubface number to a subface
2576 // number for the current FaceRefineCase
2577 return std::make_pair(face_no_guess,
2578 translate_subface_no(face_guess,
2579 subface_no,
2580 subsub_no));
2581 }
2582
2583 // if the guess was false, then we need to loop over all faces and
2584 // subfaces and find the number the hard way
2585 for (const unsigned int face_no : neighbor_cell->face_indices())
2586 {
2587 if (face_no == face_no_guess)
2588 continue;
2589
2590 const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> face =
2591 neighbor_cell->face(face_no);
2592
2593 if (!face->has_children())
2594 continue;
2595
2596 for (unsigned int subface_no = 0; subface_no < face->n_children();
2597 ++subface_no)
2598 {
2599 if (face->child_index(subface_no) == this_face_index)
2600 // call a helper function, that translates the current
2601 // subface number to a subface number for the current
2602 // FaceRefineCase
2603 return std::make_pair(face_no,
2604 translate_subface_no(face,
2605 subface_no));
2606
2607 if (face->child(subface_no)->has_children())
2608 for (unsigned int subsub_no = 0;
2609 subsub_no < face->child(subface_no)->n_children();
2610 ++subsub_no)
2611 if (face->child(subface_no)->child_index(subsub_no) ==
2612 this_face_index)
2613 // call a helper function, that translates the current
2614 // subface number and subsubface number to a subface
2615 // number for the current FaceRefineCase
2616 return std::make_pair(face_no,
2617 translate_subface_no(face,
2618 subface_no,
2619 subsub_no));
2620 }
2621 }
2622
2623 // we should never get here, since then we did not find our way
2624 // back...
2625 DEAL_II_ASSERT_UNREACHABLE();
2626 return std::make_pair(numbers::invalid_unsigned_int,
2627 numbers::invalid_unsigned_int);
2628 }
2629
2630 default:
2631 {
2632 Assert(false, ExcImpossibleInDim(1));
2633 return std::make_pair(numbers::invalid_unsigned_int,
2634 numbers::invalid_unsigned_int);
2635 }
2636 }
2637}
2638
2639
2640
2641template <int dim, int spacedim>
2642bool
2643CellAccessor<dim, spacedim>::has_periodic_neighbor(
2644 const unsigned int i_face) const
2645{
2646 /*
2647 * Implementation note: In all of the functions corresponding to periodic
2648 * faces we mainly use the Triangulation::periodic_face_map to find the
2649 * information about periodically connected faces. So, we actually search in
2650 * this std::map and return the cell_face on the other side of the periodic
2651 * boundary.
2652 *
2653 * We can not use operator[] as this would insert non-existing entries or
2654 * would require guarding with an extra std::map::find() or count().
2655 */
2656 AssertIndexRange(i_face, this->n_faces());
2657 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2658
2659 if (at_boundary(i_face) && this->tria->periodic_face_map.find(
2660 std::make_pair(cell_iterator(*this), i_face)) !=
2661 this->tria->periodic_face_map.end())
2662 return true;
2663 return false;
2664}
2665
2666
2667
2668template <int dim, int spacedim>
2669TriaIterator<CellAccessor<dim, spacedim>>
2670CellAccessor<dim, spacedim>::periodic_neighbor(const unsigned int i_face) const
2671{
2672 /*
2673 * To know, why we are using std::map::find() instead of [] operator, refer
2674 * to the implementation note in has_periodic_neighbor() function.
2675 *
2676 * my_it : the iterator to the current cell.
2677 * my_face_pair : the pair reported by periodic_face_map as its first pair
2678 * being the current cell_face.
2679 */
2680 AssertIndexRange(i_face, this->n_faces());
2681 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2682 cell_iterator current_cell(*this);
2683
2684 auto my_face_pair =
2685 this->tria->periodic_face_map.find(std::make_pair(current_cell, i_face));
2686
2687 // Make sure we are actually on a periodic boundary:
2688 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2689 TriaAccessorExceptions::ExcNoPeriodicNeighbor());
2690 return my_face_pair->second.first.first;
2691}
2692
2693
2694
2695template <int dim, int spacedim>
2696TriaIterator<CellAccessor<dim, spacedim>>
2697CellAccessor<dim, spacedim>::neighbor_or_periodic_neighbor(
2698 const unsigned int i_face) const
2699{
2700 if (!(this->face(i_face)->at_boundary()))
2701 return this->neighbor(i_face);
2702 else if (this->has_periodic_neighbor(i_face))
2703 return this->periodic_neighbor(i_face);
2704 else
2705 AssertThrow(false, TriaAccessorExceptions::ExcNoPeriodicNeighbor());
2706 // we can't come here
2707 return this->neighbor(i_face);
2708}
2709
2710
2711
2712template <int dim, int spacedim>
2713TriaIterator<CellAccessor<dim, spacedim>>
2714CellAccessor<dim, spacedim>::periodic_neighbor_child_on_subface(
2715 const unsigned int i_face,
2716 const unsigned int i_subface) const
2717{
2718 /*
2719 * To know, why we are using std::map::find() instead of [] operator, refer
2720 * to the implementation note in has_periodic_neighbor() function.
2721 *
2722 * my_it : the iterator to the current cell.
2723 * my_face_pair : the pair reported by periodic_face_map as its first pair
2724 * being the current cell_face. nb_it : the iterator to the
2725 * neighbor of current cell at i_face. face_num_of_nb : the face number of
2726 * the periodically neighboring face in the relevant element.
2727 * nb_parent_face_it: the iterator to the parent face of the periodically
2728 * neighboring face.
2729 */
2730 AssertIndexRange(i_face, this->n_faces());
2731 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2732 cell_iterator my_it(*this);
2733
2734 auto my_face_pair =
2735 this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2736 /*
2737 * There should be an assertion, which tells the user that this function
2738 * should not be used for a cell which is not located at a periodic boundary.
2739 */
2740 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2741 TriaAccessorExceptions::ExcNoPeriodicNeighbor());
2742 cell_iterator parent_nb_it = my_face_pair->second.first.first;
2743 unsigned int nb_face_num = my_face_pair->second.first.second;
2744 TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> nb_parent_face_it =
2745 parent_nb_it->face(nb_face_num);
2746 /*
2747 * We should check if the parent face of the neighbor has at least the same
2748 * number of children as i_subface.
2749 */
2750 AssertIndexRange(i_subface, nb_parent_face_it->n_children());
2751
2752 const auto [orientation, rotation, flip] =
2753 internal::split_face_orientation(my_face_pair->second.second);
2754
2755 unsigned int sub_neighbor_num =
2756 GeometryInfo<dim>::child_cell_on_face(parent_nb_it->refinement_case(),
2757 nb_face_num,
2758 i_subface,
2759 orientation,
2760 flip,
2761 rotation,
2762 nb_parent_face_it->refinement_case());
2763 return parent_nb_it->child(sub_neighbor_num);
2764}
2765
2766
2767
2768template <int dim, int spacedim>
2769std::pair<unsigned int, unsigned int>
2770CellAccessor<dim, spacedim>::periodic_neighbor_of_coarser_periodic_neighbor(
2771 const unsigned int i_face) const
2772{
2773 /*
2774 * To know, why we are using std::map::find() instead of [] operator, refer
2775 * to the implementation note in has_periodic_neighbor() function.
2776 *
2777 * my_it : the iterator to the current cell.
2778 * my_face_pair : the pair reported by periodic_face_map as its first pair
2779 * being the current cell_face. nb_it : the iterator to the periodic
2780 * neighbor. nb_face_pair : the pair reported by periodic_face_map as its
2781 * first pair being the periodic neighbor cell_face. p_nb_of_p_nb : the
2782 * iterator of the periodic neighbor of the periodic neighbor of the current
2783 * cell.
2784 */
2785 AssertIndexRange(i_face, this->n_faces());
2786 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2787 const int my_face_index = this->face_index(i_face);
2788 cell_iterator my_it(*this);
2789
2790 auto my_face_pair =
2791 this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2792 /*
2793 * There should be an assertion, which tells the user that this function
2794 * should not be used for a cell which is not located at a periodic boundary.
2795 */
2796 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2797 TriaAccessorExceptions::ExcNoPeriodicNeighbor());
2798 cell_iterator nb_it = my_face_pair->second.first.first;
2799 unsigned int face_num_of_nb = my_face_pair->second.first.second;
2800
2801 auto nb_face_pair =
2802 this->tria->periodic_face_map.find(std::make_pair(nb_it, face_num_of_nb));
2803 /*
2804 * Since, we store periodic neighbors for every cell (either active or
2805 * artificial or inactive) the nb_face_pair should also be mapped to some
2806 * cell_face pair. We assert this here.
2807 */
2808 Assert(nb_face_pair != this->tria->periodic_face_map.end(),
2809 TriaAccessorExceptions::ExcNoPeriodicNeighbor());
2810 cell_iterator p_nb_of_p_nb = nb_face_pair->second.first.first;
2811 TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> parent_face_it =
2812 p_nb_of_p_nb->face(nb_face_pair->second.first.second);
2813 for (unsigned int i_subface = 0; i_subface < parent_face_it->n_children();
2814 ++i_subface)
2815 if (parent_face_it->child_index(i_subface) == my_face_index)
2816 return std::make_pair(face_num_of_nb, i_subface);
2817 /*
2818 * Obviously, if the execution reaches to this point, some of our assumptions
2819 * should have been false. The most important one is, the user has called this
2820 * function on a face which does not have a coarser periodic neighbor.
2821 */
2822 Assert(false, TriaAccessorExceptions::ExcNeighborIsNotCoarser());
2823 return std::make_pair(numbers::invalid_unsigned_int,
2824 numbers::invalid_unsigned_int);
2825}
2826
2827
2828
2829template <int dim, int spacedim>
2830int
2831CellAccessor<dim, spacedim>::periodic_neighbor_index(
2832 const unsigned int i_face) const
2833{
2834 return periodic_neighbor(i_face)->index();
2835}
2836
2837
2838
2839template <int dim, int spacedim>
2840int
2841CellAccessor<dim, spacedim>::periodic_neighbor_level(
2842 const unsigned int i_face) const
2843{
2844 return periodic_neighbor(i_face)->level();
2845}
2846
2847
2848
2849template <int dim, int spacedim>
2850unsigned int
2851CellAccessor<dim, spacedim>::periodic_neighbor_of_periodic_neighbor(
2852 const unsigned int i_face) const
2853{
2854 return periodic_neighbor_face_no(i_face);
2855}
2856
2857
2858
2859template <int dim, int spacedim>
2860unsigned int
2861CellAccessor<dim, spacedim>::periodic_neighbor_face_no(
2862 const unsigned int i_face) const
2863{
2864 /*
2865 * To know, why we are using std::map::find() instead of [] operator, refer
2866 * to the implementation note in has_periodic_neighbor() function.
2867 *
2868 * my_it : the iterator to the current cell.
2869 * my_face_pair : the pair reported by periodic_face_map as its first pair
2870 * being the current cell_face.
2871 */
2872 AssertIndexRange(i_face, this->n_faces());
2873 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2874 cell_iterator my_it(*this);
2875
2876 auto my_face_pair =
2877 this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2878 /*
2879 * There should be an assertion, which tells the user that this function
2880 * should not be called for a cell which is not located at a periodic boundary
2881 * !
2882 */
2883 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2884 TriaAccessorExceptions::ExcNoPeriodicNeighbor());
2885 return my_face_pair->second.first.second;
2886}
2887
2888
2889
2890template <int dim, int spacedim>
2891bool
2892CellAccessor<dim, spacedim>::periodic_neighbor_is_coarser(
2893 const unsigned int i_face) const
2894{
2895 /*
2896 * To know, why we are using std::map::find() instead of [] operator, refer
2897 * to the implementation note in has_periodic_neighbor() function.
2898 *
2899 * Implementation note: Let p_nb_of_p_nb be the periodic neighbor of the
2900 * periodic neighbor of the current cell. Also, let p_face_of_p_nb_of_p_nb be
2901 * the periodic face of the p_nb_of_p_nb. If p_face_of_p_nb_of_p_nb has
2902 * children , then the periodic neighbor of the current cell is coarser than
2903 * itself. Although not tested, this implementation should work for
2904 * anisotropic refinement as well.
2905 *
2906 * my_it : the iterator to the current cell.
2907 * my_face_pair : the pair reported by periodic_face_map as its first pair
2908 * being the current cell_face. nb_it : the iterator to the periodic
2909 * neighbor. nb_face_pair : the pair reported by periodic_face_map as its
2910 * first pair being the periodic neighbor cell_face.
2911 */
2912 AssertIndexRange(i_face, this->n_faces());
2913 using cell_iterator = TriaIterator<CellAccessor<dim, spacedim>>;
2914 cell_iterator my_it(*this);
2915
2916 auto my_face_pair =
2917 this->tria->periodic_face_map.find(std::make_pair(my_it, i_face));
2918 /*
2919 * There should be an assertion, which tells the user that this function
2920 * should not be used for a cell which is not located at a periodic boundary.
2921 */
2922 Assert(my_face_pair != this->tria->periodic_face_map.end(),
2923 TriaAccessorExceptions::ExcNoPeriodicNeighbor());
2924
2925 cell_iterator nb_it = my_face_pair->second.first.first;
2926 unsigned int face_num_of_nb = my_face_pair->second.first.second;
2927
2928 auto nb_face_pair =
2929 this->tria->periodic_face_map.find(std::make_pair(nb_it, face_num_of_nb));
2930 /*
2931 * Since, we store periodic neighbors for every cell (either active or
2932 * artificial or inactive) the nb_face_pair should also be mapped to some
2933 * cell_face pair. We assert this here.
2934 */
2935 Assert(nb_face_pair != this->tria->periodic_face_map.end(),
2936 TriaAccessorExceptions::ExcNoPeriodicNeighbor());
2937 const unsigned int my_level = this->level();
2938 const unsigned int neighbor_level = nb_face_pair->second.first.first->level();
2939 Assert(my_level >= neighbor_level, ExcInternalError());
2940 return my_level > neighbor_level;
2941}
2942
2943
2944
2945template <int dim, int spacedim>
2946bool
2947CellAccessor<dim, spacedim>::at_boundary(const unsigned int i) const
2948{
2949 Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
2950 AssertIndexRange(i, this->n_faces());
2951
2952 return (neighbor_index(i) == -1);
2953}
2954
2955
2956
2957template <int dim, int spacedim>
2958bool
2959CellAccessor<dim, spacedim>::has_boundary_lines() const
2960{
2961 if (dim == 1)
2962 return at_boundary();
2963 else
2964 {
2965 for (unsigned int l = 0; l < this->n_lines(); ++l)
2966 if (this->line(l)->at_boundary())
2967 return true;
2968
2969 return false;
2970 }
2971}
2972
2973
2974
2975template <int dim, int spacedim>
2976TriaIterator<CellAccessor<dim, spacedim>>
2977CellAccessor<dim, spacedim>::neighbor_child_on_subface(
2978 const unsigned int face,
2979 const unsigned int subface) const
2980{
2981 Assert(!this->has_children(),
2982 ExcMessage("The present cell must not have children!"));
2983 Assert(!this->at_boundary(face),
2984 ExcMessage("The present cell must have a valid neighbor!"));
2985 Assert(this->neighbor(face)->has_children() == true,
2986 ExcMessage("The neighbor must have children!"));
2987
2988 switch (dim)
2989 {
2990 case 2:
2991 {
2992 Assert(this->reference_cell() == ReferenceCells::Triangle ||
2993 this->reference_cell() == ReferenceCells::Quadrilateral,
2994 ExcNotImplemented());
2995
2996 const auto neighbor = this->neighbor(face);
2997 if (neighbor->reference_cell() == ReferenceCells::Triangle)
2998 Assert(neighbor->refinement_case() ==
2999 RefinementCase<2>::isotropic_refinement,
3000 ExcNotImplemented());
3001
3002 const unsigned int neighbor_neighbor =
3003 this->neighbor_of_neighbor(face);
3004
3005 // if the refinement case is isotropic then ask the reference cell
3006 // if not we are on a quad and can ask GeometryInfo, we check for that
3007 // above
3008 const unsigned int neighbor_child_index =
3009 neighbor->refinement_case() ==
3010 RefinementCase<2>::isotropic_refinement ?
3011 neighbor->reference_cell().child_cell_on_face(
3012 neighbor_neighbor,
3013 subface,
3014 neighbor->combined_face_orientation(neighbor_neighbor)) :
3015 GeometryInfo<dim>::child_cell_on_face(neighbor->refinement_case(),
3016 neighbor_neighbor,
3017 subface,
3018 neighbor->face_orientation(
3019 neighbor_neighbor));
3020
3021 auto child = neighbor->child(neighbor_child_index);
3022 if (neighbor->reference_cell() == ReferenceCells::Triangle)
3023 Assert(!child->has_children(), ExcInternalError());
3024 while (child->has_children())
3025 {
3026 Assert((GeometryInfo<dim>::face_refinement_case(
3027 child->refinement_case(), neighbor_neighbor) ==
3028 RefinementCase<dim>::no_refinement),
3029 ExcInternalError());
3030 child = child->child(GeometryInfo<dim>::child_cell_on_face(
3031 child->refinement_case(), neighbor_neighbor, 0));
3032 }
3033
3034 return child;
3035 }
3036
3037
3038 case 3:
3039 {
3040 if (this->reference_cell() == ReferenceCells::Hexahedron)
3041 {
3042 // this function returns the neighbor's
3043 // child on a given face and
3044 // subface.
3045
3046 // we have to consider one other aspect here:
3047 // The face might be refined
3048 // anisotropically. In this case, the subface
3049 // number refers to the following, where we
3050 // look at the face from the current cell,
3051 // thus the subfaces are in standard
3052 // orientation concerning the cell
3053 //
3054 // for isotropic refinement
3055 //
3056 // *---*---*
3057 // | 2 | 3 |
3058 // *---*---*
3059 // | 0 | 1 |
3060 // *---*---*
3061 //
3062 // for 2*anisotropic refinement
3063 // (first cut_y, then cut_x)
3064 //
3065 // *---*---*
3066 // | 2 | 3 |
3067 // *---*---*
3068 // | 0 | 1 |
3069 // *---*---*
3070 //
3071 // for 2*anisotropic refinement
3072 // (first cut_x, then cut_y)
3073 //
3074 // *---*---*
3075 // | 1 | 3 |
3076 // *---*---*
3077 // | 0 | 2 |
3078 // *---*---*
3079 //
3080 // for purely anisotropic refinement:
3081 //
3082 // *---*---* *-------*
3083 // | | | | 1 |
3084 // | 0 | 1 | or *-------*
3085 // | | | | 0 |
3086 // *---*---* *-------*
3087 //
3088 // for "mixed" refinement:
3089 //
3090 // *---*---* *---*---* *---*---* *-------*
3091 // | | 2 | | 1 | | | 1 | 2 | | 2 |
3092 // | 0 *---* or *---* 2 | or *---*---* or *---*---*
3093 // | | 1 | | 0 | | | 0 | | 0 | 1 |
3094 // *---*---* *---*---* *-------* *---*---*
3095
3096 const typename Triangulation<dim, spacedim>::face_iterator
3097 mother_face = this->face(face);
3098 const unsigned int total_children =
3099 mother_face->n_active_descendants();
3100 AssertIndexRange(subface, total_children);
3101 Assert(total_children <= GeometryInfo<3>::max_children_per_face,
3102 ExcInternalError());
3103
3104 unsigned int neighbor_neighbor;
3105 TriaIterator<CellAccessor<dim, spacedim>> neighbor_child;
3106 const TriaIterator<CellAccessor<dim, spacedim>> neighbor =
3107 this->neighbor(face);
3108
3109
3110 const RefinementCase<dim - 1> mother_face_ref_case =
3111 mother_face->refinement_case();
3112 if (mother_face_ref_case ==
3113 static_cast<RefinementCase<dim - 1>>(
3114 RefinementCase<2>::cut_xy)) // total_children==4
3115 {
3116 // this case is quite easy. we are sure,
3117 // that the neighbor is not coarser.
3118
3119 // get the neighbor's number for the given
3120 // face and the neighbor
3121 neighbor_neighbor = this->neighbor_of_neighbor(face);
3122
3123 // now use the info provided by GeometryInfo
3124 // to extract the neighbors child number
3125 const unsigned int neighbor_child_index =
3126 GeometryInfo<dim>::child_cell_on_face(
3127 neighbor->refinement_case(),
3128 neighbor_neighbor,
3129 subface,
3130 neighbor->face_orientation(neighbor_neighbor),
3131 neighbor->face_flip(neighbor_neighbor),
3132 neighbor->face_rotation(neighbor_neighbor));
3133 neighbor_child = neighbor->child(neighbor_child_index);
3134
3135 // make sure that the neighbor child cell we
3136 // have found shares the desired subface.
3137 Assert((this->face(face)->child(subface) ==
3138 neighbor_child->face(neighbor_neighbor)),
3139 ExcInternalError());
3140 }
3141 else //-> the face is refined anisotropically
3142 {
3143 // first of all, we have to find the
3144 // neighbor at one of the anisotropic
3145 // children of the
3146 // mother_face. determine, which of
3147 // these we need.
3148 unsigned int first_child_to_find;
3149 unsigned int neighbor_child_index;
3150 if (total_children == 2)
3151 first_child_to_find = subface;
3152 else
3153 {
3154 first_child_to_find = subface / 2;
3155 if (total_children == 3 && subface == 1 &&
3156 !mother_face->child(0)->has_children())
3157 first_child_to_find = 1;
3158 }
3159 if (neighbor_is_coarser(face))
3160 {
3161 std::pair<unsigned int, unsigned int> indices =
3162 neighbor_of_coarser_neighbor(face);
3163 neighbor_neighbor = indices.first;
3164
3165
3166 // we have to translate our
3167 // subface_index according to the
3168 // RefineCase and subface index of
3169 // the coarser face (our face is an
3170 // anisotropic child of the coarser
3171 // face), 'a' denotes our
3172 // subface_index 0 and 'b' denotes
3173 // our subface_index 1, whereas 0...3
3174 // denote isotropic subfaces of the
3175 // coarser face
3176 //
3177 // cut_x and coarser_subface_index=0
3178 //
3179 // *---*---*
3180 // |b=2| |
3181 // | | |
3182 // |a=0| |
3183 // *---*---*
3184 //
3185 // cut_x and coarser_subface_index=1
3186 //
3187 // *---*---*
3188 // | |b=3|
3189 // | | |
3190 // | |a=1|
3191 // *---*---*
3192 //
3193 // cut_y and coarser_subface_index=0
3194 //
3195 // *-------*
3196 // | |
3197 // *-------*
3198 // |a=0 b=1|
3199 // *-------*
3200 //
3201 // cut_y and coarser_subface_index=1
3202 //
3203 // *-------*
3204 // |a=2 b=3|
3205 // *-------*
3206 // | |
3207 // *-------*
3208 unsigned int iso_subface;
3209 if (neighbor->face(neighbor_neighbor)
3210 ->refinement_case() == RefinementCase<2>::cut_x)
3211 iso_subface = 2 * first_child_to_find + indices.second;
3212 else
3213 {
3214 Assert(neighbor->face(neighbor_neighbor)
3215 ->refinement_case() ==
3216 RefinementCase<2>::cut_y,
3217 ExcInternalError());
3218 iso_subface =
3219 first_child_to_find + 2 * indices.second;
3220 }
3221 neighbor_child_index =
3222 GeometryInfo<dim>::child_cell_on_face(
3223 neighbor->refinement_case(),
3224 neighbor_neighbor,
3225 iso_subface,
3226 neighbor->face_orientation(neighbor_neighbor),
3227 neighbor->face_flip(neighbor_neighbor),
3228 neighbor->face_rotation(neighbor_neighbor));
3229 }
3230 else // neighbor is not coarser
3231 {
3232 neighbor_neighbor = neighbor_of_neighbor(face);
3233 neighbor_child_index =
3234 GeometryInfo<dim>::child_cell_on_face(
3235 neighbor->refinement_case(),
3236 neighbor_neighbor,
3237 first_child_to_find,
3238 neighbor->face_orientation(neighbor_neighbor),
3239 neighbor->face_flip(neighbor_neighbor),
3240 neighbor->face_rotation(neighbor_neighbor),
3241 mother_face_ref_case);
3242 }
3243
3244 neighbor_child = neighbor->child(neighbor_child_index);
3245 // it might be, that the neighbor_child
3246 // has children, which are not refined
3247 // along the given subface. go down that
3248 // list and deliver the last of those.
3249 while (
3250 neighbor_child->has_children() &&
3251 GeometryInfo<dim>::face_refinement_case(
3252 neighbor_child->refinement_case(), neighbor_neighbor) ==
3253 RefinementCase<2>::no_refinement)
3254 neighbor_child = neighbor_child->child(
3255 GeometryInfo<dim>::child_cell_on_face(
3256 neighbor_child->refinement_case(),
3257 neighbor_neighbor,
3258 0));
3259
3260 // if there are two total subfaces, we
3261 // are finished. if there are four we
3262 // have to get a child of our current
3263 // neighbor_child. If there are three,
3264 // we have to check which of the two
3265 // possibilities applies.
3266 if (total_children == 3)
3267 {
3268 if (mother_face->child(0)->has_children())
3269 {
3270 if (subface < 2)
3271 neighbor_child = neighbor_child->child(
3272 GeometryInfo<dim>::child_cell_on_face(
3273 neighbor_child->refinement_case(),
3274 neighbor_neighbor,
3275 subface,
3276 neighbor_child->face_orientation(
3277 neighbor_neighbor),
3278 neighbor_child->face_flip(neighbor_neighbor),
3279 neighbor_child->face_rotation(
3280 neighbor_neighbor),
3281 mother_face->child(0)->refinement_case()));
3282 }
3283 else
3284 {
3285 Assert(mother_face->child(1)->has_children(),
3286 ExcInternalError());
3287 if (subface > 0)
3288 neighbor_child = neighbor_child->child(
3289 GeometryInfo<dim>::child_cell_on_face(
3290 neighbor_child->refinement_case(),
3291 neighbor_neighbor,
3292 subface - 1,
3293 neighbor_child->face_orientation(
3294 neighbor_neighbor),
3295 neighbor_child->face_flip(neighbor_neighbor),
3296 neighbor_child->face_rotation(
3297 neighbor_neighbor),
3298 mother_face->child(1)->refinement_case()));
3299 }
3300 }
3301 else if (total_children == 4)
3302 {
3303 neighbor_child = neighbor_child->child(
3304 GeometryInfo<dim>::child_cell_on_face(
3305 neighbor_child->refinement_case(),
3306 neighbor_neighbor,
3307 subface % 2,
3308 neighbor_child->face_orientation(neighbor_neighbor),
3309 neighbor_child->face_flip(neighbor_neighbor),
3310 neighbor_child->face_rotation(neighbor_neighbor),
3311 mother_face->child(subface / 2)->refinement_case()));
3312 }
3313 }
3314
3315 // it might be, that the neighbor_child has
3316 // children, which are not refined along the
3317 // given subface. go down that list and
3318 // deliver the last of those.
3319 while (neighbor_child->has_children())
3320 neighbor_child =
3321 neighbor_child->child(GeometryInfo<dim>::child_cell_on_face(
3322 neighbor_child->refinement_case(), neighbor_neighbor, 0));
3323
3324 if constexpr (running_in_debug_mode())
3325 {
3326 // check, whether the face neighbor_child matches the
3327 // requested subface.
3328 typename Triangulation<dim, spacedim>::face_iterator
3329 requested;
3330 switch (this->subface_case(face))
3331 {
3332 case internal::SubfaceCase<3>::case_x:
3333 case internal::SubfaceCase<3>::case_y:
3334 case internal::SubfaceCase<3>::case_xy:
3335 requested = mother_face->child(subface);
3336 break;
3337 case internal::SubfaceCase<3>::case_x1y2y:
3338 case internal::SubfaceCase<3>::case_y1x2x:
3339 requested =
3340 mother_face->child(subface / 2)->child(subface % 2);
3341 break;
3342
3343 case internal::SubfaceCase<3>::case_x1y:
3344 case internal::SubfaceCase<3>::case_y1x:
3345 switch (subface)
3346 {
3347 case 0:
3348 case 1:
3349 requested = mother_face->child(0)->child(subface);
3350 break;
3351 case 2:
3352 requested = mother_face->child(1);
3353 break;
3354 default:
3355 DEAL_II_ASSERT_UNREACHABLE();
3356 }
3357 break;
3358 case internal::SubfaceCase<3>::case_x2y:
3359 case internal::SubfaceCase<3>::case_y2x:
3360 switch (subface)
3361 {
3362 case 0:
3363 requested = mother_face->child(0);
3364 break;
3365 case 1:
3366 case 2:
3367 requested =
3368 mother_face->child(1)->child(subface - 1);
3369 break;
3370 default:
3371 DEAL_II_ASSERT_UNREACHABLE();
3372 }
3373 break;
3374 default:
3375 DEAL_II_ASSERT_UNREACHABLE();
3376 break;
3377 }
3378 Assert(requested == neighbor_child->face(neighbor_neighbor),
3379 ExcInternalError());
3380 }
3381
3382 return neighbor_child;
3383 }
3384
3385 // if no reference cell type matches
3386 DEAL_II_NOT_IMPLEMENTED();
3387 return TriaIterator<CellAccessor<dim, spacedim>>();
3388 }
3389
3390 default:
3391 // if 1d or more than 3d
3392 DEAL_II_NOT_IMPLEMENTED();
3393 return TriaIterator<CellAccessor<dim, spacedim>>();
3394 }
3395}
3396
3397
3398
3399template <int structdim, int dim, int spacedim>
3400IteratorState::IteratorStates
3405
3406
3407
3408template <int structdim, int dim, int spacedim>
3409int
3414
3415
3416
3417template <int structdim, int dim, int spacedim>
3418int
3423
3424
3425// explicit instantiations
3426#include "grid/tria_accessor.inst"
3427
3428DEAL_II_NAMESPACE_CLOSE
make_array_view
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:1002
ArrayView::begin
iterator begin() const
Definition array_view.h:755
ArrayView::size
std::size_t size() const
Definition array_view.h:737
CellAccessor::recursively_set_subdomain_id
void recursively_set_subdomain_id(const types::subdomain_id new_subdomain_id) const
Definition tria_accessor.cc:2291
CellAccessor::parent
TriaIterator< CellAccessor< dim, spacedim > > parent() const
Definition tria_accessor.cc:2277
CellAccessor::set_active_cell_index
void set_active_cell_index(const unsigned int active_cell_index) const
Definition tria_accessor.cc:2244
CellAccessor::neighbor_of_neighbor_internal
unsigned int neighbor_of_neighbor_internal(const unsigned int neighbor) const
Definition tria_accessor.cc:2363
CellAccessor::periodic_neighbor
TriaIterator< CellAccessor< dim, spacedim > > periodic_neighbor(const unsigned int i) const
Definition tria_accessor.cc:2670
CellAccessor::set_direction_flag
void set_direction_flag(const bool new_direction_flag) const
Definition tria_accessor.cc:2190
CellAccessor::recursively_set_material_id
void recursively_set_material_id(const types::material_id new_material_id) const
Definition tria_accessor.cc:2150
CellAccessor::set_level_subdomain_id
void set_level_subdomain_id(const types::subdomain_id new_level_subdomain_id) const
Definition tria_accessor.cc:2178
CellAccessor::as_dof_handler_iterator
TriaActiveIterator< DoFCellAccessor< dim, spacedim, false > > as_dof_handler_iterator(const DoFHandler< dim, spacedim > &dof_handler) const
Definition tria_accessor.cc:2035
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:2977
CellAccessor::set_subdomain_id
void set_subdomain_id(const types::subdomain_id new_subdomain_id) const
Definition tria_accessor.cc:2164
CellAccessor::neighbor_is_coarser
bool neighbor_is_coarser(const unsigned int face_no) const
Definition tria_accessor.cc:2447
CellAccessor::set_global_level_cell_index
void set_global_level_cell_index(const types::global_cell_index index) const
Definition tria_accessor.cc:2266
CellAccessor::has_periodic_neighbor
bool has_periodic_neighbor(const unsigned int i) const
Definition tria_accessor.cc:2643
CellAccessor::periodic_neighbor_level
int periodic_neighbor_level(const unsigned int i) const
Definition tria_accessor.cc:2841
CellAccessor::neighbor_of_coarser_neighbor
std::pair< unsigned int, unsigned int > neighbor_of_coarser_neighbor(const unsigned int neighbor) const
Definition tria_accessor.cc:2458
CellAccessor::neighbor_of_neighbor
unsigned int neighbor_of_neighbor(const unsigned int face_no) const
Definition tria_accessor.cc:2433
CellAccessor::set_material_id
void set_material_id(const types::material_id new_material_id) const
Definition tria_accessor.cc:2136
CellAccessor::point_inside_codim
bool point_inside_codim(const Point< spacedim_ > &p) const
Definition tria_accessor.cc:2079
CellAccessor::has_boundary_lines
bool has_boundary_lines() const
Definition tria_accessor.cc:2959
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:2714
CellAccessor::periodic_neighbor_index
int periodic_neighbor_index(const unsigned int i) const
Definition tria_accessor.cc:2831
CellAccessor::periodic_neighbor_is_coarser
bool periodic_neighbor_is_coarser(const unsigned int i) const
Definition tria_accessor.cc:2892
CellAccessor::set_global_active_cell_index
void set_global_active_cell_index(const types::global_cell_index index) const
Definition tria_accessor.cc:2255
CellAccessor::get_cells_adjacent_to_line
std::set< TriaActiveIterator< CellAccessor< dim, spacedim > > > get_cells_adjacent_to_line(const unsigned int i) const
Definition tria_accessor.cc:2304
CellAccessor::set_parent
void set_parent(const unsigned int parent_index)
Definition tria_accessor.cc:2214
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:2770
CellAccessor::at_boundary
bool at_boundary() const
Definition tria_accessor.cc:2123
CellAccessor::point_inside
bool point_inside(const Point< spacedim > &p) const
CellAccessor::id
CellId id() const
Definition tria_accessor.cc:2322
CellAccessor::as_dof_handler_level_iterator
TriaIterator< DoFCellAccessor< dim, spacedim, true > > as_dof_handler_level_iterator(const DoFHandler< dim, spacedim > &dof_handler) const
Definition tria_accessor.cc:2056
CellAccessor::neighbor_or_periodic_neighbor
TriaIterator< CellAccessor< dim, spacedim > > neighbor_or_periodic_neighbor(const unsigned int i) const
Definition tria_accessor.cc:2697
CellAccessor::parent_index
int parent_index() const
Definition tria_accessor.cc:2230
CellAccessor::periodic_neighbor_of_periodic_neighbor
unsigned int periodic_neighbor_of_periodic_neighbor(const unsigned int i) const
Definition tria_accessor.cc:2851
CellAccessor::periodic_neighbor_face_no
unsigned int periodic_neighbor_face_no(const unsigned int i) const
Definition tria_accessor.cc:2861
DerivativeForm::covariant_form
DerivativeForm< 1, dim, spacedim, Number > covariant_form() const
DerivativeForm::transpose
DerivativeForm< 1, spacedim, dim, Number > transpose() const
DoFHandler::get_triangulation
const Triangulation< dim, spacedim > & get_triangulation() const
DoFHandler::level_cell_iterator
typename LevelSelector::cell_iterator level_cell_iterator
Definition dof_handler.h:511
InvalidAccessor::level
static int level()
Definition tria_accessor.cc:3410
InvalidAccessor::state
static IteratorState::IteratorStates state()
Definition tria_accessor.cc:3401
InvalidAccessor::index
static int index()
Definition tria_accessor.cc:3419
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:318
Point
Definition point.h:111
Tensor::norm
numbers::NumberTraits< Number >::real_type norm() const
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::set_user_index
void set_user_index(const unsigned int p) const
TriaAccessor::clear_user_pointer
void clear_user_pointer() const
TriaAccessor::recursively_set_user_index
void recursively_set_user_index(const unsigned int p) const
TriaAccessor::clear_user_data
void clear_user_data() 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:1670
TriaAccessor::recursively_clear_user_index
void recursively_clear_user_index() const
TriaAccessor::get_manifold
const Manifold< dim, spacedim > & get_manifold() const
TriaAccessor::recursively_set_user_pointer
void recursively_set_user_pointer(void *p) const
TriaAccessor::extent_in_direction
double extent_in_direction(const unsigned int axis) const
TriaAccessor::intermediate_point
Point< spacedim > intermediate_point(const Point< structdim > &coordinates) const
Definition tria_accessor.cc:1648
TriaAccessor::n_vertices
unsigned int n_vertices() const
TriaAccessor::recursively_clear_user_flag
void recursively_clear_user_flag() const
TriaAccessor::barycenter
Point< spacedim > barycenter() const
TriaAccessor::bounding_box
BoundingBox< spacedim > bounding_box() const
TriaAccessor::clear_user_flag
void clear_user_flag() const
TriaAccessor::recursively_set_user_flag
void recursively_set_user_flag() const
TriaAccessor::user_flag_set
bool user_flag_set() const
TriaAccessor::set_user_flag
void set_user_flag() const
TriaAccessor::user_pointer
void * user_pointer() const
TriaAccessor::center
Point< spacedim > center(const bool respect_manifold=false, const bool interpolate_from_surrounding=false) const
Definition tria_accessor.cc:1691
TriaAccessor::reference_cell
ReferenceCell< structdim > reference_cell() const
TriaAccessor::clear_user_index
void clear_user_index() const
TriaAccessor::measure
double measure() const
TriaAccessor::set_bounding_object_indices
void set_bounding_object_indices(const std::initializer_list< int > &new_indices) const
TriaAccessor::vertex
Point< spacedim > & vertex(const unsigned int i) const
TriaAccessor::user_index
unsigned int user_index() const
TriaAccessor::set_user_pointer
void set_user_pointer(void *p) const
TriaAccessor::recursively_clear_user_pointer
void recursively_clear_user_pointer() const
TriaAccessor::diameter
double diameter() const
Triangulation::get_vertices
const std::vector< Point< spacedim > > & get_vertices() const
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
running_in_debug_mode
constexpr bool running_in_debug_mode()
Definition config.h:76
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
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:497
dof_accessor.h
dof_handler.h
DEAL_II_ASSERT_UNREACHABLE
#define DEAL_II_ASSERT_UNREACHABLE()
Definition exception_macros.h:781
DEAL_II_NOT_IMPLEMENTED
#define DEAL_II_NOT_IMPLEMENTED()
Definition exception_macros.h:705
fe_q.h
geometry_info.h
level
unsigned int level
Definition grid_out.cc:4642
grid_tools_geometry.h
vertex_indices
unsigned int vertex_indices[2]
Definition grid_tools_topology.cc:941
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 exception_macros.h:543
StandardExceptions::ExcImpossibleInDim
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
TriaAccessorExceptions::ExcNeighborIsNotCoarser
static ::ExceptionBase & ExcNeighborIsNotCoarser()
TriaAccessorExceptions::ExcNeighborIsCoarser
static ::ExceptionBase & ExcNeighborIsCoarser()
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exception_macros.h:901
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition exception_macros.h:647
DoFHandler::active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition dof_handler.h:447
TriaAccessor::set_all_manifold_ids
void set_all_manifold_ids(const types::manifold_id) const
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:44
GridTools::diameter
double diameter(const Triangulation< dim, spacedim > &tria)
Definition grid_tools_geometry.cc:41
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:118
IteratorState::invalid
@ invalid
Iterator is invalid, probably due to an error.
Definition tria_iterator_base.h:39
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)
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::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
ReferenceCells::Hexahedron
constexpr ReferenceCell< 3 > Hexahedron
Definition reference_cell.h:1358
ReferenceCells::Quadrilateral
constexpr ReferenceCell< 2 > Quadrilateral
Definition reference_cell.h:1350
ReferenceCells::Triangle
constexpr ReferenceCell< 2 > Triangle
Definition reference_cell.h:1348
ReferenceCells::Tetrahedron
constexpr ReferenceCell< 3 > Tetrahedron
Definition reference_cell.h:1352
ReferenceCells::Pyramid
constexpr ReferenceCell< 3 > Pyramid
Definition reference_cell.h:1354
ReferenceCells::Wedge
constexpr ReferenceCell< 3 > Wedge
Definition reference_cell.h:1356
internal::split_face_orientation
std::tuple< bool, bool, bool > split_face_orientation(const types::geometric_orientation combined_orientation)
Definition tria_orientation.h:44
numbers::invalid_unsigned_int
constexpr unsigned int invalid_unsigned_int
Definition types.h:228
numbers::invalid_material_id
constexpr types::material_id invalid_material_id
Definition types.h:284
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6596
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6579
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6505
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6563
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.8