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