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