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