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