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