operator.h
Go to the documentation of this file.
1 /*
2  This file is part of the Ristra portage project.
3  Please see the license file at the root of this repository, or at:
4  https://github.com/laristra/portage/blob/master/LICENSE
5 */
6 #ifndef OPERATOR_H_INC_
7 #define OPERATOR_H_INC_
8 
9 #include <cassert>
10 #include <cmath>
11 #include <array>
12 #include <vector>
13 #include <stdexcept>
14 
15 #include "portage/support/basis.h"
16 #include "wonton/support/Point.h"
17 
18 namespace Portage { namespace Meshfree { namespace oper {
19 
20 using basis::Traits;
21 using basis::transfactor;
22 
26 enum Domain { Interval,
27  Quadrilateral, Triangle, Circle,
28  Hexahedron, Tetrahedron, Wedge, Sphere, LastDomain };
29 
30 enum Type { VolumeIntegral, SurfaceIntegral, LastOperator };
31 
32 template<Domain domain> class DomainTraits { public: static const size_t dimension=0; };
33 template<> class DomainTraits<Interval> { public: static const size_t dimension=1; };
34 template<> class DomainTraits<Quadrilateral>{ public: static const size_t dimension=2; };
35 template<> class DomainTraits<Triangle> { public: static const size_t dimension=2; };
36 template<> class DomainTraits<Circle> { public: static const size_t dimension=2; };
37 template<> class DomainTraits<Hexahedron> { public: static const size_t dimension=3; };
38 template<> class DomainTraits<Tetrahedron> { public: static const size_t dimension=3; };
39 template<> class DomainTraits<Wedge> { public: static const size_t dimension=3; };
40 template<> class DomainTraits<Sphere> { public: static const size_t dimension=3; };
41 
47 constexpr size_t dimension(Domain domain) {
48  switch (domain) {
49  case Interval: return DomainTraits<Interval>::dimension;
50  case Quadrilateral: return DomainTraits<Quadrilateral>::dimension;
51  case Triangle: return DomainTraits<Triangle>::dimension;
52  case Circle: return DomainTraits<Circle>::dimension;
53  case Hexahedron: return DomainTraits<Hexahedron>::dimension;
54  case Tetrahedron: return DomainTraits<Tetrahedron>::dimension;
55  case Wedge: return DomainTraits<Wedge>::dimension;
56  case Sphere: return DomainTraits<Sphere>::dimension;
57  default: return 0;
58  }
59 }
60 
67 template<int dim>
68 Domain domain_from_points(std::vector<Wonton::Point<dim>> const& points) {
69 
70  int const nb_points = points.size();
71 
72  switch(dim) {
73  case 1: return Interval;
74  case 2:
75  switch(nb_points) {
76  case 3: return Triangle;
77  case 4: return Quadrilateral;
78  default: throw std::runtime_error("invalid number of points");
79  }
80  case 3:
81  switch(nb_points) {
82  case 4: return Tetrahedron;
83  case 6: return Wedge;
84  case 8: return Hexahedron;
85  default: throw std::runtime_error("invalid number of points");
86  }
87  default: throw std::runtime_error("invalid dimension");
88  }
89 }
90 
92 // Template for integral operator base class
94 
95 template<Type type, basis::Type basis_type, Domain domain_type>
96 class OperatorBase {
97 public:
98  Type operator_type = type;
99  static constexpr basis::Type basis = basis_type;
100  static constexpr Domain domain = domain_type;
101 };
102 
104 // Template for integral operators
106 
107 template<Type type, basis::Type basis_type, Domain domain_type>
108 class Operator: public OperatorBase<type,basis_type,domain_type> {
109 public:
110  static constexpr size_t dim = dimension(domain_type);
111  static constexpr size_t operator_size = 0;
112  static constexpr size_t basis_size = 0;
113  static constexpr size_t point_size = 0;
114 
115  using result_t = std::array<std::array<double, operator_size>, basis_size>;
116  using points_t = std::array<Wonton::Point<dim>, point_size>;
117 
118  static result_t apply(const points_t p) {
119  throw std::runtime_error("invalid operator");
120  }
121 };
122 
124 // Specializations
126 
127 // 1D
128 template<>
129 class Operator<VolumeIntegral, basis::Unitary, Interval>:
130  public OperatorBase<VolumeIntegral, basis::Unitary, Interval>
131 {
132 public:
133  static constexpr size_t dim = dimension(Interval);
134  static constexpr size_t operator_size = 1;
135  static constexpr size_t basis_size =
136  basis::Traits<basis::Unitary, dim>::function_size;
137  static constexpr size_t point_size = 2;
138 
139  using result_t = std::array<std::array<double, operator_size>, basis_size>;
140  using points_t = std::array<Wonton::Point<dim>, point_size>;
141 
142  static result_t apply(const points_t p) {
143  result_t result;
144  result[0][0] = p[1][0] - p[0][0];
145  return result;
146  }
147 };
148 
149 
150 template<>
151 class Operator<VolumeIntegral, basis::Linear, Interval>:
152  public OperatorBase<VolumeIntegral, basis::Linear, Interval>
153 {
154 public:
155  static constexpr size_t dim = dimension(Interval);
156  static constexpr size_t operator_size = 1;
157  static constexpr size_t basis_size =
158  basis::Traits<basis::Linear, dim>::function_size;
159  static constexpr size_t point_size = 2;
160 
161  using result_t = std::array<std::array<double, operator_size>, basis_size>;
162  using points_t = std::array<Wonton::Point<dim>, point_size>;
163 
164  static result_t apply(const points_t p) {
165  result_t result;
166  result[0][0] = p[1][0] - p[0][0];
167  result[1][0] = .5*(pow(p[1][0],2.) - pow(p[0][0],2.));
168  return result;
169  }
170 };
171 
172 template<>
173 class Operator<VolumeIntegral, basis::Quadratic, Interval>:
174  public OperatorBase<VolumeIntegral, basis::Quadratic, Interval>
175 {
176 public:
177  static constexpr size_t dim = dimension(Interval);
178  static constexpr size_t operator_size = 1;
179  static constexpr size_t basis_size =
180  basis::Traits<basis::Quadratic, dim>::function_size;
181  static constexpr size_t point_size = 2;
182 
183  using result_t = std::array<std::array<double, operator_size>, basis_size>;
184  using points_t = std::array<Wonton::Point<dim>, point_size>;
185 
186  static result_t apply(const points_t p) {
187  result_t result;
188  result[0][0] = p[1][0] - p[0][0];
189  result[1][0] = .5*(pow(p[1][0],2.) - pow(p[0][0],2.));
190  result[2][0] = (pow(p[1][0],3.) - pow(p[0][0],3.))/6.;
191  return result;
192  }
193 };
194 
195 // 2D Triangle
196 
197 template<>
198 class Operator<VolumeIntegral, basis::Unitary, Triangle>:
199  public OperatorBase<VolumeIntegral, basis::Unitary, Interval>
200 {
201 public:
202  static constexpr size_t dim = dimension(Triangle);
203  static constexpr size_t operator_size = 1;
204  static constexpr size_t basis_size =
205  basis::Traits<basis::Unitary, dim>::function_size;
206  static constexpr size_t point_size = 3;
207 
208  using result_t = std::array<std::array<double, operator_size>, basis_size>;
209  using points_t = std::array<Wonton::Point<dim>, point_size>;
210 
211  static result_t apply(const points_t p) {
212  result_t result;
213  result[0][0] = .5*(-(p[1][1]*p[2][0]) + p[0][1]*(-p[1][0] + p[2][0]) + p[0][0]*(p[1][1] - p[2][1]) + p[1][0]*p[2][1]);
214  return result;
215  }
216 };
217 
218 template<>
219 class Operator<VolumeIntegral, basis::Linear, Triangle>:
220  public OperatorBase<VolumeIntegral, basis::Linear, Interval>
221 {
222 public:
223  static constexpr size_t dim = dimension(Triangle);
224  static constexpr size_t operator_size = 1;
225  static constexpr size_t basis_size =
226  basis::Traits<basis::Linear, dim>::function_size;
227  static constexpr size_t point_size = 3;
228 
229  using result_t = std::array<std::array<double, operator_size>, basis_size>;
230  using points_t = std::array<Wonton::Point<dim>, point_size>;
231 
232  static result_t apply(const points_t p) {
233  result_t result;
234  result[0][0] = .5*(-(p[1][1]*p[2][0]) + p[0][1]*(-p[1][0] + p[2][0]) + p[0][0]*(p[1][1] - p[2][1]) + p[1][0]*p[2][1]);
235  result[1][0] = ((p[0][0] + p[1][0] + p[2][0])*(-(p[1][1]*p[2][0]) + p[0][1]*(-p[1][0] + p[2][0]) + p[0][0]*(p[1][1] - p[2][1])+ p[1][0]*p[2][1]))/6.;
236  result[2][0] = -((p[0][1] + p[1][1] + p[2][1])*(p[0][1]*(p[1][0] - p[2][0]) + p[1][1]*p[2][0] - p[1][0]*p[2][1] + p[0][0]*(-p[1][1] + p[2][1])))/6.;
237  return result;
238  }
239 };
240 
241 template<>
242  class Operator<VolumeIntegral, basis::Quadratic, Triangle>:
243  public OperatorBase<VolumeIntegral, basis::Quadratic, Interval>
244  {
245  public:
246  static constexpr size_t dim = dimension(Triangle);
247  static constexpr size_t operator_size = 1;
248  static constexpr size_t basis_size =
249  basis::Traits<basis::Quadratic, dim>::function_size;
250  static constexpr size_t point_size = 3;
251 
252  using result_t = std::array<std::array<double, operator_size>, basis_size>;
253  using points_t = std::array<Wonton::Point<dim>, point_size>;
254 
255  static result_t apply(const points_t p) {
256  result_t result;
257  result[0][0] = .5 * (-(p[1][1] * p[2][0]) + p[0][1] * (-p[1][0] + p[2][0]) + p[0][0] * (p[1][1] - p[2][1]) + p[1][0] * p[2][1]);
258  result[1][0] = ((p[0][0] + p[1][0] + p[2][0]) * (-(p[1][1] * p[2][0]) + p[0][1] * (-p[1][0] + p[2][0]) + p[0][0] * (p[1][1] - p[2][1]) + p[1][0] * p[2][1])) / 6.;
259  result[2][0] = -((p[0][1] + p[1][1] + p[2][1]) * (p[0][1] * (p[1][0] - p[2][0]) + p[1][1] * p[2][0] - p[1][0] * p[2][1] + p[0][0] * (-p[1][1] + p[2][1]))) / 6.;
260  result[3][0] = ((pow(p[0][0], 2) + pow(p[1][0], 2) + pow(p[2][0], 2) + p[1][0] * p[2][0] + p[0][0] * (p[1][0] + p[2][0])) * (-(p[1][1] * p[2][0]) + p[0][1] * (-p[1][0] + p[2][0]) + p[0][0] * (p[1][1] - p[2][1]) + p[1][0] * p[2][1])) / 24.;
261  result[4][0] = (-(pow(p[1][1], 2) * pow(p[2][0], 2)) + pow(p[0][1], 2) * (-pow(p[1][0], 2) + pow(p[2][0], 2)) + pow(p[1][0], 2) * pow(p[2][1], 2) - 2 * pow(p[1][1], 2) * p[1][0] * p[2][0] + 2 * pow(p[2][1], 2) * p[1][0] * p[2][0] - 2 * p[0][0] * (-(pow(p[1][1], 2) * p[1][0]) + pow(p[0][1], 2) * (p[1][0] - p[2][0]) + pow(p[2][1], 2) * p[2][0]) + 2 * pow(p[1][0], 2) * p[1][1] * p[2][1] - 2 * pow(p[2][0], 2) * p[1][1] * p[2][1] + pow(p[0][0], 2) * (p[1][1] - p[2][1]) * (2 * p[0][1] + p[1][1] + p[2][1]) + p[0][1] * (-2 * pow(p[1][0], 2) * p[1][1] + 2 * pow(p[2][0], 2) * p[2][1])) / 24.;
262  result[5][0] = -((p[0][1] * (p[1][0] - p[2][0]) + p[1][1] * p[2][0] - p[1][0] * p[2][1] + p[0][0] * (-p[1][1] + p[2][1])) * (pow(p[0][1], 2) + pow(p[1][1], 2) + pow(p[2][1], 2) + p[1][1] * p[2][1] + p[0][1] * (p[1][1] + p[2][1]))) / 24.;
263  return result;
264  }
265  };
266 
267 // 2D Quadrilateral
268 
269 template<>
270 class Operator<VolumeIntegral, basis::Unitary, Quadrilateral>:
271  public OperatorBase<VolumeIntegral, basis::Unitary, Interval>
272 {
273 public:
274  static constexpr size_t dim = dimension(Quadrilateral);
275  static constexpr size_t operator_size = 1;
276  static constexpr size_t basis_size =
277  basis::Traits<basis::Unitary, dim>::function_size;
278  static constexpr size_t point_size = 4;
279 
280  using result_t = std::array<std::array<double, operator_size>, basis_size>;
281  using points_t = std::array<Wonton::Point<dim>, point_size>;
282 
283  static result_t apply(const points_t p) {
284  result_t result;
285  result[0][0] = (-(p[1][1] * p[2][0]) + p[1][0] * p[2][1] - p[2][1] * p[3][0] + p[0][1] * (-p[1][0] + p[3][0]) + p[0][0] * (p[1][1] - p[3][1]) + p[2][0] * p[3][1]) / 2;
286  return result;
287  }
288 };
289 
290 template<>
291 class Operator<VolumeIntegral, basis::Linear, Quadrilateral>:
292  public OperatorBase<VolumeIntegral, basis::Linear, Interval>
293 {
294 public:
295  static constexpr size_t dim = dimension(Quadrilateral);
296  static constexpr size_t operator_size = 1;
297  static constexpr size_t basis_size =
298  basis::Traits<basis::Linear, dim>::function_size;
299  static constexpr size_t point_size = 4;
300 
301  using result_t = std::array<std::array<double, operator_size>, basis_size>;
302  using points_t = std::array<Wonton::Point<dim>, point_size>;
303 
304  static result_t apply(const points_t p) {
305  result_t result;
306  result[0][0] = (-(p[1][1] * p[2][0]) + p[1][0] * p[2][1] - p[2][1] * p[3][0] + p[0][1] * (-p[1][0] + p[3][0]) + p[0][0] * (p[1][1] - p[3][1]) + p[2][0] * p[3][1]) / 2;
307  result[1][0] = ((-pow(p[1][0], 2) + pow(p[3][0], 2)) * p[0][1] - pow(p[2][0], 2) * p[1][1] - p[1][0] * p[1][1] * p[2][0] + pow(p[1][0], 2) * p[2][1] - pow(p[3][0], 2) * p[2][1] + p[1][0] * p[2][0] * p[2][1] - p[2][0] * p[2][1] * p[3][0] + pow(p[0][0], 2) * (p[1][1] - p[3][1]) + pow(p[2][0], 2) * p[3][1] + p[2][0] * p[3][0] * p[3][1] + p[0][0] * (p[1][0] * p[1][1] + p[0][1] * (-p[1][0] + p[3][0]) - p[3][0] * p[3][1])) / 6;
308  result[2][0] = ((pow(p[1][1], 2) - pow(p[3][1], 2)) * p[0][0] + pow(p[2][1], 2) * p[1][0] - pow(p[1][1], 2) * p[2][0] + pow(p[3][1], 2) * p[2][0] + p[1][0] * p[1][1] * p[2][1] - p[1][1] * p[2][0] * p[2][1] - pow(p[2][1], 2) * p[3][0] + pow(p[0][1], 2) * (-p[1][0] + p[3][0]) + p[2][0] * p[2][1] * p[3][1] - p[2][1] * p[3][0] * p[3][1] + p[0][1] * (-(p[1][0] * p[1][1]) + p[0][0] * (p[1][1] - p[3][1]) + p[3][0] * p[3][1])) / 6;
309  return result;
310  }
311 };
312 
313 template<>
314 class Operator<VolumeIntegral, basis::Quadratic, Quadrilateral>:
315  public OperatorBase<VolumeIntegral, basis::Quadratic, Interval>
316 {
317 public:
318  static constexpr size_t dim = dimension(Quadrilateral);
319  static constexpr size_t operator_size = 1;
320  static constexpr size_t basis_size =
321  basis::Traits<basis::Quadratic, dim>::function_size;
322  static constexpr size_t point_size = 4;
323 
324  using result_t = std::array<std::array<double, operator_size>, basis_size>;
325  using points_t = std::array<Wonton::Point<dim>, point_size>;
326 
327  static result_t apply(const points_t p) {
328  result_t result;
329  result[0][0] = (-(p[1][1]*p[2][0]) + p[1][0]*p[2][1] - p[2][1]*p[3][0] + p[0][1]*(-p[1][0] + p[3][0]) + p[0][0]*(p[1][1] - p[3][1]) + p[2][0]*p[3][1])/2;
330  result[1][0] = ((-pow(p[1][0], 2) + pow(p[3][0], 2))*p[0][1] - pow(p[2][0], 2)*p[1][1] - p[1][0]*p[1][1]*p[2][0] + pow(p[1][0], 2)*p[2][1] - pow(p[3][0], 2)*p[2][1] + p[1][0]*p[2][0]*p[2][1] - p[2][0]*p[2][1]*p[3][0] + pow(p[0][0], 2)*(p[1][1] - p[3][1]) + pow(p[2][0], 2)*p[3][1] + p[2][0]*p[3][0]*p[3][1] + p[0][0]*(p[1][0]*p[1][1] + p[0][1]*(-p[1][0] + p[3][0]) - p[3][0]*p[3][1]))/6;
331  result[2][0] = ((pow(p[1][1], 2) - pow(p[3][1], 2))*p[0][0] + pow(p[2][1], 2)*p[1][0] - pow(p[1][1], 2)*p[2][0] + pow(p[3][1], 2)*p[2][0] + p[1][0]*p[1][1]*p[2][1] - p[1][1]*p[2][0]*p[2][1] - pow(p[2][1], 2)*p[3][0] + pow(p[0][1], 2)*(-p[1][0] + p[3][0]) + p[2][0]*p[2][1]*p[3][1] - p[2][1]*p[3][0]*p[3][1] + p[0][1]*(-(p[1][0]*p[1][1]) + p[0][0]*(p[1][1] - p[3][1]) + p[3][0]*p[3][1]))/6;
332  result[3][0] = ((-pow(p[1][0], 3) + pow(p[3][0], 3))*p[0][1] - pow(p[2][0], 3)*p[1][1] - pow(p[2][0], 2)*p[1][0]*p[1][1] - pow(p[1][0], 2)*p[1][1]*p[2][0] + pow(p[1][0], 3)*p[2][1] - pow(p[3][0], 3)*p[2][1] + pow(p[2][0], 2)*p[1][0]*p[2][1] + pow(p[1][0], 2)*p[2][0]*p[2][1] - pow(p[3][0], 2)*p[2][0]*p[2][1] - pow(p[2][0], 2)*p[2][1]*p[3][0] + pow(p[0][0], 3)*(p[1][1] - p[3][1]) + pow(p[2][0], 3)*p[3][1] + pow(p[3][0], 2)*p[2][0]*p[3][1] + pow(p[2][0], 2)*p[3][0]*p[3][1] + p[0][0]*((-pow(p[1][0], 2) + pow(p[3][0], 2))*p[0][1] + pow(p[1][0], 2)*p[1][1] - pow(p[3][0], 2)*p[3][1]) + pow(p[0][0], 2)*(p[1][0]*p[1][1] + p[0][1]*(-p[1][0] + p[3][0]) - p[3][0]*p[3][1]))/24;
333  result[4][0] = (-(pow(p[1][1], 2)*pow(p[2][0], 2)) + pow(p[1][0], 2)*pow(p[2][1], 2) - pow(p[2][1], 2)*pow(p[3][0], 2) + pow(p[0][1], 2)*(-pow(p[1][0], 2) + pow(p[3][0], 2)) + pow(p[2][0], 2)*pow(p[3][1], 2) - 2*pow(p[1][1], 2)*p[1][0]*p[2][0] + 2*pow(p[2][1], 2)*p[1][0]*p[2][0] + 2*pow(p[1][0], 2)*p[1][1]*p[2][1] - 2*pow(p[2][0], 2)*p[1][1]*p[2][1] - 2*pow(p[2][1], 2)*p[2][0]*p[3][0] + 2*pow(p[3][1], 2)*p[2][0]*p[3][0] - 2*p[0][0]*(-(pow(p[1][1], 2)*p[1][0]) + pow(p[0][1], 2)*(p[1][0] - p[3][0]) + pow(p[3][1], 2)*p[3][0]) + 2*pow(p[2][0], 2)*p[2][1]*p[3][1] - 2*pow(p[3][0], 2)*p[2][1]*p[3][1] + pow(p[0][0], 2)*(p[1][1] - p[3][1])*(2*p[0][1] + p[1][1] + p[3][1]) + p[0][1]*(-2*pow(p[1][0], 2)*p[1][1] + 2*pow(p[3][0], 2)*p[3][1]))/24;
334  result[5][0] = ((pow(p[1][1], 3) - pow(p[3][1], 3))*p[0][0] + pow(p[2][1], 3)*p[1][0] + pow(p[2][1], 2)*p[1][0]*p[1][1] - pow(p[1][1], 3)*p[2][0] + pow(p[3][1], 3)*p[2][0] - pow(p[2][1], 2)*p[1][1]*p[2][0] + pow(p[1][1], 2)*p[1][0]*p[2][1] - pow(p[1][1], 2)*p[2][0]*p[2][1] + pow(p[3][1], 2)*p[2][0]*p[2][1] - pow(p[2][1], 3)*p[3][0] - pow(p[3][1], 2)*p[2][1]*p[3][0] + pow(p[0][1], 3)*(-p[1][0] + p[3][0]) + p[0][1]*((pow(p[1][1], 2) - pow(p[3][1], 2))*p[0][0] - pow(p[1][1], 2)*p[1][0] + pow(p[3][1], 2)*p[3][0]) + pow(p[2][1], 2)*p[2][0]*p[3][1] - pow(p[2][1], 2)*p[3][0]*p[3][1] + pow(p[0][1], 2)*(-(p[1][0]*p[1][1]) + p[0][0]*(p[1][1] - p[3][1]) + p[3][0]*p[3][1]))/24;
335  return result;
336  }
337 };
338 
339 // 3D Hexahedron
340 
341 template<>
342 class Operator<VolumeIntegral, basis::Unitary, Hexahedron>:
343  public OperatorBase<VolumeIntegral, basis::Unitary, Interval>
344 {
345 public:
346  static constexpr size_t dim = dimension(Hexahedron);
347  static constexpr size_t operator_size = 1;
348  static constexpr size_t basis_size =
349  basis::Traits<basis::Unitary, dim>::function_size;
350  static constexpr size_t point_size = 8;
351 
352  using result_t = std::array<std::array<double, operator_size>, basis_size>;
353  using points_t = std::array<Wonton::Point<dim>, point_size>;
354 
355  static result_t apply(const points_t p) {
356  result_t result;
357  result[0][0] = (p[0][0]*p[1][2]*p[2][1] - p[0][0]*p[1][1]*p[2][2] + p[1][2]*p[2][1]*p[3][0] - p[1][1]*p[2][2]*p[3][0] + p[0][0]*p[1][2]*p[3][1] - p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] + p[1][0]*p[2][2]*p[3][1] - p[0][0]*p[1][1]*p[3][2] + p[1][1]*p[2][0]*p[3][2] - p[0][0]*p[2][1]*p[3][2] - p[1][0]*p[2][1]*p[3][2] - p[0][0]*p[1][2]*p[4][1] + p[0][0]*p[3][2]*p[4][1] + p[0][0]*p[1][1]*p[4][2] - p[0][0]*p[3][1]*p[4][2] - p[1][2]*p[2][1]*p[5][0] + p[1][1]*p[2][2]*p[5][0] + p[1][2]*p[4][1]*p[5][0] - p[1][1]*p[4][2]*p[5][0] - p[0][0]*p[1][2]*p[5][1] + p[1][2]*p[2][0]*p[5][1] - p[1][0]*p[2][2]*p[5][1] - p[1][2]*p[4][0]*p[5][1] + p[0][0]*p[4][2]*p[5][1] + p[1][0]*p[4][2]*p[5][1] + p[0][0]*p[1][1]*p[5][2] - p[1][1]*p[2][0]*p[5][2] + p[1][0]*p[2][1]*p[5][2] + p[1][1]*p[4][0]*p[5][2] - p[0][0]*p[4][1]*p[5][2] - p[1][0]*p[4][1]*p[5][2] - p[1][2]*p[2][1]*p[6][0] + p[1][1]*p[2][2]*p[6][0] - p[2][2]*p[3][1]*p[6][0] + p[2][1]*p[3][2]*p[6][0] + p[1][2]*p[5][1]*p[6][0] + p[2][2]*p[5][1]*p[6][0] - p[4][2]*p[5][1]*p[6][0] - p[1][1]*p[5][2]*p[6][0] - p[2][1]*p[5][2]*p[6][0] + p[4][1]*p[5][2]*p[6][0] + p[1][2]*p[2][0]*p[6][1] - p[1][0]*p[2][2]*p[6][1] + p[2][2]*p[3][0]*p[6][1] - p[2][0]*p[3][2]*p[6][1] - p[1][2]*p[5][0]*p[6][1] - p[2][2]*p[5][0]*p[6][1] + p[4][2]*p[5][0]*p[6][1] + p[1][0]*p[5][2]*p[6][1] + p[2][0]*p[5][2]*p[6][1] - p[4][0]*p[5][2]*p[6][1] - p[1][1]*p[2][0]*p[6][2] + p[1][0]*p[2][1]*p[6][2] - p[2][1]*p[3][0]*p[6][2] + p[2][0]*p[3][1]*p[6][2] + p[1][1]*p[5][0]*p[6][2] + p[2][1]*p[5][0]*p[6][2] - p[4][1]*p[5][0]*p[6][2] - p[1][0]*p[5][1]*p[6][2] - p[2][0]*p[5][1]*p[6][2] + p[4][0]*p[5][1]*p[6][2] - p[2][2]*p[3][1]*p[7][0] + p[2][1]*p[3][2]*p[7][0] - p[3][2]*p[4][1]*p[7][0] + p[3][1]*p[4][2]*p[7][0] - p[4][2]*p[5][1]*p[7][0] + p[4][1]*p[5][2]*p[7][0] + p[2][2]*p[6][1]*p[7][0] + p[3][2]*p[6][1]*p[7][0] - p[4][2]*p[6][1]*p[7][0] - p[5][2]*p[6][1]*p[7][0] - p[2][1]*p[6][2]*p[7][0] - p[3][1]*p[6][2]*p[7][0] + p[4][1]*p[6][2]*p[7][0] + p[5][1]*p[6][2]*p[7][0] + p[2][2]*p[3][0]*p[7][1] + p[0][0]*p[3][2]*p[7][1] - p[2][0]*p[3][2]*p[7][1] + p[3][2]*p[4][0]*p[7][1] - p[0][0]*p[4][2]*p[7][1] - p[3][0]*p[4][2]*p[7][1] + p[4][2]*p[5][0]*p[7][1] - p[4][0]*p[5][2]*p[7][1] - p[2][2]*p[6][0]*p[7][1] - p[3][2]*p[6][0]*p[7][1] + p[4][2]*p[6][0]*p[7][1] + p[5][2]*p[6][0]*p[7][1] + p[2][0]*p[6][2]*p[7][1] + p[3][0]*p[6][2]*p[7][1] - p[4][0]*p[6][2]*p[7][1] - p[5][0]*p[6][2]*p[7][1] + p[0][2]*(p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[3][1]*p[4][0] - p[3][0]*p[4][1] + p[1][1]*(p[2][0] + p[3][0] - p[4][0] - p[5][0]) + p[4][1]*p[5][0] - p[4][0]*p[5][1] + p[1][0]*(-p[2][1] - p[3][1] + p[4][1] + p[5][1]) + p[3][1]*p[7][0] - p[4][1]*p[7][0] - p[3][0]*p[7][1] + p[4][0]*p[7][1]) - p[2][1]*p[3][0]*p[7][2] - p[0][0]*p[3][1]*p[7][2] + p[2][0]*p[3][1]*p[7][2] - p[3][1]*p[4][0]*p[7][2] + p[0][0]*p[4][1]*p[7][2] + p[3][0]*p[4][1]*p[7][2] - p[4][1]*p[5][0]*p[7][2] + p[4][0]*p[5][1]*p[7][2] + p[2][1]*p[6][0]*p[7][2] + p[3][1]*p[6][0]*p[7][2] - p[4][1]*p[6][0]*p[7][2] - p[5][1]*p[6][0]*p[7][2] - p[2][0]*p[6][1]*p[7][2] - p[3][0]*p[6][1]*p[7][2] + p[4][0]*p[6][1]*p[7][2] + p[5][0]*p[6][1]*p[7][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[2][0]*p[3][2] - p[3][2]*p[4][0] + p[3][0]*p[4][2] - p[4][2]*p[5][0] + p[1][2]*(-p[2][0] - p[3][0] + p[4][0] + p[5][0]) + p[1][0]*(p[2][2] + p[3][2] - p[4][2] - p[5][2]) + p[4][0]*p[5][2] - p[3][2]*p[7][0] + p[4][2]*p[7][0] + p[3][0]*p[7][2] - p[4][0]*p[7][2]))/12;
358  return result;
359  }
360 };
361 
362 template<>
363 class Operator<VolumeIntegral, basis::Linear, Hexahedron>:
364  public OperatorBase<VolumeIntegral, basis::Linear, Interval>
365 {
366 public:
367  static constexpr size_t dim = dimension(Hexahedron);
368  static constexpr size_t operator_size = 1;
369  static constexpr size_t basis_size =
370  basis::Traits<basis::Linear, dim>::function_size;
371  static constexpr size_t point_size = 8;
372 
373  using result_t = std::array<std::array<double, operator_size>, basis_size>;
374  using points_t = std::array<Wonton::Point<dim>, point_size>;
375 
376  static result_t apply(const points_t p) {
377  result_t result;
378  result[0][0] = (p[0][0]*p[1][2]*p[2][1] - p[0][0]*p[1][1]*p[2][2] + p[1][2]*p[2][1]*p[3][0] - p[1][1]*p[2][2]*p[3][0] + p[0][0]*p[1][2]*p[3][1] - p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] + p[1][0]*p[2][2]*p[3][1] - p[0][0]*p[1][1]*p[3][2] + p[1][1]*p[2][0]*p[3][2] - p[0][0]*p[2][1]*p[3][2] - p[1][0]*p[2][1]*p[3][2] - p[0][0]*p[1][2]*p[4][1] + p[0][0]*p[3][2]*p[4][1] + p[0][0]*p[1][1]*p[4][2] - p[0][0]*p[3][1]*p[4][2] - p[1][2]*p[2][1]*p[5][0] + p[1][1]*p[2][2]*p[5][0] + p[1][2]*p[4][1]*p[5][0] - p[1][1]*p[4][2]*p[5][0] - p[0][0]*p[1][2]*p[5][1] + p[1][2]*p[2][0]*p[5][1] - p[1][0]*p[2][2]*p[5][1] - p[1][2]*p[4][0]*p[5][1] + p[0][0]*p[4][2]*p[5][1] + p[1][0]*p[4][2]*p[5][1] + p[0][0]*p[1][1]*p[5][2] - p[1][1]*p[2][0]*p[5][2] + p[1][0]*p[2][1]*p[5][2] + p[1][1]*p[4][0]*p[5][2] - p[0][0]*p[4][1]*p[5][2] - p[1][0]*p[4][1]*p[5][2] - p[1][2]*p[2][1]*p[6][0] + p[1][1]*p[2][2]*p[6][0] - p[2][2]*p[3][1]*p[6][0] + p[2][1]*p[3][2]*p[6][0] + p[1][2]*p[5][1]*p[6][0] + p[2][2]*p[5][1]*p[6][0] - p[4][2]*p[5][1]*p[6][0] - p[1][1]*p[5][2]*p[6][0] - p[2][1]*p[5][2]*p[6][0] + p[4][1]*p[5][2]*p[6][0] + p[1][2]*p[2][0]*p[6][1] - p[1][0]*p[2][2]*p[6][1] + p[2][2]*p[3][0]*p[6][1] - p[2][0]*p[3][2]*p[6][1] - p[1][2]*p[5][0]*p[6][1] - p[2][2]*p[5][0]*p[6][1] + p[4][2]*p[5][0]*p[6][1] + p[1][0]*p[5][2]*p[6][1] + p[2][0]*p[5][2]*p[6][1] - p[4][0]*p[5][2]*p[6][1] - p[1][1]*p[2][0]*p[6][2] + p[1][0]*p[2][1]*p[6][2] - p[2][1]*p[3][0]*p[6][2] + p[2][0]*p[3][1]*p[6][2] + p[1][1]*p[5][0]*p[6][2] + p[2][1]*p[5][0]*p[6][2] - p[4][1]*p[5][0]*p[6][2] - p[1][0]*p[5][1]*p[6][2] - p[2][0]*p[5][1]*p[6][2] + p[4][0]*p[5][1]*p[6][2] - p[2][2]*p[3][1]*p[7][0] + p[2][1]*p[3][2]*p[7][0] - p[3][2]*p[4][1]*p[7][0] + p[3][1]*p[4][2]*p[7][0] - p[4][2]*p[5][1]*p[7][0] + p[4][1]*p[5][2]*p[7][0] + p[2][2]*p[6][1]*p[7][0] + p[3][2]*p[6][1]*p[7][0] - p[4][2]*p[6][1]*p[7][0] - p[5][2]*p[6][1]*p[7][0] - p[2][1]*p[6][2]*p[7][0] - p[3][1]*p[6][2]*p[7][0] + p[4][1]*p[6][2]*p[7][0] + p[5][1]*p[6][2]*p[7][0] + p[2][2]*p[3][0]*p[7][1] + p[0][0]*p[3][2]*p[7][1] - p[2][0]*p[3][2]*p[7][1] + p[3][2]*p[4][0]*p[7][1] - p[0][0]*p[4][2]*p[7][1] - p[3][0]*p[4][2]*p[7][1] + p[4][2]*p[5][0]*p[7][1] - p[4][0]*p[5][2]*p[7][1] - p[2][2]*p[6][0]*p[7][1] - p[3][2]*p[6][0]*p[7][1] + p[4][2]*p[6][0]*p[7][1] + p[5][2]*p[6][0]*p[7][1] + p[2][0]*p[6][2]*p[7][1] + p[3][0]*p[6][2]*p[7][1] - p[4][0]*p[6][2]*p[7][1] - p[5][0]*p[6][2]*p[7][1] + p[0][2]*(p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[3][1]*p[4][0] - p[3][0]*p[4][1] + p[1][1]*(p[2][0] + p[3][0] - p[4][0] - p[5][0]) + p[4][1]*p[5][0] - p[4][0]*p[5][1] + p[1][0]*(-p[2][1] - p[3][1] + p[4][1] + p[5][1]) + p[3][1]*p[7][0] - p[4][1]*p[7][0] - p[3][0]*p[7][1] + p[4][0]*p[7][1]) - p[2][1]*p[3][0]*p[7][2] - p[0][0]*p[3][1]*p[7][2] + p[2][0]*p[3][1]*p[7][2] - p[3][1]*p[4][0]*p[7][2] + p[0][0]*p[4][1]*p[7][2] + p[3][0]*p[4][1]*p[7][2] - p[4][1]*p[5][0]*p[7][2] + p[4][0]*p[5][1]*p[7][2] + p[2][1]*p[6][0]*p[7][2] + p[3][1]*p[6][0]*p[7][2] - p[4][1]*p[6][0]*p[7][2] - p[5][1]*p[6][0]*p[7][2] - p[2][0]*p[6][1]*p[7][2] - p[3][0]*p[6][1]*p[7][2] + p[4][0]*p[6][1]*p[7][2] + p[5][0]*p[6][1]*p[7][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[2][0]*p[3][2] - p[3][2]*p[4][0] + p[3][0]*p[4][2] - p[4][2]*p[5][0] + p[1][2]*(-p[2][0] - p[3][0] + p[4][0] + p[5][0]) + p[1][0]*(p[2][2] + p[3][2] - p[4][2] - p[5][2]) + p[4][0]*p[5][2] - p[3][2]*p[7][0] + p[4][2]*p[7][0] + p[3][0]*p[7][2] - p[4][0]*p[7][2]))/12;
379  result[1][0] = (-(pow(p[2][0], 2)*p[0][1]*p[1][2]) - pow(p[3][0], 2)*p[0][1]*p[1][2] + pow(p[4][0], 2)*p[0][1]*p[1][2] + pow(p[5][0], 2)*p[0][1]*p[1][2] - 2*p[0][1]*p[1][0]*p[1][2]*p[2][0] + pow(p[3][0], 2)*p[1][2]*p[2][1] - pow(p[5][0], 2)*p[1][2]*p[2][1] - pow(p[6][0], 2)*p[1][2]*p[2][1] + 2*pow(p[1][0], 2)*p[0][1]*p[2][2] - 2*pow(p[3][0], 2)*p[0][1]*p[2][2] - pow(p[3][0], 2)*p[1][1]*p[2][2] + pow(p[5][0], 2)*p[1][1]*p[2][2] + pow(p[6][0], 2)*p[1][1]*p[2][2] + p[0][1]*p[1][0]*p[2][0]*p[2][2] - p[0][1]*p[1][0]*p[1][2]*p[3][0] - p[0][1]*p[1][2]*p[2][0]*p[3][0] + p[1][0]*p[1][2]*p[2][1]*p[3][0] + 2*p[1][2]*p[2][0]*p[2][1]*p[3][0] - p[1][0]*p[1][1]*p[2][2]*p[3][0] - p[0][1]*p[2][0]*p[2][2]*p[3][0] - 2*p[1][1]*p[2][0]*p[2][2]*p[3][0] - 2*pow(p[2][0], 2)*p[1][2]*p[3][1] - p[1][0]*p[1][2]*p[2][0]*p[3][1] + pow(p[1][0], 2)*p[2][2]*p[3][1] - pow(p[6][0], 2)*p[2][2]*p[3][1] - pow(p[7][0], 2)*p[2][2]*p[3][1] + 2*p[1][0]*p[2][0]*p[2][2]*p[3][1] - p[1][2]*p[2][0]*p[3][0]*p[3][1] + p[1][0]*p[2][2]*p[3][0]*p[3][1] + pow(p[1][0], 2)*p[0][1]*p[3][2] + pow(p[2][0], 2)*p[0][1]*p[3][2] - pow(p[4][0], 2)*p[0][1]*p[3][2] - pow(p[7][0], 2)*p[0][1]*p[3][2] + 2*pow(p[2][0], 2)*p[1][1]*p[3][2] + p[0][1]*p[1][0]*p[2][0]*p[3][2] + p[1][0]*p[1][1]*p[2][0]*p[3][2] - pow(p[1][0], 2)*p[2][1]*p[3][2] + pow(p[6][0], 2)*p[2][1]*p[3][2] + pow(p[7][0], 2)*p[2][1]*p[3][2] - 2*p[1][0]*p[2][0]*p[2][1]*p[3][2] + p[0][1]*p[1][0]*p[3][0]*p[3][2] + 2*p[0][1]*p[2][0]*p[3][0]*p[3][2] + p[1][1]*p[2][0]*p[3][0]*p[3][2] - p[1][0]*p[2][1]*p[3][0]*p[3][2] + p[0][1]*p[1][0]*p[1][2]*p[4][0] - p[0][1]*p[3][0]*p[3][2]*p[4][0] + 2*pow(p[5][0], 2)*p[1][2]*p[4][1] - 2*pow(p[7][0], 2)*p[3][2]*p[4][1] - pow(p[1][0], 2)*p[0][1]*p[4][2] + pow(p[3][0], 2)*p[0][1]*p[4][2] - pow(p[5][0], 2)*p[0][1]*p[4][2] + pow(p[7][0], 2)*p[0][1]*p[4][2] - 2*pow(p[5][0], 2)*p[1][1]*p[4][2] + 2*pow(p[7][0], 2)*p[3][1]*p[4][2] - p[0][1]*p[1][0]*p[4][0]*p[4][2] + p[0][1]*p[3][0]*p[4][0]*p[4][2] + 2*p[0][1]*p[1][0]*p[1][2]*p[5][0] - 2*p[1][0]*p[1][2]*p[2][1]*p[5][0] - p[1][2]*p[2][0]*p[2][1]*p[5][0] + 2*p[1][0]*p[1][1]*p[2][2]*p[5][0] + p[1][1]*p[2][0]*p[2][2]*p[5][0] + p[0][1]*p[1][2]*p[4][0]*p[5][0] + p[1][0]*p[1][2]*p[4][1]*p[5][0] + p[1][2]*p[4][0]*p[4][1]*p[5][0] - p[0][1]*p[1][0]*p[4][2]*p[5][0] - p[1][0]*p[1][1]*p[4][2]*p[5][0] - 2*p[0][1]*p[4][0]*p[4][2]*p[5][0] - p[1][1]*p[4][0]*p[4][2]*p[5][0] + pow(p[2][0], 2)*p[1][2]*p[5][1] - pow(p[4][0], 2)*p[1][2]*p[5][1] + pow(p[6][0], 2)*p[1][2]*p[5][1] + 2*p[1][0]*p[1][2]*p[2][0]*p[5][1] - 2*pow(p[1][0], 2)*p[2][2]*p[5][1] + 2*pow(p[6][0], 2)*p[2][2]*p[5][1] - p[1][0]*p[2][0]*p[2][2]*p[5][1] - p[1][0]*p[1][2]*p[4][0]*p[5][1] + pow(p[1][0], 2)*p[4][2]*p[5][1] - pow(p[6][0], 2)*p[4][2]*p[5][1] - pow(p[7][0], 2)*p[4][2]*p[5][1] + p[1][0]*p[4][0]*p[4][2]*p[5][1] + p[1][2]*p[2][0]*p[5][0]*p[5][1] - p[1][0]*p[2][2]*p[5][0]*p[5][1] - 2*p[1][2]*p[4][0]*p[5][0]*p[5][1] + 2*p[1][0]*p[4][2]*p[5][0]*p[5][1] - 2*pow(p[1][0], 2)*p[0][1]*p[5][2] + 2*pow(p[4][0], 2)*p[0][1]*p[5][2] - pow(p[2][0], 2)*p[1][1]*p[5][2] + pow(p[4][0], 2)*p[1][1]*p[5][2] - pow(p[6][0], 2)*p[1][1]*p[5][2] - 2*p[1][0]*p[1][1]*p[2][0]*p[5][2] + 2*pow(p[1][0], 2)*p[2][1]*p[5][2] - 2*pow(p[6][0], 2)*p[2][1]*p[5][2] + p[1][0]*p[2][0]*p[2][1]*p[5][2] + p[1][0]*p[1][1]*p[4][0]*p[5][2] - pow(p[1][0], 2)*p[4][1]*p[5][2] + pow(p[6][0], 2)*p[4][1]*p[5][2] + pow(p[7][0], 2)*p[4][1]*p[5][2] - p[1][0]*p[4][0]*p[4][1]*p[5][2] - p[0][1]*p[1][0]*p[5][0]*p[5][2] - p[1][1]*p[2][0]*p[5][0]*p[5][2] + p[1][0]*p[2][1]*p[5][0]*p[5][2] + p[0][1]*p[4][0]*p[5][0]*p[5][2] + 2*p[1][1]*p[4][0]*p[5][0]*p[5][2] - 2*p[1][0]*p[4][1]*p[5][0]*p[5][2] - p[1][0]*p[1][2]*p[2][1]*p[6][0] - 2*p[1][2]*p[2][0]*p[2][1]*p[6][0] + p[1][0]*p[1][1]*p[2][2]*p[6][0] + 2*p[1][1]*p[2][0]*p[2][2]*p[6][0] - 2*p[2][0]*p[2][2]*p[3][1]*p[6][0] - p[2][2]*p[3][0]*p[3][1]*p[6][0] + 2*p[2][0]*p[2][1]*p[3][2]*p[6][0] + p[2][1]*p[3][0]*p[3][2]*p[6][0] - p[1][2]*p[2][1]*p[5][0]*p[6][0] + p[1][1]*p[2][2]*p[5][0]*p[6][0] + p[1][0]*p[1][2]*p[5][1]*p[6][0] + p[1][2]*p[2][0]*p[5][1]*p[6][0] + p[2][0]*p[2][2]*p[5][1]*p[6][0] - p[4][0]*p[4][2]*p[5][1]*p[6][0] + 2*p[1][2]*p[5][0]*p[5][1]*p[6][0] + p[2][2]*p[5][0]*p[5][1]*p[6][0] - 2*p[4][2]*p[5][0]*p[5][1]*p[6][0] - p[1][0]*p[1][1]*p[5][2]*p[6][0] - p[1][1]*p[2][0]*p[5][2]*p[6][0] - p[2][0]*p[2][1]*p[5][2]*p[6][0] + p[4][0]*p[4][1]*p[5][2]*p[6][0] - 2*p[1][1]*p[5][0]*p[5][2]*p[6][0] - p[2][1]*p[5][0]*p[5][2]*p[6][0] + 2*p[4][1]*p[5][0]*p[5][2]*p[6][0] + 2*pow(p[2][0], 2)*p[1][2]*p[6][1] - 2*pow(p[5][0], 2)*p[1][2]*p[6][1] + p[1][0]*p[1][2]*p[2][0]*p[6][1] - pow(p[1][0], 2)*p[2][2]*p[6][1] + pow(p[3][0], 2)*p[2][2]*p[6][1] - pow(p[5][0], 2)*p[2][2]*p[6][1] + pow(p[7][0], 2)*p[2][2]*p[6][1] - 2*p[1][0]*p[2][0]*p[2][2]*p[6][1] + 2*p[2][0]*p[2][2]*p[3][0]*p[6][1] - 2*pow(p[2][0], 2)*p[3][2]*p[6][1] + 2*pow(p[7][0], 2)*p[3][2]*p[6][1] - p[2][0]*p[3][0]*p[3][2]*p[6][1] + 2*pow(p[5][0], 2)*p[4][2]*p[6][1] - 2*pow(p[7][0], 2)*p[4][2]*p[6][1] - p[1][0]*p[1][2]*p[5][0]*p[6][1] - p[1][0]*p[2][2]*p[5][0]*p[6][1] - p[2][0]*p[2][2]*p[5][0]*p[6][1] + p[4][0]*p[4][2]*p[5][0]*p[6][1] + pow(p[1][0], 2)*p[5][2]*p[6][1] + pow(p[2][0], 2)*p[5][2]*p[6][1] - pow(p[4][0], 2)*p[5][2]*p[6][1] - pow(p[7][0], 2)*p[5][2]*p[6][1] + p[1][0]*p[2][0]*p[5][2]*p[6][1] + 2*p[1][0]*p[5][0]*p[5][2]*p[6][1] + p[2][0]*p[5][0]*p[5][2]*p[6][1] - 2*p[4][0]*p[5][0]*p[5][2]*p[6][1] + p[1][2]*p[2][0]*p[6][0]*p[6][1] - p[1][0]*p[2][2]*p[6][0]*p[6][1] + p[2][2]*p[3][0]*p[6][0]*p[6][1] - p[2][0]*p[3][2]*p[6][0]*p[6][1] - p[1][2]*p[5][0]*p[6][0]*p[6][1] - 2*p[2][2]*p[5][0]*p[6][0]*p[6][1] + p[4][2]*p[5][0]*p[6][0]*p[6][1] + p[1][0]*p[5][2]*p[6][0]*p[6][1] + 2*p[2][0]*p[5][2]*p[6][0]*p[6][1] - p[4][0]*p[5][2]*p[6][0]*p[6][1] - 2*pow(p[2][0], 2)*p[1][1]*p[6][2] + 2*pow(p[5][0], 2)*p[1][1]*p[6][2] - p[1][0]*p[1][1]*p[2][0]*p[6][2] + pow(p[1][0], 2)*p[2][1]*p[6][2] - pow(p[3][0], 2)*p[2][1]*p[6][2] + pow(p[5][0], 2)*p[2][1]*p[6][2] - pow(p[7][0], 2)*p[2][1]*p[6][2] + 2*p[1][0]*p[2][0]*p[2][1]*p[6][2] - 2*p[2][0]*p[2][1]*p[3][0]*p[6][2] + 2*pow(p[2][0], 2)*p[3][1]*p[6][2] - 2*pow(p[7][0], 2)*p[3][1]*p[6][2] + p[2][0]*p[3][0]*p[3][1]*p[6][2] - 2*pow(p[5][0], 2)*p[4][1]*p[6][2] + 2*pow(p[7][0], 2)*p[4][1]*p[6][2] + p[1][0]*p[1][1]*p[5][0]*p[6][2] + p[1][0]*p[2][1]*p[5][0]*p[6][2] + p[2][0]*p[2][1]*p[5][0]*p[6][2] - p[4][0]*p[4][1]*p[5][0]*p[6][2] - pow(p[1][0], 2)*p[5][1]*p[6][2] - pow(p[2][0], 2)*p[5][1]*p[6][2] + pow(p[4][0], 2)*p[5][1]*p[6][2] + pow(p[7][0], 2)*p[5][1]*p[6][2] - p[1][0]*p[2][0]*p[5][1]*p[6][2] - 2*p[1][0]*p[5][0]*p[5][1]*p[6][2] - p[2][0]*p[5][0]*p[5][1]*p[6][2] + 2*p[4][0]*p[5][0]*p[5][1]*p[6][2] - p[1][1]*p[2][0]*p[6][0]*p[6][2] + p[1][0]*p[2][1]*p[6][0]*p[6][2] - p[2][1]*p[3][0]*p[6][0]*p[6][2] + p[2][0]*p[3][1]*p[6][0]*p[6][2] + p[1][1]*p[5][0]*p[6][0]*p[6][2] + 2*p[2][1]*p[5][0]*p[6][0]*p[6][2] - p[4][1]*p[5][0]*p[6][0]*p[6][2] - p[1][0]*p[5][1]*p[6][0]*p[6][2] - 2*p[2][0]*p[5][1]*p[6][0]*p[6][2] + p[4][0]*p[5][1]*p[6][0]*p[6][2] - p[2][0]*p[2][2]*p[3][1]*p[7][0] - 2*p[2][2]*p[3][0]*p[3][1]*p[7][0] + p[2][0]*p[2][1]*p[3][2]*p[7][0] - 2*p[0][1]*p[3][0]*p[3][2]*p[7][0] + 2*p[2][1]*p[3][0]*p[3][2]*p[7][0] - p[0][1]*p[3][2]*p[4][0]*p[7][0] - p[3][0]*p[3][2]*p[4][1]*p[7][0] - p[3][2]*p[4][0]*p[4][1]*p[7][0] + p[0][1]*p[3][0]*p[4][2]*p[7][0] + p[3][0]*p[3][1]*p[4][2]*p[7][0] + 2*p[0][1]*p[4][0]*p[4][2]*p[7][0] + p[3][1]*p[4][0]*p[4][2]*p[7][0] - 2*p[4][0]*p[4][2]*p[5][1]*p[7][0] - p[4][2]*p[5][0]*p[5][1]*p[7][0] + 2*p[4][0]*p[4][1]*p[5][2]*p[7][0] + p[4][1]*p[5][0]*p[5][2]*p[7][0] - p[2][2]*p[3][1]*p[6][0]*p[7][0] + p[2][1]*p[3][2]*p[6][0]*p[7][0] - p[4][2]*p[5][1]*p[6][0]*p[7][0] + p[4][1]*p[5][2]*p[6][0]*p[7][0] + p[2][0]*p[2][2]*p[6][1]*p[7][0] + p[2][2]*p[3][0]*p[6][1]*p[7][0] + p[3][0]*p[3][2]*p[6][1]*p[7][0] - p[4][0]*p[4][2]*p[6][1]*p[7][0] - p[4][0]*p[5][2]*p[6][1]*p[7][0] - p[5][0]*p[5][2]*p[6][1]*p[7][0] + 2*p[2][2]*p[6][0]*p[6][1]*p[7][0] + p[3][2]*p[6][0]*p[6][1]*p[7][0] - p[4][2]*p[6][0]*p[6][1]*p[7][0] - 2*p[5][2]*p[6][0]*p[6][1]*p[7][0] - p[2][0]*p[2][1]*p[6][2]*p[7][0] - p[2][1]*p[3][0]*p[6][2]*p[7][0] - p[3][0]*p[3][1]*p[6][2]*p[7][0] + p[4][0]*p[4][1]*p[6][2]*p[7][0] + p[4][0]*p[5][1]*p[6][2]*p[7][0] + p[5][0]*p[5][1]*p[6][2]*p[7][0] - 2*p[2][1]*p[6][0]*p[6][2]*p[7][0] - p[3][1]*p[6][0]*p[6][2]*p[7][0] + p[4][1]*p[6][0]*p[6][2]*p[7][0] + 2*p[5][1]*p[6][0]*p[6][2]*p[7][0] + 2*pow(p[3][0], 2)*p[2][2]*p[7][1] - 2*pow(p[6][0], 2)*p[2][2]*p[7][1] + p[2][0]*p[2][2]*p[3][0]*p[7][1] - pow(p[2][0], 2)*p[3][2]*p[7][1] + pow(p[4][0], 2)*p[3][2]*p[7][1] - pow(p[6][0], 2)*p[3][2]*p[7][1] - 2*p[2][0]*p[3][0]*p[3][2]*p[7][1] + p[3][0]*p[3][2]*p[4][0]*p[7][1] - pow(p[3][0], 2)*p[4][2]*p[7][1] + pow(p[5][0], 2)*p[4][2]*p[7][1] + pow(p[6][0], 2)*p[4][2]*p[7][1] - p[3][0]*p[4][0]*p[4][2]*p[7][1] + 2*p[4][0]*p[4][2]*p[5][0]*p[7][1] - 2*pow(p[4][0], 2)*p[5][2]*p[7][1] + 2*pow(p[6][0], 2)*p[5][2]*p[7][1] - p[4][0]*p[5][0]*p[5][2]*p[7][1] - p[2][0]*p[2][2]*p[6][0]*p[7][1] - p[2][0]*p[3][2]*p[6][0]*p[7][1] - p[3][0]*p[3][2]*p[6][0]*p[7][1] + p[4][0]*p[4][2]*p[6][0]*p[7][1] + p[4][2]*p[5][0]*p[6][0]*p[7][1] + p[5][0]*p[5][2]*p[6][0]*p[7][1] + pow(p[2][0], 2)*p[6][2]*p[7][1] + pow(p[3][0], 2)*p[6][2]*p[7][1] - pow(p[4][0], 2)*p[6][2]*p[7][1] - pow(p[5][0], 2)*p[6][2]*p[7][1] + p[2][0]*p[3][0]*p[6][2]*p[7][1] - p[4][0]*p[5][0]*p[6][2]*p[7][1] + 2*p[2][0]*p[6][0]*p[6][2]*p[7][1] + p[3][0]*p[6][0]*p[6][2]*p[7][1] - p[4][0]*p[6][0]*p[6][2]*p[7][1] - 2*p[5][0]*p[6][0]*p[6][2]*p[7][1] + p[2][2]*p[3][0]*p[7][0]*p[7][1] - p[2][0]*p[3][2]*p[7][0]*p[7][1] + 2*p[3][2]*p[4][0]*p[7][0]*p[7][1] - 2*p[3][0]*p[4][2]*p[7][0]*p[7][1] + p[4][2]*p[5][0]*p[7][0]*p[7][1] - p[4][0]*p[5][2]*p[7][0]*p[7][1] - p[2][2]*p[6][0]*p[7][0]*p[7][1] - 2*p[3][2]*p[6][0]*p[7][0]*p[7][1] + 2*p[4][2]*p[6][0]*p[7][0]*p[7][1] + p[5][2]*p[6][0]*p[7][0]*p[7][1] + p[2][0]*p[6][2]*p[7][0]*p[7][1] + 2*p[3][0]*p[6][2]*p[7][0]*p[7][1] - 2*p[4][0]*p[6][2]*p[7][0]*p[7][1] - p[5][0]*p[6][2]*p[7][0]*p[7][1] + p[0][2]*(2*pow(p[3][0], 2)*p[2][1] + p[2][0]*p[2][1]*p[3][0] - pow(p[2][0], 2)*p[3][1] + pow(p[4][0], 2)*p[3][1] + pow(p[7][0], 2)*p[3][1] - 2*p[2][0]*p[3][0]*p[3][1] + p[3][0]*p[3][1]*p[4][0] - pow(p[3][0], 2)*p[4][1] + pow(p[5][0], 2)*p[4][1] - pow(p[7][0], 2)*p[4][1] - p[3][0]*p[4][0]*p[4][1] + 2*p[4][0]*p[4][1]*p[5][0] + p[1][1]*(pow(p[2][0], 2) + pow(p[3][0], 2) - pow(p[4][0], 2) - pow(p[5][0], 2) + p[2][0]*p[3][0] - p[4][0]*p[5][0]) - 2*pow(p[4][0], 2)*p[5][1] - p[4][0]*p[5][0]*p[5][1] + pow(p[1][0], 2)*(-2*p[2][1] - p[3][1] + p[4][1] + 2*p[5][1]) + p[1][0]*(-(p[3][0]*p[3][1]) - p[2][0]*(p[2][1] + p[3][1]) + p[4][0]*p[4][1] + p[1][1]*(2*p[2][0] + p[3][0] - p[4][0] - 2*p[5][0]) + p[4][1]*p[5][0] + p[5][0]*p[5][1]) + 2*p[3][0]*p[3][1]*p[7][0] + p[3][1]*p[4][0]*p[7][0] - p[3][0]*p[4][1]*p[7][0] - 2*p[4][0]*p[4][1]*p[7][0] - 2*pow(p[3][0], 2)*p[7][1] + 2*pow(p[4][0], 2)*p[7][1] - p[3][0]*p[7][0]*p[7][1] + p[4][0]*p[7][0]*p[7][1]) + 2*pow(p[3][0], 2)*p[0][1]*p[7][2] - 2*pow(p[4][0], 2)*p[0][1]*p[7][2] - 2*pow(p[3][0], 2)*p[2][1]*p[7][2] + 2*pow(p[6][0], 2)*p[2][1]*p[7][2] - p[2][0]*p[2][1]*p[3][0]*p[7][2] + pow(p[2][0], 2)*p[3][1]*p[7][2] - pow(p[4][0], 2)*p[3][1]*p[7][2] + pow(p[6][0], 2)*p[3][1]*p[7][2] + 2*p[2][0]*p[3][0]*p[3][1]*p[7][2] - p[3][0]*p[3][1]*p[4][0]*p[7][2] + pow(p[3][0], 2)*p[4][1]*p[7][2] - pow(p[5][0], 2)*p[4][1]*p[7][2] - pow(p[6][0], 2)*p[4][1]*p[7][2] + p[3][0]*p[4][0]*p[4][1]*p[7][2] - 2*p[4][0]*p[4][1]*p[5][0]*p[7][2] + 2*pow(p[4][0], 2)*p[5][1]*p[7][2] - 2*pow(p[6][0], 2)*p[5][1]*p[7][2] + p[4][0]*p[5][0]*p[5][1]*p[7][2] + p[2][0]*p[2][1]*p[6][0]*p[7][2] + p[2][0]*p[3][1]*p[6][0]*p[7][2] + p[3][0]*p[3][1]*p[6][0]*p[7][2] - p[4][0]*p[4][1]*p[6][0]*p[7][2] - p[4][1]*p[5][0]*p[6][0]*p[7][2] - p[5][0]*p[5][1]*p[6][0]*p[7][2] - pow(p[2][0], 2)*p[6][1]*p[7][2] - pow(p[3][0], 2)*p[6][1]*p[7][2] + pow(p[4][0], 2)*p[6][1]*p[7][2] + pow(p[5][0], 2)*p[6][1]*p[7][2] - p[2][0]*p[3][0]*p[6][1]*p[7][2] + p[4][0]*p[5][0]*p[6][1]*p[7][2] - 2*p[2][0]*p[6][0]*p[6][1]*p[7][2] - p[3][0]*p[6][0]*p[6][1]*p[7][2] + p[4][0]*p[6][0]*p[6][1]*p[7][2] + 2*p[5][0]*p[6][0]*p[6][1]*p[7][2] + p[0][1]*p[3][0]*p[7][0]*p[7][2] - p[2][1]*p[3][0]*p[7][0]*p[7][2] + p[2][0]*p[3][1]*p[7][0]*p[7][2] - p[0][1]*p[4][0]*p[7][0]*p[7][2] - 2*p[3][1]*p[4][0]*p[7][0]*p[7][2] + 2*p[3][0]*p[4][1]*p[7][0]*p[7][2] - p[4][1]*p[5][0]*p[7][0]*p[7][2] + p[4][0]*p[5][1]*p[7][0]*p[7][2] + p[2][1]*p[6][0]*p[7][0]*p[7][2] + 2*p[3][1]*p[6][0]*p[7][0]*p[7][2] - 2*p[4][1]*p[6][0]*p[7][0]*p[7][2] - p[5][1]*p[6][0]*p[7][0]*p[7][2] - p[2][0]*p[6][1]*p[7][0]*p[7][2] - 2*p[3][0]*p[6][1]*p[7][0]*p[7][2] + 2*p[4][0]*p[6][1]*p[7][0]*p[7][2] + p[5][0]*p[6][1]*p[7][0]*p[7][2] + pow(p[0][0], 2)*(p[2][2]*p[3][1] - p[2][1]*p[3][2] + 2*p[3][2]*p[4][1] - 2*p[3][1]*p[4][2] + p[1][2]*(p[2][1] + 2*p[3][1] - 2*p[4][1] - p[5][1]) + p[4][2]*p[5][1] - p[4][1]*p[5][2] + p[1][1]*(-p[2][2] - 2*p[3][2] + 2*p[4][2] + p[5][2]) + p[3][2]*p[7][1] - p[4][2]*p[7][1] - p[3][1]*p[7][2] + p[4][1]*p[7][2]) + p[0][0]*(2*p[1][0]*p[1][2]*p[2][1] + p[1][2]*p[2][0]*p[2][1] - 2*p[1][0]*p[1][1]*p[2][2] - p[1][1]*p[2][0]*p[2][2] + p[1][2]*p[2][1]*p[3][0] - p[1][1]*p[2][2]*p[3][0] + p[1][0]*p[1][2]*p[3][1] + p[1][0]*p[2][2]*p[3][1] + p[2][0]*p[2][2]*p[3][1] + p[1][2]*p[3][0]*p[3][1] + 2*p[2][2]*p[3][0]*p[3][1] - p[1][0]*p[1][1]*p[3][2] - p[1][0]*p[2][1]*p[3][2] - p[2][0]*p[2][1]*p[3][2] - p[1][1]*p[3][0]*p[3][2] - 2*p[2][1]*p[3][0]*p[3][2] - p[1][0]*p[1][2]*p[4][1] + p[3][0]*p[3][2]*p[4][1] - p[1][2]*p[4][0]*p[4][1] + p[3][2]*p[4][0]*p[4][1] + p[1][0]*p[1][1]*p[4][2] - p[3][0]*p[3][1]*p[4][2] + p[1][1]*p[4][0]*p[4][2] - p[3][1]*p[4][0]*p[4][2] - 2*p[1][0]*p[1][2]*p[5][1] - p[1][2]*p[4][0]*p[5][1] + p[1][0]*p[4][2]*p[5][1] + 2*p[4][0]*p[4][2]*p[5][1] - p[1][2]*p[5][0]*p[5][1] + p[4][2]*p[5][0]*p[5][1] + 2*p[1][0]*p[1][1]*p[5][2] + p[1][1]*p[4][0]*p[5][2] - p[1][0]*p[4][1]*p[5][2] - 2*p[4][0]*p[4][1]*p[5][2] + p[1][1]*p[5][0]*p[5][2] - p[4][1]*p[5][0]*p[5][2] + 2*p[3][0]*p[3][2]*p[7][1] + p[3][2]*p[4][0]*p[7][1] - p[3][0]*p[4][2]*p[7][1] - 2*p[4][0]*p[4][2]*p[7][1] + p[3][2]*p[7][0]*p[7][1] - p[4][2]*p[7][0]*p[7][1] + p[0][2]*(p[2][1]*p[3][0] - p[2][0]*p[3][1] + 2*p[3][1]*p[4][0] - 2*p[3][0]*p[4][1] + p[1][1]*(p[2][0] + 2*p[3][0] - 2*p[4][0] - p[5][0]) + p[4][1]*p[5][0] - p[4][0]*p[5][1] + p[1][0]*(-p[2][1] - 2*p[3][1] + 2*p[4][1] + p[5][1]) + p[3][1]*p[7][0] - p[4][1]*p[7][0] - p[3][0]*p[7][1] + p[4][0]*p[7][1]) - 2*p[3][0]*p[3][1]*p[7][2] - p[3][1]*p[4][0]*p[7][2] + p[3][0]*p[4][1]*p[7][2] + 2*p[4][0]*p[4][1]*p[7][2] - p[3][1]*p[7][0]*p[7][2] + p[4][1]*p[7][0]*p[7][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[2][0]*p[3][2] - 2*p[3][2]*p[4][0] + 2*p[3][0]*p[4][2] - p[4][2]*p[5][0] + p[1][2]*(-p[2][0] - 2*p[3][0] + 2*p[4][0] + p[5][0]) + p[1][0]*(p[2][2] + 2*p[3][2] - 2*p[4][2] - p[5][2]) + p[4][0]*p[5][2] - p[3][2]*p[7][0] + p[4][2]*p[7][0] + p[3][0]*p[7][2] - p[4][0]*p[7][2])))/72;
380  result[2][0] = (pow(p[2][1], 2)*p[0][0]*p[1][2] + pow(p[3][1], 2)*p[0][0]*p[1][2] - pow(p[4][1], 2)*p[0][0]*p[1][2] - pow(p[5][1], 2)*p[0][0]*p[1][2] - pow(p[3][1], 2)*p[1][2]*p[2][0] + pow(p[5][1], 2)*p[1][2]*p[2][0] + pow(p[6][1], 2)*p[1][2]*p[2][0] + 2*p[0][0]*p[1][1]*p[1][2]*p[2][1] - 2*pow(p[1][1], 2)*p[0][0]*p[2][2] + 2*pow(p[3][1], 2)*p[0][0]*p[2][2] + pow(p[3][1], 2)*p[1][0]*p[2][2] - pow(p[5][1], 2)*p[1][0]*p[2][2] - pow(p[6][1], 2)*p[1][0]*p[2][2] - p[0][0]*p[1][1]*p[2][1]*p[2][2] + 2*pow(p[2][1], 2)*p[1][2]*p[3][0] + p[1][1]*p[1][2]*p[2][1]*p[3][0] - pow(p[1][1], 2)*p[2][2]*p[3][0] + pow(p[6][1], 2)*p[2][2]*p[3][0] + pow(p[7][1], 2)*p[2][2]*p[3][0] - 2*p[1][1]*p[2][1]*p[2][2]*p[3][0] + p[0][0]*p[1][1]*p[1][2]*p[3][1] - p[1][1]*p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[1][2]*p[2][1]*p[3][1] - 2*p[1][2]*p[2][0]*p[2][1]*p[3][1] + p[1][0]*p[1][1]*p[2][2]*p[3][1] + p[0][0]*p[2][1]*p[2][2]*p[3][1] + 2*p[1][0]*p[2][1]*p[2][2]*p[3][1] + p[1][2]*p[2][1]*p[3][0]*p[3][1] - p[1][1]*p[2][2]*p[3][0]*p[3][1] - pow(p[1][1], 2)*p[0][0]*p[3][2] - pow(p[2][1], 2)*p[0][0]*p[3][2] + pow(p[4][1], 2)*p[0][0]*p[3][2] + pow(p[7][1], 2)*p[0][0]*p[3][2] - 2*pow(p[2][1], 2)*p[1][0]*p[3][2] + pow(p[1][1], 2)*p[2][0]*p[3][2] - pow(p[6][1], 2)*p[2][0]*p[3][2] - pow(p[7][1], 2)*p[2][0]*p[3][2] - p[0][0]*p[1][1]*p[2][1]*p[3][2] - p[1][0]*p[1][1]*p[2][1]*p[3][2] + 2*p[1][1]*p[2][0]*p[2][1]*p[3][2] - p[0][0]*p[1][1]*p[3][1]*p[3][2] + p[1][1]*p[2][0]*p[3][1]*p[3][2] - 2*p[0][0]*p[2][1]*p[3][1]*p[3][2] - p[1][0]*p[2][1]*p[3][1]*p[3][2] - 2*pow(p[5][1], 2)*p[1][2]*p[4][0] + 2*pow(p[7][1], 2)*p[3][2]*p[4][0] - p[0][0]*p[1][1]*p[1][2]*p[4][1] + p[0][0]*p[3][1]*p[3][2]*p[4][1] + pow(p[1][1], 2)*p[0][0]*p[4][2] - pow(p[3][1], 2)*p[0][0]*p[4][2] + pow(p[5][1], 2)*p[0][0]*p[4][2] - pow(p[7][1], 2)*p[0][0]*p[4][2] + 2*pow(p[5][1], 2)*p[1][0]*p[4][2] - 2*pow(p[7][1], 2)*p[3][0]*p[4][2] + p[0][0]*p[1][1]*p[4][1]*p[4][2] - p[0][0]*p[3][1]*p[4][1]*p[4][2] - pow(p[2][1], 2)*p[1][2]*p[5][0] + pow(p[4][1], 2)*p[1][2]*p[5][0] - pow(p[6][1], 2)*p[1][2]*p[5][0] - 2*p[1][1]*p[1][2]*p[2][1]*p[5][0] + 2*pow(p[1][1], 2)*p[2][2]*p[5][0] - 2*pow(p[6][1], 2)*p[2][2]*p[5][0] + p[1][1]*p[2][1]*p[2][2]*p[5][0] + p[1][1]*p[1][2]*p[4][1]*p[5][0] - pow(p[1][1], 2)*p[4][2]*p[5][0] + pow(p[6][1], 2)*p[4][2]*p[5][0] + pow(p[7][1], 2)*p[4][2]*p[5][0] - p[1][1]*p[4][1]*p[4][2]*p[5][0] - 2*p[0][0]*p[1][1]*p[1][2]*p[5][1] + 2*p[1][1]*p[1][2]*p[2][0]*p[5][1] + p[1][2]*p[2][0]*p[2][1]*p[5][1] - 2*p[1][0]*p[1][1]*p[2][2]*p[5][1] - p[1][0]*p[2][1]*p[2][2]*p[5][1] - p[1][1]*p[1][2]*p[4][0]*p[5][1] - p[0][0]*p[1][2]*p[4][1]*p[5][1] - p[1][2]*p[4][0]*p[4][1]*p[5][1] + p[0][0]*p[1][1]*p[4][2]*p[5][1] + p[1][0]*p[1][1]*p[4][2]*p[5][1] + 2*p[0][0]*p[4][1]*p[4][2]*p[5][1] + p[1][0]*p[4][1]*p[4][2]*p[5][1] - p[1][2]*p[2][1]*p[5][0]*p[5][1] + p[1][1]*p[2][2]*p[5][0]*p[5][1] + 2*p[1][2]*p[4][1]*p[5][0]*p[5][1] - 2*p[1][1]*p[4][2]*p[5][0]*p[5][1] + 2*pow(p[1][1], 2)*p[0][0]*p[5][2] - 2*pow(p[4][1], 2)*p[0][0]*p[5][2] + pow(p[2][1], 2)*p[1][0]*p[5][2] - pow(p[4][1], 2)*p[1][0]*p[5][2] + pow(p[6][1], 2)*p[1][0]*p[5][2] - 2*pow(p[1][1], 2)*p[2][0]*p[5][2] + 2*pow(p[6][1], 2)*p[2][0]*p[5][2] + 2*p[1][0]*p[1][1]*p[2][1]*p[5][2] - p[1][1]*p[2][0]*p[2][1]*p[5][2] + pow(p[1][1], 2)*p[4][0]*p[5][2] - pow(p[6][1], 2)*p[4][0]*p[5][2] - pow(p[7][1], 2)*p[4][0]*p[5][2] - p[1][0]*p[1][1]*p[4][1]*p[5][2] + p[1][1]*p[4][0]*p[4][1]*p[5][2] + p[0][0]*p[1][1]*p[5][1]*p[5][2] - p[1][1]*p[2][0]*p[5][1]*p[5][2] + p[1][0]*p[2][1]*p[5][1]*p[5][2] + 2*p[1][1]*p[4][0]*p[5][1]*p[5][2] - p[0][0]*p[4][1]*p[5][1]*p[5][2] - 2*p[1][0]*p[4][1]*p[5][1]*p[5][2] - 2*pow(p[2][1], 2)*p[1][2]*p[6][0] + 2*pow(p[5][1], 2)*p[1][2]*p[6][0] - p[1][1]*p[1][2]*p[2][1]*p[6][0] + pow(p[1][1], 2)*p[2][2]*p[6][0] - pow(p[3][1], 2)*p[2][2]*p[6][0] + pow(p[5][1], 2)*p[2][2]*p[6][0] - pow(p[7][1], 2)*p[2][2]*p[6][0] + 2*p[1][1]*p[2][1]*p[2][2]*p[6][0] - 2*p[2][1]*p[2][2]*p[3][1]*p[6][0] + 2*pow(p[2][1], 2)*p[3][2]*p[6][0] - 2*pow(p[7][1], 2)*p[3][2]*p[6][0] + p[2][1]*p[3][1]*p[3][2]*p[6][0] - 2*pow(p[5][1], 2)*p[4][2]*p[6][0] + 2*pow(p[7][1], 2)*p[4][2]*p[6][0] + p[1][1]*p[1][2]*p[5][1]*p[6][0] + p[1][1]*p[2][2]*p[5][1]*p[6][0] + p[2][1]*p[2][2]*p[5][1]*p[6][0] - p[4][1]*p[4][2]*p[5][1]*p[6][0] - pow(p[1][1], 2)*p[5][2]*p[6][0] - pow(p[2][1], 2)*p[5][2]*p[6][0] + pow(p[4][1], 2)*p[5][2]*p[6][0] + pow(p[7][1], 2)*p[5][2]*p[6][0] - p[1][1]*p[2][1]*p[5][2]*p[6][0] - 2*p[1][1]*p[5][1]*p[5][2]*p[6][0] - p[2][1]*p[5][1]*p[5][2]*p[6][0] + 2*p[4][1]*p[5][1]*p[5][2]*p[6][0] + p[1][1]*p[1][2]*p[2][0]*p[6][1] + 2*p[1][2]*p[2][0]*p[2][1]*p[6][1] - p[1][0]*p[1][1]*p[2][2]*p[6][1] - 2*p[1][0]*p[2][1]*p[2][2]*p[6][1] + 2*p[2][1]*p[2][2]*p[3][0]*p[6][1] + p[2][2]*p[3][0]*p[3][1]*p[6][1] - 2*p[2][0]*p[2][1]*p[3][2]*p[6][1] - p[2][0]*p[3][1]*p[3][2]*p[6][1] - p[1][1]*p[1][2]*p[5][0]*p[6][1] - p[1][2]*p[2][1]*p[5][0]*p[6][1] - p[2][1]*p[2][2]*p[5][0]*p[6][1] + p[4][1]*p[4][2]*p[5][0]*p[6][1] + p[1][2]*p[2][0]*p[5][1]*p[6][1] - p[1][0]*p[2][2]*p[5][1]*p[6][1] - 2*p[1][2]*p[5][0]*p[5][1]*p[6][1] - p[2][2]*p[5][0]*p[5][1]*p[6][1] + 2*p[4][2]*p[5][0]*p[5][1]*p[6][1] + p[1][0]*p[1][1]*p[5][2]*p[6][1] + p[1][0]*p[2][1]*p[5][2]*p[6][1] + p[2][0]*p[2][1]*p[5][2]*p[6][1] - p[4][0]*p[4][1]*p[5][2]*p[6][1] + 2*p[1][0]*p[5][1]*p[5][2]*p[6][1] + p[2][0]*p[5][1]*p[5][2]*p[6][1] - 2*p[4][0]*p[5][1]*p[5][2]*p[6][1] - p[1][2]*p[2][1]*p[6][0]*p[6][1] + p[1][1]*p[2][2]*p[6][0]*p[6][1] - p[2][2]*p[3][1]*p[6][0]*p[6][1] + p[2][1]*p[3][2]*p[6][0]*p[6][1] + p[1][2]*p[5][1]*p[6][0]*p[6][1] + 2*p[2][2]*p[5][1]*p[6][0]*p[6][1] - p[4][2]*p[5][1]*p[6][0]*p[6][1] - p[1][1]*p[5][2]*p[6][0]*p[6][1] - 2*p[2][1]*p[5][2]*p[6][0]*p[6][1] + p[4][1]*p[5][2]*p[6][0]*p[6][1] + 2*pow(p[2][1], 2)*p[1][0]*p[6][2] - 2*pow(p[5][1], 2)*p[1][0]*p[6][2] - pow(p[1][1], 2)*p[2][0]*p[6][2] + pow(p[3][1], 2)*p[2][0]*p[6][2] - pow(p[5][1], 2)*p[2][0]*p[6][2] + pow(p[7][1], 2)*p[2][0]*p[6][2] + p[1][0]*p[1][1]*p[2][1]*p[6][2] - 2*p[1][1]*p[2][0]*p[2][1]*p[6][2] - 2*pow(p[2][1], 2)*p[3][0]*p[6][2] + 2*pow(p[7][1], 2)*p[3][0]*p[6][2] + 2*p[2][0]*p[2][1]*p[3][1]*p[6][2] - p[2][1]*p[3][0]*p[3][1]*p[6][2] + 2*pow(p[5][1], 2)*p[4][0]*p[6][2] - 2*pow(p[7][1], 2)*p[4][0]*p[6][2] + pow(p[1][1], 2)*p[5][0]*p[6][2] + pow(p[2][1], 2)*p[5][0]*p[6][2] - pow(p[4][1], 2)*p[5][0]*p[6][2] - pow(p[7][1], 2)*p[5][0]*p[6][2] + p[1][1]*p[2][1]*p[5][0]*p[6][2] - p[1][0]*p[1][1]*p[5][1]*p[6][2] - p[1][1]*p[2][0]*p[5][1]*p[6][2] - p[2][0]*p[2][1]*p[5][1]*p[6][2] + p[4][0]*p[4][1]*p[5][1]*p[6][2] + 2*p[1][1]*p[5][0]*p[5][1]*p[6][2] + p[2][1]*p[5][0]*p[5][1]*p[6][2] - 2*p[4][1]*p[5][0]*p[5][1]*p[6][2] - p[1][1]*p[2][0]*p[6][1]*p[6][2] + p[1][0]*p[2][1]*p[6][1]*p[6][2] - p[2][1]*p[3][0]*p[6][1]*p[6][2] + p[2][0]*p[3][1]*p[6][1]*p[6][2] + p[1][1]*p[5][0]*p[6][1]*p[6][2] + 2*p[2][1]*p[5][0]*p[6][1]*p[6][2] - p[4][1]*p[5][0]*p[6][1]*p[6][2] - p[1][0]*p[5][1]*p[6][1]*p[6][2] - 2*p[2][0]*p[5][1]*p[6][1]*p[6][2] + p[4][0]*p[5][1]*p[6][1]*p[6][2] - 2*pow(p[3][1], 2)*p[2][2]*p[7][0] + 2*pow(p[6][1], 2)*p[2][2]*p[7][0] - p[2][1]*p[2][2]*p[3][1]*p[7][0] + pow(p[2][1], 2)*p[3][2]*p[7][0] - pow(p[4][1], 2)*p[3][2]*p[7][0] + pow(p[6][1], 2)*p[3][2]*p[7][0] + 2*p[2][1]*p[3][1]*p[3][2]*p[7][0] - p[3][1]*p[3][2]*p[4][1]*p[7][0] + pow(p[3][1], 2)*p[4][2]*p[7][0] - pow(p[5][1], 2)*p[4][2]*p[7][0] - pow(p[6][1], 2)*p[4][2]*p[7][0] + p[3][1]*p[4][1]*p[4][2]*p[7][0] - 2*p[4][1]*p[4][2]*p[5][1]*p[7][0] + 2*pow(p[4][1], 2)*p[5][2]*p[7][0] - 2*pow(p[6][1], 2)*p[5][2]*p[7][0] + p[4][1]*p[5][1]*p[5][2]*p[7][0] + p[2][1]*p[2][2]*p[6][1]*p[7][0] + p[2][1]*p[3][2]*p[6][1]*p[7][0] + p[3][1]*p[3][2]*p[6][1]*p[7][0] - p[4][1]*p[4][2]*p[6][1]*p[7][0] - p[4][2]*p[5][1]*p[6][1]*p[7][0] - p[5][1]*p[5][2]*p[6][1]*p[7][0] - pow(p[2][1], 2)*p[6][2]*p[7][0] - pow(p[3][1], 2)*p[6][2]*p[7][0] + pow(p[4][1], 2)*p[6][2]*p[7][0] + pow(p[5][1], 2)*p[6][2]*p[7][0] - p[2][1]*p[3][1]*p[6][2]*p[7][0] + p[4][1]*p[5][1]*p[6][2]*p[7][0] - 2*p[2][1]*p[6][1]*p[6][2]*p[7][0] - p[3][1]*p[6][1]*p[6][2]*p[7][0] + p[4][1]*p[6][1]*p[6][2]*p[7][0] + 2*p[5][1]*p[6][1]*p[6][2]*p[7][0] + p[2][1]*p[2][2]*p[3][0]*p[7][1] + 2*p[2][2]*p[3][0]*p[3][1]*p[7][1] - p[2][0]*p[2][1]*p[3][2]*p[7][1] + 2*p[0][0]*p[3][1]*p[3][2]*p[7][1] - 2*p[2][0]*p[3][1]*p[3][2]*p[7][1] + p[3][1]*p[3][2]*p[4][0]*p[7][1] + p[0][0]*p[3][2]*p[4][1]*p[7][1] + p[3][2]*p[4][0]*p[4][1]*p[7][1] - p[0][0]*p[3][1]*p[4][2]*p[7][1] - p[3][0]*p[3][1]*p[4][2]*p[7][1] - 2*p[0][0]*p[4][1]*p[4][2]*p[7][1] - p[3][0]*p[4][1]*p[4][2]*p[7][1] + 2*p[4][1]*p[4][2]*p[5][0]*p[7][1] + p[4][2]*p[5][0]*p[5][1]*p[7][1] - 2*p[4][0]*p[4][1]*p[5][2]*p[7][1] - p[4][0]*p[5][1]*p[5][2]*p[7][1] - p[2][1]*p[2][2]*p[6][0]*p[7][1] - p[2][2]*p[3][1]*p[6][0]*p[7][1] - p[3][1]*p[3][2]*p[6][0]*p[7][1] + p[4][1]*p[4][2]*p[6][0]*p[7][1] + p[4][1]*p[5][2]*p[6][0]*p[7][1] + p[5][1]*p[5][2]*p[6][0]*p[7][1] + p[2][2]*p[3][0]*p[6][1]*p[7][1] - p[2][0]*p[3][2]*p[6][1]*p[7][1] + p[4][2]*p[5][0]*p[6][1]*p[7][1] - p[4][0]*p[5][2]*p[6][1]*p[7][1] - 2*p[2][2]*p[6][0]*p[6][1]*p[7][1] - p[3][2]*p[6][0]*p[6][1]*p[7][1] + p[4][2]*p[6][0]*p[6][1]*p[7][1] + 2*p[5][2]*p[6][0]*p[6][1]*p[7][1] + p[2][0]*p[2][1]*p[6][2]*p[7][1] + p[2][0]*p[3][1]*p[6][2]*p[7][1] + p[3][0]*p[3][1]*p[6][2]*p[7][1] - p[4][0]*p[4][1]*p[6][2]*p[7][1] - p[4][1]*p[5][0]*p[6][2]*p[7][1] - p[5][0]*p[5][1]*p[6][2]*p[7][1] + 2*p[2][0]*p[6][1]*p[6][2]*p[7][1] + p[3][0]*p[6][1]*p[6][2]*p[7][1] - p[4][0]*p[6][1]*p[6][2]*p[7][1] - 2*p[5][0]*p[6][1]*p[6][2]*p[7][1] - p[2][2]*p[3][1]*p[7][0]*p[7][1] + p[2][1]*p[3][2]*p[7][0]*p[7][1] - 2*p[3][2]*p[4][1]*p[7][0]*p[7][1] + 2*p[3][1]*p[4][2]*p[7][0]*p[7][1] - p[4][2]*p[5][1]*p[7][0]*p[7][1] + p[4][1]*p[5][2]*p[7][0]*p[7][1] + p[2][2]*p[6][1]*p[7][0]*p[7][1] + 2*p[3][2]*p[6][1]*p[7][0]*p[7][1] - 2*p[4][2]*p[6][1]*p[7][0]*p[7][1] - p[5][2]*p[6][1]*p[7][0]*p[7][1] - p[2][1]*p[6][2]*p[7][0]*p[7][1] - 2*p[3][1]*p[6][2]*p[7][0]*p[7][1] + 2*p[4][1]*p[6][2]*p[7][0]*p[7][1] + p[5][1]*p[6][2]*p[7][0]*p[7][1] + p[0][2]*(-2*pow(p[3][1], 2)*p[2][0] + pow(p[2][1], 2)*p[3][0] - pow(p[4][1], 2)*p[3][0] - pow(p[7][1], 2)*p[3][0] - p[2][0]*p[2][1]*p[3][1] + 2*p[2][1]*p[3][0]*p[3][1] + pow(p[3][1], 2)*p[4][0] - pow(p[5][1], 2)*p[4][0] + pow(p[7][1], 2)*p[4][0] - p[3][0]*p[3][1]*p[4][1] + p[3][1]*p[4][0]*p[4][1] + pow(p[1][1], 2)*(2*p[2][0] + p[3][0] - p[4][0] - 2*p[5][0]) + 2*pow(p[4][1], 2)*p[5][0] - 2*p[4][0]*p[4][1]*p[5][1] + p[4][1]*p[5][0]*p[5][1] + p[1][0]*(-pow(p[2][1], 2) - pow(p[3][1], 2) + pow(p[4][1], 2) + pow(p[5][1], 2) - p[2][1]*p[3][1] + p[4][1]*p[5][1]) + p[1][1]*(p[2][0]*p[2][1] + p[2][1]*p[3][0] + p[3][0]*p[3][1] - p[4][0]*p[4][1] - p[4][0]*p[5][1] - p[5][0]*p[5][1] + p[1][0]*(-2*p[2][1] - p[3][1] + p[4][1] + 2*p[5][1])) + 2*pow(p[3][1], 2)*p[7][0] - 2*pow(p[4][1], 2)*p[7][0] - 2*p[3][0]*p[3][1]*p[7][1] + p[3][1]*p[4][0]*p[7][1] - p[3][0]*p[4][1]*p[7][1] + 2*p[4][0]*p[4][1]*p[7][1] + p[3][1]*p[7][0]*p[7][1] - p[4][1]*p[7][0]*p[7][1]) - 2*pow(p[3][1], 2)*p[0][0]*p[7][2] + 2*pow(p[4][1], 2)*p[0][0]*p[7][2] + 2*pow(p[3][1], 2)*p[2][0]*p[7][2] - 2*pow(p[6][1], 2)*p[2][0]*p[7][2] - pow(p[2][1], 2)*p[3][0]*p[7][2] + pow(p[4][1], 2)*p[3][0]*p[7][2] - pow(p[6][1], 2)*p[3][0]*p[7][2] + p[2][0]*p[2][1]*p[3][1]*p[7][2] - 2*p[2][1]*p[3][0]*p[3][1]*p[7][2] - pow(p[3][1], 2)*p[4][0]*p[7][2] + pow(p[5][1], 2)*p[4][0]*p[7][2] + pow(p[6][1], 2)*p[4][0]*p[7][2] + p[3][0]*p[3][1]*p[4][1]*p[7][2] - p[3][1]*p[4][0]*p[4][1]*p[7][2] - 2*pow(p[4][1], 2)*p[5][0]*p[7][2] + 2*pow(p[6][1], 2)*p[5][0]*p[7][2] + 2*p[4][0]*p[4][1]*p[5][1]*p[7][2] - p[4][1]*p[5][0]*p[5][1]*p[7][2] + pow(p[2][1], 2)*p[6][0]*p[7][2] + pow(p[3][1], 2)*p[6][0]*p[7][2] - pow(p[4][1], 2)*p[6][0]*p[7][2] - pow(p[5][1], 2)*p[6][0]*p[7][2] + p[2][1]*p[3][1]*p[6][0]*p[7][2] - p[4][1]*p[5][1]*p[6][0]*p[7][2] - p[2][0]*p[2][1]*p[6][1]*p[7][2] - p[2][1]*p[3][0]*p[6][1]*p[7][2] - p[3][0]*p[3][1]*p[6][1]*p[7][2] + p[4][0]*p[4][1]*p[6][1]*p[7][2] + p[4][0]*p[5][1]*p[6][1]*p[7][2] + p[5][0]*p[5][1]*p[6][1]*p[7][2] + 2*p[2][1]*p[6][0]*p[6][1]*p[7][2] + p[3][1]*p[6][0]*p[6][1]*p[7][2] - p[4][1]*p[6][0]*p[6][1]*p[7][2] - 2*p[5][1]*p[6][0]*p[6][1]*p[7][2] - p[2][1]*p[3][0]*p[7][1]*p[7][2] - p[0][0]*p[3][1]*p[7][1]*p[7][2] + p[2][0]*p[3][1]*p[7][1]*p[7][2] - 2*p[3][1]*p[4][0]*p[7][1]*p[7][2] + p[0][0]*p[4][1]*p[7][1]*p[7][2] + 2*p[3][0]*p[4][1]*p[7][1]*p[7][2] - p[4][1]*p[5][0]*p[7][1]*p[7][2] + p[4][0]*p[5][1]*p[7][1]*p[7][2] + p[2][1]*p[6][0]*p[7][1]*p[7][2] + 2*p[3][1]*p[6][0]*p[7][1]*p[7][2] - 2*p[4][1]*p[6][0]*p[7][1]*p[7][2] - p[5][1]*p[6][0]*p[7][1]*p[7][2] - p[2][0]*p[6][1]*p[7][1]*p[7][2] - 2*p[3][0]*p[6][1]*p[7][1]*p[7][2] + 2*p[4][0]*p[6][1]*p[7][1]*p[7][2] + p[5][0]*p[6][1]*p[7][1]*p[7][2] + pow(p[0][1], 2)*(-(p[2][2]*p[3][0]) + p[2][0]*p[3][2] - 2*p[3][2]*p[4][0] + 2*p[3][0]*p[4][2] - p[4][2]*p[5][0] + p[1][2]*(-p[2][0] - 2*p[3][0] + 2*p[4][0] + p[5][0]) + p[1][0]*(p[2][2] + 2*p[3][2] - 2*p[4][2] - p[5][2]) + p[4][0]*p[5][2] - p[3][2]*p[7][0] + p[4][2]*p[7][0] + p[3][0]*p[7][2] - p[4][0]*p[7][2]) + p[0][1]*(p[0][0]*p[1][2]*p[2][1] - p[1][2]*p[2][0]*p[2][1] + p[1][0]*p[2][1]*p[2][2] - p[2][1]*p[2][2]*p[3][0] + 2*p[0][0]*p[1][2]*p[3][1] - p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] + p[1][0]*p[2][2]*p[3][1] - p[1][2]*p[3][0]*p[3][1] - 2*p[2][2]*p[3][0]*p[3][1] - p[0][0]*p[2][1]*p[3][2] + p[2][0]*p[2][1]*p[3][2] + p[1][0]*p[3][1]*p[3][2] + 2*p[2][0]*p[3][1]*p[3][2] - p[3][1]*p[3][2]*p[4][0] - 2*p[0][0]*p[1][2]*p[4][1] + 2*p[0][0]*p[3][2]*p[4][1] + p[1][2]*p[4][0]*p[4][1] - p[3][2]*p[4][0]*p[4][1] - 2*p[0][0]*p[3][1]*p[4][2] + p[3][0]*p[3][1]*p[4][2] - p[1][0]*p[4][1]*p[4][2] + p[3][0]*p[4][1]*p[4][2] + p[1][2]*p[4][1]*p[5][0] - 2*p[4][1]*p[4][2]*p[5][0] - p[0][0]*p[1][2]*p[5][1] + p[0][0]*p[4][2]*p[5][1] + p[1][2]*p[5][0]*p[5][1] - p[4][2]*p[5][0]*p[5][1] - p[0][0]*p[4][1]*p[5][2] - p[1][0]*p[4][1]*p[5][2] + 2*p[4][0]*p[4][1]*p[5][2] - p[1][0]*p[5][1]*p[5][2] + p[4][0]*p[5][1]*p[5][2] + p[1][1]*(2*p[1][0]*p[2][2] - p[2][2]*p[3][0] + p[1][0]*p[3][2] + p[2][0]*p[3][2] - p[1][0]*p[4][2] - p[4][2]*p[5][0] + p[1][2]*(-2*p[2][0] - p[3][0] + p[4][0] + 2*p[5][0]) - 2*p[1][0]*p[5][2] + p[4][0]*p[5][2] + p[0][0]*(-p[2][2] - 2*p[3][2] + 2*p[4][2] + p[5][2])) - 2*p[3][1]*p[3][2]*p[7][0] - p[3][2]*p[4][1]*p[7][0] + p[3][1]*p[4][2]*p[7][0] + 2*p[4][1]*p[4][2]*p[7][0] + p[0][0]*p[3][2]*p[7][1] - p[0][0]*p[4][2]*p[7][1] - p[3][2]*p[7][0]*p[7][1] + p[4][2]*p[7][0]*p[7][1] + p[0][2]*(p[2][1]*p[3][0] - p[2][0]*p[3][1] + 2*p[3][1]*p[4][0] - 2*p[3][0]*p[4][1] + p[1][1]*(p[2][0] + 2*p[3][0] - 2*p[4][0] - p[5][0]) + p[4][1]*p[5][0] - p[4][0]*p[5][1] + p[1][0]*(-p[2][1] - 2*p[3][1] + 2*p[4][1] + p[5][1]) + p[3][1]*p[7][0] - p[4][1]*p[7][0] - p[3][0]*p[7][1] + p[4][0]*p[7][1]) - p[0][0]*p[3][1]*p[7][2] + 2*p[3][0]*p[3][1]*p[7][2] - p[3][1]*p[4][0]*p[7][2] + p[0][0]*p[4][1]*p[7][2] + p[3][0]*p[4][1]*p[7][2] - 2*p[4][0]*p[4][1]*p[7][2] + p[3][0]*p[7][1]*p[7][2] - p[4][0]*p[7][1]*p[7][2]))/72;
381  result[3][0] = (-(pow(p[2][2], 2)*p[0][0]*p[1][1]) - pow(p[3][2], 2)*p[0][0]*p[1][1] + pow(p[4][2], 2)*p[0][0]*p[1][1] + pow(p[5][2], 2)*p[0][0]*p[1][1] + pow(p[3][2], 2)*p[1][1]*p[2][0] - pow(p[5][2], 2)*p[1][1]*p[2][0] - pow(p[6][2], 2)*p[1][1]*p[2][0] + 2*pow(p[1][2], 2)*p[0][0]*p[2][1] - 2*pow(p[3][2], 2)*p[0][0]*p[2][1] - pow(p[3][2], 2)*p[1][0]*p[2][1] + pow(p[5][2], 2)*p[1][0]*p[2][1] + pow(p[6][2], 2)*p[1][0]*p[2][1] - 2*p[0][0]*p[1][1]*p[1][2]*p[2][2] + p[0][0]*p[1][2]*p[2][1]*p[2][2] - 2*pow(p[2][2], 2)*p[1][1]*p[3][0] + pow(p[1][2], 2)*p[2][1]*p[3][0] - pow(p[6][2], 2)*p[2][1]*p[3][0] - pow(p[7][2], 2)*p[2][1]*p[3][0] - p[1][1]*p[1][2]*p[2][2]*p[3][0] + 2*p[1][2]*p[2][1]*p[2][2]*p[3][0] + pow(p[1][2], 2)*p[0][0]*p[3][1] + pow(p[2][2], 2)*p[0][0]*p[3][1] - pow(p[4][2], 2)*p[0][0]*p[3][1] - pow(p[7][2], 2)*p[0][0]*p[3][1] + 2*pow(p[2][2], 2)*p[1][0]*p[3][1] - pow(p[1][2], 2)*p[2][0]*p[3][1] + pow(p[6][2], 2)*p[2][0]*p[3][1] + pow(p[7][2], 2)*p[2][0]*p[3][1] + p[0][0]*p[1][2]*p[2][2]*p[3][1] + p[1][0]*p[1][2]*p[2][2]*p[3][1] - 2*p[1][2]*p[2][0]*p[2][2]*p[3][1] - p[0][0]*p[1][1]*p[1][2]*p[3][2] + p[1][1]*p[1][2]*p[2][0]*p[3][2] - p[1][0]*p[1][2]*p[2][1]*p[3][2] - p[0][0]*p[1][1]*p[2][2]*p[3][2] + 2*p[1][1]*p[2][0]*p[2][2]*p[3][2] - p[0][0]*p[2][1]*p[2][2]*p[3][2] - 2*p[1][0]*p[2][1]*p[2][2]*p[3][2] + p[1][2]*p[2][1]*p[3][0]*p[3][2] - p[1][1]*p[2][2]*p[3][0]*p[3][2] + p[0][0]*p[1][2]*p[3][1]*p[3][2] - p[1][2]*p[2][0]*p[3][1]*p[3][2] + 2*p[0][0]*p[2][2]*p[3][1]*p[3][2] + p[1][0]*p[2][2]*p[3][1]*p[3][2] + 2*pow(p[5][2], 2)*p[1][1]*p[4][0] - 2*pow(p[7][2], 2)*p[3][1]*p[4][0] - pow(p[1][2], 2)*p[0][0]*p[4][1] + pow(p[3][2], 2)*p[0][0]*p[4][1] - pow(p[5][2], 2)*p[0][0]*p[4][1] + pow(p[7][2], 2)*p[0][0]*p[4][1] - 2*pow(p[5][2], 2)*p[1][0]*p[4][1] + 2*pow(p[7][2], 2)*p[3][0]*p[4][1] + p[0][0]*p[1][1]*p[1][2]*p[4][2] - p[0][0]*p[3][1]*p[3][2]*p[4][2] - p[0][0]*p[1][2]*p[4][1]*p[4][2] + p[0][0]*p[3][2]*p[4][1]*p[4][2] + pow(p[2][2], 2)*p[1][1]*p[5][0] - pow(p[4][2], 2)*p[1][1]*p[5][0] + pow(p[6][2], 2)*p[1][1]*p[5][0] - 2*pow(p[1][2], 2)*p[2][1]*p[5][0] + 2*pow(p[6][2], 2)*p[2][1]*p[5][0] + 2*p[1][1]*p[1][2]*p[2][2]*p[5][0] - p[1][2]*p[2][1]*p[2][2]*p[5][0] + pow(p[1][2], 2)*p[4][1]*p[5][0] - pow(p[6][2], 2)*p[4][1]*p[5][0] - pow(p[7][2], 2)*p[4][1]*p[5][0] - p[1][1]*p[1][2]*p[4][2]*p[5][0] + p[1][2]*p[4][1]*p[4][2]*p[5][0] - 2*pow(p[1][2], 2)*p[0][0]*p[5][1] + 2*pow(p[4][2], 2)*p[0][0]*p[5][1] - pow(p[2][2], 2)*p[1][0]*p[5][1] + pow(p[4][2], 2)*p[1][0]*p[5][1] - pow(p[6][2], 2)*p[1][0]*p[5][1] + 2*pow(p[1][2], 2)*p[2][0]*p[5][1] - 2*pow(p[6][2], 2)*p[2][0]*p[5][1] - 2*p[1][0]*p[1][2]*p[2][2]*p[5][1] + p[1][2]*p[2][0]*p[2][2]*p[5][1] - pow(p[1][2], 2)*p[4][0]*p[5][1] + pow(p[6][2], 2)*p[4][0]*p[5][1] + pow(p[7][2], 2)*p[4][0]*p[5][1] + p[1][0]*p[1][2]*p[4][2]*p[5][1] - p[1][2]*p[4][0]*p[4][2]*p[5][1] + 2*p[0][0]*p[1][1]*p[1][2]*p[5][2] - 2*p[1][1]*p[1][2]*p[2][0]*p[5][2] + 2*p[1][0]*p[1][2]*p[2][1]*p[5][2] - p[1][1]*p[2][0]*p[2][2]*p[5][2] + p[1][0]*p[2][1]*p[2][2]*p[5][2] + p[1][1]*p[1][2]*p[4][0]*p[5][2] - p[0][0]*p[1][2]*p[4][1]*p[5][2] - p[1][0]*p[1][2]*p[4][1]*p[5][2] + p[0][0]*p[1][1]*p[4][2]*p[5][2] + p[1][1]*p[4][0]*p[4][2]*p[5][2] - 2*p[0][0]*p[4][1]*p[4][2]*p[5][2] - p[1][0]*p[4][1]*p[4][2]*p[5][2] - p[1][2]*p[2][1]*p[5][0]*p[5][2] + p[1][1]*p[2][2]*p[5][0]*p[5][2] + 2*p[1][2]*p[4][1]*p[5][0]*p[5][2] - 2*p[1][1]*p[4][2]*p[5][0]*p[5][2] - p[0][0]*p[1][2]*p[5][1]*p[5][2] + p[1][2]*p[2][0]*p[5][1]*p[5][2] - p[1][0]*p[2][2]*p[5][1]*p[5][2] - 2*p[1][2]*p[4][0]*p[5][1]*p[5][2] + p[0][0]*p[4][2]*p[5][1]*p[5][2] + 2*p[1][0]*p[4][2]*p[5][1]*p[5][2] + 2*pow(p[2][2], 2)*p[1][1]*p[6][0] - 2*pow(p[5][2], 2)*p[1][1]*p[6][0] - pow(p[1][2], 2)*p[2][1]*p[6][0] + pow(p[3][2], 2)*p[2][1]*p[6][0] - pow(p[5][2], 2)*p[2][1]*p[6][0] + pow(p[7][2], 2)*p[2][1]*p[6][0] + p[1][1]*p[1][2]*p[2][2]*p[6][0] - 2*p[1][2]*p[2][1]*p[2][2]*p[6][0] - 2*pow(p[2][2], 2)*p[3][1]*p[6][0] + 2*pow(p[7][2], 2)*p[3][1]*p[6][0] + 2*p[2][1]*p[2][2]*p[3][2]*p[6][0] - p[2][2]*p[3][1]*p[3][2]*p[6][0] + 2*pow(p[5][2], 2)*p[4][1]*p[6][0] - 2*pow(p[7][2], 2)*p[4][1]*p[6][0] + pow(p[1][2], 2)*p[5][1]*p[6][0] + pow(p[2][2], 2)*p[5][1]*p[6][0] - pow(p[4][2], 2)*p[5][1]*p[6][0] - pow(p[7][2], 2)*p[5][1]*p[6][0] + p[1][2]*p[2][2]*p[5][1]*p[6][0] - p[1][1]*p[1][2]*p[5][2]*p[6][0] - p[1][2]*p[2][1]*p[5][2]*p[6][0] - p[2][1]*p[2][2]*p[5][2]*p[6][0] + p[4][1]*p[4][2]*p[5][2]*p[6][0] + 2*p[1][2]*p[5][1]*p[5][2]*p[6][0] + p[2][2]*p[5][1]*p[5][2]*p[6][0] - 2*p[4][2]*p[5][1]*p[5][2]*p[6][0] - 2*pow(p[2][2], 2)*p[1][0]*p[6][1] + 2*pow(p[5][2], 2)*p[1][0]*p[6][1] + pow(p[1][2], 2)*p[2][0]*p[6][1] - pow(p[3][2], 2)*p[2][0]*p[6][1] + pow(p[5][2], 2)*p[2][0]*p[6][1] - pow(p[7][2], 2)*p[2][0]*p[6][1] - p[1][0]*p[1][2]*p[2][2]*p[6][1] + 2*p[1][2]*p[2][0]*p[2][2]*p[6][1] + 2*pow(p[2][2], 2)*p[3][0]*p[6][1] - 2*pow(p[7][2], 2)*p[3][0]*p[6][1] - 2*p[2][0]*p[2][2]*p[3][2]*p[6][1] + p[2][2]*p[3][0]*p[3][2]*p[6][1] - 2*pow(p[5][2], 2)*p[4][0]*p[6][1] + 2*pow(p[7][2], 2)*p[4][0]*p[6][1] - pow(p[1][2], 2)*p[5][0]*p[6][1] - pow(p[2][2], 2)*p[5][0]*p[6][1] + pow(p[4][2], 2)*p[5][0]*p[6][1] + pow(p[7][2], 2)*p[5][0]*p[6][1] - p[1][2]*p[2][2]*p[5][0]*p[6][1] + p[1][0]*p[1][2]*p[5][2]*p[6][1] + p[1][2]*p[2][0]*p[5][2]*p[6][1] + p[2][0]*p[2][2]*p[5][2]*p[6][1] - p[4][0]*p[4][2]*p[5][2]*p[6][1] - 2*p[1][2]*p[5][0]*p[5][2]*p[6][1] - p[2][2]*p[5][0]*p[5][2]*p[6][1] + 2*p[4][2]*p[5][0]*p[5][2]*p[6][1] - p[1][1]*p[1][2]*p[2][0]*p[6][2] + p[1][0]*p[1][2]*p[2][1]*p[6][2] - 2*p[1][1]*p[2][0]*p[2][2]*p[6][2] + 2*p[1][0]*p[2][1]*p[2][2]*p[6][2] - 2*p[2][1]*p[2][2]*p[3][0]*p[6][2] + 2*p[2][0]*p[2][2]*p[3][1]*p[6][2] - p[2][1]*p[3][0]*p[3][2]*p[6][2] + p[2][0]*p[3][1]*p[3][2]*p[6][2] + p[1][1]*p[1][2]*p[5][0]*p[6][2] + p[1][1]*p[2][2]*p[5][0]*p[6][2] + p[2][1]*p[2][2]*p[5][0]*p[6][2] - p[4][1]*p[4][2]*p[5][0]*p[6][2] - p[1][0]*p[1][2]*p[5][1]*p[6][2] - p[1][0]*p[2][2]*p[5][1]*p[6][2] - p[2][0]*p[2][2]*p[5][1]*p[6][2] + p[4][0]*p[4][2]*p[5][1]*p[6][2] - p[1][1]*p[2][0]*p[5][2]*p[6][2] + p[1][0]*p[2][1]*p[5][2]*p[6][2] + 2*p[1][1]*p[5][0]*p[5][2]*p[6][2] + p[2][1]*p[5][0]*p[5][2]*p[6][2] - 2*p[4][1]*p[5][0]*p[5][2]*p[6][2] - 2*p[1][0]*p[5][1]*p[5][2]*p[6][2] - p[2][0]*p[5][1]*p[5][2]*p[6][2] + 2*p[4][0]*p[5][1]*p[5][2]*p[6][2] - p[1][2]*p[2][1]*p[6][0]*p[6][2] + p[1][1]*p[2][2]*p[6][0]*p[6][2] - p[2][2]*p[3][1]*p[6][0]*p[6][2] + p[2][1]*p[3][2]*p[6][0]*p[6][2] + p[1][2]*p[5][1]*p[6][0]*p[6][2] + 2*p[2][2]*p[5][1]*p[6][0]*p[6][2] - p[4][2]*p[5][1]*p[6][0]*p[6][2] - p[1][1]*p[5][2]*p[6][0]*p[6][2] - 2*p[2][1]*p[5][2]*p[6][0]*p[6][2] + p[4][1]*p[5][2]*p[6][0]*p[6][2] + p[1][2]*p[2][0]*p[6][1]*p[6][2] - p[1][0]*p[2][2]*p[6][1]*p[6][2] + p[2][2]*p[3][0]*p[6][1]*p[6][2] - p[2][0]*p[3][2]*p[6][1]*p[6][2] - p[1][2]*p[5][0]*p[6][1]*p[6][2] - 2*p[2][2]*p[5][0]*p[6][1]*p[6][2] + p[4][2]*p[5][0]*p[6][1]*p[6][2] + p[1][0]*p[5][2]*p[6][1]*p[6][2] + 2*p[2][0]*p[5][2]*p[6][1]*p[6][2] - p[4][0]*p[5][2]*p[6][1]*p[6][2] + 2*pow(p[3][2], 2)*p[2][1]*p[7][0] - 2*pow(p[6][2], 2)*p[2][1]*p[7][0] - pow(p[2][2], 2)*p[3][1]*p[7][0] + pow(p[4][2], 2)*p[3][1]*p[7][0] - pow(p[6][2], 2)*p[3][1]*p[7][0] + p[2][1]*p[2][2]*p[3][2]*p[7][0] - 2*p[2][2]*p[3][1]*p[3][2]*p[7][0] - pow(p[3][2], 2)*p[4][1]*p[7][0] + pow(p[5][2], 2)*p[4][1]*p[7][0] + pow(p[6][2], 2)*p[4][1]*p[7][0] + p[3][1]*p[3][2]*p[4][2]*p[7][0] - p[3][2]*p[4][1]*p[4][2]*p[7][0] - 2*pow(p[4][2], 2)*p[5][1]*p[7][0] + 2*pow(p[6][2], 2)*p[5][1]*p[7][0] + 2*p[4][1]*p[4][2]*p[5][2]*p[7][0] - p[4][2]*p[5][1]*p[5][2]*p[7][0] + pow(p[2][2], 2)*p[6][1]*p[7][0] + pow(p[3][2], 2)*p[6][1]*p[7][0] - pow(p[4][2], 2)*p[6][1]*p[7][0] - pow(p[5][2], 2)*p[6][1]*p[7][0] + p[2][2]*p[3][2]*p[6][1]*p[7][0] - p[4][2]*p[5][2]*p[6][1]*p[7][0] - p[2][1]*p[2][2]*p[6][2]*p[7][0] - p[2][2]*p[3][1]*p[6][2]*p[7][0] - p[3][1]*p[3][2]*p[6][2]*p[7][0] + p[4][1]*p[4][2]*p[6][2]*p[7][0] + p[4][1]*p[5][2]*p[6][2]*p[7][0] + p[5][1]*p[5][2]*p[6][2]*p[7][0] + 2*p[2][2]*p[6][1]*p[6][2]*p[7][0] + p[3][2]*p[6][1]*p[6][2]*p[7][0] - p[4][2]*p[6][1]*p[6][2]*p[7][0] - 2*p[5][2]*p[6][1]*p[6][2]*p[7][0] + 2*pow(p[3][2], 2)*p[0][0]*p[7][1] - 2*pow(p[4][2], 2)*p[0][0]*p[7][1] - 2*pow(p[3][2], 2)*p[2][0]*p[7][1] + 2*pow(p[6][2], 2)*p[2][0]*p[7][1] + pow(p[2][2], 2)*p[3][0]*p[7][1] - pow(p[4][2], 2)*p[3][0]*p[7][1] + pow(p[6][2], 2)*p[3][0]*p[7][1] - p[2][0]*p[2][2]*p[3][2]*p[7][1] + 2*p[2][2]*p[3][0]*p[3][2]*p[7][1] + pow(p[3][2], 2)*p[4][0]*p[7][1] - pow(p[5][2], 2)*p[4][0]*p[7][1] - pow(p[6][2], 2)*p[4][0]*p[7][1] - p[3][0]*p[3][2]*p[4][2]*p[7][1] + p[3][2]*p[4][0]*p[4][2]*p[7][1] + 2*pow(p[4][2], 2)*p[5][0]*p[7][1] - 2*pow(p[6][2], 2)*p[5][0]*p[7][1] - 2*p[4][0]*p[4][2]*p[5][2]*p[7][1] + p[4][2]*p[5][0]*p[5][2]*p[7][1] - pow(p[2][2], 2)*p[6][0]*p[7][1] - pow(p[3][2], 2)*p[6][0]*p[7][1] + pow(p[4][2], 2)*p[6][0]*p[7][1] + pow(p[5][2], 2)*p[6][0]*p[7][1] - p[2][2]*p[3][2]*p[6][0]*p[7][1] + p[4][2]*p[5][2]*p[6][0]*p[7][1] + p[2][0]*p[2][2]*p[6][2]*p[7][1] + p[2][2]*p[3][0]*p[6][2]*p[7][1] + p[3][0]*p[3][2]*p[6][2]*p[7][1] - p[4][0]*p[4][2]*p[6][2]*p[7][1] - p[4][0]*p[5][2]*p[6][2]*p[7][1] - p[5][0]*p[5][2]*p[6][2]*p[7][1] - 2*p[2][2]*p[6][0]*p[6][2]*p[7][1] - p[3][2]*p[6][0]*p[6][2]*p[7][1] + p[4][2]*p[6][0]*p[6][2]*p[7][1] + 2*p[5][2]*p[6][0]*p[6][2]*p[7][1] + pow(p[0][2], 2)*(p[2][1]*p[3][0] - p[2][0]*p[3][1] + 2*p[3][1]*p[4][0] - 2*p[3][0]*p[4][1] + p[1][1]*(p[2][0] + 2*p[3][0] - 2*p[4][0] - p[5][0]) + p[4][1]*p[5][0] - p[4][0]*p[5][1] + p[1][0]*(-p[2][1] - 2*p[3][1] + 2*p[4][1] + p[5][1]) + p[3][1]*p[7][0] - p[4][1]*p[7][0] - p[3][0]*p[7][1] + p[4][0]*p[7][1]) - p[2][1]*p[2][2]*p[3][0]*p[7][2] + p[2][0]*p[2][2]*p[3][1]*p[7][2] - 2*p[2][1]*p[3][0]*p[3][2]*p[7][2] - 2*p[0][0]*p[3][1]*p[3][2]*p[7][2] + 2*p[2][0]*p[3][1]*p[3][2]*p[7][2] - p[3][1]*p[3][2]*p[4][0]*p[7][2] + p[0][0]*p[3][2]*p[4][1]*p[7][2] + p[3][0]*p[3][2]*p[4][1]*p[7][2] - p[0][0]*p[3][1]*p[4][2]*p[7][2] - p[3][1]*p[4][0]*p[4][2]*p[7][2] + 2*p[0][0]*p[4][1]*p[4][2]*p[7][2] + p[3][0]*p[4][1]*p[4][2]*p[7][2] - 2*p[4][1]*p[4][2]*p[5][0]*p[7][2] + 2*p[4][0]*p[4][2]*p[5][1]*p[7][2] - p[4][1]*p[5][0]*p[5][2]*p[7][2] + p[4][0]*p[5][1]*p[5][2]*p[7][2] + p[2][1]*p[2][2]*p[6][0]*p[7][2] + p[2][1]*p[3][2]*p[6][0]*p[7][2] + p[3][1]*p[3][2]*p[6][0]*p[7][2] - p[4][1]*p[4][2]*p[6][0]*p[7][2] - p[4][2]*p[5][1]*p[6][0]*p[7][2] - p[5][1]*p[5][2]*p[6][0]*p[7][2] - p[2][0]*p[2][2]*p[6][1]*p[7][2] - p[2][0]*p[3][2]*p[6][1]*p[7][2] - p[3][0]*p[3][2]*p[6][1]*p[7][2] + p[4][0]*p[4][2]*p[6][1]*p[7][2] + p[4][2]*p[5][0]*p[6][1]*p[7][2] + p[5][0]*p[5][2]*p[6][1]*p[7][2] - p[2][1]*p[3][0]*p[6][2]*p[7][2] + p[2][0]*p[3][1]*p[6][2]*p[7][2] - p[4][1]*p[5][0]*p[6][2]*p[7][2] + p[4][0]*p[5][1]*p[6][2]*p[7][2] + 2*p[2][1]*p[6][0]*p[6][2]*p[7][2] + p[3][1]*p[6][0]*p[6][2]*p[7][2] - p[4][1]*p[6][0]*p[6][2]*p[7][2] - 2*p[5][1]*p[6][0]*p[6][2]*p[7][2] - 2*p[2][0]*p[6][1]*p[6][2]*p[7][2] - p[3][0]*p[6][1]*p[6][2]*p[7][2] + p[4][0]*p[6][1]*p[6][2]*p[7][2] + 2*p[5][0]*p[6][1]*p[6][2]*p[7][2] - p[2][2]*p[3][1]*p[7][0]*p[7][2] + p[2][1]*p[3][2]*p[7][0]*p[7][2] - 2*p[3][2]*p[4][1]*p[7][0]*p[7][2] + 2*p[3][1]*p[4][2]*p[7][0]*p[7][2] - p[4][2]*p[5][1]*p[7][0]*p[7][2] + p[4][1]*p[5][2]*p[7][0]*p[7][2] + p[2][2]*p[6][1]*p[7][0]*p[7][2] + 2*p[3][2]*p[6][1]*p[7][0]*p[7][2] - 2*p[4][2]*p[6][1]*p[7][0]*p[7][2] - p[5][2]*p[6][1]*p[7][0]*p[7][2] - p[2][1]*p[6][2]*p[7][0]*p[7][2] - 2*p[3][1]*p[6][2]*p[7][0]*p[7][2] + 2*p[4][1]*p[6][2]*p[7][0]*p[7][2] + p[5][1]*p[6][2]*p[7][0]*p[7][2] + p[2][2]*p[3][0]*p[7][1]*p[7][2] + p[0][0]*p[3][2]*p[7][1]*p[7][2] - p[2][0]*p[3][2]*p[7][1]*p[7][2] + 2*p[3][2]*p[4][0]*p[7][1]*p[7][2] - p[0][0]*p[4][2]*p[7][1]*p[7][2] - 2*p[3][0]*p[4][2]*p[7][1]*p[7][2] + p[4][2]*p[5][0]*p[7][1]*p[7][2] - p[4][0]*p[5][2]*p[7][1]*p[7][2] - p[2][2]*p[6][0]*p[7][1]*p[7][2] - 2*p[3][2]*p[6][0]*p[7][1]*p[7][2] + 2*p[4][2]*p[6][0]*p[7][1]*p[7][2] + p[5][2]*p[6][0]*p[7][1]*p[7][2] + p[2][0]*p[6][2]*p[7][1]*p[7][2] + 2*p[3][0]*p[6][2]*p[7][1]*p[7][2] - 2*p[4][0]*p[6][2]*p[7][1]*p[7][2] - p[5][0]*p[6][2]*p[7][1]*p[7][2] + p[0][1]*(2*pow(p[3][2], 2)*p[2][0] - pow(p[2][2], 2)*p[3][0] + pow(p[4][2], 2)*p[3][0] + pow(p[7][2], 2)*p[3][0] + p[2][0]*p[2][2]*p[3][2] - 2*p[2][2]*p[3][0]*p[3][2] - pow(p[3][2], 2)*p[4][0] + pow(p[5][2], 2)*p[4][0] - pow(p[7][2], 2)*p[4][0] + p[3][0]*p[3][2]*p[4][2] - p[3][2]*p[4][0]*p[4][2] - 2*pow(p[4][2], 2)*p[5][0] + pow(p[1][2], 2)*(-2*p[2][0] - p[3][0] + p[4][0] + 2*p[5][0]) + 2*p[4][0]*p[4][2]*p[5][2] - p[4][2]*p[5][0]*p[5][2] + p[1][0]*(pow(p[2][2], 2) + pow(p[3][2], 2) - pow(p[4][2], 2) - pow(p[5][2], 2) + p[2][2]*p[3][2] - p[4][2]*p[5][2]) + p[1][2]*(-(p[2][0]*p[2][2]) - p[2][2]*p[3][0] - p[3][0]*p[3][2] + p[4][0]*p[4][2] + p[1][0]*(2*p[2][2] + p[3][2] - p[4][2] - 2*p[5][2]) + p[4][0]*p[5][2] + p[5][0]*p[5][2]) - 2*pow(p[3][2], 2)*p[7][0] + 2*pow(p[4][2], 2)*p[7][0] + 2*p[3][0]*p[3][2]*p[7][2] - p[3][2]*p[4][0]*p[7][2] + p[3][0]*p[4][2]*p[7][2] - 2*p[4][0]*p[4][2]*p[7][2] - p[3][2]*p[7][0]*p[7][2] + p[4][2]*p[7][0]*p[7][2]) + p[0][2]*(p[0][0]*p[1][2]*p[2][1] - 2*p[1][0]*p[1][2]*p[2][1] - p[1][0]*p[2][1]*p[2][2] + p[1][2]*p[2][1]*p[3][0] + p[2][1]*p[2][2]*p[3][0] + 2*p[0][0]*p[1][2]*p[3][1] - p[1][0]*p[1][2]*p[3][1] - p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] - p[2][0]*p[2][2]*p[3][1] - p[0][0]*p[2][1]*p[3][2] - p[1][0]*p[2][1]*p[3][2] + 2*p[2][1]*p[3][0]*p[3][2] - p[1][0]*p[3][1]*p[3][2] - 2*p[2][0]*p[3][1]*p[3][2] + p[3][1]*p[3][2]*p[4][0] - 2*p[0][0]*p[1][2]*p[4][1] + p[1][0]*p[1][2]*p[4][1] + 2*p[0][0]*p[3][2]*p[4][1] - p[3][0]*p[3][2]*p[4][1] - 2*p[0][0]*p[3][1]*p[4][2] + p[3][1]*p[4][0]*p[4][2] + p[1][0]*p[4][1]*p[4][2] - p[3][0]*p[4][1]*p[4][2] + p[1][2]*p[4][1]*p[5][0] + 2*p[4][1]*p[4][2]*p[5][0] - p[0][0]*p[1][2]*p[5][1] + 2*p[1][0]*p[1][2]*p[5][1] - p[1][2]*p[4][0]*p[5][1] + p[0][0]*p[4][2]*p[5][1] + p[1][0]*p[4][2]*p[5][1] - 2*p[4][0]*p[4][2]*p[5][1] - p[0][0]*p[4][1]*p[5][2] + p[4][1]*p[5][0]*p[5][2] + p[1][0]*p[5][1]*p[5][2] - p[4][0]*p[5][1]*p[5][2] + p[1][1]*(p[2][0]*p[2][2] + p[2][0]*p[3][2] + p[3][0]*p[3][2] - p[4][0]*p[4][2] + p[1][2]*(2*p[2][0] + p[3][0] - p[4][0] - 2*p[5][0]) - p[4][2]*p[5][0] - p[5][0]*p[5][2] + p[0][0]*(-p[2][2] - 2*p[3][2] + 2*p[4][2] + p[5][2])) + 2*p[3][1]*p[3][2]*p[7][0] - p[3][2]*p[4][1]*p[7][0] + p[3][1]*p[4][2]*p[7][0] - 2*p[4][1]*p[4][2]*p[7][0] + p[0][0]*p[3][2]*p[7][1] - 2*p[3][0]*p[3][2]*p[7][1] + p[3][2]*p[4][0]*p[7][1] - p[0][0]*p[4][2]*p[7][1] - p[3][0]*p[4][2]*p[7][1] + 2*p[4][0]*p[4][2]*p[7][1] - p[0][0]*p[3][1]*p[7][2] + p[0][0]*p[4][1]*p[7][2] + p[3][1]*p[7][0]*p[7][2] - p[4][1]*p[7][0]*p[7][2] - p[3][0]*p[7][1]*p[7][2] + p[4][0]*p[7][1]*p[7][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[2][0]*p[3][2] - 2*p[3][2]*p[4][0] + 2*p[3][0]*p[4][2] - p[4][2]*p[5][0] + p[1][2]*(-p[2][0] - 2*p[3][0] + 2*p[4][0] + p[5][0]) + p[1][0]*(p[2][2] + 2*p[3][2] - 2*p[4][2] - p[5][2]) + p[4][0]*p[5][2] - p[3][2]*p[7][0] + p[4][2]*p[7][0] + p[3][0]*p[7][2] - p[4][0]*p[7][2])))/72;
382  return result;
383  }
384 };
385 
386 // 3D Wedge
387 
388 template<>
389 class Operator<VolumeIntegral, basis::Unitary, Wedge>:
390  public OperatorBase<VolumeIntegral, basis::Unitary, Interval>
391 {
392 public:
393  static constexpr size_t dim = dimension(Hexahedron);
394  static constexpr size_t operator_size = 1;
395  static constexpr size_t basis_size =
396  basis::Traits<basis::Unitary, dim>::function_size;
397  static constexpr size_t point_size = 6;
398 
399  using result_t = std::array<std::array<double, operator_size>, basis_size>;
400  using points_t = std::array<Wonton::Point<dim>, point_size>;
401 
402  static result_t apply(const points_t p) {
403  result_t result;
404  result[0][0] = (2*p[0][0]*p[1][2]*p[2][1] - 2*p[0][0]*p[1][1]*p[2][2] - p[0][0]*p[1][2]*p[3][1] + p[0][0]*p[2][2]*p[3][1] + p[0][0]*p[1][1]*p[3][2] - p[0][0]*p[2][1]*p[3][2] - p[1][2]*p[2][1]*p[4][0] + p[1][1]*p[2][2]*p[4][0] + p[1][2]*p[3][1]*p[4][0] - p[1][1]*p[3][2]*p[4][0] - p[0][0]*p[1][2]*p[4][1] + p[1][2]*p[2][0]*p[4][1] - p[1][0]*p[2][2]*p[4][1] - p[1][2]*p[3][0]*p[4][1] + p[0][0]*p[3][2]*p[4][1] + p[1][0]*p[3][2]*p[4][1] + p[0][0]*p[1][1]*p[4][2] - p[1][1]*p[2][0]*p[4][2] + p[1][0]*p[2][1]*p[4][2] + p[1][1]*p[3][0]*p[4][2] - p[0][0]*p[3][1]*p[4][2] - p[1][0]*p[3][1]*p[4][2] - p[1][2]*p[2][1]*p[5][0] + p[1][1]*p[2][2]*p[5][0] - p[2][2]*p[3][1]*p[5][0] + p[2][1]*p[3][2]*p[5][0] + p[1][2]*p[4][1]*p[5][0] + p[2][2]*p[4][1]*p[5][0] - 2*p[3][2]*p[4][1]*p[5][0] - p[1][1]*p[4][2]*p[5][0] - p[2][1]*p[4][2]*p[5][0] + 2*p[3][1]*p[4][2]*p[5][0] + p[1][2]*p[2][0]*p[5][1] + p[0][0]*p[2][2]*p[5][1] - p[1][0]*p[2][2]*p[5][1] + p[2][2]*p[3][0]*p[5][1] - p[0][0]*p[3][2]*p[5][1] - p[2][0]*p[3][2]*p[5][1] - p[1][2]*p[4][0]*p[5][1] - p[2][2]*p[4][0]*p[5][1] + 2*p[3][2]*p[4][0]*p[5][1] + p[1][0]*p[4][2]*p[5][1] + p[2][0]*p[4][2]*p[5][1] - 2*p[3][0]*p[4][2]*p[5][1] + p[0][2]*(p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[1][1]*(2*p[2][0] - p[3][0] - p[4][0]) + p[3][1]*p[4][0] - p[3][0]*p[4][1] + p[1][0]*(-2*p[2][1] + p[3][1] + p[4][1]) + p[2][1]*p[5][0] - p[3][1]*p[5][0] - p[2][0]*p[5][1] + p[3][0]*p[5][1]) - p[1][1]*p[2][0]*p[5][2] - p[0][0]*p[2][1]*p[5][2] + p[1][0]*p[2][1]*p[5][2] - p[2][1]*p[3][0]*p[5][2] + p[0][0]*p[3][1]*p[5][2] + p[2][0]*p[3][1]*p[5][2] + p[1][1]*p[4][0]*p[5][2] + p[2][1]*p[4][0]*p[5][2] - 2*p[3][1]*p[4][0]*p[5][2] - p[1][0]*p[4][1]*p[5][2] - p[2][0]*p[4][1]*p[5][2] + 2*p[3][0]*p[4][1]*p[5][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[2][0]*p[3][2] - p[3][2]*p[4][0] + p[1][2]*(-2*p[2][0] + p[3][0] + p[4][0]) + p[1][0]*(2*p[2][2] - p[3][2] - p[4][2]) + p[3][0]*p[4][2] - p[2][2]*p[5][0] + p[3][2]*p[5][0] + p[2][0]*p[5][2] - p[3][0]*p[5][2]))/12;
405  return result;
406  }
407 };
408 
409 template<>
410 class Operator<VolumeIntegral, basis::Linear, Wedge>:
411  public OperatorBase<VolumeIntegral, basis::Linear, Interval>
412 {
413 public:
414  static constexpr size_t dim = dimension(Wedge);
415  static constexpr size_t operator_size = 1;
416  static constexpr size_t basis_size =
417  basis::Traits<basis::Linear, dim>::function_size;
418  static constexpr size_t point_size = 6;
419 
420  using result_t = std::array<std::array<double, operator_size>, basis_size>;
421  using points_t = std::array<Wonton::Point<dim>, point_size>;
422 
423  static result_t apply(const points_t p) {
424  result_t result;
425  result[0][0] = (2*p[0][0]*p[1][2]*p[2][1] - 2*p[0][0]*p[1][1]*p[2][2] - p[0][0]*p[1][2]*p[3][1] + p[0][0]*p[2][2]*p[3][1] + p[0][0]*p[1][1]*p[3][2] - p[0][0]*p[2][1]*p[3][2] - p[1][2]*p[2][1]*p[4][0] + p[1][1]*p[2][2]*p[4][0] + p[1][2]*p[3][1]*p[4][0] - p[1][1]*p[3][2]*p[4][0] - p[0][0]*p[1][2]*p[4][1] + p[1][2]*p[2][0]*p[4][1] - p[1][0]*p[2][2]*p[4][1] - p[1][2]*p[3][0]*p[4][1] + p[0][0]*p[3][2]*p[4][1] + p[1][0]*p[3][2]*p[4][1] + p[0][0]*p[1][1]*p[4][2] - p[1][1]*p[2][0]*p[4][2] + p[1][0]*p[2][1]*p[4][2] + p[1][1]*p[3][0]*p[4][2] - p[0][0]*p[3][1]*p[4][2] - p[1][0]*p[3][1]*p[4][2] - p[1][2]*p[2][1]*p[5][0] + p[1][1]*p[2][2]*p[5][0] - p[2][2]*p[3][1]*p[5][0] + p[2][1]*p[3][2]*p[5][0] + p[1][2]*p[4][1]*p[5][0] + p[2][2]*p[4][1]*p[5][0] - 2*p[3][2]*p[4][1]*p[5][0] - p[1][1]*p[4][2]*p[5][0] - p[2][1]*p[4][2]*p[5][0] + 2*p[3][1]*p[4][2]*p[5][0] + p[1][2]*p[2][0]*p[5][1] + p[0][0]*p[2][2]*p[5][1] - p[1][0]*p[2][2]*p[5][1] + p[2][2]*p[3][0]*p[5][1] - p[0][0]*p[3][2]*p[5][1] - p[2][0]*p[3][2]*p[5][1] - p[1][2]*p[4][0]*p[5][1] - p[2][2]*p[4][0]*p[5][1] + 2*p[3][2]*p[4][0]*p[5][1] + p[1][0]*p[4][2]*p[5][1] + p[2][0]*p[4][2]*p[5][1] - 2*p[3][0]*p[4][2]*p[5][1] + p[0][2]*(p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[1][1]*(2*p[2][0] - p[3][0] - p[4][0]) + p[3][1]*p[4][0] - p[3][0]*p[4][1] + p[1][0]*(-2*p[2][1] + p[3][1] + p[4][1]) + p[2][1]*p[5][0] - p[3][1]*p[5][0] - p[2][0]*p[5][1] + p[3][0]*p[5][1]) - p[1][1]*p[2][0]*p[5][2] - p[0][0]*p[2][1]*p[5][2] + p[1][0]*p[2][1]*p[5][2] - p[2][1]*p[3][0]*p[5][2] + p[0][0]*p[3][1]*p[5][2] + p[2][0]*p[3][1]*p[5][2] + p[1][1]*p[4][0]*p[5][2] + p[2][1]*p[4][0]*p[5][2] - 2*p[3][1]*p[4][0]*p[5][2] - p[1][0]*p[4][1]*p[5][2] - p[2][0]*p[4][1]*p[5][2] + 2*p[3][0]*p[4][1]*p[5][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[2][0]*p[3][2] - p[3][2]*p[4][0] + p[1][2]*(-2*p[2][0] + p[3][0] + p[4][0]) + p[1][0]*(2*p[2][2] - p[3][2] - p[4][2]) + p[3][0]*p[4][2] - p[2][2]*p[5][0] + p[3][2]*p[5][0] + p[2][0]*p[5][2] - p[3][0]*p[5][2]))/12;
426  result[1][0] = (-3*pow(p[2][0], 2)*p[0][1]*p[1][2] + pow(p[3][0], 2)*p[0][1]*p[1][2] + pow(p[4][0], 2)*p[0][1]*p[1][2] - 3*p[0][1]*p[1][0]*p[1][2]*p[2][0] - pow(p[4][0], 2)*p[1][2]*p[2][1] - pow(p[5][0], 2)*p[1][2]*p[2][1] + 3*pow(p[1][0], 2)*p[0][1]*p[2][2] - pow(p[3][0], 2)*p[0][1]*p[2][2] - pow(p[5][0], 2)*p[0][1]*p[2][2] + pow(p[4][0], 2)*p[1][1]*p[2][2] + pow(p[5][0], 2)*p[1][1]*p[2][2] + 3*p[0][1]*p[1][0]*p[2][0]*p[2][2] + p[0][1]*p[1][0]*p[1][2]*p[3][0] - p[0][1]*p[2][0]*p[2][2]*p[3][0] + 2*pow(p[4][0], 2)*p[1][2]*p[3][1] - 2*pow(p[5][0], 2)*p[2][2]*p[3][1] - pow(p[1][0], 2)*p[0][1]*p[3][2] + pow(p[2][0], 2)*p[0][1]*p[3][2] - pow(p[4][0], 2)*p[0][1]*p[3][2] + pow(p[5][0], 2)*p[0][1]*p[3][2] - 2*pow(p[4][0], 2)*p[1][1]*p[3][2] + 2*pow(p[5][0], 2)*p[2][1]*p[3][2] - p[0][1]*p[1][0]*p[3][0]*p[3][2] + p[0][1]*p[2][0]*p[3][0]*p[3][2] + 2*p[0][1]*p[1][0]*p[1][2]*p[4][0] - 2*p[1][0]*p[1][2]*p[2][1]*p[4][0] - p[1][2]*p[2][0]*p[2][1]*p[4][0] + 2*p[1][0]*p[1][1]*p[2][2]*p[4][0] + p[1][1]*p[2][0]*p[2][2]*p[4][0] + p[0][1]*p[1][2]*p[3][0]*p[4][0] + p[1][0]*p[1][2]*p[3][1]*p[4][0] + p[1][2]*p[3][0]*p[3][1]*p[4][0] - p[0][1]*p[1][0]*p[3][2]*p[4][0] - p[1][0]*p[1][1]*p[3][2]*p[4][0] - 2*p[0][1]*p[3][0]*p[3][2]*p[4][0] - p[1][1]*p[3][0]*p[3][2]*p[4][0] + pow(p[2][0], 2)*p[1][2]*p[4][1] - pow(p[3][0], 2)*p[1][2]*p[4][1] + pow(p[5][0], 2)*p[1][2]*p[4][1] + 2*p[1][0]*p[1][2]*p[2][0]*p[4][1] - 2*pow(p[1][0], 2)*p[2][2]*p[4][1] + 2*pow(p[5][0], 2)*p[2][2]*p[4][1] - p[1][0]*p[2][0]*p[2][2]*p[4][1] - p[1][0]*p[1][2]*p[3][0]*p[4][1] + pow(p[1][0], 2)*p[3][2]*p[4][1] - 3*pow(p[5][0], 2)*p[3][2]*p[4][1] + p[1][0]*p[3][0]*p[3][2]*p[4][1] + p[1][2]*p[2][0]*p[4][0]*p[4][1] - p[1][0]*p[2][2]*p[4][0]*p[4][1] - 2*p[1][2]*p[3][0]*p[4][0]*p[4][1] + 2*p[1][0]*p[3][2]*p[4][0]*p[4][1] - 2*pow(p[1][0], 2)*p[0][1]*p[4][2] + 2*pow(p[3][0], 2)*p[0][1]*p[4][2] - pow(p[2][0], 2)*p[1][1]*p[4][2] + pow(p[3][0], 2)*p[1][1]*p[4][2] - pow(p[5][0], 2)*p[1][1]*p[4][2] - 2*p[1][0]*p[1][1]*p[2][0]*p[4][2] + 2*pow(p[1][0], 2)*p[2][1]*p[4][2] - 2*pow(p[5][0], 2)*p[2][1]*p[4][2] + p[1][0]*p[2][0]*p[2][1]*p[4][2] + p[1][0]*p[1][1]*p[3][0]*p[4][2] - pow(p[1][0], 2)*p[3][1]*p[4][2] + 3*pow(p[5][0], 2)*p[3][1]*p[4][2] - p[1][0]*p[3][0]*p[3][1]*p[4][2] - p[0][1]*p[1][0]*p[4][0]*p[4][2] - p[1][1]*p[2][0]*p[4][0]*p[4][2] + p[1][0]*p[2][1]*p[4][0]*p[4][2] + p[0][1]*p[3][0]*p[4][0]*p[4][2] + 2*p[1][1]*p[3][0]*p[4][0]*p[4][2] - 2*p[1][0]*p[3][1]*p[4][0]*p[4][2] - p[1][0]*p[1][2]*p[2][1]*p[5][0] - 2*p[1][2]*p[2][0]*p[2][1]*p[5][0] + p[1][0]*p[1][1]*p[2][2]*p[5][0] - 2*p[0][1]*p[2][0]*p[2][2]*p[5][0] + 2*p[1][1]*p[2][0]*p[2][2]*p[5][0] - p[0][1]*p[2][2]*p[3][0]*p[5][0] - p[2][0]*p[2][2]*p[3][1]*p[5][0] - p[2][2]*p[3][0]*p[3][1]*p[5][0] + p[0][1]*p[2][0]*p[3][2]*p[5][0] + p[2][0]*p[2][1]*p[3][2]*p[5][0] + 2*p[0][1]*p[3][0]*p[3][2]*p[5][0] + p[2][1]*p[3][0]*p[3][2]*p[5][0] - p[1][2]*p[2][1]*p[4][0]*p[5][0] + p[1][1]*p[2][2]*p[4][0]*p[5][0] + p[1][0]*p[1][2]*p[4][1]*p[5][0] + p[1][2]*p[2][0]*p[4][1]*p[5][0] + p[2][0]*p[2][2]*p[4][1]*p[5][0] - 3*p[3][0]*p[3][2]*p[4][1]*p[5][0] + 2*p[1][2]*p[4][0]*p[4][1]*p[5][0] + p[2][2]*p[4][0]*p[4][1]*p[5][0] - 3*p[3][2]*p[4][0]*p[4][1]*p[5][0] - p[1][0]*p[1][1]*p[4][2]*p[5][0] - p[1][1]*p[2][0]*p[4][2]*p[5][0] - p[2][0]*p[2][1]*p[4][2]*p[5][0] + 3*p[3][0]*p[3][1]*p[4][2]*p[5][0] - 2*p[1][1]*p[4][0]*p[4][2]*p[5][0] - p[2][1]*p[4][0]*p[4][2]*p[5][0] + 3*p[3][1]*p[4][0]*p[4][2]*p[5][0] + 2*pow(p[2][0], 2)*p[1][2]*p[5][1] - 2*pow(p[4][0], 2)*p[1][2]*p[5][1] + p[1][0]*p[1][2]*p[2][0]*p[5][1] - pow(p[1][0], 2)*p[2][2]*p[5][1] + pow(p[3][0], 2)*p[2][2]*p[5][1] - pow(p[4][0], 2)*p[2][2]*p[5][1] - 2*p[1][0]*p[2][0]*p[2][2]*p[5][1] + p[2][0]*p[2][2]*p[3][0]*p[5][1] - pow(p[2][0], 2)*p[3][2]*p[5][1] + 3*pow(p[4][0], 2)*p[3][2]*p[5][1] - p[2][0]*p[3][0]*p[3][2]*p[5][1] - p[1][0]*p[1][2]*p[4][0]*p[5][1] - p[1][0]*p[2][2]*p[4][0]*p[5][1] - p[2][0]*p[2][2]*p[4][0]*p[5][1] + 3*p[3][0]*p[3][2]*p[4][0]*p[5][1] + pow(p[1][0], 2)*p[4][2]*p[5][1] + pow(p[2][0], 2)*p[4][2]*p[5][1] - 3*pow(p[3][0], 2)*p[4][2]*p[5][1] + p[1][0]*p[2][0]*p[4][2]*p[5][1] + 2*p[1][0]*p[4][0]*p[4][2]*p[5][1] + p[2][0]*p[4][0]*p[4][2]*p[5][1] - 3*p[3][0]*p[4][0]*p[4][2]*p[5][1] + p[1][2]*p[2][0]*p[5][0]*p[5][1] - p[1][0]*p[2][2]*p[5][0]*p[5][1] + 2*p[2][2]*p[3][0]*p[5][0]*p[5][1] - 2*p[2][0]*p[3][2]*p[5][0]*p[5][1] - p[1][2]*p[4][0]*p[5][0]*p[5][1] - 2*p[2][2]*p[4][0]*p[5][0]*p[5][1] + 3*p[3][2]*p[4][0]*p[5][0]*p[5][1] + p[1][0]*p[4][2]*p[5][0]*p[5][1] + 2*p[2][0]*p[4][2]*p[5][0]*p[5][1] - 3*p[3][0]*p[4][2]*p[5][0]*p[5][1] + p[0][2]*(pow(p[3][0], 2)*p[2][1] + pow(p[5][0], 2)*p[2][1] + p[2][0]*p[2][1]*p[3][0] - pow(p[2][0], 2)*p[3][1] + pow(p[4][0], 2)*p[3][1] - pow(p[5][0], 2)*p[3][1] - p[2][0]*p[3][0]*p[3][1] + 2*p[3][0]*p[3][1]*p[4][0] + p[1][1]*(3*pow(p[2][0], 2) - pow(p[3][0], 2) - pow(p[4][0], 2) - p[3][0]*p[4][0]) - 2*pow(p[3][0], 2)*p[4][1] - p[3][0]*p[4][0]*p[4][1] + pow(p[1][0], 2)*(-3*p[2][1] + p[3][1] + 2*p[4][1]) + p[1][0]*(-3*p[2][0]*p[2][1] + p[3][0]*p[3][1] + p[1][1]*(3*p[2][0] - p[3][0] - 2*p[4][0]) + p[3][1]*p[4][0] + p[4][0]*p[4][1]) + 2*p[2][0]*p[2][1]*p[5][0] + p[2][1]*p[3][0]*p[5][0] - p[2][0]*p[3][1]*p[5][0] - 2*p[3][0]*p[3][1]*p[5][0] - 2*pow(p[2][0], 2)*p[5][1] + 2*pow(p[3][0], 2)*p[5][1] - p[2][0]*p[5][0]*p[5][1] + p[3][0]*p[5][0]*p[5][1]) + 2*pow(p[2][0], 2)*p[0][1]*p[5][2] - 2*pow(p[3][0], 2)*p[0][1]*p[5][2] - 2*pow(p[2][0], 2)*p[1][1]*p[5][2] + 2*pow(p[4][0], 2)*p[1][1]*p[5][2] - p[1][0]*p[1][1]*p[2][0]*p[5][2] + pow(p[1][0], 2)*p[2][1]*p[5][2] - pow(p[3][0], 2)*p[2][1]*p[5][2] + pow(p[4][0], 2)*p[2][1]*p[5][2] + 2*p[1][0]*p[2][0]*p[2][1]*p[5][2] - p[2][0]*p[2][1]*p[3][0]*p[5][2] + pow(p[2][0], 2)*p[3][1]*p[5][2] - 3*pow(p[4][0], 2)*p[3][1]*p[5][2] + p[2][0]*p[3][0]*p[3][1]*p[5][2] + p[1][0]*p[1][1]*p[4][0]*p[5][2] + p[1][0]*p[2][1]*p[4][0]*p[5][2] + p[2][0]*p[2][1]*p[4][0]*p[5][2] - 3*p[3][0]*p[3][1]*p[4][0]*p[5][2] - pow(p[1][0], 2)*p[4][1]*p[5][2] - pow(p[2][0], 2)*p[4][1]*p[5][2] + 3*pow(p[3][0], 2)*p[4][1]*p[5][2] - p[1][0]*p[2][0]*p[4][1]*p[5][2] - 2*p[1][0]*p[4][0]*p[4][1]*p[5][2] - p[2][0]*p[4][0]*p[4][1]*p[5][2] + 3*p[3][0]*p[4][0]*p[4][1]*p[5][2] + p[0][1]*p[2][0]*p[5][0]*p[5][2] - p[1][1]*p[2][0]*p[5][0]*p[5][2] + p[1][0]*p[2][1]*p[5][0]*p[5][2] - p[0][1]*p[3][0]*p[5][0]*p[5][2] - 2*p[2][1]*p[3][0]*p[5][0]*p[5][2] + 2*p[2][0]*p[3][1]*p[5][0]*p[5][2] + p[1][1]*p[4][0]*p[5][0]*p[5][2] + 2*p[2][1]*p[4][0]*p[5][0]*p[5][2] - 3*p[3][1]*p[4][0]*p[5][0]*p[5][2] - p[1][0]*p[4][1]*p[5][0]*p[5][2] - 2*p[2][0]*p[4][1]*p[5][0]*p[5][2] + 3*p[3][0]*p[4][1]*p[5][0]*p[5][2] + pow(p[0][0], 2)*(2*p[2][2]*p[3][1] - 2*p[2][1]*p[3][2] + p[1][2]*(3*p[2][1] - 2*p[3][1] - p[4][1]) + p[3][2]*p[4][1] - p[3][1]*p[4][2] + p[1][1]*(-3*p[2][2] + 2*p[3][2] + p[4][2]) + p[2][2]*p[5][1] - p[3][2]*p[5][1] - p[2][1]*p[5][2] + p[3][1]*p[5][2]) + p[0][0]*(3*p[1][0]*p[1][2]*p[2][1] + 3*p[1][2]*p[2][0]*p[2][1] - 3*p[1][0]*p[1][1]*p[2][2] - 3*p[1][1]*p[2][0]*p[2][2] - p[1][0]*p[1][2]*p[3][1] + p[2][0]*p[2][2]*p[3][1] - p[1][2]*p[3][0]*p[3][1] + p[2][2]*p[3][0]*p[3][1] + p[1][0]*p[1][1]*p[3][2] - p[2][0]*p[2][1]*p[3][2] + p[1][1]*p[3][0]*p[3][2] - p[2][1]*p[3][0]*p[3][2] - 2*p[1][0]*p[1][2]*p[4][1] - p[1][2]*p[3][0]*p[4][1] + p[1][0]*p[3][2]*p[4][1] + 2*p[3][0]*p[3][2]*p[4][1] - p[1][2]*p[4][0]*p[4][1] + p[3][2]*p[4][0]*p[4][1] + 2*p[1][0]*p[1][1]*p[4][2] + p[1][1]*p[3][0]*p[4][2] - p[1][0]*p[3][1]*p[4][2] - 2*p[3][0]*p[3][1]*p[4][2] + p[1][1]*p[4][0]*p[4][2] - p[3][1]*p[4][0]*p[4][2] + 2*p[2][0]*p[2][2]*p[5][1] + p[2][2]*p[3][0]*p[5][1] - p[2][0]*p[3][2]*p[5][1] - 2*p[3][0]*p[3][2]*p[5][1] + p[2][2]*p[5][0]*p[5][1] - p[3][2]*p[5][0]*p[5][1] + p[0][2]*(2*p[2][1]*p[3][0] - 2*p[2][0]*p[3][1] + p[1][1]*(3*p[2][0] - 2*p[3][0] - p[4][0]) + p[3][1]*p[4][0] - p[3][0]*p[4][1] + p[1][0]*(-3*p[2][1] + 2*p[3][1] + p[4][1]) + p[2][1]*p[5][0] - p[3][1]*p[5][0] - p[2][0]*p[5][1] + p[3][0]*p[5][1]) - 2*p[2][0]*p[2][1]*p[5][2] - p[2][1]*p[3][0]*p[5][2] + p[2][0]*p[3][1]*p[5][2] + 2*p[3][0]*p[3][1]*p[5][2] - p[2][1]*p[5][0]*p[5][2] + p[3][1]*p[5][0]*p[5][2] + p[0][1]*(-2*p[2][2]*p[3][0] + 2*p[2][0]*p[3][2] - p[3][2]*p[4][0] + p[1][2]*(-3*p[2][0] + 2*p[3][0] + p[4][0]) + p[1][0]*(3*p[2][2] - 2*p[3][2] - p[4][2]) + p[3][0]*p[4][2] - p[2][2]*p[5][0] + p[3][2]*p[5][0] + p[2][0]*p[5][2] - p[3][0]*p[5][2])))/72;
427  result[2][0] = (3*pow(p[2][1], 2)*p[0][0]*p[1][2] - pow(p[3][1], 2)*p[0][0]*p[1][2] - pow(p[4][1], 2)*p[0][0]*p[1][2] + pow(p[4][1], 2)*p[1][2]*p[2][0] + pow(p[5][1], 2)*p[1][2]*p[2][0] + 3*p[0][0]*p[1][1]*p[1][2]*p[2][1] - 3*pow(p[1][1], 2)*p[0][0]*p[2][2] + pow(p[3][1], 2)*p[0][0]*p[2][2] + pow(p[5][1], 2)*p[0][0]*p[2][2] - pow(p[4][1], 2)*p[1][0]*p[2][2] - pow(p[5][1], 2)*p[1][0]*p[2][2] - 3*p[0][0]*p[1][1]*p[2][1]*p[2][2] - 2*pow(p[4][1], 2)*p[1][2]*p[3][0] + 2*pow(p[5][1], 2)*p[2][2]*p[3][0] - p[0][0]*p[1][1]*p[1][2]*p[3][1] + p[0][0]*p[2][1]*p[2][2]*p[3][1] + pow(p[1][1], 2)*p[0][0]*p[3][2] - pow(p[2][1], 2)*p[0][0]*p[3][2] + pow(p[4][1], 2)*p[0][0]*p[3][2] - pow(p[5][1], 2)*p[0][0]*p[3][2] + 2*pow(p[4][1], 2)*p[1][0]*p[3][2] - 2*pow(p[5][1], 2)*p[2][0]*p[3][2] + p[0][0]*p[1][1]*p[3][1]*p[3][2] - p[0][0]*p[2][1]*p[3][1]*p[3][2] - pow(p[2][1], 2)*p[1][2]*p[4][0] + pow(p[3][1], 2)*p[1][2]*p[4][0] - pow(p[5][1], 2)*p[1][2]*p[4][0] - 2*p[1][1]*p[1][2]*p[2][1]*p[4][0] + 2*pow(p[1][1], 2)*p[2][2]*p[4][0] - 2*pow(p[5][1], 2)*p[2][2]*p[4][0] + p[1][1]*p[2][1]*p[2][2]*p[4][0] + p[1][1]*p[1][2]*p[3][1]*p[4][0] - pow(p[1][1], 2)*p[3][2]*p[4][0] + 3*pow(p[5][1], 2)*p[3][2]*p[4][0] - p[1][1]*p[3][1]*p[3][2]*p[4][0] - 2*p[0][0]*p[1][1]*p[1][2]*p[4][1] + 2*p[1][1]*p[1][2]*p[2][0]*p[4][1] + p[1][2]*p[2][0]*p[2][1]*p[4][1] - 2*p[1][0]*p[1][1]*p[2][2]*p[4][1] - p[1][0]*p[2][1]*p[2][2]*p[4][1] - p[1][1]*p[1][2]*p[3][0]*p[4][1] - p[0][0]*p[1][2]*p[3][1]*p[4][1] - p[1][2]*p[3][0]*p[3][1]*p[4][1] + p[0][0]*p[1][1]*p[3][2]*p[4][1] + p[1][0]*p[1][1]*p[3][2]*p[4][1] + 2*p[0][0]*p[3][1]*p[3][2]*p[4][1] + p[1][0]*p[3][1]*p[3][2]*p[4][1] - p[1][2]*p[2][1]*p[4][0]*p[4][1] + p[1][1]*p[2][2]*p[4][0]*p[4][1] + 2*p[1][2]*p[3][1]*p[4][0]*p[4][1] - 2*p[1][1]*p[3][2]*p[4][0]*p[4][1] + 2*pow(p[1][1], 2)*p[0][0]*p[4][2] - 2*pow(p[3][1], 2)*p[0][0]*p[4][2] + pow(p[2][1], 2)*p[1][0]*p[4][2] - pow(p[3][1], 2)*p[1][0]*p[4][2] + pow(p[5][1], 2)*p[1][0]*p[4][2] - 2*pow(p[1][1], 2)*p[2][0]*p[4][2] + 2*pow(p[5][1], 2)*p[2][0]*p[4][2] + 2*p[1][0]*p[1][1]*p[2][1]*p[4][2] - p[1][1]*p[2][0]*p[2][1]*p[4][2] + pow(p[1][1], 2)*p[3][0]*p[4][2] - 3*pow(p[5][1], 2)*p[3][0]*p[4][2] - p[1][0]*p[1][1]*p[3][1]*p[4][2] + p[1][1]*p[3][0]*p[3][1]*p[4][2] + p[0][0]*p[1][1]*p[4][1]*p[4][2] - p[1][1]*p[2][0]*p[4][1]*p[4][2] + p[1][0]*p[2][1]*p[4][1]*p[4][2] + 2*p[1][1]*p[3][0]*p[4][1]*p[4][2] - p[0][0]*p[3][1]*p[4][1]*p[4][2] - 2*p[1][0]*p[3][1]*p[4][1]*p[4][2] - 2*pow(p[2][1], 2)*p[1][2]*p[5][0] + 2*pow(p[4][1], 2)*p[1][2]*p[5][0] - p[1][1]*p[1][2]*p[2][1]*p[5][0] + pow(p[1][1], 2)*p[2][2]*p[5][0] - pow(p[3][1], 2)*p[2][2]*p[5][0] + pow(p[4][1], 2)*p[2][2]*p[5][0] + 2*p[1][1]*p[2][1]*p[2][2]*p[5][0] - p[2][1]*p[2][2]*p[3][1]*p[5][0] + pow(p[2][1], 2)*p[3][2]*p[5][0] - 3*pow(p[4][1], 2)*p[3][2]*p[5][0] + p[2][1]*p[3][1]*p[3][2]*p[5][0] + p[1][1]*p[1][2]*p[4][1]*p[5][0] + p[1][1]*p[2][2]*p[4][1]*p[5][0] + p[2][1]*p[2][2]*p[4][1]*p[5][0] - 3*p[3][1]*p[3][2]*p[4][1]*p[5][0] - pow(p[1][1], 2)*p[4][2]*p[5][0] - pow(p[2][1], 2)*p[4][2]*p[5][0] + 3*pow(p[3][1], 2)*p[4][2]*p[5][0] - p[1][1]*p[2][1]*p[4][2]*p[5][0] - 2*p[1][1]*p[4][1]*p[4][2]*p[5][0] - p[2][1]*p[4][1]*p[4][2]*p[5][0] + 3*p[3][1]*p[4][1]*p[4][2]*p[5][0] + p[1][1]*p[1][2]*p[2][0]*p[5][1] + 2*p[1][2]*p[2][0]*p[2][1]*p[5][1] - p[1][0]*p[1][1]*p[2][2]*p[5][1] + 2*p[0][0]*p[2][1]*p[2][2]*p[5][1] - 2*p[1][0]*p[2][1]*p[2][2]*p[5][1] + p[2][1]*p[2][2]*p[3][0]*p[5][1] + p[0][0]*p[2][2]*p[3][1]*p[5][1] + p[2][2]*p[3][0]*p[3][1]*p[5][1] - p[0][0]*p[2][1]*p[3][2]*p[5][1] - p[2][0]*p[2][1]*p[3][2]*p[5][1] - 2*p[0][0]*p[3][1]*p[3][2]*p[5][1] - p[2][0]*p[3][1]*p[3][2]*p[5][1] - p[1][1]*p[1][2]*p[4][0]*p[5][1] - p[1][2]*p[2][1]*p[4][0]*p[5][1] - p[2][1]*p[2][2]*p[4][0]*p[5][1] + 3*p[3][1]*p[3][2]*p[4][0]*p[5][1] + p[1][2]*p[2][0]*p[4][1]*p[5][1] - p[1][0]*p[2][2]*p[4][1]*p[5][1] - 2*p[1][2]*p[4][0]*p[4][1]*p[5][1] - p[2][2]*p[4][0]*p[4][1]*p[5][1] + 3*p[3][2]*p[4][0]*p[4][1]*p[5][1] + p[1][0]*p[1][1]*p[4][2]*p[5][1] + p[1][0]*p[2][1]*p[4][2]*p[5][1] + p[2][0]*p[2][1]*p[4][2]*p[5][1] - 3*p[3][0]*p[3][1]*p[4][2]*p[5][1] + 2*p[1][0]*p[4][1]*p[4][2]*p[5][1] + p[2][0]*p[4][1]*p[4][2]*p[5][1] - 3*p[3][0]*p[4][1]*p[4][2]*p[5][1] - p[1][2]*p[2][1]*p[5][0]*p[5][1] + p[1][1]*p[2][2]*p[5][0]*p[5][1] - 2*p[2][2]*p[3][1]*p[5][0]*p[5][1] + 2*p[2][1]*p[3][2]*p[5][0]*p[5][1] + p[1][2]*p[4][1]*p[5][0]*p[5][1] + 2*p[2][2]*p[4][1]*p[5][0]*p[5][1] - 3*p[3][2]*p[4][1]*p[5][0]*p[5][1] - p[1][1]*p[4][2]*p[5][0]*p[5][1] - 2*p[2][1]*p[4][2]*p[5][0]*p[5][1] + 3*p[3][1]*p[4][2]*p[5][0]*p[5][1] + p[0][2]*(-(pow(p[3][1], 2)*p[2][0]) - pow(p[5][1], 2)*p[2][0] + pow(p[2][1], 2)*p[3][0] - pow(p[4][1], 2)*p[3][0] + pow(p[5][1], 2)*p[3][0] - p[2][0]*p[2][1]*p[3][1] + p[2][1]*p[3][0]*p[3][1] + pow(p[1][1], 2)*(3*p[2][0] - p[3][0] - 2*p[4][0]) + 2*pow(p[3][1], 2)*p[4][0] - 2*p[3][0]*p[3][1]*p[4][1] + p[3][1]*p[4][0]*p[4][1] + p[1][0]*(-3*pow(p[2][1], 2) + pow(p[3][1], 2) + pow(p[4][1], 2) + p[3][1]*p[4][1]) - p[1][1]*(-3*p[2][0]*p[2][1] + p[3][0]*p[3][1] + p[1][0]*(3*p[2][1] - p[3][1] - 2*p[4][1]) + p[3][0]*p[4][1] + p[4][0]*p[4][1]) + 2*pow(p[2][1], 2)*p[5][0] - 2*pow(p[3][1], 2)*p[5][0] - 2*p[2][0]*p[2][1]*p[5][1] + p[2][1]*p[3][0]*p[5][1] - p[2][0]*p[3][1]*p[5][1] + 2*p[3][0]*p[3][1]*p[5][1] + p[2][1]*p[5][0]*p[5][1] - p[3][1]*p[5][0]*p[5][1]) - 2*pow(p[2][1], 2)*p[0][0]*p[5][2] + 2*pow(p[3][1], 2)*p[0][0]*p[5][2] + 2*pow(p[2][1], 2)*p[1][0]*p[5][2] - 2*pow(p[4][1], 2)*p[1][0]*p[5][2] - pow(p[1][1], 2)*p[2][0]*p[5][2] + pow(p[3][1], 2)*p[2][0]*p[5][2] - pow(p[4][1], 2)*p[2][0]*p[5][2] + p[1][0]*p[1][1]*p[2][1]*p[5][2] - 2*p[1][1]*p[2][0]*p[2][1]*p[5][2] - pow(p[2][1], 2)*p[3][0]*p[5][2] + 3*pow(p[4][1], 2)*p[3][0]*p[5][2] + p[2][0]*p[2][1]*p[3][1]*p[5][2] - p[2][1]*p[3][0]*p[3][1]*p[5][2] + pow(p[1][1], 2)*p[4][0]*p[5][2] + pow(p[2][1], 2)*p[4][0]*p[5][2] - 3*pow(p[3][1], 2)*p[4][0]*p[5][2] + p[1][1]*p[2][1]*p[4][0]*p[5][2] - p[1][0]*p[1][1]*p[4][1]*p[5][2] - p[1][1]*p[2][0]*p[4][1]*p[5][2] - p[2][0]*p[2][1]*p[4][1]*p[5][2] + 3*p[3][0]*p[3][1]*p[4][1]*p[5][2] + 2*p[1][1]*p[4][0]*p[4][1]*p[5][2] + p[2][1]*p[4][0]*p[4][1]*p[5][2] - 3*p[3][1]*p[4][0]*p[4][1]*p[5][2] - p[1][1]*p[2][0]*p[5][1]*p[5][2] - p[0][0]*p[2][1]*p[5][1]*p[5][2] + p[1][0]*p[2][1]*p[5][1]*p[5][2] - 2*p[2][1]*p[3][0]*p[5][1]*p[5][2] + p[0][0]*p[3][1]*p[5][1]*p[5][2] + 2*p[2][0]*p[3][1]*p[5][1]*p[5][2] + p[1][1]*p[4][0]*p[5][1]*p[5][2] + 2*p[2][1]*p[4][0]*p[5][1]*p[5][2] - 3*p[3][1]*p[4][0]*p[5][1]*p[5][2] - p[1][0]*p[4][1]*p[5][1]*p[5][2] - 2*p[2][0]*p[4][1]*p[5][1]*p[5][2] + 3*p[3][0]*p[4][1]*p[5][1]*p[5][2] + pow(p[0][1], 2)*(-2*p[2][2]*p[3][0] + 2*p[2][0]*p[3][2] - p[3][2]*p[4][0] + p[1][2]*(-3*p[2][0] + 2*p[3][0] + p[4][0]) + p[1][0]*(3*p[2][2] - 2*p[3][2] - p[4][2]) + p[3][0]*p[4][2] - p[2][2]*p[5][0] + p[3][2]*p[5][0] + p[2][0]*p[5][2] - p[3][0]*p[5][2]) + p[0][1]*(3*p[0][0]*p[1][2]*p[2][1] - 3*p[1][2]*p[2][0]*p[2][1] + 3*p[1][0]*p[2][1]*p[2][2] - p[2][1]*p[2][2]*p[3][0] - 2*p[0][0]*p[1][2]*p[3][1] + 2*p[0][0]*p[2][2]*p[3][1] + p[1][2]*p[3][0]*p[3][1] - p[2][2]*p[3][0]*p[3][1] - 2*p[0][0]*p[2][1]*p[3][2] + p[2][0]*p[2][1]*p[3][2] - p[1][0]*p[3][1]*p[3][2] + p[2][0]*p[3][1]*p[3][2] + p[1][2]*p[3][1]*p[4][0] - 2*p[3][1]*p[3][2]*p[4][0] - p[0][0]*p[1][2]*p[4][1] + p[0][0]*p[3][2]*p[4][1] + p[1][2]*p[4][0]*p[4][1] - p[3][2]*p[4][0]*p[4][1] - p[0][0]*p[3][1]*p[4][2] - p[1][0]*p[3][1]*p[4][2] + 2*p[3][0]*p[3][1]*p[4][2] - p[1][0]*p[4][1]*p[4][2] + p[3][0]*p[4][1]*p[4][2] + p[1][1]*(3*p[1][0]*p[2][2] - p[1][0]*p[3][2] - p[3][2]*p[4][0] + p[1][2]*(-3*p[2][0] + p[3][0] + 2*p[4][0]) - 2*p[1][0]*p[4][2] + p[3][0]*p[4][2] + p[0][0]*(-3*p[2][2] + 2*p[3][2] + p[4][2])) - 2*p[2][1]*p[2][2]*p[5][0] - p[2][2]*p[3][1]*p[5][0] + p[2][1]*p[3][2]*p[5][0] + 2*p[3][1]*p[3][2]*p[5][0] + p[0][0]*p[2][2]*p[5][1] - p[0][0]*p[3][2]*p[5][1] - p[2][2]*p[5][0]*p[5][1] + p[3][2]*p[5][0]*p[5][1] + p[0][2]*(2*p[2][1]*p[3][0] - 2*p[2][0]*p[3][1] + p[1][1]*(3*p[2][0] - 2*p[3][0] - p[4][0]) + p[3][1]*p[4][0] - p[3][0]*p[4][1] + p[1][0]*(-3*p[2][1] + 2*p[3][1] + p[4][1]) + p[2][1]*p[5][0] - p[3][1]*p[5][0] - p[2][0]*p[5][1] + p[3][0]*p[5][1]) - p[0][0]*p[2][1]*p[5][2] + 2*p[2][0]*p[2][1]*p[5][2] - p[2][1]*p[3][0]*p[5][2] + p[0][0]*p[3][1]*p[5][2] + p[2][0]*p[3][1]*p[5][2] - 2*p[3][0]*p[3][1]*p[5][2] + p[2][0]*p[5][1]*p[5][2] - p[3][0]*p[5][1]*p[5][2]))/72;
428  result[3][0] = (-3*pow(p[2][2], 2)*p[0][0]*p[1][1] + pow(p[3][2], 2)*p[0][0]*p[1][1] + pow(p[4][2], 2)*p[0][0]*p[1][1] - pow(p[4][2], 2)*p[1][1]*p[2][0] - pow(p[5][2], 2)*p[1][1]*p[2][0] + 3*pow(p[1][2], 2)*p[0][0]*p[2][1] - pow(p[3][2], 2)*p[0][0]*p[2][1] - pow(p[5][2], 2)*p[0][0]*p[2][1] + pow(p[4][2], 2)*p[1][0]*p[2][1] + pow(p[5][2], 2)*p[1][0]*p[2][1] - 3*p[0][0]*p[1][1]*p[1][2]*p[2][2] + 3*p[0][0]*p[1][2]*p[2][1]*p[2][2] + 2*pow(p[4][2], 2)*p[1][1]*p[3][0] - 2*pow(p[5][2], 2)*p[2][1]*p[3][0] - pow(p[1][2], 2)*p[0][0]*p[3][1] + pow(p[2][2], 2)*p[0][0]*p[3][1] - pow(p[4][2], 2)*p[0][0]*p[3][1] + pow(p[5][2], 2)*p[0][0]*p[3][1] - 2*pow(p[4][2], 2)*p[1][0]*p[3][1] + 2*pow(p[5][2], 2)*p[2][0]*p[3][1] + p[0][0]*p[1][1]*p[1][2]*p[3][2] - p[0][0]*p[2][1]*p[2][2]*p[3][2] - p[0][0]*p[1][2]*p[3][1]*p[3][2] + p[0][0]*p[2][2]*p[3][1]*p[3][2] + pow(p[2][2], 2)*p[1][1]*p[4][0] - pow(p[3][2], 2)*p[1][1]*p[4][0] + pow(p[5][2], 2)*p[1][1]*p[4][0] - 2*pow(p[1][2], 2)*p[2][1]*p[4][0] + 2*pow(p[5][2], 2)*p[2][1]*p[4][0] + 2*p[1][1]*p[1][2]*p[2][2]*p[4][0] - p[1][2]*p[2][1]*p[2][2]*p[4][0] + pow(p[1][2], 2)*p[3][1]*p[4][0] - 3*pow(p[5][2], 2)*p[3][1]*p[4][0] - p[1][1]*p[1][2]*p[3][2]*p[4][0] + p[1][2]*p[3][1]*p[3][2]*p[4][0] - 2*pow(p[1][2], 2)*p[0][0]*p[4][1] + 2*pow(p[3][2], 2)*p[0][0]*p[4][1] - pow(p[2][2], 2)*p[1][0]*p[4][1] + pow(p[3][2], 2)*p[1][0]*p[4][1] - pow(p[5][2], 2)*p[1][0]*p[4][1] + 2*pow(p[1][2], 2)*p[2][0]*p[4][1] - 2*pow(p[5][2], 2)*p[2][0]*p[4][1] - 2*p[1][0]*p[1][2]*p[2][2]*p[4][1] + p[1][2]*p[2][0]*p[2][2]*p[4][1] - pow(p[1][2], 2)*p[3][0]*p[4][1] + 3*pow(p[5][2], 2)*p[3][0]*p[4][1] + p[1][0]*p[1][2]*p[3][2]*p[4][1] - p[1][2]*p[3][0]*p[3][2]*p[4][1] + 2*p[0][0]*p[1][1]*p[1][2]*p[4][2] - 2*p[1][1]*p[1][2]*p[2][0]*p[4][2] + 2*p[1][0]*p[1][2]*p[2][1]*p[4][2] - p[1][1]*p[2][0]*p[2][2]*p[4][2] + p[1][0]*p[2][1]*p[2][2]*p[4][2] + p[1][1]*p[1][2]*p[3][0]*p[4][2] - p[0][0]*p[1][2]*p[3][1]*p[4][2] - p[1][0]*p[1][2]*p[3][1]*p[4][2] + p[0][0]*p[1][1]*p[3][2]*p[4][2] + p[1][1]*p[3][0]*p[3][2]*p[4][2] - 2*p[0][0]*p[3][1]*p[3][2]*p[4][2] - p[1][0]*p[3][1]*p[3][2]*p[4][2] - p[1][2]*p[2][1]*p[4][0]*p[4][2] + p[1][1]*p[2][2]*p[4][0]*p[4][2] + 2*p[1][2]*p[3][1]*p[4][0]*p[4][2] - 2*p[1][1]*p[3][2]*p[4][0]*p[4][2] - p[0][0]*p[1][2]*p[4][1]*p[4][2] + p[1][2]*p[2][0]*p[4][1]*p[4][2] - p[1][0]*p[2][2]*p[4][1]*p[4][2] - 2*p[1][2]*p[3][0]*p[4][1]*p[4][2] + p[0][0]*p[3][2]*p[4][1]*p[4][2] + 2*p[1][0]*p[3][2]*p[4][1]*p[4][2] + 2*pow(p[2][2], 2)*p[1][1]*p[5][0] - 2*pow(p[4][2], 2)*p[1][1]*p[5][0] - pow(p[1][2], 2)*p[2][1]*p[5][0] + pow(p[3][2], 2)*p[2][1]*p[5][0] - pow(p[4][2], 2)*p[2][1]*p[5][0] + p[1][1]*p[1][2]*p[2][2]*p[5][0] - 2*p[1][2]*p[2][1]*p[2][2]*p[5][0] - pow(p[2][2], 2)*p[3][1]*p[5][0] + 3*pow(p[4][2], 2)*p[3][1]*p[5][0] + p[2][1]*p[2][2]*p[3][2]*p[5][0] - p[2][2]*p[3][1]*p[3][2]*p[5][0] + pow(p[1][2], 2)*p[4][1]*p[5][0] + pow(p[2][2], 2)*p[4][1]*p[5][0] - 3*pow(p[3][2], 2)*p[4][1]*p[5][0] + p[1][2]*p[2][2]*p[4][1]*p[5][0] - p[1][1]*p[1][2]*p[4][2]*p[5][0] - p[1][2]*p[2][1]*p[4][2]*p[5][0] - p[2][1]*p[2][2]*p[4][2]*p[5][0] + 3*p[3][1]*p[3][2]*p[4][2]*p[5][0] + 2*p[1][2]*p[4][1]*p[4][2]*p[5][0] + p[2][2]*p[4][1]*p[4][2]*p[5][0] - 3*p[3][2]*p[4][1]*p[4][2]*p[5][0] + 2*pow(p[2][2], 2)*p[0][0]*p[5][1] - 2*pow(p[3][2], 2)*p[0][0]*p[5][1] - 2*pow(p[2][2], 2)*p[1][0]*p[5][1] + 2*pow(p[4][2], 2)*p[1][0]*p[5][1] + pow(p[1][2], 2)*p[2][0]*p[5][1] - pow(p[3][2], 2)*p[2][0]*p[5][1] + pow(p[4][2], 2)*p[2][0]*p[5][1] - p[1][0]*p[1][2]*p[2][2]*p[5][1] + 2*p[1][2]*p[2][0]*p[2][2]*p[5][1] + pow(p[2][2], 2)*p[3][0]*p[5][1] - 3*pow(p[4][2], 2)*p[3][0]*p[5][1] - p[2][0]*p[2][2]*p[3][2]*p[5][1] + p[2][2]*p[3][0]*p[3][2]*p[5][1] - pow(p[1][2], 2)*p[4][0]*p[5][1] - pow(p[2][2], 2)*p[4][0]*p[5][1] + 3*pow(p[3][2], 2)*p[4][0]*p[5][1] - p[1][2]*p[2][2]*p[4][0]*p[5][1] + p[1][0]*p[1][2]*p[4][2]*p[5][1] + p[1][2]*p[2][0]*p[4][2]*p[5][1] + p[2][0]*p[2][2]*p[4][2]*p[5][1] - 3*p[3][0]*p[3][2]*p[4][2]*p[5][1] - 2*p[1][2]*p[4][0]*p[4][2]*p[5][1] - p[2][2]*p[4][0]*p[4][2]*p[5][1] + 3*p[3][2]*p[4][0]*p[4][2]*p[5][1] + pow(p[0][2], 2)*(2*p[2][1]*p[3][0] - 2*p[2][0]*p[3][1] + p[1][1]*(3*p[2][0] - 2*p[3][0] - p[4][0]) + p[3][1]*p[4][0] - p[3][0]*p[4][1] + p[1][0]*(-3*p[2][1] + 2*p[3][1] + p[4][1]) + p[2][1]*p[5][0] - p[3][1]*p[5][0] - p[2][0]*p[5][1] + p[3][0]*p[5][1]) - p[1][1]*p[1][2]*p[2][0]*p[5][2] + p[1][0]*p[1][2]*p[2][1]*p[5][2] - 2*p[1][1]*p[2][0]*p[2][2]*p[5][2] - 2*p[0][0]*p[2][1]*p[2][2]*p[5][2] + 2*p[1][0]*p[2][1]*p[2][2]*p[5][2] - p[2][1]*p[2][2]*p[3][0]*p[5][2] + p[0][0]*p[2][2]*p[3][1]*p[5][2] + p[2][0]*p[2][2]*p[3][1]*p[5][2] - p[0][0]*p[2][1]*p[3][2]*p[5][2] - p[2][1]*p[3][0]*p[3][2]*p[5][2] + 2*p[0][0]*p[3][1]*p[3][2]*p[5][2] + p[2][0]*p[3][1]*p[3][2]*p[5][2] + p[1][1]*p[1][2]*p[4][0]*p[5][2] + p[1][1]*p[2][2]*p[4][0]*p[5][2] + p[2][1]*p[2][2]*p[4][0]*p[5][2] - 3*p[3][1]*p[3][2]*p[4][0]*p[5][2] - p[1][0]*p[1][2]*p[4][1]*p[5][2] - p[1][0]*p[2][2]*p[4][1]*p[5][2] - p[2][0]*p[2][2]*p[4][1]*p[5][2] + 3*p[3][0]*p[3][2]*p[4][1]*p[5][2] - p[1][1]*p[2][0]*p[4][2]*p[5][2] + p[1][0]*p[2][1]*p[4][2]*p[5][2] + 2*p[1][1]*p[4][0]*p[4][2]*p[5][2] + p[2][1]*p[4][0]*p[4][2]*p[5][2] - 3*p[3][1]*p[4][0]*p[4][2]*p[5][2] - 2*p[1][0]*p[4][1]*p[4][2]*p[5][2] - p[2][0]*p[4][1]*p[4][2]*p[5][2] + 3*p[3][0]*p[4][1]*p[4][2]*p[5][2] - p[1][2]*p[2][1]*p[5][0]*p[5][2] + p[1][1]*p[2][2]*p[5][0]*p[5][2] - 2*p[2][2]*p[3][1]*p[5][0]*p[5][2] + 2*p[2][1]*p[3][2]*p[5][0]*p[5][2] + p[1][2]*p[4][1]*p[5][0]*p[5][2] + 2*p[2][2]*p[4][1]*p[5][0]*p[5][2] - 3*p[3][2]*p[4][1]*p[5][0]*p[5][2] - p[1][1]*p[4][2]*p[5][0]*p[5][2] - 2*p[2][1]*p[4][2]*p[5][0]*p[5][2] + 3*p[3][1]*p[4][2]*p[5][0]*p[5][2] + p[1][2]*p[2][0]*p[5][1]*p[5][2] + p[0][0]*p[2][2]*p[5][1]*p[5][2] - p[1][0]*p[2][2]*p[5][1]*p[5][2] + 2*p[2][2]*p[3][0]*p[5][1]*p[5][2] - p[0][0]*p[3][2]*p[5][1]*p[5][2] - 2*p[2][0]*p[3][2]*p[5][1]*p[5][2] - p[1][2]*p[4][0]*p[5][1]*p[5][2] - 2*p[2][2]*p[4][0]*p[5][1]*p[5][2] + 3*p[3][2]*p[4][0]*p[5][1]*p[5][2] + p[1][0]*p[4][2]*p[5][1]*p[5][2] + 2*p[2][0]*p[4][2]*p[5][1]*p[5][2] - 3*p[3][0]*p[4][2]*p[5][1]*p[5][2] + p[0][1]*(pow(p[3][2], 2)*p[2][0] + pow(p[5][2], 2)*p[2][0] - pow(p[2][2], 2)*p[3][0] + pow(p[4][2], 2)*p[3][0] - pow(p[5][2], 2)*p[3][0] + p[2][0]*p[2][2]*p[3][2] - p[2][2]*p[3][0]*p[3][2] - 2*pow(p[3][2], 2)*p[4][0] + pow(p[1][2], 2)*(-3*p[2][0] + p[3][0] + 2*p[4][0]) + 2*p[3][0]*p[3][2]*p[4][2] - p[3][2]*p[4][0]*p[4][2] + p[1][0]*(3*pow(p[2][2], 2) - pow(p[3][2], 2) - pow(p[4][2], 2) - p[3][2]*p[4][2]) + p[1][2]*(-3*p[2][0]*p[2][2] + p[3][0]*p[3][2] + p[1][0]*(3*p[2][2] - p[3][2] - 2*p[4][2]) + p[3][0]*p[4][2] + p[4][0]*p[4][2]) - 2*pow(p[2][2], 2)*p[5][0] + 2*pow(p[3][2], 2)*p[5][0] + 2*p[2][0]*p[2][2]*p[5][2] - p[2][2]*p[3][0]*p[5][2] + p[2][0]*p[3][2]*p[5][2] - 2*p[3][0]*p[3][2]*p[5][2] - p[2][2]*p[5][0]*p[5][2] + p[3][2]*p[5][0]*p[5][2]) + p[0][2]*(3*p[0][0]*p[1][2]*p[2][1] - 3*p[1][0]*p[1][2]*p[2][1] - 3*p[1][0]*p[2][1]*p[2][2] + p[2][1]*p[2][2]*p[3][0] - 2*p[0][0]*p[1][2]*p[3][1] + p[1][0]*p[1][2]*p[3][1] + 2*p[0][0]*p[2][2]*p[3][1] - p[2][0]*p[2][2]*p[3][1] - 2*p[0][0]*p[2][1]*p[3][2] + p[2][1]*p[3][0]*p[3][2] + p[1][0]*p[3][1]*p[3][2] - p[2][0]*p[3][1]*p[3][2] + p[1][2]*p[3][1]*p[4][0] + 2*p[3][1]*p[3][2]*p[4][0] - p[0][0]*p[1][2]*p[4][1] + 2*p[1][0]*p[1][2]*p[4][1] - p[1][2]*p[3][0]*p[4][1] + p[0][0]*p[3][2]*p[4][1] + p[1][0]*p[3][2]*p[4][1] - 2*p[3][0]*p[3][2]*p[4][1] - p[0][0]*p[3][1]*p[4][2] + p[3][1]*p[4][0]*p[4][2] + p[1][0]*p[4][1]*p[4][2] - p[3][0]*p[4][1]*p[4][2] + p[1][1]*(3*p[2][0]*p[2][2] - p[3][0]*p[3][2] + p[1][2]*(3*p[2][0] - p[3][0] - 2*p[4][0]) - p[3][2]*p[4][0] - p[4][0]*p[4][2] + p[0][0]*(-3*p[2][2] + 2*p[3][2] + p[4][2])) + 2*p[2][1]*p[2][2]*p[5][0] - p[2][2]*p[3][1]*p[5][0] + p[2][1]*p[3][2]*p[5][0] - 2*p[3][1]*p[3][2]*p[5][0] + p[0][0]*p[2][2]*p[5][1] - 2*p[2][0]*p[2][2]*p[5][1] + p[2][2]*p[3][0]*p[5][1] - p[0][0]*p[3][2]*p[5][1] - p[2][0]*p[3][2]*p[5][1] + 2*p[3][0]*p[3][2]*p[5][1] - p[0][0]*p[2][1]*p[5][2] + p[0][0]*p[3][1]*p[5][2] + p[2][1]*p[5][0]*p[5][2] - p[3][1]*p[5][0]*p[5][2] - p[2][0]*p[5][1]*p[5][2] + p[3][0]*p[5][1]*p[5][2] + p[0][1]*(-2*p[2][2]*p[3][0] + 2*p[2][0]*p[3][2] - p[3][2]*p[4][0] + p[1][2]*(-3*p[2][0] + 2*p[3][0] + p[4][0]) + p[1][0]*(3*p[2][2] - 2*p[3][2] - p[4][2]) + p[3][0]*p[4][2] - p[2][2]*p[5][0] + p[3][2]*p[5][0] + p[2][0]*p[5][2] - p[3][0]*p[5][2])))/72;
429  return result;
430  }
431 };
432 
433 // 3D Tetrahedron
434 
435 template<>
436 class Operator<VolumeIntegral, basis::Unitary, Tetrahedron>:
437  public OperatorBase<VolumeIntegral, basis::Unitary, Interval>
438 {
439 public:
440  static constexpr size_t dim = dimension(Tetrahedron);
441  static constexpr size_t operator_size = 1;
442  static constexpr size_t basis_size =
443  basis::Traits<basis::Unitary, dim>::function_size;
444  static constexpr size_t point_size = 4;
445 
446  using result_t = std::array<std::array<double, operator_size>, basis_size>;
447  using points_t = std::array<Wonton::Point<dim>, point_size>;
448 
449  static result_t apply(const points_t p) {
450  result_t result;
451  result[0][0] = (p[0][0]*p[1][2]*p[2][1] - p[0][0]*p[1][1]*p[2][2] - p[1][2]*p[2][1]*p[3][0] + p[1][1]*p[2][2]*p[3][0] - p[0][0]*p[1][2]*p[3][1] + p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] - p[1][0]*p[2][2]*p[3][1] + p[0][2]*(p[1][1]*(p[2][0] - p[3][0]) + p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[1][0]*(-p[2][1] + p[3][1])) + p[0][0]*p[1][1]*p[3][2] - p[1][1]*p[2][0]*p[3][2] - p[0][0]*p[2][1]*p[3][2] + p[1][0]*p[2][1]*p[3][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[1][2]*(-p[2][0] + p[3][0]) + p[1][0]*(p[2][2] - p[3][2]) + p[2][0]*p[3][2]))/6;
452  return result;
453  }
454 };
455 
456 template<>
457 class Operator<VolumeIntegral, basis::Linear, Tetrahedron>:
458  public OperatorBase<VolumeIntegral, basis::Linear, Interval>
459 {
460 public:
461  static constexpr size_t dim = dimension(Tetrahedron);
462  static constexpr size_t operator_size = 1;
463  static constexpr size_t basis_size =
464  basis::Traits<basis::Linear, dim>::function_size;
465  static constexpr size_t point_size = 4;
466 
467  using result_t = std::array<std::array<double, operator_size>, basis_size>;
468  using points_t = std::array<Wonton::Point<dim>, point_size>;
469 
470  static result_t apply(const points_t p) {
471  result_t result;
472  result[0][0] = (p[0][0]*p[1][2]*p[2][1] - p[0][0]*p[1][1]*p[2][2] - p[1][2]*p[2][1]*p[3][0] + p[1][1]*p[2][2]*p[3][0] - p[0][0]*p[1][2]*p[3][1] + p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] - p[1][0]*p[2][2]*p[3][1] + p[0][2]*(p[1][1]*(p[2][0] - p[3][0]) + p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[1][0]*(-p[2][1] + p[3][1])) + p[0][0]*p[1][1]*p[3][2] - p[1][1]*p[2][0]*p[3][2] - p[0][0]*p[2][1]*p[3][2] + p[1][0]*p[2][1]*p[3][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[1][2]*(-p[2][0] + p[3][0]) + p[1][0]*(p[2][2] - p[3][2]) + p[2][0]*p[3][2]))/6;
473  result[1][0] = ((p[0][0] + p[1][0] + p[2][0] + p[3][0])*(p[0][0]*p[1][2]*p[2][1] - p[0][0]*p[1][1]*p[2][2] - p[1][2]*p[2][1]*p[3][0] + p[1][1]*p[2][2]*p[3][0] - p[0][0]*p[1][2]*p[3][1] + p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] - p[1][0]*p[2][2]*p[3][1] + p[0][2]*(p[1][1]*(p[2][0] - p[3][0]) + p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[1][0]*(-p[2][1] + p[3][1])) + p[0][0]*p[1][1]*p[3][2] - p[1][1]*p[2][0]*p[3][2] - p[0][0]*p[2][1]*p[3][2] + p[1][0]*p[2][1]*p[3][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[1][2]*(-p[2][0] + p[3][0]) + p[1][0]*(p[2][2] - p[3][2]) + p[2][0]*p[3][2])))/24;
474  result[2][0] = ((p[0][1] + p[1][1] + p[2][1] + p[3][1])*(p[0][0]*p[1][2]*p[2][1] - p[0][0]*p[1][1]*p[2][2] - p[1][2]*p[2][1]*p[3][0] + p[1][1]*p[2][2]*p[3][0] - p[0][0]*p[1][2]*p[3][1] + p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] - p[1][0]*p[2][2]*p[3][1] + p[0][2]*(p[1][1]*(p[2][0] - p[3][0]) + p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[1][0]*(-p[2][1] + p[3][1])) + p[0][0]*p[1][1]*p[3][2] - p[1][1]*p[2][0]*p[3][2] - p[0][0]*p[2][1]*p[3][2] + p[1][0]*p[2][1]*p[3][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[1][2]*(-p[2][0] + p[3][0]) + p[1][0]*(p[2][2] - p[3][2]) + p[2][0]*p[3][2])))/24;
475  result[3][0] = ((p[0][2] + p[1][2] + p[2][2] + p[3][2])*(p[0][0]*p[1][2]*p[2][1] - p[0][0]*p[1][1]*p[2][2] - p[1][2]*p[2][1]*p[3][0] + p[1][1]*p[2][2]*p[3][0] - p[0][0]*p[1][2]*p[3][1] + p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] - p[1][0]*p[2][2]*p[3][1] + p[0][2]*(p[1][1]*(p[2][0] - p[3][0]) + p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[1][0]*(-p[2][1] + p[3][1])) + p[0][0]*p[1][1]*p[3][2] - p[1][1]*p[2][0]*p[3][2] - p[0][0]*p[2][1]*p[3][2] + p[1][0]*p[2][1]*p[3][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[1][2]*(-p[2][0] + p[3][0]) + p[1][0]*(p[2][2] - p[3][2]) + p[2][0]*p[3][2])))/24;
476  return result;
477  }
478 };
479 
480 template<>
481 class Operator<VolumeIntegral, basis::Quadratic, Tetrahedron>:
482  public OperatorBase<VolumeIntegral, basis::Quadratic, Interval>
483 {
484 public:
485  static constexpr size_t dim = dimension(Tetrahedron);
486  static constexpr size_t operator_size = 1;
487  static constexpr size_t basis_size =
488  basis::Traits<basis::Quadratic, dim>::function_size;
489  static constexpr size_t point_size = 4;
490 
491  using result_t = std::array<std::array<double, operator_size>, basis_size>;
492  using points_t = std::array<Wonton::Point<dim>, point_size>;
493 
494  static result_t apply(const points_t p) {
495  result_t result;
496  result[0][0] = (p[0][0]*p[1][2]*p[2][1] - p[0][0]*p[1][1]*p[2][2] - p[1][2]*p[2][1]*p[3][0] + p[1][1]*p[2][2]*p[3][0] - p[0][0]*p[1][2]*p[3][1] + p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] - p[1][0]*p[2][2]*p[3][1] + p[0][2]*(p[1][1]*(p[2][0] - p[3][0]) + p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[1][0]*(-p[2][1] + p[3][1])) + p[0][0]*p[1][1]*p[3][2] - p[1][1]*p[2][0]*p[3][2] - p[0][0]*p[2][1]*p[3][2] + p[1][0]*p[2][1]*p[3][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[1][2]*(-p[2][0] + p[3][0]) + p[1][0]*(p[2][2] - p[3][2]) + p[2][0]*p[3][2]))/6;
497  result[1][0] = ((p[0][0] + p[1][0] + p[2][0] + p[3][0])*(p[0][0]*p[1][2]*p[2][1] - p[0][0]*p[1][1]*p[2][2] - p[1][2]*p[2][1]*p[3][0] + p[1][1]*p[2][2]*p[3][0] - p[0][0]*p[1][2]*p[3][1] + p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] - p[1][0]*p[2][2]*p[3][1] + p[0][2]*(p[1][1]*(p[2][0] - p[3][0]) + p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[1][0]*(-p[2][1] + p[3][1])) + p[0][0]*p[1][1]*p[3][2] - p[1][1]*p[2][0]*p[3][2] - p[0][0]*p[2][1]*p[3][2] + p[1][0]*p[2][1]*p[3][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[1][2]*(-p[2][0] + p[3][0]) + p[1][0]*(p[2][2] - p[3][2]) + p[2][0]*p[3][2])))/24;
498  result[2][0] = ((p[0][1] + p[1][1] + p[2][1] + p[3][1])*(p[0][0]*p[1][2]*p[2][1] - p[0][0]*p[1][1]*p[2][2] - p[1][2]*p[2][1]*p[3][0] + p[1][1]*p[2][2]*p[3][0] - p[0][0]*p[1][2]*p[3][1] + p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] - p[1][0]*p[2][2]*p[3][1] + p[0][2]*(p[1][1]*(p[2][0] - p[3][0]) + p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[1][0]*(-p[2][1] + p[3][1])) + p[0][0]*p[1][1]*p[3][2] - p[1][1]*p[2][0]*p[3][2] - p[0][0]*p[2][1]*p[3][2] + p[1][0]*p[2][1]*p[3][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[1][2]*(-p[2][0] + p[3][0]) + p[1][0]*(p[2][2] - p[3][2]) + p[2][0]*p[3][2])))/24;
499  result[3][0] = ((p[0][2] + p[1][2] + p[2][2] + p[3][2])*(p[0][0]*p[1][2]*p[2][1] - p[0][0]*p[1][1]*p[2][2] - p[1][2]*p[2][1]*p[3][0] + p[1][1]*p[2][2]*p[3][0] - p[0][0]*p[1][2]*p[3][1] + p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] - p[1][0]*p[2][2]*p[3][1] + p[0][2]*(p[1][1]*(p[2][0] - p[3][0]) + p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[1][0]*(-p[2][1] + p[3][1])) + p[0][0]*p[1][1]*p[3][2] - p[1][1]*p[2][0]*p[3][2] - p[0][0]*p[2][1]*p[3][2] + p[1][0]*p[2][1]*p[3][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[1][2]*(-p[2][0] + p[3][0]) + p[1][0]*(p[2][2] - p[3][2]) + p[2][0]*p[3][2])))/24;
500  result[4][0] = ((pow(p[0][0], 2) + pow(p[1][0], 2) + pow(p[2][0], 2) + pow(p[3][0], 2) + p[2][0]*p[3][0] + p[1][0]*(p[2][0] + p[3][0]) + p[0][0]*(p[1][0] + p[2][0] + p[3][0]))*(p[0][0]*p[1][2]*p[2][1] - p[0][0]*p[1][1]*p[2][2] - p[1][2]*p[2][1]*p[3][0] + p[1][1]*p[2][2]*p[3][0] - p[0][0]*p[1][2]*p[3][1] + p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] - p[1][0]*p[2][2]*p[3][1] + p[0][2]*(p[1][1]*(p[2][0] - p[3][0]) + p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[1][0]*(-p[2][1] + p[3][1])) + p[0][0]*p[1][1]*p[3][2] - p[1][1]*p[2][0]*p[3][2] - p[0][0]*p[2][1]*p[3][2] + p[1][0]*p[2][1]*p[3][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[1][2]*(-p[2][0] + p[3][0]) + p[1][0]*(p[2][2] - p[3][2]) + p[2][0]*p[3][2])))/120;
501  result[5][0] = -((2*p[1][0]*p[1][1] + p[1][1]*p[2][0] + p[1][0]*p[2][1] + 2*p[2][0]*p[2][1] + p[1][1]*p[3][0] + p[2][1]*p[3][0] + p[0][1]*(p[1][0] + p[2][0] + p[3][0]) + p[1][0]*p[3][1] + p[2][0]*p[3][1] + 2*p[3][0]*p[3][1] + p[0][0]*(2*p[0][1] + p[1][1] + p[2][1] + p[3][1]))*(-(p[0][0]*p[1][2]*p[2][1]) + p[0][0]*p[1][1]*p[2][2] + p[1][2]*p[2][1]*p[3][0] - p[1][1]*p[2][2]*p[3][0] + p[0][0]*p[1][2]*p[3][1] - p[1][2]*p[2][0]*p[3][1] - p[0][0]*p[2][2]*p[3][1] + p[1][0]*p[2][2]*p[3][1] + p[0][2]*(-(p[2][1]*p[3][0]) + p[1][1]*(-p[2][0] + p[3][0]) + p[1][0]*(p[2][1] - p[3][1]) + p[2][0]*p[3][1]) - p[0][0]*p[1][1]*p[3][2] + p[1][1]*p[2][0]*p[3][2] + p[0][0]*p[2][1]*p[3][2] - p[1][0]*p[2][1]*p[3][2] + p[0][1]*(p[1][2]*(p[2][0] - p[3][0]) + p[2][2]*p[3][0] - p[2][0]*p[3][2] + p[1][0]*(-p[2][2] + p[3][2]))))/120;
502  result[6][0] = ((2*p[1][0]*p[1][2] + p[1][2]*p[2][0] + p[1][0]*p[2][2] + 2*p[2][0]*p[2][2] + p[1][2]*p[3][0] + p[2][2]*p[3][0] + p[0][2]*(p[1][0] + p[2][0] + p[3][0]) + p[1][0]*p[3][2] + p[2][0]*p[3][2] + 2*p[3][0]*p[3][2] + p[0][0]*(2*p[0][2] + p[1][2] + p[2][2] + p[3][2]))*(p[0][0]*p[1][2]*p[2][1] - p[0][0]*p[1][1]*p[2][2] - p[1][2]*p[2][1]*p[3][0] + p[1][1]*p[2][2]*p[3][0] - p[0][0]*p[1][2]*p[3][1] + p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] - p[1][0]*p[2][2]*p[3][1] + p[0][2]*(p[1][1]*(p[2][0] - p[3][0]) + p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[1][0]*(-p[2][1] + p[3][1])) + p[0][0]*p[1][1]*p[3][2] - p[1][1]*p[2][0]*p[3][2] - p[0][0]*p[2][1]*p[3][2] + p[1][0]*p[2][1]*p[3][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[1][2]*(-p[2][0] + p[3][0]) + p[1][0]*(p[2][2] - p[3][2]) + p[2][0]*p[3][2])))/120;
503  result[7][0] = ((pow(p[0][1], 2) + pow(p[1][1], 2) + pow(p[2][1], 2) + pow(p[3][1], 2) + p[2][1]*p[3][1] + p[1][1]*(p[2][1] + p[3][1]) + p[0][1]*(p[1][1] + p[2][1] + p[3][1]))*(p[0][0]*p[1][2]*p[2][1] - p[0][0]*p[1][1]*p[2][2] - p[1][2]*p[2][1]*p[3][0] + p[1][1]*p[2][2]*p[3][0] - p[0][0]*p[1][2]*p[3][1] + p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] - p[1][0]*p[2][2]*p[3][1] + p[0][2]*(p[1][1]*(p[2][0] - p[3][0]) + p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[1][0]*(-p[2][1] + p[3][1])) + p[0][0]*p[1][1]*p[3][2] - p[1][1]*p[2][0]*p[3][2] - p[0][0]*p[2][1]*p[3][2] + p[1][0]*p[2][1]*p[3][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[1][2]*(-p[2][0] + p[3][0]) + p[1][0]*(p[2][2] - p[3][2]) + p[2][0]*p[3][2])))/120;
504  result[8][0] = ((2*p[1][1]*p[1][2] + p[1][2]*p[2][1] + p[1][1]*p[2][2] + 2*p[2][1]*p[2][2] + p[1][2]*p[3][1] + p[2][2]*p[3][1] + p[0][2]*(p[1][1] + p[2][1] + p[3][1]) + p[1][1]*p[3][2] + p[2][1]*p[3][2] + 2*p[3][1]*p[3][2] + p[0][1]*(2*p[0][2] + p[1][2] + p[2][2] + p[3][2]))*(p[0][0]*p[1][2]*p[2][1] - p[0][0]*p[1][1]*p[2][2] - p[1][2]*p[2][1]*p[3][0] + p[1][1]*p[2][2]*p[3][0] - p[0][0]*p[1][2]*p[3][1] + p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] - p[1][0]*p[2][2]*p[3][1] + p[0][2]*(p[1][1]*(p[2][0] - p[3][0]) + p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[1][0]*(-p[2][1] + p[3][1])) + p[0][0]*p[1][1]*p[3][2] - p[1][1]*p[2][0]*p[3][2] - p[0][0]*p[2][1]*p[3][2] + p[1][0]*p[2][1]*p[3][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[1][2]*(-p[2][0] + p[3][0]) + p[1][0]*(p[2][2] - p[3][2]) + p[2][0]*p[3][2])))/120;
505  result[9][0] = ((pow(p[0][2], 2) + pow(p[1][2], 2) + pow(p[2][2], 2) + pow(p[3][2], 2) + p[2][2]*p[3][2] + p[1][2]*(p[2][2] + p[3][2]) + p[0][2]*(p[1][2] + p[2][2] + p[3][2]))*(p[0][0]*p[1][2]*p[2][1] - p[0][0]*p[1][1]*p[2][2] - p[1][2]*p[2][1]*p[3][0] + p[1][1]*p[2][2]*p[3][0] - p[0][0]*p[1][2]*p[3][1] + p[1][2]*p[2][0]*p[3][1] + p[0][0]*p[2][2]*p[3][1] - p[1][0]*p[2][2]*p[3][1] + p[0][2]*(p[1][1]*(p[2][0] - p[3][0]) + p[2][1]*p[3][0] - p[2][0]*p[3][1] + p[1][0]*(-p[2][1] + p[3][1])) + p[0][0]*p[1][1]*p[3][2] - p[1][1]*p[2][0]*p[3][2] - p[0][0]*p[2][1]*p[3][2] + p[1][0]*p[2][1]*p[3][2] + p[0][1]*(-(p[2][2]*p[3][0]) + p[1][2]*(-p[2][0] + p[3][0]) + p[1][0]*(p[2][2] - p[3][2]) + p[2][0]*p[3][2])))/120;
506  return result;
507  }
508 };
509 
511 // Dynamic size information
513 inline std::array<size_t,3> size_info(Type type, basis::Type basis, Domain domain) {
514 
515  if (type == VolumeIntegral) {
516  switch (basis) {
517  case basis::Unitary:
518  switch (domain) {
519  case Interval: {
520  using OP = Operator<VolumeIntegral, basis::Unitary, Interval>;
521  return { OP::operator_size, OP::basis_size, OP::point_size };
522  }
523  case Quadrilateral: {
524  using OP = Operator<VolumeIntegral, basis::Unitary, Quadrilateral>;
525  return { OP::operator_size, OP::basis_size, OP::point_size };
526  }
527  case Triangle: {
528  using OP = Operator<VolumeIntegral, basis::Unitary, Triangle>;
529  return { OP::operator_size, OP::basis_size, OP::point_size };
530  }
531  case Hexahedron: {
532  using OP = Operator<VolumeIntegral, basis::Unitary, Hexahedron>;
533  return { OP::operator_size, OP::basis_size, OP::point_size };
534  }
535  case Wedge: {
536  using OP = Operator<VolumeIntegral, basis::Unitary, Wedge>;
537  return { OP::operator_size, OP::basis_size, OP::point_size };
538  }
539  case Tetrahedron: {
540  using OP = Operator<VolumeIntegral, basis::Unitary, Tetrahedron>;
541  return { OP::operator_size, OP::basis_size, OP::point_size };
542  }
543  default: break;
544  }
545  break;
546  case basis::Linear:
547  switch (domain) {
548  case Interval: {
549  using OP = Operator<VolumeIntegral, basis::Linear, Interval>;
550  return { OP::operator_size, OP::basis_size, OP::point_size };
551  }
552  case Quadrilateral: {
553  using OP = Operator<VolumeIntegral, basis::Linear, Quadrilateral>;
554  return { OP::operator_size, OP::basis_size, OP::point_size };
555  }
556  case Triangle: {
557  using OP = Operator<VolumeIntegral, basis::Linear, Triangle>;
558  return { OP::operator_size, OP::basis_size, OP::point_size };
559  }
560  case Hexahedron: {
561  using OP = Operator<VolumeIntegral, basis::Linear, Hexahedron>;
562  return { OP::operator_size, OP::basis_size, OP::point_size };
563  }
564  case Wedge: {
565  using OP = Operator<VolumeIntegral, basis::Linear, Wedge>;
566  return { OP::operator_size, OP::basis_size, OP::point_size };
567  }
568  case Tetrahedron: {
569  using OP = Operator<VolumeIntegral, basis::Linear, Tetrahedron>;
570  return { OP::operator_size, OP::basis_size, OP::point_size };
571  }
572  default: break;
573  }
574  break;
575  case basis::Quadratic:
576  switch (domain) {
577  case Interval: {
578  using OP = Operator<VolumeIntegral, basis::Quadratic, Interval>;
579  return { OP::operator_size, OP::basis_size, OP::point_size };
580  }
581  case Quadrilateral: {
582  using OP = Operator<VolumeIntegral, basis::Quadratic, Quadrilateral>;
583  return { OP::operator_size, OP::basis_size, OP::point_size };
584  }
585  case Triangle: {
586  using OP = Operator<VolumeIntegral, basis::Quadratic, Triangle>;
587  return { OP::operator_size, OP::basis_size, OP::point_size };
588  }
589  case Tetrahedron: {
590  using OP = Operator<VolumeIntegral, basis::Quadratic, Tetrahedron>;
591  return { OP::operator_size, OP::basis_size, OP::point_size };
592  }
593  default: break;
594  }
595  break;
596  default: break;
597  }
598  }
599  return { 0, 0, 0 };
600 }
601 
603 // Helper template functions
605 
606 template<class OP>
607 inline void copy_points(const std::vector<Wonton::Point<OP::dim>>& points,
608  typename OP::points_t& apts) {
609  int const num_points = OP::point_size;
610  for (int i = 0; i < num_points; i++) {
611  apts[i] = points[i];
612  }
613 }
614 
615 template<class OP>
616 inline Wonton::Point<OP::dim> centroid(const typename OP::points_t& apts) {
617 
618  int const dim = OP::dim;
619  int const num_points = OP::point_size;
620  Wonton::Point<dim> cent {};
621 
622  for (int i = 0; i < num_points; i++) {
623  for (int j = 0; j < dim; j++) {
624  cent[j] += apts[i][j];
625  }
626  }
627 
628  for (int j = 0; j < dim; j++) {
629  cent[j] /= num_points;
630  }
631  return cent;
632 }
633 
634 template<class OP>
635 inline void shift_points(const Wonton::Point<OP::dim> c, typename OP::points_t& apts) {
636 
637  int const num_points = OP::point_size;
638  int const dim = OP::dim;
639 
640  for (int i = 0; i < num_points; i++) {
641  for (int j = 0; j < dim; j++) {
642  apts[i][j] = apts[i][j] - c[j];
643  }
644  }
645 }
646 
647 template<class OP>
648 inline void resize_result(std::vector<std::vector<double>>& result) {
649  result.resize(OP::basis_size);
650  for (auto iter = result.begin(); iter != result.end(); iter++) {
651  iter->resize(OP::operator_size);
652  }
653 }
654 
655 template<class OP>
656 inline void copy_result(const typename OP::result_t& ares,
657  std::vector<std::vector<double>>& result) {
658  int const num_basis = OP::basis_size;
659  int const num_operators = OP::operator_size;
660 
661  for (int i = 0; i < num_basis; i++) {
662  for (int j = 0; j < num_operators; j++) {
663  result[i][j] = ares[i][j];
664  }
665  }
666 }
667 
668 template<class OP>
669 inline void get_result(const std::vector<Wonton::Point<OP::dim>>& points,
670  std::vector<std::vector<double>>& result, const bool center = true) {
671  resize_result<OP>(result);
672  typename OP::points_t apts;
673  copy_points<OP>(points, apts);
674 
675  if (center) {
676  Wonton::Point<OP::dim> c = centroid<OP>(apts);
677  shift_points<OP>(c, apts);
678  auto tf = basis::transfactor<OP::dim>(OP::basis, c);
679  typename OP::result_t ares = OP::apply(apts);
680 
681  int const num_basis = OP::basis_size;
682  int const num_operators = OP::operator_size;
683 
684  for (int i = 0; i < num_basis; i++) {
685  for (int j = 0; j < num_operators; j++) {
686  result[i][j] = 0.;
687  }
688  }
689 
690  for (int j = 0; j < num_operators; j++) {
691  for (int i = 0; i < num_basis; i++) {
692  for (int k = 0; k < num_basis; k++) {
693  result[i][j] += tf[i][k] * ares[k][j];
694  }
695  }
696  }
697  } else {
698  typename OP::result_t ares = OP::apply(apts);
699  copy_result<OP>(ares, result);
700  }
701 }
702 
703 template<class OP>
704 inline std::vector<std::vector<double>>
705 get_result(const std::vector<Wonton::Point<OP::dim>>& points, bool center = true) {
706  std::vector<std::vector<double>> result;
707  get_result<OP>(points, result, center);
708  return result;
709 }
710 
712 // Dynamic Accessor
714 
715 template<size_t dim>
716 void apply(const Type type, const basis::Type basis_type,
717  const Domain domain_type, const std::vector<Wonton::Point<dim>> &points,
718  std::vector<std::vector<double>> &result);
719 
720 template<>
721 inline
722 void apply<1>(const Type type, const basis::Type basis_type,
723  const Domain domain_type, const std::vector<Wonton::Point<1>>& points,
724  std::vector<std::vector<double>>& result) {
725  bool center = true;
726  switch (type) {
727  case VolumeIntegral:
728  switch (domain_type) {
729  case Interval:
730  switch (basis_type) {
731  case basis::Unitary: {
732  using OP = Operator<VolumeIntegral, basis::Unitary, Interval>;
733  get_result<OP>(points, result, center);
734  }
735  break;
736  case basis::Linear: {
737  using OP = Operator<VolumeIntegral, basis::Linear, Interval>;
738  get_result<OP>(points, result, center);
739  }
740  break;
741  case basis::Quadratic: {
742  using OP = Operator<VolumeIntegral, basis::Quadratic, Interval>;
743  get_result<OP>(points, result, center);
744  }
745  break;
746  default:
747  throw (std::runtime_error("invalid basis"));
748  }
749  break;
750  default:
751  throw (std::runtime_error("invalid domain"));
752  }
753  break;
754  default:
755  throw (std::runtime_error("invalid operator"));
756  }
757 }
758 
759 template<>
760 inline
761 void apply<2>(const Type type, const basis::Type basis_type,
762  const Domain domain_type, const std::vector<Wonton::Point<2>>& points,
763  std::vector<std::vector<double>>& result) {
764  bool center = true;
765  switch (type) {
766  case VolumeIntegral:
767  switch (domain_type) {
768  case Quadrilateral:
769  switch (basis_type) {
770  case basis::Unitary: {
771  using OP = Operator<VolumeIntegral, basis::Unitary, Quadrilateral>;
772  get_result<OP>(points, result, center);
773  }
774  break;
775  case basis::Linear: {
776  using OP = Operator<VolumeIntegral, basis::Linear, Quadrilateral>;
777  get_result<OP>(points, result, center);
778  }
779  break;
780  case basis::Quadratic: {
781  using OP = Operator<VolumeIntegral, basis::Quadratic, Quadrilateral>;
782  get_result<OP>(points, result, center);
783  }
784  break;
785  default:
786  throw (std::runtime_error("invalid basis"));
787  }
788  break;
789  case Triangle:
790  switch (basis_type) {
791  case basis::Unitary: {
792  using OP = Operator<VolumeIntegral, basis::Unitary, Triangle>;
793  get_result<OP>(points, result, center);
794  }
795  break;
796  case basis::Linear: {
797  using OP = Operator<VolumeIntegral, basis::Linear, Triangle>;
798  get_result<OP>(points, result, center);
799  }
800  break;
801  case basis::Quadratic: {
802  using OP = Operator<VolumeIntegral, basis::Quadratic, Triangle>;
803  get_result<OP>(points, result, center);
804  }
805  break;
806  default:
807  throw (std::runtime_error("invalid basis"));
808  }
809  break;
810  default:
811  throw (std::runtime_error("invalid domain"));
812  }
813  break;
814  default:
815  throw (std::runtime_error("invalid operator"));
816  }
817 }
818 
819 template<>
820 inline
821 void apply<3>(const Type type, const basis::Type basis_type,
822  const Domain domain_type, const std::vector<Wonton::Point<3>>& points,
823  std::vector<std::vector<double>>& result) {
824  bool center = true;
825  switch (type) {
826  case VolumeIntegral:
827  switch (domain_type) {
828  case Hexahedron:
829  switch (basis_type) {
830  case basis::Unitary: {
831  using OP = Operator<VolumeIntegral, basis::Unitary, Hexahedron>;
832  get_result<OP>(points, result, center);
833  }
834  break;
835  case basis::Linear: {
836  using OP = Operator<VolumeIntegral, basis::Linear, Hexahedron>;
837  get_result<OP>(points, result, center);
838  }
839  break;
840  default:
841  throw (std::runtime_error("invalid basis"));
842  }
843  break;
844  case Wedge:
845  switch (basis_type) {
846  case basis::Unitary: {
847  using OP = Operator<VolumeIntegral, basis::Unitary, Wedge>;
848  get_result<OP>(points, result, center);
849  }
850  break;
851  case basis::Linear: {
852  using OP = Operator<VolumeIntegral, basis::Linear, Wedge>;
853  get_result<OP>(points, result, center);
854  }
855  break;
856  default:
857  throw (std::runtime_error("invalid basis"));
858  }
859  break;
860  case Tetrahedron:
861  switch (basis_type) {
862  case basis::Unitary: {
863  using OP = Operator<VolumeIntegral, basis::Unitary, Tetrahedron>;
864  get_result<OP>(points, result, center);
865  }
866  break;
867  case basis::Linear: {
868  using OP = Operator<VolumeIntegral, basis::Linear, Tetrahedron>;
869  get_result<OP>(points, result, center);
870  }
871  break;
872  case basis::Quadratic: {
873  using OP = Operator<VolumeIntegral, basis::Quadratic, Tetrahedron>;
874  get_result<OP>(points, result, center);
875  }
876  break;
877  default:
878  throw (std::runtime_error("invalid basis"));
879  }
880  break;
881  default:
882  throw (std::runtime_error("invalid domain"));
883  }
884  break;
885  default:
886  throw (std::runtime_error("invalid operator"));
887  }
888 }
889 
890 }}} // namespace Portage::Meshfree::Operator
891 
892 #endif // OPERATOR_H_INC_
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Wedge >::apply
static result_t apply(const points_t p)
Definition: operator.h:423
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Tetrahedron >::result_t
std::array< std::array< double, operator_size >, basis_size > result_t
Definition: operator.h:467
Portage::Meshfree::oper::Interval
Definition: operator.h:26
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Hexahedron >::apply
static result_t apply(const points_t p)
Definition: operator.h:355
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Quadratic, Tetrahedron >::points_t
std::array< Wonton::Point< dim >, point_size > points_t
Definition: operator.h:492
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Hexahedron >::result_t
std::array< std::array< double, operator_size >, basis_size > result_t
Definition: operator.h:373
Portage::Meshfree::oper::DomainTraits
Definition: operator.h:32
Portage::Meshfree::oper::LastDomain
Definition: operator.h:28
Portage::Meshfree::oper::DomainTraits::dimension
static const size_t dimension
Definition: operator.h:32
Portage::Meshfree::basis::Traits
Definition: basis.h:29
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Interval >::points_t
std::array< Wonton::Point< dim >, point_size > points_t
Definition: operator.h:162
Portage::vector
std::vector< T > vector
Definition: portage.h:238
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Quadrilateral >::points_t
std::array< Wonton::Point< dim >, point_size > points_t
Definition: operator.h:302
Portage::Meshfree::oper::apply< 3 >
void apply< 3 >(const Type type, const basis::Type basis_type, const Domain domain_type, const std::vector< Wonton::Point< 3 >> &points, std::vector< std::vector< double >> &result)
Definition: operator.h:821
Portage::Meshfree::oper::Operator::points_t
std::array< Wonton::Point< dim >, point_size > points_t
Definition: operator.h:116
Portage::Meshfree::oper::Circle
Definition: operator.h:27
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Wedge >::result_t
std::array< std::array< double, operator_size >, basis_size > result_t
Definition: operator.h:399
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Wedge >::apply
static result_t apply(const points_t p)
Definition: operator.h:402
Portage::Meshfree::oper::centroid
Wonton::Point< OP::dim > centroid(const typename OP::points_t &apts)
Definition: operator.h:616
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Quadrilateral >::points_t
std::array< Wonton::Point< dim >, point_size > points_t
Definition: operator.h:281
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Quadratic, Quadrilateral >::result_t
std::array< std::array< double, operator_size >, basis_size > result_t
Definition: operator.h:324
Portage::Meshfree::oper::Triangle
Definition: operator.h:27
Portage::Meshfree::oper::copy_result
void copy_result(const typename OP::result_t &ares, std::vector< std::vector< double >> &result)
Definition: operator.h:656
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Interval >::result_t
std::array< std::array< double, operator_size >, basis_size > result_t
Definition: operator.h:161
Portage::Meshfree::oper::Operator
Definition: operator.h:108
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Quadratic, Tetrahedron >::result_t
std::array< std::array< double, operator_size >, basis_size > result_t
Definition: operator.h:491
Portage::Meshfree::oper::Domain
Domain
Definition: operator.h:26
Portage::Meshfree::basis::Linear
Definition: basis.h:20
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Quadratic, Quadrilateral >::points_t
std::array< Wonton::Point< dim >, point_size > points_t
Definition: operator.h:325
Portage::Meshfree::oper::domain_from_points
Domain domain_from_points(std::vector< Wonton::Point< dim >> const &points)
Definition: operator.h:68
Portage::Meshfree::oper::Operator::apply
static result_t apply(const points_t p)
Definition: operator.h:118
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Tetrahedron >::apply
static result_t apply(const points_t p)
Definition: operator.h:449
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Quadrilateral >::apply
static result_t apply(const points_t p)
Definition: operator.h:283
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Tetrahedron >::points_t
std::array< Wonton::Point< dim >, point_size > points_t
Definition: operator.h:468
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Hexahedron >::points_t
std::array< Wonton::Point< dim >, point_size > points_t
Definition: operator.h:353
Portage::Meshfree::oper::copy_points
void copy_points(const std::vector< Wonton::Point< OP::dim >> &points, typename OP::points_t &apts)
Definition: operator.h:607
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Triangle >::result_t
std::array< std::array< double, operator_size >, basis_size > result_t
Definition: operator.h:229
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Tetrahedron >::apply
static result_t apply(const points_t p)
Definition: operator.h:470
Portage::Meshfree::oper::apply
void apply(const Type type, const basis::Type basis_type, const Domain domain_type, const std::vector< Wonton::Point< dim >> &points, std::vector< std::vector< double >> &result)
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Wedge >::result_t
std::array< std::array< double, operator_size >, basis_size > result_t
Definition: operator.h:420
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Quadrilateral >::result_t
std::array< std::array< double, operator_size >, basis_size > result_t
Definition: operator.h:301
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Interval >::apply
static result_t apply(const points_t p)
Definition: operator.h:142
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Quadratic, Interval >::apply
static result_t apply(const points_t p)
Definition: operator.h:186
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Quadratic, Triangle >::result_t
std::array< std::array< double, operator_size >, basis_size > result_t
Definition: operator.h:252
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Interval >::points_t
std::array< Wonton::Point< dim >, point_size > points_t
Definition: operator.h:140
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Tetrahedron >::points_t
std::array< Wonton::Point< dim >, point_size > points_t
Definition: operator.h:447
Portage::Meshfree::oper::LastOperator
Definition: operator.h:30
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Quadratic, Triangle >::apply
static result_t apply(const points_t p)
Definition: operator.h:255
Portage::Meshfree::oper::get_result
void get_result(const std::vector< Wonton::Point< OP::dim >> &points, std::vector< std::vector< double >> &result, const bool center=true)
Definition: operator.h:669
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Quadratic, Interval >::result_t
std::array< std::array< double, operator_size >, basis_size > result_t
Definition: operator.h:183
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Wedge >::points_t
std::array< Wonton::Point< dim >, point_size > points_t
Definition: operator.h:400
Portage::Meshfree::oper::apply< 1 >
void apply< 1 >(const Type type, const basis::Type basis_type, const Domain domain_type, const std::vector< Wonton::Point< 1 >> &points, std::vector< std::vector< double >> &result)
Definition: operator.h:722
Portage::Meshfree::basis::Type
Type
Definition: basis.h:20
Portage::Meshfree::basis::Unitary
Definition: basis.h:20
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Quadratic, Triangle >::points_t
std::array< Wonton::Point< dim >, point_size > points_t
Definition: operator.h:253
Portage::Meshfree::oper::OperatorBase
Definition: operator.h:96
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Triangle >::apply
static result_t apply(const points_t p)
Definition: operator.h:232
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Interval >::result_t
std::array< std::array< double, operator_size >, basis_size > result_t
Definition: operator.h:139
Portage::Meshfree::oper::Type
Type
Definition: operator.h:30
Portage::Meshfree::basis::transfactor
Traits< type, dim >::matrix_t transfactor(const Wonton::Point< dim > &c)
Definition: basis.h:133
Portage::Meshfree::reference::OP
oper::Operator< type, basis_type, domain_type > OP
Definition: operator_references.h:25
Portage::Meshfree::oper::apply< 2 >
void apply< 2 >(const Type type, const basis::Type basis_type, const Domain domain_type, const std::vector< Wonton::Point< 2 >> &points, std::vector< std::vector< double >> &result)
Definition: operator.h:761
Portage::Meshfree::oper::Tetrahedron
Definition: operator.h:28
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Quadrilateral >::result_t
std::array< std::array< double, operator_size >, basis_size > result_t
Definition: operator.h:280
Portage::Meshfree::oper::resize_result
void resize_result(std::vector< std::vector< double >> &result)
Definition: operator.h:648
Portage
Definition: coredriver.h:42
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Quadratic, Quadrilateral >::apply
static result_t apply(const points_t p)
Definition: operator.h:327
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Quadratic, Interval >::points_t
std::array< Wonton::Point< dim >, point_size > points_t
Definition: operator.h:184
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Hexahedron >::apply
static result_t apply(const points_t p)
Definition: operator.h:376
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Triangle >::result_t
std::array< std::array< double, operator_size >, basis_size > result_t
Definition: operator.h:208
Portage::Meshfree::oper::Wedge
Definition: operator.h:28
Portage::Meshfree::oper::shift_points
void shift_points(const Wonton::Point< OP::dim > c, typename OP::points_t &apts)
Definition: operator.h:635
Portage::Meshfree::reference::points
constexpr std::vector< Point< oper::dimension(domain)> > points()
Definition: operator_references.h:29
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Quadratic, Tetrahedron >::apply
static result_t apply(const points_t p)
Definition: operator.h:494
Portage::Meshfree::oper::Operator::result_t
std::array< std::array< double, operator_size >, basis_size > result_t
Definition: operator.h:115
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Triangle >::points_t
std::array< Wonton::Point< dim >, point_size > points_t
Definition: operator.h:209
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Interval >::apply
static result_t apply(const points_t p)
Definition: operator.h:164
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Hexahedron >::points_t
std::array< Wonton::Point< dim >, point_size > points_t
Definition: operator.h:374
Portage::Meshfree::oper::Sphere
Definition: operator.h:28
Portage::Meshfree::oper::Hexahedron
Definition: operator.h:28
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Triangle >::points_t
std::array< Wonton::Point< dim >, point_size > points_t
Definition: operator.h:230
Portage::Meshfree::oper::size_info
std::array< size_t, 3 > size_info(Type type, basis::Type basis, Domain domain)
Definition: operator.h:513
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Quadrilateral >::apply
static result_t apply(const points_t p)
Definition: operator.h:304
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Triangle >::apply
static result_t apply(const points_t p)
Definition: operator.h:211
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Linear, Wedge >::points_t
std::array< Wonton::Point< dim >, point_size > points_t
Definition: operator.h:421
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Hexahedron >::result_t
std::array< std::array< double, operator_size >, basis_size > result_t
Definition: operator.h:352
Portage::Meshfree::oper::Operator< VolumeIntegral, basis::Unitary, Tetrahedron >::result_t
std::array< std::array< double, operator_size >, basis_size > result_t
Definition: operator.h:446