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