Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2004 Sandia Corporation and Argonne National
5 : : Laboratory. Under the terms of Contract DE-AC04-94AL85000
6 : : with Sandia Corporation, the U.S. Government retains certain
7 : : rights in this software.
8 : :
9 : : This library is free software; you can redistribute it and/or
10 : : modify it under the terms of the GNU Lesser General Public
11 : : License as published by the Free Software Foundation; either
12 : : version 2.1 of the License, or (at your option) any later version.
13 : :
14 : : This library is distributed in the hope that it will be useful,
15 : : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 : : Lesser General Public License for more details.
18 : :
19 : : You should have received a copy of the GNU Lesser General Public License
20 : : (lgpl.txt) along with this library; if not, write to the Free Software
21 : : Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 : :
23 : : [email protected], [email protected], [email protected],
24 : : [email protected], [email protected], [email protected]
25 : :
26 : : ***************************************************************** */
27 : : // -*- Mode : c++; tab-width: 3; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 3
28 : : // -*-
29 : :
30 : : /*! \file MeanRatioFunctions.hpp
31 : :
32 : : Header that defines generalized Mean Ratio function, gradient, and hessian
33 : : evaluations.
34 : :
35 : : \author Todd Munson
36 : : \date 2003-06-11
37 : : */
38 : :
39 : : #ifndef MeanRatioFunctions_hpp
40 : : #define MeanRatioFunctions_hpp
41 : :
42 : : #include <math.h>
43 : : #include "Mesquite.hpp"
44 : : #include "Vector3D.hpp"
45 : : #include "Matrix3D.hpp"
46 : : #include "Exponent.hpp"
47 : :
48 : : namespace MBMesquite
49 : : {
50 : : /***************************************************************************/
51 : : /* Gradient calculation courtesy of Paul Hovland. The code was modified */
52 : : /* to reduce the number of flops and intermediate variables, and improve */
53 : : /* the locality of reference. */
54 : : /***************************************************************************/
55 : : /* The Hessian calculation is computes everyting in (x,y,z) blocks and */
56 : : /* stores only the upper triangular blocks. The results are put into */
57 : : /* (c1,c2,c3,c4) order prior to returning the Hessian matrix. */
58 : : /***************************************************************************/
59 : : /* The form of the function, gradient, and Hessian follows: */
60 : : /* o(x) = a * pow(f(A(x)), b) * pow(g(A(x)), c) */
61 : : /* where A(x) is the incidence matrix generated by: */
62 : : /* [x1-x0 x2-x0 x3-x0] */
63 : : /* A(x) = [y1-y0 y2-y0 y3-y0] * inv(W) */
64 : : /* [z1-z0 z2-z0 z3-z0] */
65 : : /* and f() is the squared Frobenius norm of A(x), and g() is the */
66 : : /* determinant of A(x). */
67 : : /* */
68 : : /* The gradient is calculated as follows: */
69 : : /* alpha := a*b*pow(f(A(x)),b-1)*pow(g(A(x)),c) */
70 : : /* beta := a*c*pow(f(A(x)),b)*pow(g(A(x)),c-1) */
71 : : /* */
72 : : /* do/dx = (alpha * (df/dA) + beta * (dg/dA)) (dA/dx) */
73 : : /* */
74 : : /* (df/dA)_i = 2*A_i */
75 : : /* (dg/dA)_i = A_j*A_k - A_l*A_m for some {j,k,l,m} */
76 : : /* */
77 : : /* d^2o/dx^2 = (dA/dx)' * ((d alpha/dA) * (df/dA) + */
78 : : /* (d beta/dA) * (dg/dA) */
79 : : /* alpha * (d^2f/dA^2) */
80 : : /* beta * (d^2g/dA^2)) * (dA/dx) */
81 : : /* */
82 : : /* Note: since A(x) is a linear function, there are no terms involving */
83 : : /* d^2A/dx^2 since this matrix is zero. */
84 : : /* */
85 : : /* gamma := a*b*c*pow(f(A(x)),b-1)*pow(g(A(x)),c-1) */
86 : : /* delta := a*c*(c-1)*pow(f(A(x)),b)*pow(g(A(x)),c-2) */
87 : : /* psi := a*b*(b-1)*pow(f(A(x)),b-2)*pow(g(A(x)),c) */
88 : : /* */
89 : : /* d^2o/dx^2 = (dA/dx)' * (gamma*((dg/dA)'*(df/dA) + (df/dA)'*(dg/dA)) + */
90 : : /* delta* (dg/dA)'*(dg/dA) + */
91 : : /* psi* (df/dA)'*(df/dA) + */
92 : : /* alpha*(d^2f/dA^2) + */
93 : : /* beta*(d^2g/dA^2)) * (dA/dx) */
94 : : /* */
95 : : /* Note: (df/dA) and (dg/dA) are row vectors and we only calculate the */
96 : : /* upper triangular part of the inner matrices. */
97 : : /* */
98 : : /* For regular tetrahedral elements, we have the following: */
99 : : /* */
100 : : /* [-1 1 0 0 ] */
101 : : /* M = [-sqrt(3) -sqrt(3) 2*sqrt(3) 0 ] */
102 : : /* [-sqrt(6) -sqrt(6) -sqrt(6) 3*sqrt(6) ] */
103 : : /* */
104 : : /* [M 0 0] */
105 : : /* dA/dx = [0 M 0] */
106 : : /* [0 0 M] */
107 : : /***************************************************************************/
108 : :
109 : : /*****************************************************************************/
110 : : /* Not all compilers substitute out constants (especially the square root). */
111 : : /* Therefore, they are substituted out manually. The values below were */
112 : : /* calculated on a solaris machine using long doubles. I believe they are */
113 : : /* accurate. */
114 : : /*****************************************************************************/
115 : :
116 : : #define isqrt3 5.77350269189625797959429519858e-01 /* 1.0/sqrt(3.0)*/
117 : : #define tisqrt3 1.15470053837925159591885903972e+00 /* 2.0/sqrt(3.0)*/
118 : : #define isqrt6 4.08248290463863052509822647505e-01 /* 1.0/sqrt(6.0)*/
119 : : #define tisqrt6 1.22474487139158915752946794252e+00 /* 3.0/sqrt(6.0)*/
120 : :
121 : : /*****************************************************************************/
122 : : /* The following set of functions reference triangular elements to an */
123 : : /* equilateral triangle in the plane defined by the normal. They are */
124 : : /* used when assessing the quality of a triangular element. A zero */
125 : : /* return value indicates success, while a nonzero value indicates failure. */
126 : : /*****************************************************************************/
127 : :
128 : : /*****************************************************************************/
129 : : /* Function evaluation requires 44 flops. */
130 : : /* Reductions possible when b == 1 or c == 1 */
131 : : /*****************************************************************************/
132 : 367 : inline bool m_fcn_2e( double& obj, const Vector3D x[3], const Vector3D& n, const double a, const Exponent& b,
133 : : const Exponent& c )
134 : : {
135 : : double matr[9], f;
136 : : double g;
137 : :
138 : : /* Calculate M = [A*inv(W) n] */
139 [ + - ][ + - ]: 367 : matr[0] = x[1][0] - x[0][0];
140 [ + - ][ + - ]: 367 : matr[1] = ( 2.0 * x[2][0] - x[1][0] - x[0][0] ) * isqrt3;
[ + - ]
141 [ + - ]: 367 : matr[2] = n[0];
142 : :
143 [ + - ][ + - ]: 367 : matr[3] = x[1][1] - x[0][1];
144 [ + - ][ + - ]: 367 : matr[4] = ( 2.0 * x[2][1] - x[1][1] - x[0][1] ) * isqrt3;
[ + - ]
145 [ + - ]: 367 : matr[5] = n[1];
146 : :
147 [ + - ][ + - ]: 367 : matr[6] = x[1][2] - x[0][2];
148 [ + - ][ + - ]: 367 : matr[7] = ( 2.0 * x[2][2] - x[1][2] - x[0][2] ) * isqrt3;
[ + - ]
149 [ + - ]: 367 : matr[8] = n[2];
150 : :
151 : : /* Calculate det(M). */
152 : 367 : g = matr[0] * ( matr[4] * matr[8] - matr[5] * matr[7] ) + matr[3] * ( matr[2] * matr[7] - matr[1] * matr[8] ) +
153 : 367 : matr[6] * ( matr[1] * matr[5] - matr[2] * matr[4] );
154 [ + + ]: 367 : if( g < MSQ_MIN )
155 : : {
156 : 3 : obj = g;
157 : 3 : return false;
158 : : }
159 : :
160 : : /* Calculate norm(M). */
161 : 364 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[3] * matr[3] + matr[4] * matr[4] + matr[6] * matr[6] +
162 : 364 : matr[7] * matr[7];
163 : :
164 : : /* Calculate objective function. */
165 [ + - ][ + - ]: 364 : obj = a * b.raise( f ) * c.raise( g );
166 : 367 : return true;
167 : : }
168 : :
169 : : /*****************************************************************************/
170 : : /* Gradient evaluation requires 88 flops. */
171 : : /* Reductions possible when b == 1 or c == 1 */
172 : : /*****************************************************************************/
173 : 101 : inline bool g_fcn_2e( double& obj, Vector3D g_obj[3], const Vector3D x[3], const Vector3D& n, const double a,
174 : : const Exponent& b, const Exponent& c )
175 : : {
176 : : double matr[9], f;
177 : : double adj_m[9], g; // adj_m[2,5,8] not used
178 : : double loc1, loc2, loc3;
179 : :
180 : : /* Calculate M = [A*inv(W) n] */
181 [ + - ][ + - ]: 101 : matr[0] = x[1][0] - x[0][0];
182 [ + - ][ + - ]: 101 : matr[1] = ( 2.0 * x[2][0] - x[1][0] - x[0][0] ) * isqrt3;
[ + - ]
183 [ + - ]: 101 : matr[2] = n[0];
184 : :
185 [ + - ][ + - ]: 101 : matr[3] = x[1][1] - x[0][1];
186 [ + - ][ + - ]: 101 : matr[4] = ( 2.0 * x[2][1] - x[1][1] - x[0][1] ) * isqrt3;
[ + - ]
187 [ + - ]: 101 : matr[5] = n[1];
188 : :
189 [ + - ][ + - ]: 101 : matr[6] = x[1][2] - x[0][2];
190 [ + - ][ + - ]: 101 : matr[7] = ( 2.0 * x[2][2] - x[1][2] - x[0][2] ) * isqrt3;
[ + - ]
191 [ + - ]: 101 : matr[8] = n[2];
192 : :
193 : : /* Calculate det([n M]). */
194 : 101 : loc1 = matr[4] * matr[8] - matr[5] * matr[7];
195 : 101 : loc2 = matr[2] * matr[7] - matr[1] * matr[8];
196 : 101 : loc3 = matr[1] * matr[5] - matr[2] * matr[4];
197 : 101 : g = matr[0] * loc1 + matr[3] * loc2 + matr[6] * loc3;
198 [ - + ]: 101 : if( g < MSQ_MIN )
199 : : {
200 : 0 : obj = g;
201 : 0 : return false;
202 : : }
203 : :
204 : : /* Calculate norm(M). */
205 : 101 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[3] * matr[3] + matr[4] * matr[4] + matr[6] * matr[6] +
206 : 101 : matr[7] * matr[7];
207 : :
208 : : /* Calculate objective function. */
209 [ + - ][ + - ]: 101 : obj = a * b.raise( f ) * c.raise( g );
210 : :
211 : : /* Calculate the derivative of the objective function. */
212 [ + - ]: 101 : f = b.value() * obj / f; /* Constant on nabla f */
213 [ + - ]: 101 : g = c.value() * obj / g; /* Constant on nable g */
214 : 101 : f *= 2.0; /* Modification for nabla f */
215 : :
216 : 101 : adj_m[0] = matr[0] * f + loc1 * g;
217 : 101 : adj_m[3] = matr[3] * f + loc2 * g;
218 : 101 : adj_m[6] = matr[6] * f + loc3 * g;
219 : :
220 : 101 : loc1 = matr[0] * g;
221 : 101 : loc2 = matr[3] * g;
222 : 101 : loc3 = matr[6] * g;
223 : :
224 : 101 : adj_m[1] = matr[1] * f + loc3 * matr[5] - loc2 * matr[8];
225 : 101 : adj_m[4] = matr[4] * f + loc1 * matr[8] - loc3 * matr[2];
226 : 101 : adj_m[7] = matr[7] * f + loc2 * matr[2] - loc1 * matr[5];
227 : :
228 : 101 : loc1 = isqrt3 * adj_m[1];
229 [ + - ]: 101 : g_obj[0][0] = -adj_m[0] - loc1;
230 [ + - ]: 101 : g_obj[1][0] = adj_m[0] - loc1;
231 [ + - ]: 101 : g_obj[2][0] = 2.0 * loc1;
232 : :
233 : 101 : loc1 = isqrt3 * adj_m[4];
234 [ + - ]: 101 : g_obj[0][1] = -adj_m[3] - loc1;
235 [ + - ]: 101 : g_obj[1][1] = adj_m[3] - loc1;
236 [ + - ]: 101 : g_obj[2][1] = 2.0 * loc1;
237 : :
238 : 101 : loc1 = isqrt3 * adj_m[7];
239 [ + - ]: 101 : g_obj[0][2] = -adj_m[6] - loc1;
240 [ + - ]: 101 : g_obj[1][2] = adj_m[6] - loc1;
241 [ + - ]: 101 : g_obj[2][2] = 2.0 * loc1;
242 : 101 : return true;
243 : : }
244 : :
245 : : /*****************************************************************************/
246 : : /* Hessian evaluation requires 316 flops. */
247 : : /* Reductions possible when b == 1 or c == 1 */
248 : : /*****************************************************************************/
249 : 0 : inline bool h_fcn_2e( double& obj, Vector3D g_obj[3], Matrix3D h_obj[6], const Vector3D x[3], const Vector3D& n,
250 : : const double a, const Exponent& b, const Exponent& c )
251 : : {
252 : : double matr[9], f;
253 : : double adj_m[9], g; // adj_m[2,5,8] not used
254 : : double dg[9]; // dg[2,5,8] not used
255 : : double loc0, loc1, loc2, loc3, loc4;
256 : : double A[12], J_A[6], J_B[9], J_C[9], cross; // only 2x2 corners used
257 : :
258 : : /* Calculate M = [A*inv(W) n] */
259 [ # # ][ # # ]: 0 : matr[0] = x[1][0] - x[0][0];
260 [ # # ][ # # ]: 0 : matr[1] = ( 2.0 * x[2][0] - x[1][0] - x[0][0] ) * isqrt3;
[ # # ]
261 [ # # ]: 0 : matr[2] = n[0];
262 : :
263 [ # # ][ # # ]: 0 : matr[3] = x[1][1] - x[0][1];
264 [ # # ][ # # ]: 0 : matr[4] = ( 2.0 * x[2][1] - x[1][1] - x[0][1] ) * isqrt3;
[ # # ]
265 [ # # ]: 0 : matr[5] = n[1];
266 : :
267 [ # # ][ # # ]: 0 : matr[6] = x[1][2] - x[0][2];
268 [ # # ][ # # ]: 0 : matr[7] = ( 2.0 * x[2][2] - x[1][2] - x[0][2] ) * isqrt3;
[ # # ]
269 [ # # ]: 0 : matr[8] = n[2];
270 : :
271 : : /* Calculate det([n M]). */
272 : 0 : dg[0] = matr[4] * matr[8] - matr[5] * matr[7];
273 : 0 : dg[3] = matr[2] * matr[7] - matr[1] * matr[8];
274 : 0 : dg[6] = matr[1] * matr[5] - matr[2] * matr[4];
275 : 0 : g = matr[0] * dg[0] + matr[3] * dg[3] + matr[6] * dg[6];
276 [ # # ]: 0 : if( g < MSQ_MIN )
277 : : {
278 : 0 : obj = g;
279 : 0 : return false;
280 : : }
281 : :
282 : : /* Calculate norm(M). */
283 : 0 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[3] * matr[3] + matr[4] * matr[4] + matr[6] * matr[6] +
284 : 0 : matr[7] * matr[7];
285 : :
286 : 0 : loc3 = f;
287 : 0 : loc4 = g;
288 : :
289 : : /* Calculate objective function. */
290 [ # # ][ # # ]: 0 : obj = a * b.raise( f ) * c.raise( g );
291 : :
292 : : /* Calculate the derivative of the objective function. */
293 [ # # ]: 0 : f = b.value() * obj / f; /* Constant on nabla f */
294 [ # # ]: 0 : g = c.value() * obj / g; /* Constant on nable g */
295 : 0 : f *= 2.0; /* Modification for nabla f */
296 : :
297 : 0 : dg[1] = matr[5] * matr[6] - matr[3] * matr[8];
298 : 0 : dg[4] = matr[0] * matr[8] - matr[2] * matr[6];
299 : 0 : dg[7] = matr[2] * matr[3] - matr[0] * matr[5];
300 : :
301 : 0 : adj_m[0] = matr[0] * f + dg[0] * g;
302 : 0 : adj_m[1] = matr[1] * f + dg[1] * g;
303 : 0 : adj_m[3] = matr[3] * f + dg[3] * g;
304 : 0 : adj_m[4] = matr[4] * f + dg[4] * g;
305 : 0 : adj_m[6] = matr[6] * f + dg[6] * g;
306 : 0 : adj_m[7] = matr[7] * f + dg[7] * g;
307 : :
308 : 0 : loc1 = isqrt3 * adj_m[1];
309 [ # # ]: 0 : g_obj[0][0] = -adj_m[0] - loc1;
310 [ # # ]: 0 : g_obj[1][0] = adj_m[0] - loc1;
311 [ # # ]: 0 : g_obj[2][0] = 2.0 * loc1;
312 : :
313 : 0 : loc1 = isqrt3 * adj_m[4];
314 [ # # ]: 0 : g_obj[0][1] = -adj_m[3] - loc1;
315 [ # # ]: 0 : g_obj[1][1] = adj_m[3] - loc1;
316 [ # # ]: 0 : g_obj[2][1] = 2.0 * loc1;
317 : :
318 : 0 : loc1 = isqrt3 * adj_m[7];
319 [ # # ]: 0 : g_obj[0][2] = -adj_m[6] - loc1;
320 [ # # ]: 0 : g_obj[1][2] = adj_m[6] - loc1;
321 [ # # ]: 0 : g_obj[2][2] = 2.0 * loc1;
322 : :
323 : : /* Calculate the hessian of the objective. */
324 : 0 : loc0 = f; /* Constant on nabla^2 f */
325 : 0 : loc1 = g; /* Constant on nabla^2 g */
326 [ # # ]: 0 : cross = f * c.value() / loc4; /* Constant on nabla g nabla f */
327 [ # # ]: 0 : f = f * ( b.value() - 1 ) / loc3; /* Constant on nabla f nabla f */
328 [ # # ]: 0 : g = g * ( c.value() - 1 ) / loc4; /* Constant on nabla g nabla g */
329 : 0 : f *= 2.0; /* Modification for nabla f */
330 : :
331 : : /* First block of rows */
332 : 0 : loc3 = matr[0] * f + dg[0] * cross;
333 : 0 : loc4 = dg[0] * g + matr[0] * cross;
334 : :
335 : 0 : J_A[0] = loc0 + loc3 * matr[0] + loc4 * dg[0];
336 : 0 : J_A[1] = loc3 * matr[1] + loc4 * dg[1];
337 : 0 : J_B[0] = loc3 * matr[3] + loc4 * dg[3];
338 : 0 : J_B[1] = loc3 * matr[4] + loc4 * dg[4];
339 : 0 : J_C[0] = loc3 * matr[6] + loc4 * dg[6];
340 : 0 : J_C[1] = loc3 * matr[7] + loc4 * dg[7];
341 : :
342 : 0 : loc3 = matr[1] * f + dg[1] * cross;
343 : 0 : loc4 = dg[1] * g + matr[1] * cross;
344 : :
345 : 0 : J_A[3] = loc0 + loc3 * matr[1] + loc4 * dg[1];
346 : 0 : J_B[3] = loc3 * matr[3] + loc4 * dg[3];
347 : 0 : J_B[4] = loc3 * matr[4] + loc4 * dg[4];
348 : 0 : J_C[3] = loc3 * matr[6] + loc4 * dg[6];
349 : 0 : J_C[4] = loc3 * matr[7] + loc4 * dg[7];
350 : :
351 : : /* First diagonal block */
352 : 0 : loc2 = isqrt3 * J_A[1];
353 : 0 : A[0] = -J_A[0] - loc2;
354 : 0 : A[1] = J_A[0] - loc2;
355 : :
356 : 0 : loc2 = isqrt3 * J_A[3];
357 : 0 : A[4] = -J_A[1] - loc2;
358 : 0 : A[5] = J_A[1] - loc2;
359 : 0 : A[6] = 2.0 * loc2;
360 : :
361 : 0 : loc2 = isqrt3 * A[4];
362 [ # # ]: 0 : h_obj[0][0][0] = -A[0] - loc2;
363 [ # # ]: 0 : h_obj[1][0][0] = A[0] - loc2;
364 [ # # ]: 0 : h_obj[2][0][0] = 2.0 * loc2;
365 : :
366 : 0 : loc2 = isqrt3 * A[5];
367 [ # # ]: 0 : h_obj[3][0][0] = A[1] - loc2;
368 [ # # ]: 0 : h_obj[4][0][0] = 2.0 * loc2;
369 : :
370 [ # # ]: 0 : h_obj[5][0][0] = tisqrt3 * A[6];
371 : :
372 : : /* First off-diagonal block */
373 : 0 : loc2 = matr[8] * loc1;
374 : 0 : J_B[1] += loc2;
375 : 0 : J_B[3] -= loc2;
376 : :
377 : 0 : loc2 = isqrt3 * J_B[3];
378 : 0 : A[0] = -J_B[0] - loc2;
379 : 0 : A[1] = J_B[0] - loc2;
380 : 0 : A[2] = 2.0 * loc2;
381 : :
382 : 0 : loc2 = isqrt3 * J_B[4];
383 : 0 : A[4] = -J_B[1] - loc2;
384 : 0 : A[5] = J_B[1] - loc2;
385 : 0 : A[6] = 2.0 * loc2;
386 : :
387 : 0 : loc2 = isqrt3 * A[4];
388 [ # # ]: 0 : h_obj[0][0][1] = -A[0] - loc2;
389 [ # # ]: 0 : h_obj[1][0][1] = A[0] - loc2;
390 [ # # ]: 0 : h_obj[2][0][1] = 2.0 * loc2;
391 : :
392 : 0 : loc2 = isqrt3 * A[5];
393 [ # # ]: 0 : h_obj[1][1][0] = -A[1] - loc2;
394 [ # # ]: 0 : h_obj[3][0][1] = A[1] - loc2;
395 [ # # ]: 0 : h_obj[4][0][1] = 2.0 * loc2;
396 : :
397 : 0 : loc2 = isqrt3 * A[6];
398 [ # # ]: 0 : h_obj[2][1][0] = -A[2] - loc2;
399 [ # # ]: 0 : h_obj[4][1][0] = A[2] - loc2;
400 [ # # ]: 0 : h_obj[5][0][1] = 2.0 * loc2;
401 : :
402 : : /* Second off-diagonal block */
403 : 0 : loc2 = matr[5] * loc1;
404 : 0 : J_C[1] -= loc2;
405 : 0 : J_C[3] += loc2;
406 : :
407 : 0 : loc2 = isqrt3 * J_C[3];
408 : 0 : A[0] = -J_C[0] - loc2;
409 : 0 : A[1] = J_C[0] - loc2;
410 : 0 : A[2] = 2.0 * loc2;
411 : :
412 : 0 : loc2 = isqrt3 * J_C[4];
413 : 0 : A[4] = -J_C[1] - loc2;
414 : 0 : A[5] = J_C[1] - loc2;
415 : 0 : A[6] = 2.0 * loc2;
416 : :
417 : 0 : loc2 = isqrt3 * A[4];
418 [ # # ]: 0 : h_obj[0][0][2] = -A[0] - loc2;
419 [ # # ]: 0 : h_obj[1][0][2] = A[0] - loc2;
420 [ # # ]: 0 : h_obj[2][0][2] = 2.0 * loc2;
421 : :
422 : 0 : loc2 = isqrt3 * A[5];
423 [ # # ]: 0 : h_obj[1][2][0] = -A[1] - loc2;
424 [ # # ]: 0 : h_obj[3][0][2] = A[1] - loc2;
425 [ # # ]: 0 : h_obj[4][0][2] = 2.0 * loc2;
426 : :
427 : 0 : loc2 = isqrt3 * A[6];
428 [ # # ]: 0 : h_obj[2][2][0] = -A[2] - loc2;
429 [ # # ]: 0 : h_obj[4][2][0] = A[2] - loc2;
430 [ # # ]: 0 : h_obj[5][0][2] = 2.0 * loc2;
431 : :
432 : : /* Second block of rows */
433 : 0 : loc3 = matr[3] * f + dg[3] * cross;
434 : 0 : loc4 = dg[3] * g + matr[3] * cross;
435 : :
436 : 0 : J_A[0] = loc0 + loc3 * matr[3] + loc4 * dg[3];
437 : 0 : J_A[1] = loc3 * matr[4] + loc4 * dg[4];
438 : 0 : J_B[0] = loc3 * matr[6] + loc4 * dg[6];
439 : 0 : J_B[1] = loc3 * matr[7] + loc4 * dg[7];
440 : :
441 : 0 : loc3 = matr[4] * f + dg[4] * cross;
442 : 0 : loc4 = dg[4] * g + matr[4] * cross;
443 : :
444 : 0 : J_A[3] = loc0 + loc3 * matr[4] + loc4 * dg[4];
445 : 0 : J_B[3] = loc3 * matr[6] + loc4 * dg[6];
446 : 0 : J_B[4] = loc3 * matr[7] + loc4 * dg[7];
447 : :
448 : : /* Second diagonal block */
449 : 0 : loc2 = isqrt3 * J_A[1];
450 : 0 : A[0] = -J_A[0] - loc2;
451 : 0 : A[1] = J_A[0] - loc2;
452 : :
453 : 0 : loc2 = isqrt3 * J_A[3];
454 : 0 : A[4] = -J_A[1] - loc2;
455 : 0 : A[5] = J_A[1] - loc2;
456 : 0 : A[6] = 2.0 * loc2;
457 : :
458 : 0 : loc2 = isqrt3 * A[4];
459 [ # # ]: 0 : h_obj[0][1][1] = -A[0] - loc2;
460 [ # # ]: 0 : h_obj[1][1][1] = A[0] - loc2;
461 [ # # ]: 0 : h_obj[2][1][1] = 2.0 * loc2;
462 : :
463 : 0 : loc2 = isqrt3 * A[5];
464 [ # # ]: 0 : h_obj[3][1][1] = A[1] - loc2;
465 [ # # ]: 0 : h_obj[4][1][1] = 2.0 * loc2;
466 : :
467 [ # # ]: 0 : h_obj[5][1][1] = tisqrt3 * A[6];
468 : :
469 : : /* Third off-diagonal block */
470 : 0 : loc2 = matr[2] * loc1;
471 : 0 : J_B[1] += loc2;
472 : 0 : J_B[3] -= loc2;
473 : :
474 : 0 : loc2 = isqrt3 * J_B[3];
475 : 0 : A[0] = -J_B[0] - loc2;
476 : 0 : A[1] = J_B[0] - loc2;
477 : 0 : A[2] = 2.0 * loc2;
478 : :
479 : 0 : loc2 = isqrt3 * J_B[4];
480 : 0 : A[4] = -J_B[1] - loc2;
481 : 0 : A[5] = J_B[1] - loc2;
482 : 0 : A[6] = 2.0 * loc2;
483 : :
484 : 0 : loc2 = isqrt3 * A[4];
485 [ # # ]: 0 : h_obj[0][1][2] = -A[0] - loc2;
486 [ # # ]: 0 : h_obj[1][1][2] = A[0] - loc2;
487 [ # # ]: 0 : h_obj[2][1][2] = 2.0 * loc2;
488 : :
489 : 0 : loc2 = isqrt3 * A[5];
490 [ # # ]: 0 : h_obj[1][2][1] = -A[1] - loc2;
491 [ # # ]: 0 : h_obj[3][1][2] = A[1] - loc2;
492 [ # # ]: 0 : h_obj[4][1][2] = 2.0 * loc2;
493 : :
494 : 0 : loc2 = isqrt3 * A[6];
495 [ # # ]: 0 : h_obj[2][2][1] = -A[2] - loc2;
496 [ # # ]: 0 : h_obj[4][2][1] = A[2] - loc2;
497 [ # # ]: 0 : h_obj[5][1][2] = 2.0 * loc2;
498 : :
499 : : /* Third block of rows */
500 : 0 : loc3 = matr[6] * f + dg[6] * cross;
501 : 0 : loc4 = dg[6] * g + matr[6] * cross;
502 : :
503 : 0 : J_A[0] = loc0 + loc3 * matr[6] + loc4 * dg[6];
504 : 0 : J_A[1] = loc3 * matr[7] + loc4 * dg[7];
505 : :
506 : 0 : loc3 = matr[7] * f + dg[7] * cross;
507 : 0 : loc4 = dg[7] * g + matr[7] * cross;
508 : :
509 : 0 : J_A[3] = loc0 + loc3 * matr[7] + loc4 * dg[7];
510 : :
511 : : /* Third diagonal block */
512 : 0 : loc2 = isqrt3 * J_A[1];
513 : 0 : A[0] = -J_A[0] - loc2;
514 : 0 : A[1] = J_A[0] - loc2;
515 : :
516 : 0 : loc2 = isqrt3 * J_A[3];
517 : 0 : A[4] = -J_A[1] - loc2;
518 : 0 : A[5] = J_A[1] - loc2;
519 : 0 : A[6] = 2.0 * loc2;
520 : :
521 : 0 : loc2 = isqrt3 * A[4];
522 [ # # ]: 0 : h_obj[0][2][2] = -A[0] - loc2;
523 [ # # ]: 0 : h_obj[1][2][2] = A[0] - loc2;
524 [ # # ]: 0 : h_obj[2][2][2] = 2.0 * loc2;
525 : :
526 : 0 : loc2 = isqrt3 * A[5];
527 [ # # ]: 0 : h_obj[3][2][2] = A[1] - loc2;
528 [ # # ]: 0 : h_obj[4][2][2] = 2.0 * loc2;
529 : :
530 [ # # ]: 0 : h_obj[5][2][2] = tisqrt3 * A[6];
531 : :
532 : : // completes diagonal blocks.
533 [ # # ]: 0 : h_obj[0].fill_lower_triangle();
534 [ # # ]: 0 : h_obj[3].fill_lower_triangle();
535 [ # # ]: 0 : h_obj[5].fill_lower_triangle();
536 : 0 : return true;
537 : : }
538 : :
539 : : /*****************************************************************************/
540 : : /* The following set of functions reference triangular elements to an */
541 : : /* right triangle in the plane defined by the normal. They are used when */
542 : : /* assessing the quality of a quadrilateral elements. A zero return value */
543 : : /* indicates success, while a nonzero value indicates failure. */
544 : : /*****************************************************************************/
545 : :
546 : : /*****************************************************************************/
547 : : /* Function evaluation -- requires 41 flops. */
548 : : /* Reductions possible when b == 1, c == 1, or d == 1 */
549 : : /*****************************************************************************/
550 : 71052 : inline bool m_fcn_2i( double& obj, const Vector3D x[3], const Vector3D& n, const double a, const Exponent& b,
551 : : const Exponent& c, const Vector3D& d )
552 : : {
553 : : double matr[9];
554 : : double f;
555 : : double g;
556 : :
557 : : /* Calculate M = A*inv(W). */
558 [ + - ][ + - ]: 71052 : matr[0] = d[0] * ( x[1][0] - x[0][0] );
[ + - ]
559 [ + - ][ + - ]: 71052 : matr[1] = d[1] * ( x[2][0] - x[0][0] );
[ + - ]
560 [ + - ]: 71052 : matr[2] = n[0];
561 : :
562 [ + - ][ + - ]: 71052 : matr[3] = d[0] * ( x[1][1] - x[0][1] );
[ + - ]
563 [ + - ][ + - ]: 71052 : matr[4] = d[1] * ( x[2][1] - x[0][1] );
[ + - ]
564 [ + - ]: 71052 : matr[5] = n[1];
565 : :
566 [ + - ][ + - ]: 71052 : matr[6] = d[0] * ( x[1][2] - x[0][2] );
[ + - ]
567 [ + - ][ + - ]: 71052 : matr[7] = d[1] * ( x[2][2] - x[0][2] );
[ + - ]
568 [ + - ]: 71052 : matr[8] = n[2];
569 : :
570 : : /* Calculate det(M). */
571 : 71052 : g = matr[0] * ( matr[4] * matr[8] - matr[5] * matr[7] ) + matr[3] * ( matr[2] * matr[7] - matr[1] * matr[8] ) +
572 : 71052 : matr[6] * ( matr[1] * matr[5] - matr[2] * matr[4] );
573 [ + + ]: 71052 : if( g < MSQ_MIN )
574 : : {
575 : 28 : obj = g;
576 : 28 : return false;
577 : : }
578 : :
579 : : /* Calculate norm(M). */
580 : 71024 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[3] * matr[3] + matr[4] * matr[4] + matr[6] * matr[6] +
581 : 71024 : matr[7] * matr[7];
582 : :
583 : : /* Calculate objective function. */
584 [ + - ][ + - ]: 71024 : obj = a * b.raise( f ) * c.raise( g );
585 : 71052 : return true;
586 : : }
587 : :
588 : : /*****************************************************************************/
589 : : /* Gradient requires 82 flops. */
590 : : /* Reductions possible when b == 1, c == 1, or d == 1 */
591 : : /*****************************************************************************/
592 : 2128 : inline bool g_fcn_2i( double& obj, Vector3D g_obj[3], const Vector3D x[3], const Vector3D& n, const double a,
593 : : const Exponent& b, const Exponent& c, const Vector3D& d )
594 : : {
595 : : double matr[9], f;
596 : : double adj_m[9], g; // adj_m[2,5,8] not used
597 : : double loc1, loc2, loc3;
598 : :
599 : : /* Calculate M = [A*inv(W) n] */
600 [ + - ][ + - ]: 2128 : matr[0] = d[0] * ( x[1][0] - x[0][0] );
[ + - ]
601 [ + - ][ + - ]: 2128 : matr[1] = d[1] * ( x[2][0] - x[0][0] );
[ + - ]
602 [ + - ]: 2128 : matr[2] = n[0];
603 : :
604 [ + - ][ + - ]: 2128 : matr[3] = d[0] * ( x[1][1] - x[0][1] );
[ + - ]
605 [ + - ][ + - ]: 2128 : matr[4] = d[1] * ( x[2][1] - x[0][1] );
[ + - ]
606 [ + - ]: 2128 : matr[5] = n[1];
607 : :
608 [ + - ][ + - ]: 2128 : matr[6] = d[0] * ( x[1][2] - x[0][2] );
[ + - ]
609 [ + - ][ + - ]: 2128 : matr[7] = d[1] * ( x[2][2] - x[0][2] );
[ + - ]
610 [ + - ]: 2128 : matr[8] = n[2];
611 : :
612 : : /* Calculate det([n M]). */
613 : 2128 : loc1 = matr[4] * matr[8] - matr[5] * matr[7];
614 : 2128 : loc2 = matr[2] * matr[7] - matr[1] * matr[8];
615 : 2128 : loc3 = matr[1] * matr[5] - matr[2] * matr[4];
616 : 2128 : g = matr[0] * loc1 + matr[3] * loc2 + matr[6] * loc3;
617 [ - + ]: 2128 : if( g < MSQ_MIN )
618 : : {
619 : 0 : obj = g;
620 : 0 : return false;
621 : : }
622 : :
623 : : /* Calculate norm(M). */
624 : 2128 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[3] * matr[3] + matr[4] * matr[4] + matr[6] * matr[6] +
625 : 2128 : matr[7] * matr[7];
626 : :
627 : : /* Calculate objective function. */
628 [ + - ][ + - ]: 2128 : obj = a * b.raise( f ) * c.raise( g );
629 : :
630 : : /* Calculate the derivative of the objective function. */
631 [ + - ]: 2128 : f = b.value() * obj / f; /* Constant on nabla f */
632 [ + - ]: 2128 : g = c.value() * obj / g; /* Constant on nable g */
633 : 2128 : f *= 2.0; /* Modification for nabla f */
634 : :
635 [ + - ]: 2128 : adj_m[0] = d[0] * ( matr[0] * f + loc1 * g );
636 [ + - ]: 2128 : adj_m[3] = d[0] * ( matr[3] * f + loc2 * g );
637 [ + - ]: 2128 : adj_m[6] = d[0] * ( matr[6] * f + loc3 * g );
638 : :
639 : 2128 : loc1 = matr[0] * g;
640 : 2128 : loc2 = matr[3] * g;
641 : 2128 : loc3 = matr[6] * g;
642 : :
643 [ + - ]: 2128 : adj_m[1] = d[1] * ( matr[1] * f + loc3 * matr[5] - loc2 * matr[8] );
644 [ + - ]: 2128 : adj_m[4] = d[1] * ( matr[4] * f + loc1 * matr[8] - loc3 * matr[2] );
645 [ + - ]: 2128 : adj_m[7] = d[1] * ( matr[7] * f + loc2 * matr[2] - loc1 * matr[5] );
646 : :
647 [ + - ]: 2128 : g_obj[0][0] = -adj_m[0] - adj_m[1];
648 [ + - ]: 2128 : g_obj[1][0] = adj_m[0];
649 [ + - ]: 2128 : g_obj[2][0] = adj_m[1];
650 : :
651 [ + - ]: 2128 : g_obj[0][1] = -adj_m[3] - adj_m[4];
652 [ + - ]: 2128 : g_obj[1][1] = adj_m[3];
653 [ + - ]: 2128 : g_obj[2][1] = adj_m[4];
654 : :
655 [ + - ]: 2128 : g_obj[0][2] = -adj_m[6] - adj_m[7];
656 [ + - ]: 2128 : g_obj[1][2] = adj_m[6];
657 [ + - ]: 2128 : g_obj[2][2] = adj_m[7];
658 : 2128 : return true;
659 : : }
660 : :
661 : : /*****************************************************************************/
662 : : /* Hessian requires 253 flops. */
663 : : /* Reductions possible when b == 1, c == 1, or d == 1 */
664 : : /*****************************************************************************/
665 : 0 : inline bool h_fcn_2i( double& obj, Vector3D g_obj[3], Matrix3D h_obj[6], const Vector3D x[3], const Vector3D& n,
666 : : const double a, const Exponent& b, const Exponent& c, const Vector3D& d )
667 : : {
668 : : double matr[9], f;
669 : : double adj_m[9], g; // adj_m[2,5,8] not used
670 : : double dg[9]; // dg[2,5,8] not used
671 : : double loc0, loc1, loc2, loc3, loc4;
672 : : double A[12], J_A[6], J_B[9], J_C[9], cross; // only 2x2 corners used
673 : :
674 [ # # ][ # # ]: 0 : const double scale[3] = { d[0] * d[0], d[0] * d[1], d[1] * d[1] };
[ # # ][ # # ]
[ # # ][ # # ]
675 : :
676 : : /* Calculate M = [A*inv(W) n] */
677 [ # # ][ # # ]: 0 : matr[0] = d[0] * ( x[1][0] - x[0][0] );
[ # # ]
678 [ # # ][ # # ]: 0 : matr[1] = d[1] * ( x[2][0] - x[0][0] );
[ # # ]
679 [ # # ]: 0 : matr[2] = n[0];
680 : :
681 [ # # ][ # # ]: 0 : matr[3] = d[0] * ( x[1][1] - x[0][1] );
[ # # ]
682 [ # # ][ # # ]: 0 : matr[4] = d[1] * ( x[2][1] - x[0][1] );
[ # # ]
683 [ # # ]: 0 : matr[5] = n[1];
684 : :
685 [ # # ][ # # ]: 0 : matr[6] = d[0] * ( x[1][2] - x[0][2] );
[ # # ]
686 [ # # ][ # # ]: 0 : matr[7] = d[1] * ( x[2][2] - x[0][2] );
[ # # ]
687 [ # # ]: 0 : matr[8] = n[2];
688 : :
689 : : /* Calculate det([n M]). */
690 : 0 : dg[0] = matr[4] * matr[8] - matr[5] * matr[7];
691 : 0 : dg[3] = matr[2] * matr[7] - matr[1] * matr[8];
692 : 0 : dg[6] = matr[1] * matr[5] - matr[2] * matr[4];
693 : 0 : g = matr[0] * dg[0] + matr[3] * dg[3] + matr[6] * dg[6];
694 [ # # ]: 0 : if( g < MSQ_MIN )
695 : : {
696 : 0 : obj = g;
697 : 0 : return false;
698 : : }
699 : :
700 : : /* Calculate norm(M). */
701 : 0 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[3] * matr[3] + matr[4] * matr[4] + matr[6] * matr[6] +
702 : 0 : matr[7] * matr[7];
703 : :
704 : 0 : loc3 = f;
705 : 0 : loc4 = g;
706 : :
707 : : /* Calculate objective function. */
708 [ # # ][ # # ]: 0 : obj = a * b.raise( f ) * c.raise( g );
709 : :
710 : : /* Calculate the derivative of the objective function. */
711 [ # # ]: 0 : f = b.value() * obj / f; /* Constant on nabla f */
712 [ # # ]: 0 : g = c.value() * obj / g; /* Constant on nable g */
713 : 0 : f *= 2.0; /* Modification for nabla f */
714 : :
715 : 0 : dg[1] = matr[5] * matr[6] - matr[3] * matr[8];
716 : 0 : dg[4] = matr[0] * matr[8] - matr[2] * matr[6];
717 : 0 : dg[7] = matr[2] * matr[3] - matr[0] * matr[5];
718 : :
719 [ # # ]: 0 : adj_m[0] = d[0] * ( matr[0] * f + dg[0] * g );
720 [ # # ]: 0 : adj_m[1] = d[1] * ( matr[1] * f + dg[1] * g );
721 [ # # ]: 0 : adj_m[3] = d[0] * ( matr[3] * f + dg[3] * g );
722 [ # # ]: 0 : adj_m[4] = d[1] * ( matr[4] * f + dg[4] * g );
723 [ # # ]: 0 : adj_m[6] = d[0] * ( matr[6] * f + dg[6] * g );
724 [ # # ]: 0 : adj_m[7] = d[1] * ( matr[7] * f + dg[7] * g );
725 : :
726 [ # # ]: 0 : g_obj[0][0] = -adj_m[0] - adj_m[1];
727 [ # # ]: 0 : g_obj[1][0] = adj_m[0];
728 [ # # ]: 0 : g_obj[2][0] = adj_m[1];
729 : :
730 [ # # ]: 0 : g_obj[0][1] = -adj_m[3] - adj_m[4];
731 [ # # ]: 0 : g_obj[1][1] = adj_m[3];
732 [ # # ]: 0 : g_obj[2][1] = adj_m[4];
733 : :
734 [ # # ]: 0 : g_obj[0][2] = -adj_m[6] - adj_m[7];
735 [ # # ]: 0 : g_obj[1][2] = adj_m[6];
736 [ # # ]: 0 : g_obj[2][2] = adj_m[7];
737 : :
738 : : /* Calculate the hessian of the objective. */
739 : 0 : loc0 = f; /* Constant on nabla^2 f */
740 : 0 : loc1 = g; /* Constant on nabla^2 g */
741 [ # # ]: 0 : cross = f * c.value() / loc4; /* Constant on nabla g nabla f */
742 [ # # ]: 0 : f = f * ( b.value() - 1 ) / loc3; /* Constant on nabla f nabla f */
743 [ # # ]: 0 : g = g * ( c.value() - 1 ) / loc4; /* Constant on nabla g nabla g */
744 : 0 : f *= 2.0; /* Modification for nabla f */
745 : :
746 : : /* First block of rows */
747 : 0 : loc3 = matr[0] * f + dg[0] * cross;
748 : 0 : loc4 = dg[0] * g + matr[0] * cross;
749 : :
750 : 0 : J_A[0] = loc0 + loc3 * matr[0] + loc4 * dg[0];
751 : 0 : J_A[1] = loc3 * matr[1] + loc4 * dg[1];
752 : 0 : J_B[0] = loc3 * matr[3] + loc4 * dg[3];
753 : 0 : J_B[1] = loc3 * matr[4] + loc4 * dg[4];
754 : 0 : J_C[0] = loc3 * matr[6] + loc4 * dg[6];
755 : 0 : J_C[1] = loc3 * matr[7] + loc4 * dg[7];
756 : :
757 : 0 : loc3 = matr[1] * f + dg[1] * cross;
758 : 0 : loc4 = dg[1] * g + matr[1] * cross;
759 : :
760 : 0 : J_A[3] = loc0 + loc3 * matr[1] + loc4 * dg[1];
761 : 0 : J_B[3] = loc3 * matr[3] + loc4 * dg[3];
762 : 0 : J_B[4] = loc3 * matr[4] + loc4 * dg[4];
763 : 0 : J_C[3] = loc3 * matr[6] + loc4 * dg[6];
764 : 0 : J_C[4] = loc3 * matr[7] + loc4 * dg[7];
765 : :
766 : : /* First diagonal block */
767 : 0 : J_A[0] *= scale[0];
768 : 0 : J_A[1] *= scale[1];
769 : 0 : J_A[3] *= scale[2];
770 : :
771 : 0 : A[0] = -J_A[0] - J_A[1];
772 : 0 : A[4] = -J_A[1] - J_A[3];
773 : :
774 [ # # ]: 0 : h_obj[0][0][0] = -A[0] - A[4];
775 [ # # ]: 0 : h_obj[1][0][0] = A[0];
776 [ # # ]: 0 : h_obj[2][0][0] = A[4];
777 : :
778 [ # # ]: 0 : h_obj[3][0][0] = J_A[0];
779 [ # # ]: 0 : h_obj[4][0][0] = J_A[1];
780 : :
781 [ # # ]: 0 : h_obj[5][0][0] = J_A[3];
782 : :
783 : : /* First off-diagonal block */
784 : 0 : loc2 = matr[8] * loc1;
785 : 0 : J_B[1] += loc2;
786 : 0 : J_B[3] -= loc2;
787 : :
788 : 0 : J_B[0] *= scale[0];
789 : 0 : J_B[1] *= scale[1];
790 : 0 : J_B[3] *= scale[1];
791 : 0 : J_B[4] *= scale[2];
792 : :
793 : 0 : A[0] = -J_B[0] - J_B[3];
794 : 0 : A[4] = -J_B[1] - J_B[4];
795 : :
796 [ # # ]: 0 : h_obj[0][0][1] = -A[0] - A[4];
797 [ # # ]: 0 : h_obj[1][0][1] = A[0];
798 [ # # ]: 0 : h_obj[2][0][1] = A[4];
799 : :
800 [ # # ]: 0 : h_obj[1][1][0] = -J_B[0] - J_B[1];
801 [ # # ]: 0 : h_obj[3][0][1] = J_B[0];
802 [ # # ]: 0 : h_obj[4][0][1] = J_B[1];
803 : :
804 [ # # ]: 0 : h_obj[2][1][0] = -J_B[3] - J_B[4];
805 [ # # ]: 0 : h_obj[4][1][0] = J_B[3];
806 [ # # ]: 0 : h_obj[5][0][1] = J_B[4];
807 : :
808 : : /* Second off-diagonal block */
809 : 0 : loc2 = matr[5] * loc1;
810 : 0 : J_C[1] -= loc2;
811 : 0 : J_C[3] += loc2;
812 : :
813 : 0 : J_C[0] *= scale[0];
814 : 0 : J_C[1] *= scale[1];
815 : 0 : J_C[3] *= scale[1];
816 : 0 : J_C[4] *= scale[2];
817 : :
818 : 0 : A[0] = -J_C[0] - J_C[3];
819 : 0 : A[4] = -J_C[1] - J_C[4];
820 : :
821 [ # # ]: 0 : h_obj[0][0][2] = -A[0] - A[4];
822 [ # # ]: 0 : h_obj[1][0][2] = A[0];
823 [ # # ]: 0 : h_obj[2][0][2] = A[4];
824 : :
825 [ # # ]: 0 : h_obj[1][2][0] = -J_C[0] - J_C[1];
826 [ # # ]: 0 : h_obj[3][0][2] = J_C[0];
827 [ # # ]: 0 : h_obj[4][0][2] = J_C[1];
828 : :
829 [ # # ]: 0 : h_obj[2][2][0] = -J_C[3] - J_C[4];
830 [ # # ]: 0 : h_obj[4][2][0] = J_C[3];
831 [ # # ]: 0 : h_obj[5][0][2] = J_C[4];
832 : :
833 : : /* Second block of rows */
834 : 0 : loc3 = matr[3] * f + dg[3] * cross;
835 : 0 : loc4 = dg[3] * g + matr[3] * cross;
836 : :
837 : 0 : J_A[0] = loc0 + loc3 * matr[3] + loc4 * dg[3];
838 : 0 : J_A[1] = loc3 * matr[4] + loc4 * dg[4];
839 : 0 : J_B[0] = loc3 * matr[6] + loc4 * dg[6];
840 : 0 : J_B[1] = loc3 * matr[7] + loc4 * dg[7];
841 : :
842 : 0 : loc3 = matr[4] * f + dg[4] * cross;
843 : 0 : loc4 = dg[4] * g + matr[4] * cross;
844 : :
845 : 0 : J_A[3] = loc0 + loc3 * matr[4] + loc4 * dg[4];
846 : 0 : J_B[3] = loc3 * matr[6] + loc4 * dg[6];
847 : 0 : J_B[4] = loc3 * matr[7] + loc4 * dg[7];
848 : :
849 : : /* Second diagonal block */
850 : 0 : J_A[0] *= scale[0];
851 : 0 : J_A[1] *= scale[1];
852 : 0 : J_A[3] *= scale[2];
853 : :
854 : 0 : A[0] = -J_A[0] - J_A[1];
855 : 0 : A[4] = -J_A[1] - J_A[3];
856 : :
857 [ # # ]: 0 : h_obj[0][1][1] = -A[0] - A[4];
858 [ # # ]: 0 : h_obj[1][1][1] = A[0];
859 [ # # ]: 0 : h_obj[2][1][1] = A[4];
860 : :
861 [ # # ]: 0 : h_obj[3][1][1] = J_A[0];
862 [ # # ]: 0 : h_obj[4][1][1] = J_A[1];
863 : :
864 [ # # ]: 0 : h_obj[5][1][1] = J_A[3];
865 : :
866 : : /* Third off-diagonal block */
867 : 0 : loc2 = matr[2] * loc1;
868 : 0 : J_B[1] += loc2;
869 : 0 : J_B[3] -= loc2;
870 : :
871 : 0 : J_B[0] *= scale[0];
872 : 0 : J_B[1] *= scale[1];
873 : 0 : J_B[3] *= scale[1];
874 : 0 : J_B[4] *= scale[2];
875 : :
876 : 0 : A[0] = -J_B[0] - J_B[3];
877 : 0 : A[4] = -J_B[1] - J_B[4];
878 : :
879 [ # # ]: 0 : h_obj[0][1][2] = -A[0] - A[4];
880 [ # # ]: 0 : h_obj[1][1][2] = A[0];
881 [ # # ]: 0 : h_obj[2][1][2] = A[4];
882 : :
883 [ # # ]: 0 : h_obj[1][2][1] = -J_B[0] - J_B[1];
884 [ # # ]: 0 : h_obj[3][1][2] = J_B[0];
885 [ # # ]: 0 : h_obj[4][1][2] = J_B[1];
886 : :
887 [ # # ]: 0 : h_obj[2][2][1] = -J_B[3] - J_B[4];
888 [ # # ]: 0 : h_obj[4][2][1] = J_B[3];
889 [ # # ]: 0 : h_obj[5][1][2] = J_B[4];
890 : :
891 : : /* Third block of rows */
892 : 0 : loc3 = matr[6] * f + dg[6] * cross;
893 : 0 : loc4 = dg[6] * g + matr[6] * cross;
894 : :
895 : 0 : J_A[0] = loc0 + loc3 * matr[6] + loc4 * dg[6];
896 : 0 : J_A[1] = loc3 * matr[7] + loc4 * dg[7];
897 : :
898 : 0 : loc3 = matr[7] * f + dg[7] * cross;
899 : 0 : loc4 = dg[7] * g + matr[7] * cross;
900 : :
901 : 0 : J_A[3] = loc0 + loc3 * matr[7] + loc4 * dg[7];
902 : :
903 : : /* Third diagonal block */
904 : 0 : J_A[0] *= scale[0];
905 : 0 : J_A[1] *= scale[1];
906 : 0 : J_A[3] *= scale[2];
907 : :
908 : 0 : A[0] = -J_A[0] - J_A[1];
909 : 0 : A[4] = -J_A[1] - J_A[3];
910 : :
911 [ # # ]: 0 : h_obj[0][2][2] = -A[0] - A[4];
912 [ # # ]: 0 : h_obj[1][2][2] = A[0];
913 [ # # ]: 0 : h_obj[2][2][2] = A[4];
914 : :
915 [ # # ]: 0 : h_obj[3][2][2] = J_A[0];
916 [ # # ]: 0 : h_obj[4][2][2] = J_A[1];
917 : :
918 [ # # ]: 0 : h_obj[5][2][2] = J_A[3];
919 : :
920 : : // completes diagonal blocks.
921 [ # # ]: 0 : h_obj[0].fill_lower_triangle();
922 [ # # ]: 0 : h_obj[3].fill_lower_triangle();
923 [ # # ]: 0 : h_obj[5].fill_lower_triangle();
924 : 0 : return true;
925 : : }
926 : :
927 : : /*****************************************************************************/
928 : : /* The following set of functions reference tetrahedral elements to a */
929 : : /* regular tetrahedron. They are used when assessing the quality of a */
930 : : /* tetrahedral element. A zero return value indicates success, while */
931 : : /* a nonzero value indicates failure. */
932 : : /*****************************************************************************/
933 : :
934 : : /*****************************************************************************/
935 : : /* Function evaluation requires 62 flops. */
936 : : /* Reductions possible when b == 1 or c == 1 */
937 : : /*****************************************************************************/
938 : 57270 : inline bool m_fcn_3e( double& obj, const Vector3D x[4], const double a, const Exponent& b, const Exponent& c )
939 : : {
940 : : double matr[9], f;
941 : : double g;
942 : :
943 : : /* Calculate M = A*inv(W). */
944 [ + - ][ + - ]: 57270 : f = x[1][0] + x[0][0];
945 [ + - ][ + - ]: 57270 : matr[0] = x[1][0] - x[0][0];
946 [ + - ]: 57270 : matr[1] = ( 2.0 * x[2][0] - f ) * isqrt3;
947 [ + - ][ + - ]: 57270 : matr[2] = ( 3.0 * x[3][0] - x[2][0] - f ) * isqrt6;
948 : :
949 [ + - ][ + - ]: 57270 : f = x[1][1] + x[0][1];
950 [ + - ][ + - ]: 57270 : matr[3] = x[1][1] - x[0][1];
951 [ + - ]: 57270 : matr[4] = ( 2.0 * x[2][1] - f ) * isqrt3;
952 [ + - ][ + - ]: 57270 : matr[5] = ( 3.0 * x[3][1] - x[2][1] - f ) * isqrt6;
953 : :
954 [ + - ][ + - ]: 57270 : f = x[1][2] + x[0][2];
955 [ + - ][ + - ]: 57270 : matr[6] = x[1][2] - x[0][2];
956 [ + - ]: 57270 : matr[7] = ( 2.0 * x[2][2] - f ) * isqrt3;
957 [ + - ][ + - ]: 57270 : matr[8] = ( 3.0 * x[3][2] - x[2][2] - f ) * isqrt6;
958 : :
959 : : /* Calculate det(M). */
960 : 57270 : g = matr[0] * ( matr[4] * matr[8] - matr[5] * matr[7] ) + matr[1] * ( matr[5] * matr[6] - matr[3] * matr[8] ) +
961 : 57270 : matr[2] * ( matr[3] * matr[7] - matr[4] * matr[6] );
962 [ - + ]: 57270 : if( g < MSQ_MIN )
963 : : {
964 : 0 : obj = g;
965 : 0 : return false;
966 : : }
967 : :
968 : : /* Calculate norm(M). */
969 : 114540 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[2] * matr[2] + matr[3] * matr[3] + matr[4] * matr[4] +
970 : 114540 : matr[5] * matr[5] + matr[6] * matr[6] + matr[7] * matr[7] + matr[8] * matr[8];
971 : :
972 : : /* Calculate objective function. */
973 [ + - ][ + - ]: 57270 : obj = a * b.raise( f ) * c.raise( g );
974 : 57270 : return true;
975 : : }
976 : :
977 : : /*****************************************************************************/
978 : : /* Gradient evaluation requires 133 flops. */
979 : : /* Reductions possible when b == 1 or c == 1 */
980 : : /*****************************************************************************/
981 : 136059 : inline bool g_fcn_3e( double& obj, Vector3D g_obj[4], const Vector3D x[4], const double a, const Exponent& b,
982 : : const Exponent& c )
983 : : {
984 : : double matr[9], f;
985 : : double adj_m[9], g;
986 : : double loc1, loc2, loc3;
987 : :
988 : : /* Calculate M = A*inv(W). */
989 [ + - ][ + - ]: 136059 : f = x[1][0] + x[0][0];
990 [ + - ][ + - ]: 136059 : matr[0] = x[1][0] - x[0][0];
991 [ + - ]: 136059 : matr[1] = ( 2.0 * x[2][0] - f ) * isqrt3;
992 [ + - ][ + - ]: 136059 : matr[2] = ( 3.0 * x[3][0] - x[2][0] - f ) * isqrt6;
993 : :
994 [ + - ][ + - ]: 136059 : f = x[1][1] + x[0][1];
995 [ + - ][ + - ]: 136059 : matr[3] = x[1][1] - x[0][1];
996 [ + - ]: 136059 : matr[4] = ( 2.0 * x[2][1] - f ) * isqrt3;
997 [ + - ][ + - ]: 136059 : matr[5] = ( 3.0 * x[3][1] - x[2][1] - f ) * isqrt6;
998 : :
999 [ + - ][ + - ]: 136059 : f = x[1][2] + x[0][2];
1000 [ + - ][ + - ]: 136059 : matr[6] = x[1][2] - x[0][2];
1001 [ + - ]: 136059 : matr[7] = ( 2.0 * x[2][2] - f ) * isqrt3;
1002 [ + - ][ + - ]: 136059 : matr[8] = ( 3.0 * x[3][2] - x[2][2] - f ) * isqrt6;
1003 : :
1004 : : /* Calculate det(M). */
1005 : 136059 : loc1 = matr[4] * matr[8] - matr[5] * matr[7];
1006 : 136059 : loc2 = matr[5] * matr[6] - matr[3] * matr[8];
1007 : 136059 : loc3 = matr[3] * matr[7] - matr[4] * matr[6];
1008 : 136059 : g = matr[0] * loc1 + matr[1] * loc2 + matr[2] * loc3;
1009 [ - + ]: 136059 : if( g < MSQ_MIN )
1010 : : {
1011 : 0 : obj = g;
1012 : 0 : return false;
1013 : : }
1014 : :
1015 : : /* Calculate norm(M). */
1016 : 272118 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[2] * matr[2] + matr[3] * matr[3] + matr[4] * matr[4] +
1017 : 272118 : matr[5] * matr[5] + matr[6] * matr[6] + matr[7] * matr[7] + matr[8] * matr[8];
1018 : :
1019 : : /* Calculate objective function. */
1020 [ + - ][ + - ]: 136059 : obj = a * b.raise( f ) * c.raise( g );
1021 : :
1022 : : /* Calculate the derivative of the objective function. */
1023 [ + - ]: 136059 : f = b.value() * obj / f; /* Constant on nabla f */
1024 [ + - ]: 136059 : g = c.value() * obj / g; /* Constant on nable g */
1025 : 136059 : f *= 2.0; /* Modification for nabla f */
1026 : :
1027 : 136059 : adj_m[0] = matr[0] * f + loc1 * g;
1028 : 136059 : adj_m[1] = matr[1] * f + loc2 * g;
1029 : 136059 : adj_m[2] = matr[2] * f + loc3 * g;
1030 : :
1031 : 136059 : loc1 = matr[0] * g;
1032 : 136059 : loc2 = matr[1] * g;
1033 : 136059 : loc3 = matr[2] * g;
1034 : :
1035 : 136059 : adj_m[3] = matr[3] * f + loc3 * matr[7] - loc2 * matr[8];
1036 : 136059 : adj_m[4] = matr[4] * f + loc1 * matr[8] - loc3 * matr[6];
1037 : 136059 : adj_m[5] = matr[5] * f + loc2 * matr[6] - loc1 * matr[7];
1038 : :
1039 : 136059 : adj_m[6] = matr[6] * f + loc2 * matr[5] - loc3 * matr[4];
1040 : 136059 : adj_m[7] = matr[7] * f + loc3 * matr[3] - loc1 * matr[5];
1041 : 136059 : adj_m[8] = matr[8] * f + loc1 * matr[4] - loc2 * matr[3];
1042 : :
1043 : 136059 : loc1 = isqrt3 * adj_m[1];
1044 : 136059 : loc2 = isqrt6 * adj_m[2];
1045 : 136059 : loc3 = loc1 + loc2;
1046 [ + - ]: 136059 : g_obj[0][0] = -adj_m[0] - loc3;
1047 [ + - ]: 136059 : g_obj[1][0] = adj_m[0] - loc3;
1048 [ + - ]: 136059 : g_obj[2][0] = 2.0 * loc1 - loc2;
1049 [ + - ]: 136059 : g_obj[3][0] = 3.0 * loc2;
1050 : :
1051 : 136059 : loc1 = isqrt3 * adj_m[4];
1052 : 136059 : loc2 = isqrt6 * adj_m[5];
1053 : 136059 : loc3 = loc1 + loc2;
1054 [ + - ]: 136059 : g_obj[0][1] = -adj_m[3] - loc3;
1055 [ + - ]: 136059 : g_obj[1][1] = adj_m[3] - loc3;
1056 [ + - ]: 136059 : g_obj[2][1] = 2.0 * loc1 - loc2;
1057 [ + - ]: 136059 : g_obj[3][1] = 3.0 * loc2;
1058 : :
1059 : 136059 : loc1 = isqrt3 * adj_m[7];
1060 : 136059 : loc2 = isqrt6 * adj_m[8];
1061 : 136059 : loc3 = loc1 + loc2;
1062 [ + - ]: 136059 : g_obj[0][2] = -adj_m[6] - loc3;
1063 [ + - ]: 136059 : g_obj[1][2] = adj_m[6] - loc3;
1064 [ + - ]: 136059 : g_obj[2][2] = 2.0 * loc1 - loc2;
1065 [ + - ]: 136059 : g_obj[3][2] = 3.0 * loc2;
1066 : 136059 : return true;
1067 : : }
1068 : :
1069 : : /*****************************************************************************/
1070 : : /* Hessian evaluation requires 634 flops. */
1071 : : /* Reductions possible when b == 1 or c == 1 */
1072 : : /*****************************************************************************/
1073 : 125015 : inline bool h_fcn_3e( double& obj, Vector3D g_obj[4], Matrix3D h_obj[10], const Vector3D x[4], const double a,
1074 : : const Exponent& b, const Exponent& c )
1075 : : {
1076 : : double matr[9], f;
1077 : : double adj_m[9], g;
1078 : : double dg[9], loc0, loc1, loc2, loc3, loc4;
1079 : : double A[12], J_A[6], J_B[9], J_C[9], cross;
1080 : :
1081 : : /* Calculate M = A*inv(W). */
1082 [ + - ][ + - ]: 125015 : f = x[1][0] + x[0][0];
1083 [ + - ][ + - ]: 125015 : matr[0] = x[1][0] - x[0][0];
1084 [ + - ]: 125015 : matr[1] = ( 2.0 * x[2][0] - f ) * isqrt3;
1085 [ + - ][ + - ]: 125015 : matr[2] = ( 3.0 * x[3][0] - x[2][0] - f ) * isqrt6;
1086 : :
1087 [ + - ][ + - ]: 125015 : f = x[1][1] + x[0][1];
1088 [ + - ][ + - ]: 125015 : matr[3] = x[1][1] - x[0][1];
1089 [ + - ]: 125015 : matr[4] = ( 2.0 * x[2][1] - f ) * isqrt3;
1090 [ + - ][ + - ]: 125015 : matr[5] = ( 3.0 * x[3][1] - x[2][1] - f ) * isqrt6;
1091 : :
1092 [ + - ][ + - ]: 125015 : f = x[1][2] + x[0][2];
1093 [ + - ][ + - ]: 125015 : matr[6] = x[1][2] - x[0][2];
1094 [ + - ]: 125015 : matr[7] = ( 2.0 * x[2][2] - f ) * isqrt3;
1095 [ + - ][ + - ]: 125015 : matr[8] = ( 3.0 * x[3][2] - x[2][2] - f ) * isqrt6;
1096 : :
1097 : : /* Calculate det(M). */
1098 : 125015 : dg[0] = matr[4] * matr[8] - matr[5] * matr[7];
1099 : 125015 : dg[1] = matr[5] * matr[6] - matr[3] * matr[8];
1100 : 125015 : dg[2] = matr[3] * matr[7] - matr[4] * matr[6];
1101 : 125015 : g = matr[0] * dg[0] + matr[1] * dg[1] + matr[2] * dg[2];
1102 [ - + ]: 125015 : if( g < MSQ_MIN )
1103 : : {
1104 : 0 : obj = g;
1105 : 0 : return false;
1106 : : }
1107 : :
1108 : : /* Calculate norm(M). */
1109 : 250030 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[2] * matr[2] + matr[3] * matr[3] + matr[4] * matr[4] +
1110 : 250030 : matr[5] * matr[5] + matr[6] * matr[6] + matr[7] * matr[7] + matr[8] * matr[8];
1111 : :
1112 : 125015 : loc3 = f;
1113 : 125015 : loc4 = g;
1114 : :
1115 : : /* Calculate objective function. */
1116 [ + - ][ + - ]: 125015 : obj = a * b.raise( f ) * c.raise( g );
1117 : :
1118 : : /* Calculate the derivative of the objective function. */
1119 [ + - ]: 125015 : f = b.value() * obj / f; /* Constant on nabla f */
1120 [ + - ]: 125015 : g = c.value() * obj / g; /* Constant on nable g */
1121 : 125015 : f *= 2.0; /* Modification for nabla f */
1122 : :
1123 : 125015 : dg[3] = matr[2] * matr[7] - matr[1] * matr[8];
1124 : 125015 : dg[4] = matr[0] * matr[8] - matr[2] * matr[6];
1125 : 125015 : dg[5] = matr[1] * matr[6] - matr[0] * matr[7];
1126 : 125015 : dg[6] = matr[1] * matr[5] - matr[2] * matr[4];
1127 : 125015 : dg[7] = matr[2] * matr[3] - matr[0] * matr[5];
1128 : 125015 : dg[8] = matr[0] * matr[4] - matr[1] * matr[3];
1129 : :
1130 : 125015 : adj_m[0] = matr[0] * f + dg[0] * g;
1131 : 125015 : adj_m[1] = matr[1] * f + dg[1] * g;
1132 : 125015 : adj_m[2] = matr[2] * f + dg[2] * g;
1133 : 125015 : adj_m[3] = matr[3] * f + dg[3] * g;
1134 : 125015 : adj_m[4] = matr[4] * f + dg[4] * g;
1135 : 125015 : adj_m[5] = matr[5] * f + dg[5] * g;
1136 : 125015 : adj_m[6] = matr[6] * f + dg[6] * g;
1137 : 125015 : adj_m[7] = matr[7] * f + dg[7] * g;
1138 : 125015 : adj_m[8] = matr[8] * f + dg[8] * g;
1139 : :
1140 : 125015 : loc0 = isqrt3 * adj_m[1];
1141 : 125015 : loc1 = isqrt6 * adj_m[2];
1142 : 125015 : loc2 = loc0 + loc1;
1143 [ + - ]: 125015 : g_obj[0][0] = -adj_m[0] - loc2;
1144 [ + - ]: 125015 : g_obj[1][0] = adj_m[0] - loc2;
1145 [ + - ]: 125015 : g_obj[2][0] = 2.0 * loc0 - loc1;
1146 [ + - ]: 125015 : g_obj[3][0] = 3.0 * loc1;
1147 : :
1148 : 125015 : loc0 = isqrt3 * adj_m[4];
1149 : 125015 : loc1 = isqrt6 * adj_m[5];
1150 : 125015 : loc2 = loc0 + loc1;
1151 [ + - ]: 125015 : g_obj[0][1] = -adj_m[3] - loc2;
1152 [ + - ]: 125015 : g_obj[1][1] = adj_m[3] - loc2;
1153 [ + - ]: 125015 : g_obj[2][1] = 2.0 * loc0 - loc1;
1154 [ + - ]: 125015 : g_obj[3][1] = 3.0 * loc1;
1155 : :
1156 : 125015 : loc0 = isqrt3 * adj_m[7];
1157 : 125015 : loc1 = isqrt6 * adj_m[8];
1158 : 125015 : loc2 = loc0 + loc1;
1159 [ + - ]: 125015 : g_obj[0][2] = -adj_m[6] - loc2;
1160 [ + - ]: 125015 : g_obj[1][2] = adj_m[6] - loc2;
1161 [ + - ]: 125015 : g_obj[2][2] = 2.0 * loc0 - loc1;
1162 [ + - ]: 125015 : g_obj[3][2] = 3.0 * loc1;
1163 : :
1164 : : /* Calculate the hessian of the objective. */
1165 : 125015 : loc0 = f; /* Constant on nabla^2 f */
1166 : 125015 : loc1 = g; /* Constant on nabla^2 g */
1167 [ + - ]: 125015 : cross = f * c.value() / loc4; /* Constant on nabla g nabla f */
1168 [ + - ]: 125015 : f = f * ( b.value() - 1 ) / loc3; /* Constant on nabla f nabla f */
1169 [ + - ]: 125015 : g = g * ( c.value() - 1 ) / loc4; /* Constant on nabla g nabla g */
1170 : 125015 : f *= 2.0; /* Modification for nabla f */
1171 : :
1172 : : /* First block of rows */
1173 : 125015 : loc3 = matr[0] * f + dg[0] * cross;
1174 : 125015 : loc4 = dg[0] * g + matr[0] * cross;
1175 : :
1176 : 125015 : J_A[0] = loc0 + loc3 * matr[0] + loc4 * dg[0];
1177 : 125015 : J_A[1] = loc3 * matr[1] + loc4 * dg[1];
1178 : 125015 : J_A[2] = loc3 * matr[2] + loc4 * dg[2];
1179 : 125015 : J_B[0] = loc3 * matr[3] + loc4 * dg[3];
1180 : 125015 : J_B[1] = loc3 * matr[4] + loc4 * dg[4];
1181 : 125015 : J_B[2] = loc3 * matr[5] + loc4 * dg[5];
1182 : 125015 : J_C[0] = loc3 * matr[6] + loc4 * dg[6];
1183 : 125015 : J_C[1] = loc3 * matr[7] + loc4 * dg[7];
1184 : 125015 : J_C[2] = loc3 * matr[8] + loc4 * dg[8];
1185 : :
1186 : 125015 : loc3 = matr[1] * f + dg[1] * cross;
1187 : 125015 : loc4 = dg[1] * g + matr[1] * cross;
1188 : :
1189 : 125015 : J_A[3] = loc0 + loc3 * matr[1] + loc4 * dg[1];
1190 : 125015 : J_A[4] = loc3 * matr[2] + loc4 * dg[2];
1191 : 125015 : J_B[3] = loc3 * matr[3] + loc4 * dg[3];
1192 : 125015 : J_B[4] = loc3 * matr[4] + loc4 * dg[4];
1193 : 125015 : J_B[5] = loc3 * matr[5] + loc4 * dg[5];
1194 : 125015 : J_C[3] = loc3 * matr[6] + loc4 * dg[6];
1195 : 125015 : J_C[4] = loc3 * matr[7] + loc4 * dg[7];
1196 : 125015 : J_C[5] = loc3 * matr[8] + loc4 * dg[8];
1197 : :
1198 : 125015 : loc3 = matr[2] * f + dg[2] * cross;
1199 : 125015 : loc4 = dg[2] * g + matr[2] * cross;
1200 : :
1201 : 125015 : J_A[5] = loc0 + loc3 * matr[2] + loc4 * dg[2];
1202 : 125015 : J_B[6] = loc3 * matr[3] + loc4 * dg[3];
1203 : 125015 : J_B[7] = loc3 * matr[4] + loc4 * dg[4];
1204 : 125015 : J_B[8] = loc3 * matr[5] + loc4 * dg[5];
1205 : 125015 : J_C[6] = loc3 * matr[6] + loc4 * dg[6];
1206 : 125015 : J_C[7] = loc3 * matr[7] + loc4 * dg[7];
1207 : 125015 : J_C[8] = loc3 * matr[8] + loc4 * dg[8];
1208 : :
1209 : : /* First diagonal block */
1210 : 125015 : loc2 = isqrt3 * J_A[1];
1211 : 125015 : loc3 = isqrt6 * J_A[2];
1212 : 125015 : loc4 = loc2 + loc3;
1213 : :
1214 : 125015 : A[0] = -J_A[0] - loc4;
1215 : 125015 : A[1] = J_A[0] - loc4;
1216 : :
1217 : 125015 : loc2 = isqrt3 * J_A[3];
1218 : 125015 : loc3 = isqrt6 * J_A[4];
1219 : 125015 : loc4 = loc2 + loc3;
1220 : :
1221 : 125015 : A[4] = -J_A[1] - loc4;
1222 : 125015 : A[5] = J_A[1] - loc4;
1223 : 125015 : A[6] = 2.0 * loc2 - loc3;
1224 : :
1225 : 125015 : loc2 = isqrt3 * J_A[4];
1226 : 125015 : loc3 = isqrt6 * J_A[5];
1227 : 125015 : loc4 = loc2 + loc3;
1228 : :
1229 : 125015 : A[8] = -J_A[2] - loc4;
1230 : 125015 : A[9] = J_A[2] - loc4;
1231 : 125015 : A[10] = 2.0 * loc2 - loc3;
1232 : 125015 : A[11] = 3.0 * loc3;
1233 : :
1234 : 125015 : loc2 = isqrt3 * A[4];
1235 : 125015 : loc3 = isqrt6 * A[8];
1236 : 125015 : loc4 = loc2 + loc3;
1237 : :
1238 [ + - ]: 125015 : h_obj[0][0][0] = -A[0] - loc4;
1239 [ + - ]: 125015 : h_obj[1][0][0] = A[0] - loc4;
1240 [ + - ]: 125015 : h_obj[2][0][0] = 2.0 * loc2 - loc3;
1241 [ + - ]: 125015 : h_obj[3][0][0] = 3.0 * loc3;
1242 : :
1243 : 125015 : loc2 = isqrt3 * A[5];
1244 : 125015 : loc3 = isqrt6 * A[9];
1245 : :
1246 [ + - ]: 125015 : h_obj[4][0][0] = A[1] - loc2 - loc3;
1247 [ + - ]: 125015 : h_obj[5][0][0] = 2.0 * loc2 - loc3;
1248 [ + - ]: 125015 : h_obj[6][0][0] = 3.0 * loc3;
1249 : :
1250 : 125015 : loc3 = isqrt6 * A[10];
1251 [ + - ]: 125015 : h_obj[7][0][0] = tisqrt3 * A[6] - loc3;
1252 [ + - ]: 125015 : h_obj[8][0][0] = 3.0 * loc3;
1253 : :
1254 [ + - ]: 125015 : h_obj[9][0][0] = tisqrt6 * A[11];
1255 : :
1256 : : /* First off-diagonal block */
1257 : 125015 : loc2 = matr[8] * loc1;
1258 : 125015 : J_B[1] += loc2;
1259 : 125015 : J_B[3] -= loc2;
1260 : :
1261 : 125015 : loc2 = matr[7] * loc1;
1262 : 125015 : J_B[2] -= loc2;
1263 : 125015 : J_B[6] += loc2;
1264 : :
1265 : 125015 : loc2 = matr[6] * loc1;
1266 : 125015 : J_B[5] += loc2;
1267 : 125015 : J_B[7] -= loc2;
1268 : :
1269 : 125015 : loc2 = isqrt3 * J_B[3];
1270 : 125015 : loc3 = isqrt6 * J_B[6];
1271 : 125015 : loc4 = loc2 + loc3;
1272 : :
1273 : 125015 : A[0] = -J_B[0] - loc4;
1274 : 125015 : A[1] = J_B[0] - loc4;
1275 : 125015 : A[2] = 2.0 * loc2 - loc3;
1276 : 125015 : A[3] = 3.0 * loc3;
1277 : :
1278 : 125015 : loc2 = isqrt3 * J_B[4];
1279 : 125015 : loc3 = isqrt6 * J_B[7];
1280 : 125015 : loc4 = loc2 + loc3;
1281 : :
1282 : 125015 : A[4] = -J_B[1] - loc4;
1283 : 125015 : A[5] = J_B[1] - loc4;
1284 : 125015 : A[6] = 2.0 * loc2 - loc3;
1285 : 125015 : A[7] = 3.0 * loc3;
1286 : :
1287 : 125015 : loc2 = isqrt3 * J_B[5];
1288 : 125015 : loc3 = isqrt6 * J_B[8];
1289 : 125015 : loc4 = loc2 + loc3;
1290 : :
1291 : 125015 : A[8] = -J_B[2] - loc4;
1292 : 125015 : A[9] = J_B[2] - loc4;
1293 : 125015 : A[10] = 2.0 * loc2 - loc3;
1294 : 125015 : A[11] = 3.0 * loc3;
1295 : :
1296 : 125015 : loc2 = isqrt3 * A[4];
1297 : 125015 : loc3 = isqrt6 * A[8];
1298 : 125015 : loc4 = loc2 + loc3;
1299 : :
1300 [ + - ]: 125015 : h_obj[0][0][1] = -A[0] - loc4;
1301 [ + - ]: 125015 : h_obj[1][0][1] = A[0] - loc4;
1302 [ + - ]: 125015 : h_obj[2][0][1] = 2.0 * loc2 - loc3;
1303 [ + - ]: 125015 : h_obj[3][0][1] = 3.0 * loc3;
1304 : :
1305 : 125015 : loc2 = isqrt3 * A[5];
1306 : 125015 : loc3 = isqrt6 * A[9];
1307 : 125015 : loc4 = loc2 + loc3;
1308 : :
1309 [ + - ]: 125015 : h_obj[1][1][0] = -A[1] - loc4;
1310 [ + - ]: 125015 : h_obj[4][0][1] = A[1] - loc4;
1311 [ + - ]: 125015 : h_obj[5][0][1] = 2.0 * loc2 - loc3;
1312 [ + - ]: 125015 : h_obj[6][0][1] = 3.0 * loc3;
1313 : :
1314 : 125015 : loc2 = isqrt3 * A[6];
1315 : 125015 : loc3 = isqrt6 * A[10];
1316 : 125015 : loc4 = loc2 + loc3;
1317 : :
1318 [ + - ]: 125015 : h_obj[2][1][0] = -A[2] - loc4;
1319 [ + - ]: 125015 : h_obj[5][1][0] = A[2] - loc4;
1320 [ + - ]: 125015 : h_obj[7][0][1] = 2.0 * loc2 - loc3;
1321 [ + - ]: 125015 : h_obj[8][0][1] = 3.0 * loc3;
1322 : :
1323 : 125015 : loc2 = isqrt3 * A[7];
1324 : 125015 : loc3 = isqrt6 * A[11];
1325 : 125015 : loc4 = loc2 + loc3;
1326 : :
1327 [ + - ]: 125015 : h_obj[3][1][0] = -A[3] - loc4;
1328 [ + - ]: 125015 : h_obj[6][1][0] = A[3] - loc4;
1329 [ + - ]: 125015 : h_obj[8][1][0] = 2.0 * loc2 - loc3;
1330 [ + - ]: 125015 : h_obj[9][0][1] = 3.0 * loc3;
1331 : :
1332 : : /* Second off-diagonal block */
1333 : 125015 : loc2 = matr[5] * loc1;
1334 : 125015 : J_C[1] -= loc2;
1335 : 125015 : J_C[3] += loc2;
1336 : :
1337 : 125015 : loc2 = matr[4] * loc1;
1338 : 125015 : J_C[2] += loc2;
1339 : 125015 : J_C[6] -= loc2;
1340 : :
1341 : 125015 : loc2 = matr[3] * loc1;
1342 : 125015 : J_C[5] -= loc2;
1343 : 125015 : J_C[7] += loc2;
1344 : :
1345 : 125015 : loc2 = isqrt3 * J_C[3];
1346 : 125015 : loc3 = isqrt6 * J_C[6];
1347 : 125015 : loc4 = loc2 + loc3;
1348 : :
1349 : 125015 : A[0] = -J_C[0] - loc4;
1350 : 125015 : A[1] = J_C[0] - loc4;
1351 : 125015 : A[2] = 2.0 * loc2 - loc3;
1352 : 125015 : A[3] = 3.0 * loc3;
1353 : :
1354 : 125015 : loc2 = isqrt3 * J_C[4];
1355 : 125015 : loc3 = isqrt6 * J_C[7];
1356 : 125015 : loc4 = loc2 + loc3;
1357 : :
1358 : 125015 : A[4] = -J_C[1] - loc4;
1359 : 125015 : A[5] = J_C[1] - loc4;
1360 : 125015 : A[6] = 2.0 * loc2 - loc3;
1361 : 125015 : A[7] = 3.0 * loc3;
1362 : :
1363 : 125015 : loc2 = isqrt3 * J_C[5];
1364 : 125015 : loc3 = isqrt6 * J_C[8];
1365 : 125015 : loc4 = loc2 + loc3;
1366 : :
1367 : 125015 : A[8] = -J_C[2] - loc4;
1368 : 125015 : A[9] = J_C[2] - loc4;
1369 : 125015 : A[10] = 2.0 * loc2 - loc3;
1370 : 125015 : A[11] = 3.0 * loc3;
1371 : :
1372 : 125015 : loc2 = isqrt3 * A[4];
1373 : 125015 : loc3 = isqrt6 * A[8];
1374 : 125015 : loc4 = loc2 + loc3;
1375 : :
1376 [ + - ]: 125015 : h_obj[0][0][2] = -A[0] - loc4;
1377 [ + - ]: 125015 : h_obj[1][0][2] = A[0] - loc4;
1378 [ + - ]: 125015 : h_obj[2][0][2] = 2.0 * loc2 - loc3;
1379 [ + - ]: 125015 : h_obj[3][0][2] = 3.0 * loc3;
1380 : :
1381 : 125015 : loc2 = isqrt3 * A[5];
1382 : 125015 : loc3 = isqrt6 * A[9];
1383 : 125015 : loc4 = loc2 + loc3;
1384 : :
1385 [ + - ]: 125015 : h_obj[1][2][0] = -A[1] - loc4;
1386 [ + - ]: 125015 : h_obj[4][0][2] = A[1] - loc4;
1387 [ + - ]: 125015 : h_obj[5][0][2] = 2.0 * loc2 - loc3;
1388 [ + - ]: 125015 : h_obj[6][0][2] = 3.0 * loc3;
1389 : :
1390 : 125015 : loc2 = isqrt3 * A[6];
1391 : 125015 : loc3 = isqrt6 * A[10];
1392 : 125015 : loc4 = loc2 + loc3;
1393 : :
1394 [ + - ]: 125015 : h_obj[2][2][0] = -A[2] - loc4;
1395 [ + - ]: 125015 : h_obj[5][2][0] = A[2] - loc4;
1396 [ + - ]: 125015 : h_obj[7][0][2] = 2.0 * loc2 - loc3;
1397 [ + - ]: 125015 : h_obj[8][0][2] = 3.0 * loc3;
1398 : :
1399 : 125015 : loc2 = isqrt3 * A[7];
1400 : 125015 : loc3 = isqrt6 * A[11];
1401 : 125015 : loc4 = loc2 + loc3;
1402 : :
1403 [ + - ]: 125015 : h_obj[3][2][0] = -A[3] - loc4;
1404 [ + - ]: 125015 : h_obj[6][2][0] = A[3] - loc4;
1405 [ + - ]: 125015 : h_obj[8][2][0] = 2.0 * loc2 - loc3;
1406 [ + - ]: 125015 : h_obj[9][0][2] = 3.0 * loc3;
1407 : :
1408 : : /* Second block of rows */
1409 : 125015 : loc3 = matr[3] * f + dg[3] * cross;
1410 : 125015 : loc4 = dg[3] * g + matr[3] * cross;
1411 : :
1412 : 125015 : J_A[0] = loc0 + loc3 * matr[3] + loc4 * dg[3];
1413 : 125015 : J_A[1] = loc3 * matr[4] + loc4 * dg[4];
1414 : 125015 : J_A[2] = loc3 * matr[5] + loc4 * dg[5];
1415 : 125015 : J_B[0] = loc3 * matr[6] + loc4 * dg[6];
1416 : 125015 : J_B[1] = loc3 * matr[7] + loc4 * dg[7];
1417 : 125015 : J_B[2] = loc3 * matr[8] + loc4 * dg[8];
1418 : :
1419 : 125015 : loc3 = matr[4] * f + dg[4] * cross;
1420 : 125015 : loc4 = dg[4] * g + matr[4] * cross;
1421 : :
1422 : 125015 : J_A[3] = loc0 + loc3 * matr[4] + loc4 * dg[4];
1423 : 125015 : J_A[4] = loc3 * matr[5] + loc4 * dg[5];
1424 : 125015 : J_B[3] = loc3 * matr[6] + loc4 * dg[6];
1425 : 125015 : J_B[4] = loc3 * matr[7] + loc4 * dg[7];
1426 : 125015 : J_B[5] = loc3 * matr[8] + loc4 * dg[8];
1427 : :
1428 : 125015 : loc3 = matr[5] * f + dg[5] * cross;
1429 : 125015 : loc4 = dg[5] * g + matr[5] * cross;
1430 : :
1431 : 125015 : J_A[5] = loc0 + loc3 * matr[5] + loc4 * dg[5];
1432 : 125015 : J_B[6] = loc3 * matr[6] + loc4 * dg[6];
1433 : 125015 : J_B[7] = loc3 * matr[7] + loc4 * dg[7];
1434 : 125015 : J_B[8] = loc3 * matr[8] + loc4 * dg[8];
1435 : :
1436 : : /* Second diagonal block */
1437 : 125015 : loc2 = isqrt3 * J_A[1];
1438 : 125015 : loc3 = isqrt6 * J_A[2];
1439 : 125015 : loc4 = loc2 + loc3;
1440 : :
1441 : 125015 : A[0] = -J_A[0] - loc4;
1442 : 125015 : A[1] = J_A[0] - loc4;
1443 : :
1444 : 125015 : loc2 = isqrt3 * J_A[3];
1445 : 125015 : loc3 = isqrt6 * J_A[4];
1446 : 125015 : loc4 = loc2 + loc3;
1447 : :
1448 : 125015 : A[4] = -J_A[1] - loc4;
1449 : 125015 : A[5] = J_A[1] - loc4;
1450 : 125015 : A[6] = 2.0 * loc2 - loc3;
1451 : :
1452 : 125015 : loc2 = isqrt3 * J_A[4];
1453 : 125015 : loc3 = isqrt6 * J_A[5];
1454 : 125015 : loc4 = loc2 + loc3;
1455 : :
1456 : 125015 : A[8] = -J_A[2] - loc4;
1457 : 125015 : A[9] = J_A[2] - loc4;
1458 : 125015 : A[10] = 2.0 * loc2 - loc3;
1459 : 125015 : A[11] = 3.0 * loc3;
1460 : :
1461 : 125015 : loc2 = isqrt3 * A[4];
1462 : 125015 : loc3 = isqrt6 * A[8];
1463 : 125015 : loc4 = loc2 + loc3;
1464 : :
1465 [ + - ]: 125015 : h_obj[0][1][1] = -A[0] - loc4;
1466 [ + - ]: 125015 : h_obj[1][1][1] = A[0] - loc4;
1467 [ + - ]: 125015 : h_obj[2][1][1] = 2.0 * loc2 - loc3;
1468 [ + - ]: 125015 : h_obj[3][1][1] = 3.0 * loc3;
1469 : :
1470 : 125015 : loc2 = isqrt3 * A[5];
1471 : 125015 : loc3 = isqrt6 * A[9];
1472 : :
1473 [ + - ]: 125015 : h_obj[4][1][1] = A[1] - loc2 - loc3;
1474 [ + - ]: 125015 : h_obj[5][1][1] = 2.0 * loc2 - loc3;
1475 [ + - ]: 125015 : h_obj[6][1][1] = 3.0 * loc3;
1476 : :
1477 : 125015 : loc3 = isqrt6 * A[10];
1478 [ + - ]: 125015 : h_obj[7][1][1] = tisqrt3 * A[6] - loc3;
1479 [ + - ]: 125015 : h_obj[8][1][1] = 3.0 * loc3;
1480 : :
1481 [ + - ]: 125015 : h_obj[9][1][1] = tisqrt6 * A[11];
1482 : :
1483 : : /* Third off-diagonal block */
1484 : 125015 : loc2 = matr[2] * loc1;
1485 : 125015 : J_B[1] += loc2;
1486 : 125015 : J_B[3] -= loc2;
1487 : :
1488 : 125015 : loc2 = matr[1] * loc1;
1489 : 125015 : J_B[2] -= loc2;
1490 : 125015 : J_B[6] += loc2;
1491 : :
1492 : 125015 : loc2 = matr[0] * loc1;
1493 : 125015 : J_B[5] += loc2;
1494 : 125015 : J_B[7] -= loc2;
1495 : :
1496 : 125015 : loc2 = isqrt3 * J_B[3];
1497 : 125015 : loc3 = isqrt6 * J_B[6];
1498 : 125015 : loc4 = loc2 + loc3;
1499 : :
1500 : 125015 : A[0] = -J_B[0] - loc4;
1501 : 125015 : A[1] = J_B[0] - loc4;
1502 : 125015 : A[2] = 2.0 * loc2 - loc3;
1503 : 125015 : A[3] = 3.0 * loc3;
1504 : :
1505 : 125015 : loc2 = isqrt3 * J_B[4];
1506 : 125015 : loc3 = isqrt6 * J_B[7];
1507 : 125015 : loc4 = loc2 + loc3;
1508 : :
1509 : 125015 : A[4] = -J_B[1] - loc4;
1510 : 125015 : A[5] = J_B[1] - loc4;
1511 : 125015 : A[6] = 2.0 * loc2 - loc3;
1512 : 125015 : A[7] = 3.0 * loc3;
1513 : :
1514 : 125015 : loc2 = isqrt3 * J_B[5];
1515 : 125015 : loc3 = isqrt6 * J_B[8];
1516 : 125015 : loc4 = loc2 + loc3;
1517 : :
1518 : 125015 : A[8] = -J_B[2] - loc4;
1519 : 125015 : A[9] = J_B[2] - loc4;
1520 : 125015 : A[10] = 2.0 * loc2 - loc3;
1521 : 125015 : A[11] = 3.0 * loc3;
1522 : :
1523 : 125015 : loc2 = isqrt3 * A[4];
1524 : 125015 : loc3 = isqrt6 * A[8];
1525 : 125015 : loc4 = loc2 + loc3;
1526 : :
1527 [ + - ]: 125015 : h_obj[0][1][2] = -A[0] - loc4;
1528 [ + - ]: 125015 : h_obj[1][1][2] = A[0] - loc4;
1529 [ + - ]: 125015 : h_obj[2][1][2] = 2.0 * loc2 - loc3;
1530 [ + - ]: 125015 : h_obj[3][1][2] = 3.0 * loc3;
1531 : :
1532 : 125015 : loc2 = isqrt3 * A[5];
1533 : 125015 : loc3 = isqrt6 * A[9];
1534 : 125015 : loc4 = loc2 + loc3;
1535 : :
1536 [ + - ]: 125015 : h_obj[1][2][1] = -A[1] - loc4;
1537 [ + - ]: 125015 : h_obj[4][1][2] = A[1] - loc4;
1538 [ + - ]: 125015 : h_obj[5][1][2] = 2.0 * loc2 - loc3;
1539 [ + - ]: 125015 : h_obj[6][1][2] = 3.0 * loc3;
1540 : :
1541 : 125015 : loc2 = isqrt3 * A[6];
1542 : 125015 : loc3 = isqrt6 * A[10];
1543 : 125015 : loc4 = loc2 + loc3;
1544 : :
1545 [ + - ]: 125015 : h_obj[2][2][1] = -A[2] - loc4;
1546 [ + - ]: 125015 : h_obj[5][2][1] = A[2] - loc4;
1547 [ + - ]: 125015 : h_obj[7][1][2] = 2.0 * loc2 - loc3;
1548 [ + - ]: 125015 : h_obj[8][1][2] = 3.0 * loc3;
1549 : :
1550 : 125015 : loc2 = isqrt3 * A[7];
1551 : 125015 : loc3 = isqrt6 * A[11];
1552 : 125015 : loc4 = loc2 + loc3;
1553 : :
1554 [ + - ]: 125015 : h_obj[3][2][1] = -A[3] - loc4;
1555 [ + - ]: 125015 : h_obj[6][2][1] = A[3] - loc4;
1556 [ + - ]: 125015 : h_obj[8][2][1] = 2.0 * loc2 - loc3;
1557 [ + - ]: 125015 : h_obj[9][1][2] = 3.0 * loc3;
1558 : :
1559 : : /* Third block of rows */
1560 : 125015 : loc3 = matr[6] * f + dg[6] * cross;
1561 : 125015 : loc4 = dg[6] * g + matr[6] * cross;
1562 : :
1563 : 125015 : J_A[0] = loc0 + loc3 * matr[6] + loc4 * dg[6];
1564 : 125015 : J_A[1] = loc3 * matr[7] + loc4 * dg[7];
1565 : 125015 : J_A[2] = loc3 * matr[8] + loc4 * dg[8];
1566 : :
1567 : 125015 : loc3 = matr[7] * f + dg[7] * cross;
1568 : 125015 : loc4 = dg[7] * g + matr[7] * cross;
1569 : :
1570 : 125015 : J_A[3] = loc0 + loc3 * matr[7] + loc4 * dg[7];
1571 : 125015 : J_A[4] = loc3 * matr[8] + loc4 * dg[8];
1572 : :
1573 : 125015 : loc3 = matr[8] * f + dg[8] * cross;
1574 : 125015 : loc4 = dg[8] * g + matr[8] * cross;
1575 : :
1576 : 125015 : J_A[5] = loc0 + loc3 * matr[8] + loc4 * dg[8];
1577 : :
1578 : : /* Third diagonal block */
1579 : 125015 : loc2 = isqrt3 * J_A[1];
1580 : 125015 : loc3 = isqrt6 * J_A[2];
1581 : 125015 : loc4 = loc2 + loc3;
1582 : :
1583 : 125015 : A[0] = -J_A[0] - loc4;
1584 : 125015 : A[1] = J_A[0] - loc4;
1585 : :
1586 : 125015 : loc2 = isqrt3 * J_A[3];
1587 : 125015 : loc3 = isqrt6 * J_A[4];
1588 : 125015 : loc4 = loc2 + loc3;
1589 : :
1590 : 125015 : A[4] = -J_A[1] - loc4;
1591 : 125015 : A[5] = J_A[1] - loc4;
1592 : 125015 : A[6] = 2.0 * loc2 - loc3;
1593 : :
1594 : 125015 : loc2 = isqrt3 * J_A[4];
1595 : 125015 : loc3 = isqrt6 * J_A[5];
1596 : 125015 : loc4 = loc2 + loc3;
1597 : :
1598 : 125015 : A[8] = -J_A[2] - loc4;
1599 : 125015 : A[9] = J_A[2] - loc4;
1600 : 125015 : A[10] = 2.0 * loc2 - loc3;
1601 : 125015 : A[11] = 3.0 * loc3;
1602 : :
1603 : 125015 : loc2 = isqrt3 * A[4];
1604 : 125015 : loc3 = isqrt6 * A[8];
1605 : 125015 : loc4 = loc2 + loc3;
1606 : :
1607 [ + - ]: 125015 : h_obj[0][2][2] = -A[0] - loc4;
1608 [ + - ]: 125015 : h_obj[1][2][2] = A[0] - loc4;
1609 [ + - ]: 125015 : h_obj[2][2][2] = 2.0 * loc2 - loc3;
1610 [ + - ]: 125015 : h_obj[3][2][2] = 3.0 * loc3;
1611 : :
1612 : 125015 : loc2 = isqrt3 * A[5];
1613 : 125015 : loc3 = isqrt6 * A[9];
1614 : :
1615 [ + - ]: 125015 : h_obj[4][2][2] = A[1] - loc2 - loc3;
1616 [ + - ]: 125015 : h_obj[5][2][2] = 2.0 * loc2 - loc3;
1617 [ + - ]: 125015 : h_obj[6][2][2] = 3.0 * loc3;
1618 : :
1619 : 125015 : loc3 = isqrt6 * A[10];
1620 [ + - ]: 125015 : h_obj[7][2][2] = tisqrt3 * A[6] - loc3;
1621 [ + - ]: 125015 : h_obj[8][2][2] = 3.0 * loc3;
1622 : :
1623 [ + - ]: 125015 : h_obj[9][2][2] = tisqrt6 * A[11];
1624 : :
1625 : : // completes diagonal blocks.
1626 [ + - ]: 125015 : h_obj[0].fill_lower_triangle();
1627 [ + - ]: 125015 : h_obj[4].fill_lower_triangle();
1628 [ + - ]: 125015 : h_obj[7].fill_lower_triangle();
1629 [ + - ]: 125015 : h_obj[9].fill_lower_triangle();
1630 : 125015 : return true;
1631 : : }
1632 : :
1633 : : /*****************************************************************************/
1634 : : /* The following set of functions reference tetrahedral elements to a */
1635 : : /* equilateral tetrahedron. These functions are specialized toward */
1636 : : /* obtaining the gradient and Hessian with respect to a single vertex. */
1637 : : /*****************************************************************************/
1638 : :
1639 : : inline bool g_fcn_3e_v3( double& obj, Vector3D& g_obj, const Vector3D x[4], const double a, const Exponent& b,
1640 : : const Exponent& c )
1641 : : {
1642 : : double matr[9], f, g;
1643 : : double loc1, loc2, loc3;
1644 : :
1645 : : /* Calculate M = A*inv(W). */
1646 : : f = x[1][0] + x[0][0];
1647 : : matr[0] = x[1][0] - x[0][0];
1648 : : matr[1] = ( 2.0 * x[2][0] - f ) * isqrt3;
1649 : : matr[2] = ( 3.0 * x[3][0] - x[2][0] - f ) * isqrt6;
1650 : :
1651 : : f = x[1][1] + x[0][1];
1652 : : matr[3] = x[1][1] - x[0][1];
1653 : : matr[4] = ( 2.0 * x[2][1] - f ) * isqrt3;
1654 : : matr[5] = ( 3.0 * x[3][1] - x[2][1] - f ) * isqrt6;
1655 : :
1656 : : f = x[1][2] + x[0][2];
1657 : : matr[6] = x[1][2] - x[0][2];
1658 : : matr[7] = ( 2.0 * x[2][2] - f ) * isqrt3;
1659 : : matr[8] = ( 3.0 * x[3][2] - x[2][2] - f ) * isqrt6;
1660 : :
1661 : : /* Calculate det(M). */
1662 : : loc1 = matr[4] * matr[8] - matr[5] * matr[7];
1663 : : loc2 = matr[5] * matr[6] - matr[3] * matr[8];
1664 : : loc3 = matr[3] * matr[7] - matr[4] * matr[6];
1665 : : g = matr[0] * loc1 + matr[1] * loc2 + matr[2] * loc3;
1666 : : if( g < MSQ_MIN )
1667 : : {
1668 : : obj = g;
1669 : : return false;
1670 : : }
1671 : :
1672 : : /* Calculate norm(M). */
1673 : : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[2] * matr[2] + matr[3] * matr[3] + matr[4] * matr[4] +
1674 : : matr[5] * matr[5] + matr[6] * matr[6] + matr[7] * matr[7] + matr[8] * matr[8];
1675 : :
1676 : : /* Calculate objective function. */
1677 : : obj = a * b.raise( f ) * c.raise( g );
1678 : :
1679 : : /* Calculate the derivative of the objective function. */
1680 : : f = b.value() * obj / f; /* Constant on nabla f */
1681 : : g = c.value() * obj / g; /* Constant on nable g */
1682 : : f *= 2.0; /* Modification for nabla f */
1683 : :
1684 : : g_obj[0] = tisqrt6 * ( matr[2] * f + loc3 * g );
1685 : :
1686 : : loc1 = matr[0] * g;
1687 : : loc2 = matr[1] * g;
1688 : :
1689 : : g_obj[1] = tisqrt6 * ( matr[5] * f + loc2 * matr[6] - loc1 * matr[7] );
1690 : : g_obj[2] = tisqrt6 * ( matr[8] * f + loc1 * matr[4] - loc2 * matr[3] );
1691 : : return true;
1692 : : }
1693 : :
1694 : : inline bool g_fcn_3e_v0( double& obj, Vector3D& g_obj, const Vector3D x[4], const double a, const Exponent& b,
1695 : : const Exponent& c )
1696 : : {
1697 : : static Vector3D my_x[4];
1698 : :
1699 : : my_x[0] = x[1];
1700 : : my_x[1] = x[3];
1701 : : my_x[2] = x[2];
1702 : : my_x[3] = x[0];
1703 : : return g_fcn_3e_v3( obj, g_obj, my_x, a, b, c );
1704 : : }
1705 : :
1706 : : inline bool g_fcn_3e_v1( double& obj, Vector3D& g_obj, const Vector3D x[4], const double a, const Exponent& b,
1707 : : const Exponent& c )
1708 : : {
1709 : : static Vector3D my_x[4];
1710 : :
1711 : : my_x[0] = x[0];
1712 : : my_x[1] = x[2];
1713 : : my_x[2] = x[3];
1714 : : my_x[3] = x[1];
1715 : : return g_fcn_3e_v3( obj, g_obj, my_x, a, b, c );
1716 : : }
1717 : :
1718 : : inline bool g_fcn_3e_v2( double& obj, Vector3D& g_obj, const Vector3D x[4], const double a, const Exponent& b,
1719 : : const Exponent& c )
1720 : : {
1721 : : static Vector3D my_x[4];
1722 : :
1723 : : my_x[0] = x[1];
1724 : : my_x[1] = x[0];
1725 : : my_x[2] = x[3];
1726 : : my_x[3] = x[2];
1727 : : return g_fcn_3e_v3( obj, g_obj, my_x, a, b, c );
1728 : : }
1729 : :
1730 : : inline bool h_fcn_3e_v3( double& obj, Vector3D& g_obj, Matrix3D& h_obj, const Vector3D x[4], const double a,
1731 : : const Exponent& b, const Exponent& c )
1732 : : {
1733 : : double matr[9], f, g;
1734 : : double dg[9], loc0, /*loc1,*/ loc3, loc4;
1735 : : double cross;
1736 : :
1737 : : /* Calculate M = A*inv(W). */
1738 : : f = x[1][0] + x[0][0];
1739 : : matr[0] = x[1][0] - x[0][0];
1740 : : matr[1] = ( 2.0 * x[2][0] - f ) * isqrt3;
1741 : : matr[2] = ( 3.0 * x[3][0] - x[2][0] - f ) * isqrt6;
1742 : :
1743 : : f = x[1][1] + x[0][1];
1744 : : matr[3] = x[1][1] - x[0][1];
1745 : : matr[4] = ( 2.0 * x[2][1] - f ) * isqrt3;
1746 : : matr[5] = ( 3.0 * x[3][1] - x[2][1] - f ) * isqrt6;
1747 : :
1748 : : f = x[1][2] + x[0][2];
1749 : : matr[6] = x[1][2] - x[0][2];
1750 : : matr[7] = ( 2.0 * x[2][2] - f ) * isqrt3;
1751 : : matr[8] = ( 3.0 * x[3][2] - x[2][2] - f ) * isqrt6;
1752 : :
1753 : : /* Calculate det(M). */
1754 : : dg[0] = matr[4] * matr[8] - matr[5] * matr[7];
1755 : : dg[1] = matr[5] * matr[6] - matr[3] * matr[8];
1756 : : dg[2] = matr[3] * matr[7] - matr[4] * matr[6];
1757 : : g = matr[0] * dg[0] + matr[1] * dg[1] + matr[2] * dg[2];
1758 : : if( g < MSQ_MIN )
1759 : : {
1760 : : obj = g;
1761 : : return false;
1762 : : }
1763 : :
1764 : : /* Calculate norm(M). */
1765 : : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[2] * matr[2] + matr[3] * matr[3] + matr[4] * matr[4] +
1766 : : matr[5] * matr[5] + matr[6] * matr[6] + matr[7] * matr[7] + matr[8] * matr[8];
1767 : :
1768 : : loc3 = f;
1769 : : loc4 = g;
1770 : :
1771 : : /* Calculate objective function. */
1772 : : obj = a * b.raise( f ) * c.raise( g );
1773 : :
1774 : : /* Calculate the derivative of the objective function. */
1775 : : f = b.value() * obj / f; /* Constant on nabla f */
1776 : : g = c.value() * obj / g; /* Constant on nable g */
1777 : : f *= 2.0; /* Modification for nabla f */
1778 : :
1779 : : dg[5] = matr[1] * matr[6] - matr[0] * matr[7];
1780 : : dg[8] = matr[0] * matr[4] - matr[1] * matr[3];
1781 : :
1782 : : g_obj[0] = tisqrt6 * ( matr[2] * f + dg[2] * g );
1783 : : g_obj[1] = tisqrt6 * ( matr[5] * f + dg[5] * g );
1784 : : g_obj[2] = tisqrt6 * ( matr[8] * f + dg[8] * g );
1785 : :
1786 : : /* Calculate the hessian of the objective. */
1787 : : loc0 = f; /* Constant on nabla^2 f */
1788 : : // loc1 = g; /* Constant on nabla^2 g */
1789 : : cross = f * c.value() / loc4; /* Constant on nabla g nabla f */
1790 : : f = f * ( b.value() - 1 ) / loc3; /* Constant on nabla f nabla f */
1791 : : g = g * ( c.value() - 1 ) / loc4; /* Constant on nabla g nabla g */
1792 : : f *= 2.0; /* Modification for nabla f */
1793 : :
1794 : : /* First block of rows */
1795 : : loc3 = matr[2] * f + dg[2] * cross;
1796 : : loc4 = dg[2] * g + matr[2] * cross;
1797 : :
1798 : : h_obj[0][0] = 1.5 * ( loc0 + loc3 * matr[2] + loc4 * dg[2] );
1799 : : h_obj[0][1] = 1.5 * ( loc3 * matr[5] + loc4 * dg[5] );
1800 : : h_obj[0][2] = 1.5 * ( loc3 * matr[8] + loc4 * dg[8] );
1801 : :
1802 : : /* Second block of rows */
1803 : : loc3 = matr[5] * f + dg[5] * cross;
1804 : : loc4 = dg[5] * g + matr[5] * cross;
1805 : :
1806 : : h_obj[1][1] = 1.5 * ( loc0 + loc3 * matr[5] + loc4 * dg[5] );
1807 : : h_obj[1][2] = 1.5 * ( loc3 * matr[8] + loc4 * dg[8] );
1808 : :
1809 : : /* Third block of rows */
1810 : : loc3 = matr[8] * f + dg[8] * cross;
1811 : : loc4 = dg[8] * g + matr[8] * cross;
1812 : :
1813 : : h_obj[2][2] = 1.5 * ( loc0 + loc3 * matr[8] + loc4 * dg[8] );
1814 : :
1815 : : // completes diagonal block.
1816 : : h_obj.fill_lower_triangle();
1817 : : return true;
1818 : : }
1819 : :
1820 : : inline bool h_fcn_3e_v0( double& obj, Vector3D& g_obj, Matrix3D& h_obj, const Vector3D x[4], const double a,
1821 : : const Exponent& b, const Exponent& c )
1822 : : {
1823 : : static Vector3D my_x[4];
1824 : :
1825 : : my_x[0] = x[1];
1826 : : my_x[1] = x[3];
1827 : : my_x[2] = x[2];
1828 : : my_x[3] = x[0];
1829 : : return h_fcn_3e_v3( obj, g_obj, h_obj, my_x, a, b, c );
1830 : : }
1831 : :
1832 : : inline bool h_fcn_3e_v1( double& obj, Vector3D& g_obj, Matrix3D& h_obj, const Vector3D x[4], const double a,
1833 : : const Exponent& b, const Exponent& c )
1834 : : {
1835 : : static Vector3D my_x[4];
1836 : :
1837 : : my_x[0] = x[0];
1838 : : my_x[1] = x[2];
1839 : : my_x[2] = x[3];
1840 : : my_x[3] = x[1];
1841 : : return h_fcn_3e_v3( obj, g_obj, h_obj, my_x, a, b, c );
1842 : : }
1843 : :
1844 : : inline bool h_fcn_3e_v2( double& obj, Vector3D& g_obj, Matrix3D& h_obj, const Vector3D x[4], const double a,
1845 : : const Exponent& b, const Exponent& c )
1846 : : {
1847 : : static Vector3D my_x[4];
1848 : :
1849 : : my_x[0] = x[1];
1850 : : my_x[1] = x[0];
1851 : : my_x[2] = x[3];
1852 : : my_x[3] = x[2];
1853 : : return h_fcn_3e_v3( obj, g_obj, h_obj, my_x, a, b, c );
1854 : : }
1855 : :
1856 : : /*****************************************************************************/
1857 : : /* The following set of functions reference tetrahedral elements to a */
1858 : : /* right tetrahedron. They are used when assessing the quality of a */
1859 : : /* hexahedral element. A zero return value indicates success, while */
1860 : : /* a nonzero value indicates failure. */
1861 : : /*****************************************************************************/
1862 : :
1863 : : /*****************************************************************************/
1864 : : /* Function evaluation requires 53 flops. */
1865 : : /* Reductions possible when b == 1, c == 1, or d == 1 */
1866 : : /*****************************************************************************/
1867 : 63219 : inline bool m_fcn_3i( double& obj, const Vector3D x[4], const double a, const Exponent& b, const Exponent& c,
1868 : : const Vector3D& d )
1869 : : {
1870 : : double matr[9], f;
1871 : : double g;
1872 : :
1873 : : /* Calculate M = A*inv(W). */
1874 [ + - ][ + - ]: 63219 : matr[0] = d[0] * ( x[1][0] - x[0][0] );
[ + - ]
1875 [ + - ][ + - ]: 63219 : matr[1] = d[1] * ( x[2][0] - x[0][0] );
[ + - ]
1876 [ + - ][ + - ]: 63219 : matr[2] = d[2] * ( x[3][0] - x[0][0] );
[ + - ]
1877 : :
1878 [ + - ][ + - ]: 63219 : matr[3] = d[0] * ( x[1][1] - x[0][1] );
[ + - ]
1879 [ + - ][ + - ]: 63219 : matr[4] = d[1] * ( x[2][1] - x[0][1] );
[ + - ]
1880 [ + - ][ + - ]: 63219 : matr[5] = d[2] * ( x[3][1] - x[0][1] );
[ + - ]
1881 : :
1882 [ + - ][ + - ]: 63219 : matr[6] = d[0] * ( x[1][2] - x[0][2] );
[ + - ]
1883 [ + - ][ + - ]: 63219 : matr[7] = d[1] * ( x[2][2] - x[0][2] );
[ + - ]
1884 [ + - ][ + - ]: 63219 : matr[8] = d[2] * ( x[3][2] - x[0][2] );
[ + - ]
1885 : :
1886 : : /* Calculate det(M). */
1887 : 63219 : g = matr[0] * ( matr[4] * matr[8] - matr[5] * matr[7] ) + matr[1] * ( matr[5] * matr[6] - matr[3] * matr[8] ) +
1888 : 63219 : matr[2] * ( matr[3] * matr[7] - matr[4] * matr[6] );
1889 [ + + ]: 63219 : if( g < MSQ_MIN )
1890 : : {
1891 : 1 : obj = g;
1892 : 1 : return false;
1893 : : }
1894 : :
1895 : : /* Calculate norm(M). */
1896 : 126436 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[2] * matr[2] + matr[3] * matr[3] + matr[4] * matr[4] +
1897 : 126436 : matr[5] * matr[5] + matr[6] * matr[6] + matr[7] * matr[7] + matr[8] * matr[8];
1898 : :
1899 : : /* Calculate objective function. */
1900 [ + - ][ + - ]: 63218 : obj = a * b.raise( f ) * c.raise( g );
1901 : 63219 : return true;
1902 : : }
1903 : :
1904 : : /*****************************************************************************/
1905 : : /* Gradient evaluation requires 115 flops. */
1906 : : /* Reductions possible when b == 1, c == 1, or d == 1 */
1907 : : /*****************************************************************************/
1908 : 108648 : inline bool g_fcn_3i( double& obj, Vector3D g_obj[4], const Vector3D x[4], const double a, const Exponent& b,
1909 : : const Exponent& c, const Vector3D& d )
1910 : : {
1911 : : double matr[9], f;
1912 : : double adj_m[9], g;
1913 : : double loc1, loc2, loc3;
1914 : :
1915 : : /* Calculate M = A*inv(W). */
1916 [ + - ][ + - ]: 108648 : matr[0] = d[0] * ( x[1][0] - x[0][0] );
[ + - ]
1917 [ + - ][ + - ]: 108648 : matr[1] = d[1] * ( x[2][0] - x[0][0] );
[ + - ]
1918 [ + - ][ + - ]: 108648 : matr[2] = d[2] * ( x[3][0] - x[0][0] );
[ + - ]
1919 : :
1920 [ + - ][ + - ]: 108648 : matr[3] = d[0] * ( x[1][1] - x[0][1] );
[ + - ]
1921 [ + - ][ + - ]: 108648 : matr[4] = d[1] * ( x[2][1] - x[0][1] );
[ + - ]
1922 [ + - ][ + - ]: 108648 : matr[5] = d[2] * ( x[3][1] - x[0][1] );
[ + - ]
1923 : :
1924 [ + - ][ + - ]: 108648 : matr[6] = d[0] * ( x[1][2] - x[0][2] );
[ + - ]
1925 [ + - ][ + - ]: 108648 : matr[7] = d[1] * ( x[2][2] - x[0][2] );
[ + - ]
1926 [ + - ][ + - ]: 108648 : matr[8] = d[2] * ( x[3][2] - x[0][2] );
[ + - ]
1927 : :
1928 : : /* Calculate det(M). */
1929 : 108648 : loc1 = matr[4] * matr[8] - matr[5] * matr[7];
1930 : 108648 : loc2 = matr[5] * matr[6] - matr[3] * matr[8];
1931 : 108648 : loc3 = matr[3] * matr[7] - matr[4] * matr[6];
1932 : 108648 : g = matr[0] * loc1 + matr[1] * loc2 + matr[2] * loc3;
1933 [ - + ]: 108648 : if( g < MSQ_MIN )
1934 : : {
1935 : 0 : obj = g;
1936 : 0 : return false;
1937 : : }
1938 : :
1939 : : /* Calculate norm(M). */
1940 : 217296 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[2] * matr[2] + matr[3] * matr[3] + matr[4] * matr[4] +
1941 : 217296 : matr[5] * matr[5] + matr[6] * matr[6] + matr[7] * matr[7] + matr[8] * matr[8];
1942 : :
1943 : : /* Calculate objective function. */
1944 [ + - ][ + - ]: 108648 : obj = a * b.raise( f ) * c.raise( g );
1945 : :
1946 : : /* Calculate the derivative of the objective function. */
1947 [ + - ]: 108648 : f = b.value() * obj / f; /* Constant on nabla f */
1948 [ + - ]: 108648 : g = c.value() * obj / g; /* Constant on nable g */
1949 : 108648 : f *= 2.0; /* Modification for nabla f */
1950 : :
1951 [ + - ]: 108648 : adj_m[0] = d[0] * ( matr[0] * f + loc1 * g );
1952 [ + - ]: 108648 : adj_m[1] = d[1] * ( matr[1] * f + loc2 * g );
1953 [ + - ]: 108648 : adj_m[2] = d[2] * ( matr[2] * f + loc3 * g );
1954 : :
1955 : 108648 : loc1 = matr[0] * g;
1956 : 108648 : loc2 = matr[1] * g;
1957 : 108648 : loc3 = matr[2] * g;
1958 : :
1959 [ + - ]: 108648 : adj_m[3] = d[0] * ( matr[3] * f + loc3 * matr[7] - loc2 * matr[8] );
1960 [ + - ]: 108648 : adj_m[4] = d[1] * ( matr[4] * f + loc1 * matr[8] - loc3 * matr[6] );
1961 [ + - ]: 108648 : adj_m[5] = d[2] * ( matr[5] * f + loc2 * matr[6] - loc1 * matr[7] );
1962 : :
1963 [ + - ]: 108648 : adj_m[6] = d[0] * ( matr[6] * f + loc2 * matr[5] - loc3 * matr[4] );
1964 [ + - ]: 108648 : adj_m[7] = d[1] * ( matr[7] * f + loc3 * matr[3] - loc1 * matr[5] );
1965 [ + - ]: 108648 : adj_m[8] = d[2] * ( matr[8] * f + loc1 * matr[4] - loc2 * matr[3] );
1966 : :
1967 [ + - ]: 108648 : g_obj[0][0] = -adj_m[0] - adj_m[1] - adj_m[2];
1968 [ + - ]: 108648 : g_obj[1][0] = adj_m[0];
1969 [ + - ]: 108648 : g_obj[2][0] = adj_m[1];
1970 [ + - ]: 108648 : g_obj[3][0] = adj_m[2];
1971 : :
1972 [ + - ]: 108648 : g_obj[0][1] = -adj_m[3] - adj_m[4] - adj_m[5];
1973 [ + - ]: 108648 : g_obj[1][1] = adj_m[3];
1974 [ + - ]: 108648 : g_obj[2][1] = adj_m[4];
1975 [ + - ]: 108648 : g_obj[3][1] = adj_m[5];
1976 : :
1977 [ + - ]: 108648 : g_obj[0][2] = -adj_m[6] - adj_m[7] - adj_m[8];
1978 [ + - ]: 108648 : g_obj[1][2] = adj_m[6];
1979 [ + - ]: 108648 : g_obj[2][2] = adj_m[7];
1980 [ + - ]: 108648 : g_obj[3][2] = adj_m[8];
1981 : 108648 : return true;
1982 : : }
1983 : :
1984 : : /*****************************************************************************/
1985 : : /* Hessian evaluation requires 469 flops. */
1986 : : /* Reductions possible when b == 1, c == 1, or d == 1 */
1987 : : /*****************************************************************************/
1988 : 98840 : inline int h_fcn_3i( double& obj, Vector3D g_obj[4], Matrix3D h_obj[10], const Vector3D x[4], const double a,
1989 : : const Exponent& b, const Exponent& c, const Vector3D& d )
1990 : : {
1991 : : double matr[9], f;
1992 : : double adj_m[9], g;
1993 : : double dg[9], loc0, loc1, loc2, loc3, loc4;
1994 : : double A[3], J_A[6], J_B[9], J_C[9], cross;
1995 : :
1996 [ + - ][ + - ]: 98840 : const double scale[6] = { d[0] * d[0], d[0] * d[1], d[0] * d[2], d[1] * d[1], d[1] * d[2], d[2] * d[2] };
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
1997 : :
1998 : : /* Calculate M = A*inv(W). */
1999 [ + - ][ + - ]: 98840 : matr[0] = d[0] * ( x[1][0] - x[0][0] );
[ + - ]
2000 [ + - ][ + - ]: 98840 : matr[1] = d[1] * ( x[2][0] - x[0][0] );
[ + - ]
2001 [ + - ][ + - ]: 98840 : matr[2] = d[2] * ( x[3][0] - x[0][0] );
[ + - ]
2002 : :
2003 [ + - ][ + - ]: 98840 : matr[3] = d[0] * ( x[1][1] - x[0][1] );
[ + - ]
2004 [ + - ][ + - ]: 98840 : matr[4] = d[1] * ( x[2][1] - x[0][1] );
[ + - ]
2005 [ + - ][ + - ]: 98840 : matr[5] = d[2] * ( x[3][1] - x[0][1] );
[ + - ]
2006 : :
2007 [ + - ][ + - ]: 98840 : matr[6] = d[0] * ( x[1][2] - x[0][2] );
[ + - ]
2008 [ + - ][ + - ]: 98840 : matr[7] = d[1] * ( x[2][2] - x[0][2] );
[ + - ]
2009 [ + - ][ + - ]: 98840 : matr[8] = d[2] * ( x[3][2] - x[0][2] );
[ + - ]
2010 : :
2011 : : /* Calculate det(M). */
2012 : 98840 : dg[0] = matr[4] * matr[8] - matr[5] * matr[7];
2013 : 98840 : dg[1] = matr[5] * matr[6] - matr[3] * matr[8];
2014 : 98840 : dg[2] = matr[3] * matr[7] - matr[4] * matr[6];
2015 : 98840 : g = matr[0] * dg[0] + matr[1] * dg[1] + matr[2] * dg[2];
2016 [ - + ]: 98840 : if( g < MSQ_MIN )
2017 : : {
2018 : 0 : obj = g;
2019 : 0 : return false;
2020 : : }
2021 : :
2022 : : /* Calculate norm(M). */
2023 : 197680 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[2] * matr[2] + matr[3] * matr[3] + matr[4] * matr[4] +
2024 : 197680 : matr[5] * matr[5] + matr[6] * matr[6] + matr[7] * matr[7] + matr[8] * matr[8];
2025 : :
2026 : 98840 : loc3 = f;
2027 : 98840 : loc4 = g;
2028 : :
2029 : : /* Calculate objective function. */
2030 [ + - ][ + - ]: 98840 : obj = a * b.raise( f ) * c.raise( g );
2031 : :
2032 : : /* Calculate the derivative of the objective function. */
2033 [ + - ]: 98840 : f = b.value() * obj / f; /* Constant on nabla f */
2034 [ + - ]: 98840 : g = c.value() * obj / g; /* Constant on nable g */
2035 : 98840 : f *= 2.0; /* Modification for nabla f */
2036 : :
2037 : 98840 : dg[3] = matr[2] * matr[7] - matr[1] * matr[8];
2038 : 98840 : dg[4] = matr[0] * matr[8] - matr[2] * matr[6];
2039 : 98840 : dg[5] = matr[1] * matr[6] - matr[0] * matr[7];
2040 : 98840 : dg[6] = matr[1] * matr[5] - matr[2] * matr[4];
2041 : 98840 : dg[7] = matr[2] * matr[3] - matr[0] * matr[5];
2042 : 98840 : dg[8] = matr[0] * matr[4] - matr[1] * matr[3];
2043 : :
2044 [ + - ]: 98840 : adj_m[0] = d[0] * ( matr[0] * f + dg[0] * g );
2045 [ + - ]: 98840 : adj_m[1] = d[1] * ( matr[1] * f + dg[1] * g );
2046 [ + - ]: 98840 : adj_m[2] = d[2] * ( matr[2] * f + dg[2] * g );
2047 [ + - ]: 98840 : adj_m[3] = d[0] * ( matr[3] * f + dg[3] * g );
2048 [ + - ]: 98840 : adj_m[4] = d[1] * ( matr[4] * f + dg[4] * g );
2049 [ + - ]: 98840 : adj_m[5] = d[2] * ( matr[5] * f + dg[5] * g );
2050 [ + - ]: 98840 : adj_m[6] = d[0] * ( matr[6] * f + dg[6] * g );
2051 [ + - ]: 98840 : adj_m[7] = d[1] * ( matr[7] * f + dg[7] * g );
2052 [ + - ]: 98840 : adj_m[8] = d[2] * ( matr[8] * f + dg[8] * g );
2053 : :
2054 [ + - ]: 98840 : g_obj[0][0] = -adj_m[0] - adj_m[1] - adj_m[2];
2055 [ + - ]: 98840 : g_obj[1][0] = adj_m[0];
2056 [ + - ]: 98840 : g_obj[2][0] = adj_m[1];
2057 [ + - ]: 98840 : g_obj[3][0] = adj_m[2];
2058 : :
2059 [ + - ]: 98840 : g_obj[0][1] = -adj_m[3] - adj_m[4] - adj_m[5];
2060 [ + - ]: 98840 : g_obj[1][1] = adj_m[3];
2061 [ + - ]: 98840 : g_obj[2][1] = adj_m[4];
2062 [ + - ]: 98840 : g_obj[3][1] = adj_m[5];
2063 : :
2064 [ + - ]: 98840 : g_obj[0][2] = -adj_m[6] - adj_m[7] - adj_m[8];
2065 [ + - ]: 98840 : g_obj[1][2] = adj_m[6];
2066 [ + - ]: 98840 : g_obj[2][2] = adj_m[7];
2067 [ + - ]: 98840 : g_obj[3][2] = adj_m[8];
2068 : :
2069 : : /* Calculate the hessian of the objective. */
2070 : 98840 : loc0 = f; /* Constant on nabla^2 f */
2071 : 98840 : loc1 = g; /* Constant on nabla^2 g */
2072 [ + - ]: 98840 : cross = f * c.value() / loc4; /* Constant on nabla g nabla f */
2073 [ + - ]: 98840 : f = f * ( b.value() - 1 ) / loc3; /* Constant on nabla f nabla f */
2074 [ + - ]: 98840 : g = g * ( c.value() - 1 ) / loc4; /* Constant on nabla g nabla g */
2075 : 98840 : f *= 2.0; /* Modification for nabla f */
2076 : :
2077 : : /* First block of rows */
2078 : 98840 : loc3 = matr[0] * f + dg[0] * cross;
2079 : 98840 : loc4 = dg[0] * g + matr[0] * cross;
2080 : :
2081 : 98840 : J_A[0] = loc0 + loc3 * matr[0] + loc4 * dg[0];
2082 : 98840 : J_A[1] = loc3 * matr[1] + loc4 * dg[1];
2083 : 98840 : J_A[2] = loc3 * matr[2] + loc4 * dg[2];
2084 : 98840 : J_B[0] = loc3 * matr[3] + loc4 * dg[3];
2085 : 98840 : J_B[1] = loc3 * matr[4] + loc4 * dg[4];
2086 : 98840 : J_B[2] = loc3 * matr[5] + loc4 * dg[5];
2087 : 98840 : J_C[0] = loc3 * matr[6] + loc4 * dg[6];
2088 : 98840 : J_C[1] = loc3 * matr[7] + loc4 * dg[7];
2089 : 98840 : J_C[2] = loc3 * matr[8] + loc4 * dg[8];
2090 : :
2091 : 98840 : loc3 = matr[1] * f + dg[1] * cross;
2092 : 98840 : loc4 = dg[1] * g + matr[1] * cross;
2093 : :
2094 : 98840 : J_A[3] = loc0 + loc3 * matr[1] + loc4 * dg[1];
2095 : 98840 : J_A[4] = loc3 * matr[2] + loc4 * dg[2];
2096 : 98840 : J_B[3] = loc3 * matr[3] + loc4 * dg[3];
2097 : 98840 : J_B[4] = loc3 * matr[4] + loc4 * dg[4];
2098 : 98840 : J_B[5] = loc3 * matr[5] + loc4 * dg[5];
2099 : 98840 : J_C[3] = loc3 * matr[6] + loc4 * dg[6];
2100 : 98840 : J_C[4] = loc3 * matr[7] + loc4 * dg[7];
2101 : 98840 : J_C[5] = loc3 * matr[8] + loc4 * dg[8];
2102 : :
2103 : 98840 : loc3 = matr[2] * f + dg[2] * cross;
2104 : 98840 : loc4 = dg[2] * g + matr[2] * cross;
2105 : :
2106 : 98840 : J_A[5] = loc0 + loc3 * matr[2] + loc4 * dg[2];
2107 : 98840 : J_B[6] = loc3 * matr[3] + loc4 * dg[3];
2108 : 98840 : J_B[7] = loc3 * matr[4] + loc4 * dg[4];
2109 : 98840 : J_B[8] = loc3 * matr[5] + loc4 * dg[5];
2110 : 98840 : J_C[6] = loc3 * matr[6] + loc4 * dg[6];
2111 : 98840 : J_C[7] = loc3 * matr[7] + loc4 * dg[7];
2112 : 98840 : J_C[8] = loc3 * matr[8] + loc4 * dg[8];
2113 : :
2114 : : /* First diagonal block */
2115 : 98840 : J_A[0] *= scale[0];
2116 : 98840 : J_A[1] *= scale[1];
2117 : 98840 : J_A[2] *= scale[2];
2118 : 98840 : J_A[3] *= scale[3];
2119 : 98840 : J_A[4] *= scale[4];
2120 : 98840 : J_A[5] *= scale[5];
2121 : :
2122 : 98840 : A[0] = -J_A[0] - J_A[1] - J_A[2];
2123 : 98840 : A[1] = -J_A[1] - J_A[3] - J_A[4];
2124 : 98840 : A[2] = -J_A[2] - J_A[4] - J_A[5];
2125 : :
2126 [ + - ]: 98840 : h_obj[0][0][0] = -A[0] - A[1] - A[2];
2127 [ + - ]: 98840 : h_obj[1][0][0] = A[0];
2128 [ + - ]: 98840 : h_obj[2][0][0] = A[1];
2129 [ + - ]: 98840 : h_obj[3][0][0] = A[2];
2130 : :
2131 [ + - ]: 98840 : h_obj[4][0][0] = J_A[0];
2132 [ + - ]: 98840 : h_obj[5][0][0] = J_A[1];
2133 [ + - ]: 98840 : h_obj[6][0][0] = J_A[2];
2134 : :
2135 [ + - ]: 98840 : h_obj[7][0][0] = J_A[3];
2136 [ + - ]: 98840 : h_obj[8][0][0] = J_A[4];
2137 : :
2138 [ + - ]: 98840 : h_obj[9][0][0] = J_A[5];
2139 : :
2140 : : /* First off-diagonal block */
2141 : 98840 : loc2 = matr[8] * loc1;
2142 : 98840 : J_B[1] += loc2;
2143 : 98840 : J_B[3] -= loc2;
2144 : :
2145 : 98840 : loc2 = matr[7] * loc1;
2146 : 98840 : J_B[2] -= loc2;
2147 : 98840 : J_B[6] += loc2;
2148 : :
2149 : 98840 : loc2 = matr[6] * loc1;
2150 : 98840 : J_B[5] += loc2;
2151 : 98840 : J_B[7] -= loc2;
2152 : :
2153 : 98840 : J_B[0] *= scale[0];
2154 : 98840 : J_B[1] *= scale[1];
2155 : 98840 : J_B[2] *= scale[2];
2156 : 98840 : J_B[3] *= scale[1];
2157 : 98840 : J_B[4] *= scale[3];
2158 : 98840 : J_B[5] *= scale[4];
2159 : 98840 : J_B[6] *= scale[2];
2160 : 98840 : J_B[7] *= scale[4];
2161 : 98840 : J_B[8] *= scale[5];
2162 : :
2163 : 98840 : A[0] = -J_B[0] - J_B[3] - J_B[6];
2164 : 98840 : A[1] = -J_B[1] - J_B[4] - J_B[7];
2165 : 98840 : A[2] = -J_B[2] - J_B[5] - J_B[8];
2166 : :
2167 [ + - ]: 98840 : h_obj[0][0][1] = -A[0] - A[1] - A[2];
2168 [ + - ]: 98840 : h_obj[1][0][1] = A[0];
2169 [ + - ]: 98840 : h_obj[2][0][1] = A[1];
2170 [ + - ]: 98840 : h_obj[3][0][1] = A[2];
2171 : :
2172 [ + - ]: 98840 : h_obj[1][1][0] = -J_B[0] - J_B[1] - J_B[2];
2173 [ + - ]: 98840 : h_obj[4][0][1] = J_B[0];
2174 [ + - ]: 98840 : h_obj[5][0][1] = J_B[1];
2175 [ + - ]: 98840 : h_obj[6][0][1] = J_B[2];
2176 : :
2177 [ + - ]: 98840 : h_obj[2][1][0] = -J_B[3] - J_B[4] - J_B[5];
2178 [ + - ]: 98840 : h_obj[5][1][0] = J_B[3];
2179 [ + - ]: 98840 : h_obj[7][0][1] = J_B[4];
2180 [ + - ]: 98840 : h_obj[8][0][1] = J_B[5];
2181 : :
2182 [ + - ]: 98840 : h_obj[3][1][0] = -J_B[6] - J_B[7] - J_B[8];
2183 [ + - ]: 98840 : h_obj[6][1][0] = J_B[6];
2184 [ + - ]: 98840 : h_obj[8][1][0] = J_B[7];
2185 [ + - ]: 98840 : h_obj[9][0][1] = J_B[8];
2186 : :
2187 : : /* Second off-diagonal block */
2188 : 98840 : loc2 = matr[5] * loc1;
2189 : 98840 : J_C[1] -= loc2;
2190 : 98840 : J_C[3] += loc2;
2191 : :
2192 : 98840 : loc2 = matr[4] * loc1;
2193 : 98840 : J_C[2] += loc2;
2194 : 98840 : J_C[6] -= loc2;
2195 : :
2196 : 98840 : loc2 = matr[3] * loc1;
2197 : 98840 : J_C[5] -= loc2;
2198 : 98840 : J_C[7] += loc2;
2199 : :
2200 : 98840 : J_C[0] *= scale[0];
2201 : 98840 : J_C[1] *= scale[1];
2202 : 98840 : J_C[2] *= scale[2];
2203 : 98840 : J_C[3] *= scale[1];
2204 : 98840 : J_C[4] *= scale[3];
2205 : 98840 : J_C[5] *= scale[4];
2206 : 98840 : J_C[6] *= scale[2];
2207 : 98840 : J_C[7] *= scale[4];
2208 : 98840 : J_C[8] *= scale[5];
2209 : :
2210 : 98840 : A[0] = -J_C[0] - J_C[3] - J_C[6];
2211 : 98840 : A[1] = -J_C[1] - J_C[4] - J_C[7];
2212 : 98840 : A[2] = -J_C[2] - J_C[5] - J_C[8];
2213 : :
2214 [ + - ]: 98840 : h_obj[0][0][2] = -A[0] - A[1] - A[2];
2215 [ + - ]: 98840 : h_obj[1][0][2] = A[0];
2216 [ + - ]: 98840 : h_obj[2][0][2] = A[1];
2217 [ + - ]: 98840 : h_obj[3][0][2] = A[2];
2218 : :
2219 [ + - ]: 98840 : h_obj[1][2][0] = -J_C[0] - J_C[1] - J_C[2];
2220 [ + - ]: 98840 : h_obj[4][0][2] = J_C[0];
2221 [ + - ]: 98840 : h_obj[5][0][2] = J_C[1];
2222 [ + - ]: 98840 : h_obj[6][0][2] = J_C[2];
2223 : :
2224 [ + - ]: 98840 : h_obj[2][2][0] = -J_C[3] - J_C[4] - J_C[5];
2225 [ + - ]: 98840 : h_obj[5][2][0] = J_C[3];
2226 [ + - ]: 98840 : h_obj[7][0][2] = J_C[4];
2227 [ + - ]: 98840 : h_obj[8][0][2] = J_C[5];
2228 : :
2229 [ + - ]: 98840 : h_obj[3][2][0] = -J_C[6] - J_C[7] - J_C[8];
2230 [ + - ]: 98840 : h_obj[6][2][0] = J_C[6];
2231 [ + - ]: 98840 : h_obj[8][2][0] = J_C[7];
2232 [ + - ]: 98840 : h_obj[9][0][2] = J_C[8];
2233 : :
2234 : : /* Second block of rows */
2235 : 98840 : loc3 = matr[3] * f + dg[3] * cross;
2236 : 98840 : loc4 = dg[3] * g + matr[3] * cross;
2237 : :
2238 : 98840 : J_A[0] = loc0 + loc3 * matr[3] + loc4 * dg[3];
2239 : 98840 : J_A[1] = loc3 * matr[4] + loc4 * dg[4];
2240 : 98840 : J_A[2] = loc3 * matr[5] + loc4 * dg[5];
2241 : 98840 : J_B[0] = loc3 * matr[6] + loc4 * dg[6];
2242 : 98840 : J_B[1] = loc3 * matr[7] + loc4 * dg[7];
2243 : 98840 : J_B[2] = loc3 * matr[8] + loc4 * dg[8];
2244 : :
2245 : 98840 : loc3 = matr[4] * f + dg[4] * cross;
2246 : 98840 : loc4 = dg[4] * g + matr[4] * cross;
2247 : :
2248 : 98840 : J_A[3] = loc0 + loc3 * matr[4] + loc4 * dg[4];
2249 : 98840 : J_A[4] = loc3 * matr[5] + loc4 * dg[5];
2250 : 98840 : J_B[3] = loc3 * matr[6] + loc4 * dg[6];
2251 : 98840 : J_B[4] = loc3 * matr[7] + loc4 * dg[7];
2252 : 98840 : J_B[5] = loc3 * matr[8] + loc4 * dg[8];
2253 : :
2254 : 98840 : loc3 = matr[5] * f + dg[5] * cross;
2255 : 98840 : loc4 = dg[5] * g + matr[5] * cross;
2256 : :
2257 : 98840 : J_A[5] = loc0 + loc3 * matr[5] + loc4 * dg[5];
2258 : 98840 : J_B[6] = loc3 * matr[6] + loc4 * dg[6];
2259 : 98840 : J_B[7] = loc3 * matr[7] + loc4 * dg[7];
2260 : 98840 : J_B[8] = loc3 * matr[8] + loc4 * dg[8];
2261 : :
2262 : : /* Second diagonal block */
2263 : 98840 : J_A[0] *= scale[0];
2264 : 98840 : J_A[1] *= scale[1];
2265 : 98840 : J_A[2] *= scale[2];
2266 : 98840 : J_A[3] *= scale[3];
2267 : 98840 : J_A[4] *= scale[4];
2268 : 98840 : J_A[5] *= scale[5];
2269 : :
2270 : 98840 : A[0] = -J_A[0] - J_A[1] - J_A[2];
2271 : 98840 : A[1] = -J_A[1] - J_A[3] - J_A[4];
2272 : 98840 : A[2] = -J_A[2] - J_A[4] - J_A[5];
2273 : :
2274 [ + - ]: 98840 : h_obj[0][1][1] = -A[0] - A[1] - A[2];
2275 [ + - ]: 98840 : h_obj[1][1][1] = A[0];
2276 [ + - ]: 98840 : h_obj[2][1][1] = A[1];
2277 [ + - ]: 98840 : h_obj[3][1][1] = A[2];
2278 : :
2279 [ + - ]: 98840 : h_obj[4][1][1] = J_A[0];
2280 [ + - ]: 98840 : h_obj[5][1][1] = J_A[1];
2281 [ + - ]: 98840 : h_obj[6][1][1] = J_A[2];
2282 : :
2283 [ + - ]: 98840 : h_obj[7][1][1] = J_A[3];
2284 [ + - ]: 98840 : h_obj[8][1][1] = J_A[4];
2285 : :
2286 [ + - ]: 98840 : h_obj[9][1][1] = J_A[5];
2287 : :
2288 : : /* Third off-diagonal block */
2289 : 98840 : loc2 = matr[2] * loc1;
2290 : 98840 : J_B[1] += loc2;
2291 : 98840 : J_B[3] -= loc2;
2292 : :
2293 : 98840 : loc2 = matr[1] * loc1;
2294 : 98840 : J_B[2] -= loc2;
2295 : 98840 : J_B[6] += loc2;
2296 : :
2297 : 98840 : loc2 = matr[0] * loc1;
2298 : 98840 : J_B[5] += loc2;
2299 : 98840 : J_B[7] -= loc2;
2300 : :
2301 : 98840 : J_B[0] *= scale[0];
2302 : 98840 : J_B[1] *= scale[1];
2303 : 98840 : J_B[2] *= scale[2];
2304 : 98840 : J_B[3] *= scale[1];
2305 : 98840 : J_B[4] *= scale[3];
2306 : 98840 : J_B[5] *= scale[4];
2307 : 98840 : J_B[6] *= scale[2];
2308 : 98840 : J_B[7] *= scale[4];
2309 : 98840 : J_B[8] *= scale[5];
2310 : :
2311 : 98840 : A[0] = -J_B[0] - J_B[3] - J_B[6];
2312 : 98840 : A[1] = -J_B[1] - J_B[4] - J_B[7];
2313 : 98840 : A[2] = -J_B[2] - J_B[5] - J_B[8];
2314 : :
2315 [ + - ]: 98840 : h_obj[0][1][2] = -A[0] - A[1] - A[2];
2316 [ + - ]: 98840 : h_obj[1][1][2] = A[0];
2317 [ + - ]: 98840 : h_obj[2][1][2] = A[1];
2318 [ + - ]: 98840 : h_obj[3][1][2] = A[2];
2319 : :
2320 [ + - ]: 98840 : h_obj[1][2][1] = -J_B[0] - J_B[1] - J_B[2];
2321 [ + - ]: 98840 : h_obj[4][1][2] = J_B[0];
2322 [ + - ]: 98840 : h_obj[5][1][2] = J_B[1];
2323 [ + - ]: 98840 : h_obj[6][1][2] = J_B[2];
2324 : :
2325 [ + - ]: 98840 : h_obj[2][2][1] = -J_B[3] - J_B[4] - J_B[5];
2326 [ + - ]: 98840 : h_obj[5][2][1] = J_B[3];
2327 [ + - ]: 98840 : h_obj[7][1][2] = J_B[4];
2328 [ + - ]: 98840 : h_obj[8][1][2] = J_B[5];
2329 : :
2330 [ + - ]: 98840 : h_obj[3][2][1] = -J_B[6] - J_B[7] - J_B[8];
2331 [ + - ]: 98840 : h_obj[6][2][1] = J_B[6];
2332 [ + - ]: 98840 : h_obj[8][2][1] = J_B[7];
2333 [ + - ]: 98840 : h_obj[9][1][2] = J_B[8];
2334 : :
2335 : : /* Third block of rows */
2336 : 98840 : loc3 = matr[6] * f + dg[6] * cross;
2337 : 98840 : loc4 = dg[6] * g + matr[6] * cross;
2338 : :
2339 : 98840 : J_A[0] = loc0 + loc3 * matr[6] + loc4 * dg[6];
2340 : 98840 : J_A[1] = loc3 * matr[7] + loc4 * dg[7];
2341 : 98840 : J_A[2] = loc3 * matr[8] + loc4 * dg[8];
2342 : :
2343 : 98840 : loc3 = matr[7] * f + dg[7] * cross;
2344 : 98840 : loc4 = dg[7] * g + matr[7] * cross;
2345 : :
2346 : 98840 : J_A[3] = loc0 + loc3 * matr[7] + loc4 * dg[7];
2347 : 98840 : J_A[4] = loc3 * matr[8] + loc4 * dg[8];
2348 : :
2349 : 98840 : loc3 = matr[8] * f + dg[8] * cross;
2350 : 98840 : loc4 = dg[8] * g + matr[8] * cross;
2351 : :
2352 : 98840 : J_A[5] = loc0 + loc3 * matr[8] + loc4 * dg[8];
2353 : :
2354 : : /* Third diagonal block */
2355 : 98840 : J_A[0] *= scale[0];
2356 : 98840 : J_A[1] *= scale[1];
2357 : 98840 : J_A[2] *= scale[2];
2358 : 98840 : J_A[3] *= scale[3];
2359 : 98840 : J_A[4] *= scale[4];
2360 : 98840 : J_A[5] *= scale[5];
2361 : :
2362 : 98840 : A[0] = -J_A[0] - J_A[1] - J_A[2];
2363 : 98840 : A[1] = -J_A[1] - J_A[3] - J_A[4];
2364 : 98840 : A[2] = -J_A[2] - J_A[4] - J_A[5];
2365 : :
2366 [ + - ]: 98840 : h_obj[0][2][2] = -A[0] - A[1] - A[2];
2367 [ + - ]: 98840 : h_obj[1][2][2] = A[0];
2368 [ + - ]: 98840 : h_obj[2][2][2] = A[1];
2369 [ + - ]: 98840 : h_obj[3][2][2] = A[2];
2370 : :
2371 [ + - ]: 98840 : h_obj[4][2][2] = J_A[0];
2372 [ + - ]: 98840 : h_obj[5][2][2] = J_A[1];
2373 [ + - ]: 98840 : h_obj[6][2][2] = J_A[2];
2374 : :
2375 [ + - ]: 98840 : h_obj[7][2][2] = J_A[3];
2376 [ + - ]: 98840 : h_obj[8][2][2] = J_A[4];
2377 : :
2378 [ + - ]: 98840 : h_obj[9][2][2] = J_A[5];
2379 : :
2380 : : // completes diagonal blocks.
2381 [ + - ]: 98840 : h_obj[0].fill_lower_triangle();
2382 [ + - ]: 98840 : h_obj[4].fill_lower_triangle();
2383 [ + - ]: 98840 : h_obj[7].fill_lower_triangle();
2384 [ + - ]: 98840 : h_obj[9].fill_lower_triangle();
2385 : 98840 : return true;
2386 : : }
2387 : :
2388 : : /*********************************************************************
2389 : : * Reference tetrahedral elements to halves of an ideal pyramid.
2390 : : * Vertices should be ordered such that the edge between the 2nd and
2391 : : * 3rd vertices is the one longer edge in the reference tetrahedron
2392 : : * (the diagonal of the base of the pyramid of which the tetrahedron
2393 : : * is one half of).
2394 : : * 1 0 1/2 1 0 -1/sqrt(2)
2395 : : * W = 0 1 1/2 inv(W) = 0 1 -1/sqrt(2)
2396 : : * 0 0 1/sqrt(2) 0 0 sqrt(2)
2397 : : *
2398 : : *********************************************************************/
2399 : 4553 : inline bool m_fcn_3p( double& obj, const Vector3D x[4], const double a, const Exponent& b, const Exponent& c )
2400 : : {
2401 : 4553 : const double h = 0.5; /* h = 1 / (2*height) */
2402 : :
2403 : : double matr[9], f;
2404 : : double g;
2405 : :
2406 : : /* Calculate M = A*inv(W). */
2407 [ + - ][ + - ]: 4553 : matr[0] = x[1][0] - x[0][0];
2408 [ + - ][ + - ]: 4553 : matr[1] = x[2][0] - x[0][0];
2409 [ + - ][ + - ]: 4553 : matr[2] = ( 2.0 * x[3][0] - x[1][0] - x[2][0] ) * h;
[ + - ]
2410 : :
2411 [ + - ][ + - ]: 4553 : matr[3] = x[1][1] - x[0][1];
2412 [ + - ][ + - ]: 4553 : matr[4] = x[2][1] - x[0][1];
2413 [ + - ][ + - ]: 4553 : matr[5] = ( 2.0 * x[3][1] - x[1][1] - x[2][1] ) * h;
[ + - ]
2414 : :
2415 [ + - ][ + - ]: 4553 : matr[6] = x[1][2] - x[0][2];
2416 [ + - ][ + - ]: 4553 : matr[7] = x[2][2] - x[0][2];
2417 [ + - ][ + - ]: 4553 : matr[8] = ( 2.0 * x[3][2] - x[1][2] - x[2][2] ) * h;
[ + - ]
2418 : :
2419 : : /* Calculate det(M). */
2420 : 4553 : g = matr[0] * ( matr[4] * matr[8] - matr[5] * matr[7] ) + matr[1] * ( matr[5] * matr[6] - matr[3] * matr[8] ) +
2421 : 4553 : matr[2] * ( matr[3] * matr[7] - matr[4] * matr[6] );
2422 [ + + ]: 4553 : if( g < MSQ_MIN )
2423 : : {
2424 : 3 : obj = g;
2425 : 3 : return false;
2426 : : }
2427 : :
2428 : : /* Calculate norm(M). */
2429 : 9100 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[2] * matr[2] + matr[3] * matr[3] + matr[4] * matr[4] +
2430 : 9100 : matr[5] * matr[5] + matr[6] * matr[6] + matr[7] * matr[7] + matr[8] * matr[8];
2431 : :
2432 : : /* Calculate objective function. */
2433 [ + - ][ + - ]: 4550 : obj = a * b.raise( f ) * c.raise( g );
2434 : 4553 : return true;
2435 : : }
2436 : :
2437 : 9490 : inline bool g_fcn_3p( double& obj, Vector3D g_obj[4], const Vector3D x[4], const double a, const Exponent& b,
2438 : : const Exponent& c )
2439 : : {
2440 : 9490 : const double h = 0.5; /* h = 1 / (2*height) */
2441 : 9490 : const double th = 1.0; /* h = 1 / (height) */
2442 : :
2443 : : double matr[9], f;
2444 : : double adj_m[9], g;
2445 : : double loc1, loc2, loc3;
2446 : :
2447 : : /* Calculate M = A*inv(W). */
2448 [ + - ][ + - ]: 9490 : matr[0] = x[1][0] - x[0][0];
2449 [ + - ][ + - ]: 9490 : matr[1] = x[2][0] - x[0][0];
2450 [ + - ][ + - ]: 9490 : matr[2] = ( 2.0 * x[3][0] - x[1][0] - x[2][0] ) * h;
[ + - ]
2451 : :
2452 [ + - ][ + - ]: 9490 : matr[3] = x[1][1] - x[0][1];
2453 [ + - ][ + - ]: 9490 : matr[4] = x[2][1] - x[0][1];
2454 [ + - ][ + - ]: 9490 : matr[5] = ( 2.0 * x[3][1] - x[1][1] - x[2][1] ) * h;
[ + - ]
2455 : :
2456 [ + - ][ + - ]: 9490 : matr[6] = x[1][2] - x[0][2];
2457 [ + - ][ + - ]: 9490 : matr[7] = x[2][2] - x[0][2];
2458 [ + - ][ + - ]: 9490 : matr[8] = ( 2.0 * x[3][2] - x[1][2] - x[2][2] ) * h;
[ + - ]
2459 : :
2460 : : /* Calculate det(M). */
2461 : 9490 : loc1 = matr[4] * matr[8] - matr[5] * matr[7];
2462 : 9490 : loc2 = matr[5] * matr[6] - matr[3] * matr[8];
2463 : 9490 : loc3 = matr[3] * matr[7] - matr[4] * matr[6];
2464 : 9490 : g = matr[0] * loc1 + matr[1] * loc2 + matr[2] * loc3;
2465 [ + + ]: 9490 : if( g < MSQ_MIN )
2466 : : {
2467 : 2 : obj = g;
2468 : 2 : return false;
2469 : : }
2470 : :
2471 : : /* Calculate norm(M). */
2472 : 18976 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[2] * matr[2] + matr[3] * matr[3] + matr[4] * matr[4] +
2473 : 18976 : matr[5] * matr[5] + matr[6] * matr[6] + matr[7] * matr[7] + matr[8] * matr[8];
2474 : :
2475 : : /* Calculate objective function. */
2476 [ + - ][ + - ]: 9488 : obj = a * b.raise( f ) * c.raise( g );
2477 : :
2478 : : /* Calculate the derivative of the objective function. */
2479 [ + - ]: 9488 : f = b.value() * obj / f; /* Constant on nabla f */
2480 [ + - ]: 9488 : g = c.value() * obj / g; /* Constant on nable g */
2481 : 9488 : f *= 2.0; /* Modification for nabla f */
2482 : :
2483 : 9488 : adj_m[0] = matr[0] * f + loc1 * g;
2484 : 9488 : adj_m[1] = matr[1] * f + loc2 * g;
2485 : 9488 : adj_m[2] = matr[2] * f + loc3 * g;
2486 : :
2487 : 9488 : loc1 = matr[0] * g;
2488 : 9488 : loc2 = matr[1] * g;
2489 : 9488 : loc3 = matr[2] * g;
2490 : :
2491 : 9488 : adj_m[3] = matr[3] * f + loc3 * matr[7] - loc2 * matr[8];
2492 : 9488 : adj_m[4] = matr[4] * f + loc1 * matr[8] - loc3 * matr[6];
2493 : 9488 : adj_m[5] = matr[5] * f + loc2 * matr[6] - loc1 * matr[7];
2494 : :
2495 : 9488 : adj_m[6] = matr[6] * f + loc2 * matr[5] - loc3 * matr[4];
2496 : 9488 : adj_m[7] = matr[7] * f + loc3 * matr[3] - loc1 * matr[5];
2497 : 9488 : adj_m[8] = matr[8] * f + loc1 * matr[4] - loc2 * matr[3];
2498 : :
2499 [ + - ]: 9488 : g_obj[0][0] = -adj_m[0] - adj_m[1];
2500 [ + - ]: 9488 : g_obj[1][0] = adj_m[0] - h * adj_m[2];
2501 [ + - ]: 9488 : g_obj[2][0] = adj_m[1] - h * adj_m[2];
2502 [ + - ]: 9488 : g_obj[3][0] = th * adj_m[2];
2503 : :
2504 [ + - ]: 9488 : g_obj[0][1] = -adj_m[3] - adj_m[4];
2505 [ + - ]: 9488 : g_obj[1][1] = adj_m[3] - h * adj_m[5];
2506 [ + - ]: 9488 : g_obj[2][1] = adj_m[4] - h * adj_m[5];
2507 [ + - ]: 9488 : g_obj[3][1] = th * adj_m[5];
2508 : :
2509 [ + - ]: 9488 : g_obj[0][2] = -adj_m[6] - adj_m[7];
2510 [ + - ]: 9488 : g_obj[1][2] = adj_m[6] - h * adj_m[8];
2511 [ + - ]: 9488 : g_obj[2][2] = adj_m[7] - h * adj_m[8];
2512 [ + - ]: 9488 : g_obj[3][2] = th * adj_m[8];
2513 : 9490 : return true;
2514 : : }
2515 : :
2516 : 9224 : inline bool h_fcn_3p( double& obj, Vector3D g_obj[4], Matrix3D h_obj[10], const Vector3D x[4], const double a,
2517 : : const Exponent& b, const Exponent& c )
2518 : : {
2519 : 9224 : const double h = 0.5; /* h = 1 / (2*height) */
2520 : 9224 : const double th = 1.0; /* h = 1 / (height) */
2521 : :
2522 : : double matr[9], f;
2523 : : double adj_m[9], g;
2524 : : double dg[9], loc0, loc1, loc2, loc3, loc4;
2525 : : double A[12], J_A[6], J_B[9], J_C[9], cross;
2526 : :
2527 : : /* Calculate M = A*inv(W). */
2528 [ + - ][ + - ]: 9224 : matr[0] = x[1][0] - x[0][0];
2529 [ + - ][ + - ]: 9224 : matr[1] = x[2][0] - x[0][0];
2530 [ + - ][ + - ]: 9224 : matr[2] = ( 2.0 * x[3][0] - x[1][0] - x[2][0] ) * h;
[ + - ]
2531 : :
2532 [ + - ][ + - ]: 9224 : matr[3] = x[1][1] - x[0][1];
2533 [ + - ][ + - ]: 9224 : matr[4] = x[2][1] - x[0][1];
2534 [ + - ][ + - ]: 9224 : matr[5] = ( 2.0 * x[3][1] - x[1][1] - x[2][1] ) * h;
[ + - ]
2535 : :
2536 [ + - ][ + - ]: 9224 : matr[6] = x[1][2] - x[0][2];
2537 [ + - ][ + - ]: 9224 : matr[7] = x[2][2] - x[0][2];
2538 [ + - ][ + - ]: 9224 : matr[8] = ( 2.0 * x[3][2] - x[1][2] - x[2][2] ) * h;
[ + - ]
2539 : :
2540 : : /* Calculate det(M). */
2541 : 9224 : dg[0] = matr[4] * matr[8] - matr[5] * matr[7];
2542 : 9224 : dg[1] = matr[5] * matr[6] - matr[3] * matr[8];
2543 : 9224 : dg[2] = matr[3] * matr[7] - matr[4] * matr[6];
2544 : 9224 : g = matr[0] * dg[0] + matr[1] * dg[1] + matr[2] * dg[2];
2545 [ - + ]: 9224 : if( g < MSQ_MIN )
2546 : : {
2547 : 0 : obj = g;
2548 : 0 : return false;
2549 : : }
2550 : :
2551 : 9224 : dg[3] = matr[2] * matr[7] - matr[1] * matr[8];
2552 : 9224 : dg[4] = matr[0] * matr[8] - matr[2] * matr[6];
2553 : 9224 : dg[5] = matr[1] * matr[6] - matr[0] * matr[7];
2554 : 9224 : dg[6] = matr[1] * matr[5] - matr[2] * matr[4];
2555 : 9224 : dg[7] = matr[2] * matr[3] - matr[0] * matr[5];
2556 : 9224 : dg[8] = matr[0] * matr[4] - matr[1] * matr[3];
2557 : :
2558 : : /* Calculate norm(M). */
2559 : 18448 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[2] * matr[2] + matr[3] * matr[3] + matr[4] * matr[4] +
2560 : 18448 : matr[5] * matr[5] + matr[6] * matr[6] + matr[7] * matr[7] + matr[8] * matr[8];
2561 : :
2562 : 9224 : loc3 = f;
2563 : 9224 : loc4 = g;
2564 : :
2565 : : /* Calculate objective function. */
2566 [ + - ][ + - ]: 9224 : obj = a * b.raise( f ) * c.raise( g );
2567 : :
2568 : : /* Calculate the derivative of the objective function. */
2569 [ + - ]: 9224 : f = b.value() * obj / f; /* Constant on nabla f */
2570 [ + - ]: 9224 : g = c.value() * obj / g; /* Constant on nable g */
2571 : 9224 : f *= 2.0; /* Modification for nabla f */
2572 : :
2573 : 9224 : adj_m[0] = matr[0] * f + dg[0] * g;
2574 : 9224 : adj_m[1] = matr[1] * f + dg[1] * g;
2575 : 9224 : adj_m[2] = matr[2] * f + dg[2] * g;
2576 : 9224 : adj_m[3] = matr[3] * f + dg[3] * g;
2577 : 9224 : adj_m[4] = matr[4] * f + dg[4] * g;
2578 : 9224 : adj_m[5] = matr[5] * f + dg[5] * g;
2579 : 9224 : adj_m[6] = matr[6] * f + dg[6] * g;
2580 : 9224 : adj_m[7] = matr[7] * f + dg[7] * g;
2581 : 9224 : adj_m[8] = matr[8] * f + dg[8] * g;
2582 : :
2583 [ + - ]: 9224 : g_obj[0][0] = -adj_m[0] - adj_m[1];
2584 [ + - ]: 9224 : g_obj[1][0] = adj_m[0] - h * adj_m[2];
2585 [ + - ]: 9224 : g_obj[2][0] = adj_m[1] - h * adj_m[2];
2586 [ + - ]: 9224 : g_obj[3][0] = th * adj_m[2];
2587 : :
2588 [ + - ]: 9224 : g_obj[0][1] = -adj_m[3] - adj_m[4];
2589 [ + - ]: 9224 : g_obj[1][1] = adj_m[3] - h * adj_m[5];
2590 [ + - ]: 9224 : g_obj[2][1] = adj_m[4] - h * adj_m[5];
2591 [ + - ]: 9224 : g_obj[3][1] = th * adj_m[5];
2592 : :
2593 [ + - ]: 9224 : g_obj[0][2] = -adj_m[6] - adj_m[7];
2594 [ + - ]: 9224 : g_obj[1][2] = adj_m[6] - h * adj_m[8];
2595 [ + - ]: 9224 : g_obj[2][2] = adj_m[7] - h * adj_m[8];
2596 [ + - ]: 9224 : g_obj[3][2] = th * adj_m[8];
2597 : :
2598 : : /* Calculate the hessian of the objective. */
2599 : 9224 : loc0 = f; /* Constant on nabla^2 f */
2600 : 9224 : loc1 = g; /* Constant on nabla^2 g */
2601 [ + - ]: 9224 : cross = f * c.value() / loc4; /* Constant on nabla g nabla f */
2602 [ + - ]: 9224 : f = f * ( b.value() - 1 ) / loc3; /* Constant on nabla f nabla f */
2603 [ + - ]: 9224 : g = g * ( c.value() - 1 ) / loc4; /* Constant on nabla g nabla g */
2604 : 9224 : f *= 2.0; /* Modification for nabla f */
2605 : :
2606 : : /* First block of rows */
2607 : 9224 : loc3 = matr[0] * f + dg[0] * cross;
2608 : 9224 : loc4 = dg[0] * g + matr[0] * cross;
2609 : :
2610 : 9224 : J_A[0] = loc0 + loc3 * matr[0] + loc4 * dg[0];
2611 : 9224 : J_A[1] = loc3 * matr[1] + loc4 * dg[1];
2612 : 9224 : J_A[2] = loc3 * matr[2] + loc4 * dg[2];
2613 : 9224 : J_B[0] = loc3 * matr[3] + loc4 * dg[3];
2614 : 9224 : J_B[1] = loc3 * matr[4] + loc4 * dg[4];
2615 : 9224 : J_B[2] = loc3 * matr[5] + loc4 * dg[5];
2616 : 9224 : J_C[0] = loc3 * matr[6] + loc4 * dg[6];
2617 : 9224 : J_C[1] = loc3 * matr[7] + loc4 * dg[7];
2618 : 9224 : J_C[2] = loc3 * matr[8] + loc4 * dg[8];
2619 : :
2620 : 9224 : loc3 = matr[1] * f + dg[1] * cross;
2621 : 9224 : loc4 = dg[1] * g + matr[1] * cross;
2622 : :
2623 : 9224 : J_A[3] = loc0 + loc3 * matr[1] + loc4 * dg[1];
2624 : 9224 : J_A[4] = loc3 * matr[2] + loc4 * dg[2];
2625 : 9224 : J_B[3] = loc3 * matr[3] + loc4 * dg[3];
2626 : 9224 : J_B[4] = loc3 * matr[4] + loc4 * dg[4];
2627 : 9224 : J_B[5] = loc3 * matr[5] + loc4 * dg[5];
2628 : 9224 : J_C[3] = loc3 * matr[6] + loc4 * dg[6];
2629 : 9224 : J_C[4] = loc3 * matr[7] + loc4 * dg[7];
2630 : 9224 : J_C[5] = loc3 * matr[8] + loc4 * dg[8];
2631 : :
2632 : 9224 : loc3 = matr[2] * f + dg[2] * cross;
2633 : 9224 : loc4 = dg[2] * g + matr[2] * cross;
2634 : :
2635 : 9224 : J_A[5] = loc0 + loc3 * matr[2] + loc4 * dg[2];
2636 : 9224 : J_B[6] = loc3 * matr[3] + loc4 * dg[3];
2637 : 9224 : J_B[7] = loc3 * matr[4] + loc4 * dg[4];
2638 : 9224 : J_B[8] = loc3 * matr[5] + loc4 * dg[5];
2639 : 9224 : J_C[6] = loc3 * matr[6] + loc4 * dg[6];
2640 : 9224 : J_C[7] = loc3 * matr[7] + loc4 * dg[7];
2641 : 9224 : J_C[8] = loc3 * matr[8] + loc4 * dg[8];
2642 : :
2643 : : /* First diagonal block */
2644 : 9224 : A[0] = -J_A[0] - J_A[1];
2645 : 9224 : A[1] = J_A[0] - h * J_A[2];
2646 : 9224 : A[2] = J_A[1] - h * J_A[2];
2647 : 9224 : A[3] = th * J_A[2];
2648 : :
2649 : 9224 : A[4] = -J_A[1] - J_A[3];
2650 : 9224 : A[5] = J_A[1] - h * J_A[4];
2651 : 9224 : A[6] = J_A[3] - h * J_A[4];
2652 : 9224 : A[7] = th * J_A[4];
2653 : :
2654 : 9224 : A[8] = -J_A[2] - J_A[4];
2655 : 9224 : A[9] = J_A[2] - h * J_A[5];
2656 : 9224 : A[10] = J_A[4] - h * J_A[5];
2657 : 9224 : A[11] = th * J_A[5];
2658 : :
2659 [ + - ]: 9224 : h_obj[0][0][0] = -A[0] - A[4];
2660 [ + - ]: 9224 : h_obj[1][0][0] = A[0] - h * A[8];
2661 [ + - ]: 9224 : h_obj[2][0][0] = A[4] - h * A[8];
2662 [ + - ]: 9224 : h_obj[3][0][0] = th * A[8];
2663 : :
2664 [ + - ]: 9224 : h_obj[4][0][0] = A[1] - h * A[9];
2665 [ + - ]: 9224 : h_obj[5][0][0] = A[5] - h * A[9];
2666 [ + - ]: 9224 : h_obj[6][0][0] = th * A[9];
2667 : :
2668 [ + - ]: 9224 : h_obj[7][0][0] = A[6] - h * A[10];
2669 [ + - ]: 9224 : h_obj[8][0][0] = th * A[10];
2670 : :
2671 [ + - ]: 9224 : h_obj[9][0][0] = th * A[11];
2672 : :
2673 : : /* First off-diagonal block */
2674 : 9224 : loc2 = matr[8] * loc1;
2675 : 9224 : J_B[1] += loc2;
2676 : 9224 : J_B[3] -= loc2;
2677 : :
2678 : 9224 : loc2 = matr[7] * loc1;
2679 : 9224 : J_B[2] -= loc2;
2680 : 9224 : J_B[6] += loc2;
2681 : :
2682 : 9224 : loc2 = matr[6] * loc1;
2683 : 9224 : J_B[5] += loc2;
2684 : 9224 : J_B[7] -= loc2;
2685 : :
2686 : 9224 : A[0] = -J_B[0] - J_B[1];
2687 : 9224 : A[1] = J_B[0] - h * J_B[2];
2688 : 9224 : A[2] = J_B[1] - h * J_B[2];
2689 : 9224 : A[3] = th * J_B[2];
2690 : :
2691 : 9224 : A[4] = -J_B[3] - J_B[4];
2692 : 9224 : A[5] = J_B[3] - h * J_B[5];
2693 : 9224 : A[6] = J_B[4] - h * J_B[5];
2694 : 9224 : A[7] = th * J_B[5];
2695 : :
2696 : 9224 : A[8] = -J_B[6] - J_B[7];
2697 : 9224 : A[9] = J_B[6] - h * J_B[8];
2698 : 9224 : A[10] = J_B[7] - h * J_B[8];
2699 : 9224 : A[11] = th * J_B[8];
2700 : :
2701 [ + - ]: 9224 : h_obj[0][0][1] = -A[0] - A[4];
2702 [ + - ]: 9224 : h_obj[1][1][0] = A[0] - h * A[8];
2703 [ + - ]: 9224 : h_obj[2][1][0] = A[4] - h * A[8];
2704 [ + - ]: 9224 : h_obj[3][1][0] = th * A[8];
2705 : :
2706 [ + - ]: 9224 : h_obj[1][0][1] = -A[1] - A[5];
2707 [ + - ]: 9224 : h_obj[4][0][1] = A[1] - h * A[9];
2708 [ + - ]: 9224 : h_obj[5][1][0] = A[5] - h * A[9];
2709 [ + - ]: 9224 : h_obj[6][1][0] = th * A[9];
2710 : :
2711 [ + - ]: 9224 : h_obj[2][0][1] = -A[2] - A[6];
2712 [ + - ]: 9224 : h_obj[5][0][1] = A[2] - h * A[10];
2713 [ + - ]: 9224 : h_obj[7][0][1] = A[6] - h * A[10];
2714 [ + - ]: 9224 : h_obj[8][1][0] = th * A[10];
2715 : :
2716 [ + - ]: 9224 : h_obj[3][0][1] = -A[3] - A[7];
2717 [ + - ]: 9224 : h_obj[6][0][1] = A[3] - h * A[11];
2718 [ + - ]: 9224 : h_obj[8][0][1] = A[7] - h * A[11];
2719 [ + - ]: 9224 : h_obj[9][0][1] = th * A[11];
2720 : :
2721 : : /* Second off-diagonal block */
2722 : 9224 : loc2 = matr[5] * loc1;
2723 : 9224 : J_C[1] -= loc2;
2724 : 9224 : J_C[3] += loc2;
2725 : :
2726 : 9224 : loc2 = matr[4] * loc1;
2727 : 9224 : J_C[2] += loc2;
2728 : 9224 : J_C[6] -= loc2;
2729 : :
2730 : 9224 : loc2 = matr[3] * loc1;
2731 : 9224 : J_C[5] -= loc2;
2732 : 9224 : J_C[7] += loc2;
2733 : :
2734 : 9224 : A[0] = -J_C[0] - J_C[1];
2735 : 9224 : A[1] = J_C[0] - h * J_C[2];
2736 : 9224 : A[2] = J_C[1] - h * J_C[2];
2737 : 9224 : A[3] = th * J_C[2];
2738 : :
2739 : 9224 : A[4] = -J_C[3] - J_C[4];
2740 : 9224 : A[5] = J_C[3] - h * J_C[5];
2741 : 9224 : A[6] = J_C[4] - h * J_C[5];
2742 : 9224 : A[7] = th * J_C[5];
2743 : :
2744 : 9224 : A[8] = -J_C[6] - J_C[7];
2745 : 9224 : A[9] = J_C[6] - h * J_C[8];
2746 : 9224 : A[10] = J_C[7] - h * J_C[8];
2747 : 9224 : A[11] = th * J_C[8];
2748 : :
2749 [ + - ]: 9224 : h_obj[0][0][2] = -A[0] - A[4];
2750 [ + - ]: 9224 : h_obj[1][2][0] = A[0] - h * A[8];
2751 [ + - ]: 9224 : h_obj[2][2][0] = A[4] - h * A[8];
2752 [ + - ]: 9224 : h_obj[3][2][0] = th * A[8];
2753 : :
2754 [ + - ]: 9224 : h_obj[1][0][2] = -A[1] - A[5];
2755 [ + - ]: 9224 : h_obj[4][0][2] = A[1] - h * A[9];
2756 [ + - ]: 9224 : h_obj[5][2][0] = A[5] - h * A[9];
2757 [ + - ]: 9224 : h_obj[6][2][0] = th * A[9];
2758 : :
2759 [ + - ]: 9224 : h_obj[2][0][2] = -A[2] - A[6];
2760 [ + - ]: 9224 : h_obj[5][0][2] = A[2] - h * A[10];
2761 [ + - ]: 9224 : h_obj[7][0][2] = A[6] - h * A[10];
2762 [ + - ]: 9224 : h_obj[8][2][0] = th * A[10];
2763 : :
2764 [ + - ]: 9224 : h_obj[3][0][2] = -A[3] - A[7];
2765 [ + - ]: 9224 : h_obj[6][0][2] = A[3] - h * A[11];
2766 [ + - ]: 9224 : h_obj[8][0][2] = A[7] - h * A[11];
2767 [ + - ]: 9224 : h_obj[9][0][2] = th * A[11];
2768 : :
2769 : : /* Second block of rows */
2770 : 9224 : loc3 = matr[3] * f + dg[3] * cross;
2771 : 9224 : loc4 = dg[3] * g + matr[3] * cross;
2772 : :
2773 : 9224 : J_A[0] = loc0 + loc3 * matr[3] + loc4 * dg[3];
2774 : 9224 : J_A[1] = loc3 * matr[4] + loc4 * dg[4];
2775 : 9224 : J_A[2] = loc3 * matr[5] + loc4 * dg[5];
2776 : 9224 : J_B[0] = loc3 * matr[6] + loc4 * dg[6];
2777 : 9224 : J_B[1] = loc3 * matr[7] + loc4 * dg[7];
2778 : 9224 : J_B[2] = loc3 * matr[8] + loc4 * dg[8];
2779 : :
2780 : 9224 : loc3 = matr[4] * f + dg[4] * cross;
2781 : 9224 : loc4 = dg[4] * g + matr[4] * cross;
2782 : :
2783 : 9224 : J_A[3] = loc0 + loc3 * matr[4] + loc4 * dg[4];
2784 : 9224 : J_A[4] = loc3 * matr[5] + loc4 * dg[5];
2785 : 9224 : J_B[3] = loc3 * matr[6] + loc4 * dg[6];
2786 : 9224 : J_B[4] = loc3 * matr[7] + loc4 * dg[7];
2787 : 9224 : J_B[5] = loc3 * matr[8] + loc4 * dg[8];
2788 : :
2789 : 9224 : loc3 = matr[5] * f + dg[5] * cross;
2790 : 9224 : loc4 = dg[5] * g + matr[5] * cross;
2791 : :
2792 : 9224 : J_A[5] = loc0 + loc3 * matr[5] + loc4 * dg[5];
2793 : 9224 : J_B[6] = loc3 * matr[6] + loc4 * dg[6];
2794 : 9224 : J_B[7] = loc3 * matr[7] + loc4 * dg[7];
2795 : 9224 : J_B[8] = loc3 * matr[8] + loc4 * dg[8];
2796 : :
2797 : : /* Second diagonal block */
2798 : 9224 : A[0] = -J_A[0] - J_A[1];
2799 : 9224 : A[1] = J_A[0] - h * J_A[2];
2800 : 9224 : A[2] = J_A[1] - h * J_A[2];
2801 : 9224 : A[3] = th * J_A[2];
2802 : :
2803 : 9224 : A[4] = -J_A[1] - J_A[3];
2804 : 9224 : A[5] = J_A[1] - h * J_A[4];
2805 : 9224 : A[6] = J_A[3] - h * J_A[4];
2806 : 9224 : A[7] = th * J_A[4];
2807 : :
2808 : 9224 : A[8] = -J_A[2] - J_A[4];
2809 : 9224 : A[9] = J_A[2] - h * J_A[5];
2810 : 9224 : A[10] = J_A[4] - h * J_A[5];
2811 : 9224 : A[11] = th * J_A[5];
2812 : :
2813 [ + - ]: 9224 : h_obj[0][1][1] = -A[0] - A[4];
2814 [ + - ]: 9224 : h_obj[1][1][1] = A[0] - h * A[8];
2815 [ + - ]: 9224 : h_obj[2][1][1] = A[4] - h * A[8];
2816 [ + - ]: 9224 : h_obj[3][1][1] = th * A[8];
2817 : :
2818 [ + - ]: 9224 : h_obj[4][1][1] = A[1] - h * A[9];
2819 [ + - ]: 9224 : h_obj[5][1][1] = A[5] - h * A[9];
2820 [ + - ]: 9224 : h_obj[6][1][1] = th * A[9];
2821 : :
2822 [ + - ]: 9224 : h_obj[7][1][1] = A[6] - h * A[10];
2823 [ + - ]: 9224 : h_obj[8][1][1] = th * A[10];
2824 : :
2825 [ + - ]: 9224 : h_obj[9][1][1] = th * A[11];
2826 : :
2827 : : /* Third off-diagonal block */
2828 : 9224 : loc2 = matr[2] * loc1;
2829 : 9224 : J_B[1] += loc2;
2830 : 9224 : J_B[3] -= loc2;
2831 : :
2832 : 9224 : loc2 = matr[1] * loc1;
2833 : 9224 : J_B[2] -= loc2;
2834 : 9224 : J_B[6] += loc2;
2835 : :
2836 : 9224 : loc2 = matr[0] * loc1;
2837 : 9224 : J_B[5] += loc2;
2838 : 9224 : J_B[7] -= loc2;
2839 : :
2840 : 9224 : A[0] = -J_B[0] - J_B[1];
2841 : 9224 : A[1] = J_B[0] - h * J_B[2];
2842 : 9224 : A[2] = J_B[1] - h * J_B[2];
2843 : 9224 : A[3] = th * J_B[2];
2844 : :
2845 : 9224 : A[4] = -J_B[3] - J_B[4];
2846 : 9224 : A[5] = J_B[3] - h * J_B[5];
2847 : 9224 : A[6] = J_B[4] - h * J_B[5];
2848 : 9224 : A[7] = th * J_B[5];
2849 : :
2850 : 9224 : A[8] = -J_B[6] - J_B[7];
2851 : 9224 : A[9] = J_B[6] - h * J_B[8];
2852 : 9224 : A[10] = J_B[7] - h * J_B[8];
2853 : 9224 : A[11] = th * J_B[8];
2854 : :
2855 [ + - ]: 9224 : h_obj[0][1][2] = -A[0] - A[4];
2856 [ + - ]: 9224 : h_obj[1][2][1] = A[0] - h * A[8];
2857 [ + - ]: 9224 : h_obj[2][2][1] = A[4] - h * A[8];
2858 [ + - ]: 9224 : h_obj[3][2][1] = th * A[8];
2859 : :
2860 [ + - ]: 9224 : h_obj[1][1][2] = -A[1] - A[5];
2861 [ + - ]: 9224 : h_obj[4][1][2] = A[1] - h * A[9];
2862 [ + - ]: 9224 : h_obj[5][2][1] = A[5] - h * A[9];
2863 [ + - ]: 9224 : h_obj[6][2][1] = th * A[9];
2864 : :
2865 [ + - ]: 9224 : h_obj[2][1][2] = -A[2] - A[6];
2866 [ + - ]: 9224 : h_obj[5][1][2] = A[2] - h * A[10];
2867 [ + - ]: 9224 : h_obj[7][1][2] = A[6] - h * A[10];
2868 [ + - ]: 9224 : h_obj[8][2][1] = th * A[10];
2869 : :
2870 [ + - ]: 9224 : h_obj[3][1][2] = -A[3] - A[7];
2871 [ + - ]: 9224 : h_obj[6][1][2] = A[3] - h * A[11];
2872 [ + - ]: 9224 : h_obj[8][1][2] = A[7] - h * A[11];
2873 [ + - ]: 9224 : h_obj[9][1][2] = th * A[11];
2874 : :
2875 : : /* Third block of rows */
2876 : 9224 : loc3 = matr[6] * f + dg[6] * cross;
2877 : 9224 : loc4 = dg[6] * g + matr[6] * cross;
2878 : :
2879 : 9224 : J_A[0] = loc0 + loc3 * matr[6] + loc4 * dg[6];
2880 : 9224 : J_A[1] = loc3 * matr[7] + loc4 * dg[7];
2881 : 9224 : J_A[2] = loc3 * matr[8] + loc4 * dg[8];
2882 : :
2883 : 9224 : loc3 = matr[7] * f + dg[7] * cross;
2884 : 9224 : loc4 = dg[7] * g + matr[7] * cross;
2885 : :
2886 : 9224 : J_A[3] = loc0 + loc3 * matr[7] + loc4 * dg[7];
2887 : 9224 : J_A[4] = loc3 * matr[8] + loc4 * dg[8];
2888 : :
2889 : 9224 : loc3 = matr[8] * f + dg[8] * cross;
2890 : 9224 : loc4 = dg[8] * g + matr[8] * cross;
2891 : :
2892 : 9224 : J_A[5] = loc0 + loc3 * matr[8] + loc4 * dg[8];
2893 : :
2894 : : /* Third diagonal block */
2895 : 9224 : A[0] = -J_A[0] - J_A[1];
2896 : 9224 : A[1] = J_A[0] - h * J_A[2];
2897 : 9224 : A[2] = J_A[1] - h * J_A[2];
2898 : 9224 : A[3] = th * J_A[2];
2899 : :
2900 : 9224 : A[4] = -J_A[1] - J_A[3];
2901 : 9224 : A[5] = J_A[1] - h * J_A[4];
2902 : 9224 : A[6] = J_A[3] - h * J_A[4];
2903 : 9224 : A[7] = th * J_A[4];
2904 : :
2905 : 9224 : A[8] = -J_A[2] - J_A[4];
2906 : 9224 : A[9] = J_A[2] - h * J_A[5];
2907 : 9224 : A[10] = J_A[4] - h * J_A[5];
2908 : 9224 : A[11] = th * J_A[5];
2909 : :
2910 [ + - ]: 9224 : h_obj[0][2][2] = -A[0] - A[4];
2911 [ + - ]: 9224 : h_obj[1][2][2] = A[0] - h * A[8];
2912 [ + - ]: 9224 : h_obj[2][2][2] = A[4] - h * A[8];
2913 [ + - ]: 9224 : h_obj[3][2][2] = th * A[8];
2914 : :
2915 [ + - ]: 9224 : h_obj[4][2][2] = A[1] - h * A[9];
2916 [ + - ]: 9224 : h_obj[5][2][2] = A[5] - h * A[9];
2917 [ + - ]: 9224 : h_obj[6][2][2] = th * A[9];
2918 : :
2919 [ + - ]: 9224 : h_obj[7][2][2] = A[6] - h * A[10];
2920 [ + - ]: 9224 : h_obj[8][2][2] = th * A[10];
2921 : :
2922 [ + - ]: 9224 : h_obj[9][2][2] = th * A[11];
2923 : :
2924 : : #if 0
2925 : : int i, j, k, l;
2926 : : double nobj;
2927 : : Vector3D nx[4];
2928 : : Vector3D ng_obj[4];
2929 : :
2930 : : for (i = 0; i < 4; ++i) {
2931 : : nx[i] = x[i];
2932 : : ng_obj[i] = 0.0;
2933 : : }
2934 : :
2935 : : m_fcn_3p(obj, x, a, b, c);
2936 : : for (i = 0; i < 4; ++i) {
2937 : : for (j = 0; j < 3; ++j) {
2938 : : nx[i][j] += 1e-6;
2939 : : m_fcn_3p(nobj, nx, a, b, c);
2940 : : nx[i][j] = x[i][j];
2941 : :
2942 : : ng_obj[i][j] = (nobj - obj) / 1e-6;
2943 : : std::cout << i << " " << j << " " << g_obj[i][j] << " " << ng_obj[i][j] << std::endl;
2944 : : }
2945 : : }
2946 : : std::cout << std::endl;
2947 : :
2948 : : for (i = 0; i < 4; ++i) {
2949 : : for (j = 0; j < 3; ++j) {
2950 : : nx[i][j] += 1e-6;
2951 : : g_fcn_3p(nobj, ng_obj, nx, a, b, c);
2952 : : nx[i][j] = x[i][j];
2953 : :
2954 : : printf("%d %d", i, j);
2955 : : for (k = 0; k < 4; ++k) {
2956 : : for (l = 0; l < 3; ++l) {
2957 : : printf(" % 5.4e", (ng_obj[k][l] - g_obj[k][l]) / 1e-6);
2958 : : }
2959 : : }
2960 : : printf("\n");
2961 : : }
2962 : : }
2963 : :
2964 : : for (i = 0; i < 10; ++i) {
2965 : : std::cout << i << std::endl << h_obj[i] << std::endl << std::endl;
2966 : : }
2967 : : #endif
2968 : :
2969 : : // completes diagonal blocks.
2970 [ + - ]: 9224 : h_obj[0].fill_lower_triangle();
2971 [ + - ]: 9224 : h_obj[4].fill_lower_triangle();
2972 [ + - ]: 9224 : h_obj[7].fill_lower_triangle();
2973 [ + - ]: 9224 : h_obj[9].fill_lower_triangle();
2974 : 9224 : return true;
2975 : : }
2976 : :
2977 : : /*********************************************************************
2978 : : * Reference tetrahedral elements to corners of an ideal wedge element.
2979 : : * Vertices should be ordered such that the first three vertices form
2980 : : * the ideally-equaliteral face of the tetrahedron (the end of the
2981 : : * wedge) and the first and fourth vertices form the tetrahedral edge
2982 : : * orthogonal to the ideally-equaliteral face (the lateral edge of the
2983 : : * edge.)
2984 : : * 1 1/2 0 1 -1/sqrt(3) 0
2985 : : * W = 0 sqrt(3)/2 0 inv(W) = 0 2/sqrt(3) 0
2986 : : * 0 0 1 0 0 1
2987 : : *
2988 : : *********************************************************************/
2989 : 36 : inline bool m_fcn_3w( double& obj, const Vector3D x[4], const double a, const Exponent& b, const Exponent& c )
2990 : : {
2991 : : double matr[9], f, g;
2992 : : double loc1, loc2, loc3;
2993 : :
2994 : : /* Calculate M = A*inv(W). */
2995 [ + - ][ + - ]: 36 : matr[0] = x[1][0] - x[0][0];
2996 [ + - ][ + - ]: 36 : matr[1] = isqrt3 * ( 2 * x[2][0] - x[1][0] - x[0][0] );
[ + - ]
2997 [ + - ][ + - ]: 36 : matr[2] = x[3][0] - x[0][0];
2998 : :
2999 [ + - ][ + - ]: 36 : matr[3] = x[1][1] - x[0][1];
3000 [ + - ][ + - ]: 36 : matr[4] = isqrt3 * ( 2 * x[2][1] - x[1][1] - x[0][1] );
[ + - ]
3001 [ + - ][ + - ]: 36 : matr[5] = x[3][1] - x[0][1];
3002 : :
3003 [ + - ][ + - ]: 36 : matr[6] = x[1][2] - x[0][2];
3004 [ + - ][ + - ]: 36 : matr[7] = isqrt3 * ( 2 * x[2][2] - x[1][2] - x[0][2] );
[ + - ]
3005 [ + - ][ + - ]: 36 : matr[8] = x[3][2] - x[0][2];
3006 : :
3007 : : /* Calculate det(M). */
3008 : 36 : loc1 = matr[4] * matr[8] - matr[5] * matr[7];
3009 : 36 : loc2 = matr[5] * matr[6] - matr[3] * matr[8];
3010 : 36 : loc3 = matr[3] * matr[7] - matr[4] * matr[6];
3011 : 36 : g = matr[0] * loc1 + matr[1] * loc2 + matr[2] * loc3;
3012 [ - + ]: 36 : if( g < MSQ_MIN )
3013 : : {
3014 : 0 : obj = g;
3015 : 0 : return false;
3016 : : }
3017 : :
3018 : : /* Calculate norm(M). */
3019 : 72 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[2] * matr[2] + matr[3] * matr[3] + matr[4] * matr[4] +
3020 : 72 : matr[5] * matr[5] + matr[6] * matr[6] + matr[7] * matr[7] + matr[8] * matr[8];
3021 : :
3022 : : /* Calculate objective function. */
3023 [ + - ][ + - ]: 36 : obj = a * b.raise( f ) * c.raise( g );
3024 : 36 : return true;
3025 : : }
3026 : :
3027 : 812 : inline bool g_fcn_3w( double& obj, Vector3D g_obj[4], const Vector3D x[4], const double a, const Exponent& b,
3028 : : const Exponent& c )
3029 : : {
3030 : : double matr[9], f;
3031 : : double adj_m[9], g;
3032 : : double loc1, loc2, loc3;
3033 : :
3034 : : /* Calculate M = A*inv(W). */
3035 [ + - ][ + - ]: 812 : matr[0] = x[1][0] - x[0][0];
3036 [ + - ][ + - ]: 812 : matr[1] = isqrt3 * ( 2 * x[2][0] - x[1][0] - x[0][0] );
[ + - ]
3037 [ + - ][ + - ]: 812 : matr[2] = x[3][0] - x[0][0];
3038 : :
3039 [ + - ][ + - ]: 812 : matr[3] = x[1][1] - x[0][1];
3040 [ + - ][ + - ]: 812 : matr[4] = isqrt3 * ( 2 * x[2][1] - x[1][1] - x[0][1] );
[ + - ]
3041 [ + - ][ + - ]: 812 : matr[5] = x[3][1] - x[0][1];
3042 : :
3043 [ + - ][ + - ]: 812 : matr[6] = x[1][2] - x[0][2];
3044 [ + - ][ + - ]: 812 : matr[7] = isqrt3 * ( 2 * x[2][2] - x[1][2] - x[0][2] );
[ + - ]
3045 [ + - ][ + - ]: 812 : matr[8] = x[3][2] - x[0][2];
3046 : :
3047 : : /* Calculate det(M). */
3048 : 812 : loc1 = matr[4] * matr[8] - matr[5] * matr[7];
3049 : 812 : loc2 = matr[5] * matr[6] - matr[3] * matr[8];
3050 : 812 : loc3 = matr[3] * matr[7] - matr[4] * matr[6];
3051 : 812 : g = matr[0] * loc1 + matr[1] * loc2 + matr[2] * loc3;
3052 [ + + ]: 812 : if( g < MSQ_MIN )
3053 : : {
3054 : 1 : obj = g;
3055 : 1 : return false;
3056 : : }
3057 : :
3058 : : /* Calculate norm(M). */
3059 : 1622 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[2] * matr[2] + matr[3] * matr[3] + matr[4] * matr[4] +
3060 : 1622 : matr[5] * matr[5] + matr[6] * matr[6] + matr[7] * matr[7] + matr[8] * matr[8];
3061 : :
3062 : : /* Calculate objective function. */
3063 [ + - ][ + - ]: 811 : obj = a * b.raise( f ) * c.raise( g );
3064 : :
3065 : : /* Calculate the derivative of the objective function. */
3066 [ + - ]: 811 : f = b.value() * obj / f; /* Constant on nabla f */
3067 [ + - ]: 811 : g = c.value() * obj / g; /* Constant on nable g */
3068 : 811 : f *= 2.0; /* Modification for nabla f */
3069 : :
3070 : 811 : adj_m[0] = matr[0] * f + loc1 * g;
3071 : 811 : adj_m[1] = matr[1] * f + loc2 * g;
3072 : 811 : adj_m[2] = matr[2] * f + loc3 * g;
3073 : :
3074 : 811 : loc1 = matr[0] * g;
3075 : 811 : loc2 = matr[1] * g;
3076 : 811 : loc3 = matr[2] * g;
3077 : :
3078 : 811 : adj_m[3] = matr[3] * f + loc3 * matr[7] - loc2 * matr[8];
3079 : 811 : adj_m[4] = matr[4] * f + loc1 * matr[8] - loc3 * matr[6];
3080 : 811 : adj_m[5] = matr[5] * f + loc2 * matr[6] - loc1 * matr[7];
3081 : :
3082 : 811 : adj_m[6] = matr[6] * f + loc2 * matr[5] - loc3 * matr[4];
3083 : 811 : adj_m[7] = matr[7] * f + loc3 * matr[3] - loc1 * matr[5];
3084 : 811 : adj_m[8] = matr[8] * f + loc1 * matr[4] - loc2 * matr[3];
3085 : :
3086 [ + - ]: 811 : g_obj[0][0] = -adj_m[0] - isqrt3 * adj_m[1] - adj_m[2];
3087 [ + - ]: 811 : g_obj[1][0] = adj_m[0] - isqrt3 * adj_m[1];
3088 [ + - ]: 811 : g_obj[2][0] = tisqrt3 * adj_m[1];
3089 [ + - ]: 811 : g_obj[3][0] = adj_m[2];
3090 : :
3091 [ + - ]: 811 : g_obj[0][1] = -adj_m[3] - isqrt3 * adj_m[4] - adj_m[5];
3092 [ + - ]: 811 : g_obj[1][1] = adj_m[3] - isqrt3 * adj_m[4];
3093 [ + - ]: 811 : g_obj[2][1] = tisqrt3 * adj_m[4];
3094 [ + - ]: 811 : g_obj[3][1] = adj_m[5];
3095 : :
3096 [ + - ]: 811 : g_obj[0][2] = -adj_m[6] - isqrt3 * adj_m[7] - adj_m[8];
3097 [ + - ]: 811 : g_obj[1][2] = adj_m[6] - isqrt3 * adj_m[7];
3098 [ + - ]: 811 : g_obj[2][2] = tisqrt3 * adj_m[7];
3099 [ + - ]: 811 : g_obj[3][2] = adj_m[8];
3100 : :
3101 : 812 : return true;
3102 : : }
3103 : :
3104 : 792 : inline bool h_fcn_3w( double& obj, Vector3D g_obj[4], Matrix3D h_obj[10], const Vector3D x[4], const double a,
3105 : : const Exponent& b, const Exponent& c )
3106 : : {
3107 : : double matr[9], f;
3108 : : double adj_m[9], g;
3109 : : double dg[9], loc0, loc1, loc2, loc3, loc4;
3110 : : double A[12], J_A[6], J_B[9], J_C[9], cross;
3111 : :
3112 : : /* Calculate M = A*inv(W). */
3113 [ + - ][ + - ]: 792 : matr[0] = x[1][0] - x[0][0];
3114 [ + - ][ + - ]: 792 : matr[1] = isqrt3 * ( 2 * x[2][0] - x[1][0] - x[0][0] );
[ + - ]
3115 [ + - ][ + - ]: 792 : matr[2] = x[3][0] - x[0][0];
3116 : :
3117 [ + - ][ + - ]: 792 : matr[3] = x[1][1] - x[0][1];
3118 [ + - ][ + - ]: 792 : matr[4] = isqrt3 * ( 2 * x[2][1] - x[1][1] - x[0][1] );
[ + - ]
3119 [ + - ][ + - ]: 792 : matr[5] = x[3][1] - x[0][1];
3120 : :
3121 [ + - ][ + - ]: 792 : matr[6] = x[1][2] - x[0][2];
3122 [ + - ][ + - ]: 792 : matr[7] = isqrt3 * ( 2 * x[2][2] - x[1][2] - x[0][2] );
[ + - ]
3123 [ + - ][ + - ]: 792 : matr[8] = x[3][2] - x[0][2];
3124 : :
3125 : : /* Calculate det(M). */
3126 : 792 : dg[0] = matr[4] * matr[8] - matr[5] * matr[7];
3127 : 792 : dg[1] = matr[5] * matr[6] - matr[3] * matr[8];
3128 : 792 : dg[2] = matr[3] * matr[7] - matr[4] * matr[6];
3129 : 792 : g = matr[0] * dg[0] + matr[1] * dg[1] + matr[2] * dg[2];
3130 [ - + ]: 792 : if( g < MSQ_MIN )
3131 : : {
3132 : 0 : obj = g;
3133 : 0 : return false;
3134 : : }
3135 : :
3136 : 792 : dg[3] = matr[2] * matr[7] - matr[1] * matr[8];
3137 : 792 : dg[4] = matr[0] * matr[8] - matr[2] * matr[6];
3138 : 792 : dg[5] = matr[1] * matr[6] - matr[0] * matr[7];
3139 : 792 : dg[6] = matr[1] * matr[5] - matr[2] * matr[4];
3140 : 792 : dg[7] = matr[2] * matr[3] - matr[0] * matr[5];
3141 : 792 : dg[8] = matr[0] * matr[4] - matr[1] * matr[3];
3142 : :
3143 : : /* Calculate norm(M). */
3144 : 1584 : f = matr[0] * matr[0] + matr[1] * matr[1] + matr[2] * matr[2] + matr[3] * matr[3] + matr[4] * matr[4] +
3145 : 1584 : matr[5] * matr[5] + matr[6] * matr[6] + matr[7] * matr[7] + matr[8] * matr[8];
3146 : :
3147 : 792 : loc3 = f;
3148 : 792 : loc4 = g;
3149 : :
3150 : : /* Calculate objective function. */
3151 [ + - ][ + - ]: 792 : obj = a * b.raise( f ) * c.raise( g );
3152 : :
3153 : : /* Calculate the derivative of the objective function. */
3154 [ + - ]: 792 : f = b.value() * obj / f; /* Constant on nabla f */
3155 [ + - ]: 792 : g = c.value() * obj / g; /* Constant on nable g */
3156 : 792 : f *= 2.0; /* Modification for nabla f */
3157 : :
3158 : 792 : adj_m[0] = matr[0] * f + dg[0] * g;
3159 : 792 : adj_m[1] = matr[1] * f + dg[1] * g;
3160 : 792 : adj_m[2] = matr[2] * f + dg[2] * g;
3161 : 792 : adj_m[3] = matr[3] * f + dg[3] * g;
3162 : 792 : adj_m[4] = matr[4] * f + dg[4] * g;
3163 : 792 : adj_m[5] = matr[5] * f + dg[5] * g;
3164 : 792 : adj_m[6] = matr[6] * f + dg[6] * g;
3165 : 792 : adj_m[7] = matr[7] * f + dg[7] * g;
3166 : 792 : adj_m[8] = matr[8] * f + dg[8] * g;
3167 : :
3168 [ + - ]: 792 : g_obj[0][0] = -adj_m[0] - isqrt3 * adj_m[1] - adj_m[2];
3169 [ + - ]: 792 : g_obj[1][0] = adj_m[0] - isqrt3 * adj_m[1];
3170 [ + - ]: 792 : g_obj[2][0] = tisqrt3 * adj_m[1];
3171 [ + - ]: 792 : g_obj[3][0] = adj_m[2];
3172 : :
3173 [ + - ]: 792 : g_obj[0][1] = -adj_m[3] - isqrt3 * adj_m[4] - adj_m[5];
3174 [ + - ]: 792 : g_obj[1][1] = adj_m[3] - isqrt3 * adj_m[4];
3175 [ + - ]: 792 : g_obj[2][1] = tisqrt3 * adj_m[4];
3176 [ + - ]: 792 : g_obj[3][1] = adj_m[5];
3177 : :
3178 [ + - ]: 792 : g_obj[0][2] = -adj_m[6] - isqrt3 * adj_m[7] - adj_m[8];
3179 [ + - ]: 792 : g_obj[1][2] = adj_m[6] - isqrt3 * adj_m[7];
3180 [ + - ]: 792 : g_obj[2][2] = tisqrt3 * adj_m[7];
3181 [ + - ]: 792 : g_obj[3][2] = adj_m[8];
3182 : :
3183 : : /* Calculate the hessian of the objective. */
3184 : 792 : loc0 = f; /* Constant on nabla^2 f */
3185 : 792 : loc1 = g; /* Constant on nabla^2 g */
3186 [ + - ]: 792 : cross = f * c.value() / loc4; /* Constant on nabla g nabla f */
3187 [ + - ]: 792 : f = f * ( b.value() - 1 ) / loc3; /* Constant on nabla f nabla f */
3188 [ + - ]: 792 : g = g * ( c.value() - 1 ) / loc4; /* Constant on nabla g nabla g */
3189 : 792 : f *= 2.0; /* Modification for nabla f */
3190 : :
3191 : : /* First block of rows */
3192 : 792 : loc3 = matr[0] * f + dg[0] * cross;
3193 : 792 : loc4 = dg[0] * g + matr[0] * cross;
3194 : :
3195 : 792 : J_A[0] = loc0 + loc3 * matr[0] + loc4 * dg[0];
3196 : 792 : J_A[1] = loc3 * matr[1] + loc4 * dg[1];
3197 : 792 : J_A[2] = loc3 * matr[2] + loc4 * dg[2];
3198 : 792 : J_B[0] = loc3 * matr[3] + loc4 * dg[3];
3199 : 792 : J_B[1] = loc3 * matr[4] + loc4 * dg[4];
3200 : 792 : J_B[2] = loc3 * matr[5] + loc4 * dg[5];
3201 : 792 : J_C[0] = loc3 * matr[6] + loc4 * dg[6];
3202 : 792 : J_C[1] = loc3 * matr[7] + loc4 * dg[7];
3203 : 792 : J_C[2] = loc3 * matr[8] + loc4 * dg[8];
3204 : :
3205 : 792 : loc3 = matr[1] * f + dg[1] * cross;
3206 : 792 : loc4 = dg[1] * g + matr[1] * cross;
3207 : :
3208 : 792 : J_A[3] = loc0 + loc3 * matr[1] + loc4 * dg[1];
3209 : 792 : J_A[4] = loc3 * matr[2] + loc4 * dg[2];
3210 : 792 : J_B[3] = loc3 * matr[3] + loc4 * dg[3];
3211 : 792 : J_B[4] = loc3 * matr[4] + loc4 * dg[4];
3212 : 792 : J_B[5] = loc3 * matr[5] + loc4 * dg[5];
3213 : 792 : J_C[3] = loc3 * matr[6] + loc4 * dg[6];
3214 : 792 : J_C[4] = loc3 * matr[7] + loc4 * dg[7];
3215 : 792 : J_C[5] = loc3 * matr[8] + loc4 * dg[8];
3216 : :
3217 : 792 : loc3 = matr[2] * f + dg[2] * cross;
3218 : 792 : loc4 = dg[2] * g + matr[2] * cross;
3219 : :
3220 : 792 : J_A[5] = loc0 + loc3 * matr[2] + loc4 * dg[2];
3221 : 792 : J_B[6] = loc3 * matr[3] + loc4 * dg[3];
3222 : 792 : J_B[7] = loc3 * matr[4] + loc4 * dg[4];
3223 : 792 : J_B[8] = loc3 * matr[5] + loc4 * dg[5];
3224 : 792 : J_C[6] = loc3 * matr[6] + loc4 * dg[6];
3225 : 792 : J_C[7] = loc3 * matr[7] + loc4 * dg[7];
3226 : 792 : J_C[8] = loc3 * matr[8] + loc4 * dg[8];
3227 : :
3228 : : /* First diagonal block */
3229 : 792 : A[0] = -J_A[0] - isqrt3 * J_A[1] - J_A[2];
3230 : 792 : A[1] = J_A[0] - isqrt3 * J_A[1];
3231 : : // A[2] = tisqrt3*J_A[1];
3232 : : // A[3] = J_A[2];
3233 : :
3234 : 792 : A[4] = -J_A[1] - isqrt3 * J_A[3] - J_A[4];
3235 : 792 : A[5] = J_A[1] - isqrt3 * J_A[3];
3236 : 792 : A[6] = tisqrt3 * J_A[3];
3237 : : // A[7] = J_A[4];
3238 : :
3239 : 792 : A[8] = -J_A[2] - isqrt3 * J_A[4] - J_A[5];
3240 : 792 : A[9] = J_A[2] - isqrt3 * J_A[4];
3241 : 792 : A[10] = tisqrt3 * J_A[4];
3242 : 792 : A[11] = J_A[5];
3243 : :
3244 [ + - ]: 792 : h_obj[0][0][0] = -A[0] - isqrt3 * A[4] - A[8];
3245 [ + - ]: 792 : h_obj[1][0][0] = A[0] - isqrt3 * A[4];
3246 [ + - ]: 792 : h_obj[2][0][0] = tisqrt3 * A[4];
3247 [ + - ]: 792 : h_obj[3][0][0] = A[8];
3248 : :
3249 [ + - ]: 792 : h_obj[4][0][0] = A[1] - isqrt3 * A[5];
3250 [ + - ]: 792 : h_obj[5][0][0] = tisqrt3 * A[5];
3251 [ + - ]: 792 : h_obj[6][0][0] = A[9];
3252 : :
3253 [ + - ]: 792 : h_obj[7][0][0] = tisqrt3 * A[6];
3254 [ + - ]: 792 : h_obj[8][0][0] = A[10];
3255 : :
3256 [ + - ]: 792 : h_obj[9][0][0] = A[11];
3257 : :
3258 : : /* First off-diagonal block */
3259 : 792 : loc2 = matr[8] * loc1;
3260 : 792 : J_B[1] += loc2;
3261 : 792 : J_B[3] -= loc2;
3262 : :
3263 : 792 : loc2 = matr[7] * loc1;
3264 : 792 : J_B[2] -= loc2;
3265 : 792 : J_B[6] += loc2;
3266 : :
3267 : 792 : loc2 = matr[6] * loc1;
3268 : 792 : J_B[5] += loc2;
3269 : 792 : J_B[7] -= loc2;
3270 : :
3271 : 792 : A[0] = -J_B[0] - isqrt3 * J_B[1] - J_B[2];
3272 : 792 : A[1] = J_B[0] - isqrt3 * J_B[1];
3273 : 792 : A[2] = tisqrt3 * J_B[1];
3274 : 792 : A[3] = J_B[2];
3275 : :
3276 : 792 : A[4] = -J_B[3] - isqrt3 * J_B[4] - J_B[5];
3277 : 792 : A[5] = J_B[3] - isqrt3 * J_B[4];
3278 : 792 : A[6] = tisqrt3 * J_B[4];
3279 : 792 : A[7] = J_B[5];
3280 : :
3281 : 792 : A[8] = -J_B[6] - isqrt3 * J_B[7] - J_B[8];
3282 : 792 : A[9] = J_B[6] - isqrt3 * J_B[7];
3283 : 792 : A[10] = tisqrt3 * J_B[7];
3284 : 792 : A[11] = J_B[8];
3285 : :
3286 [ + - ]: 792 : h_obj[0][0][1] = -A[0] - isqrt3 * A[4] - A[8];
3287 [ + - ]: 792 : h_obj[1][1][0] = A[0] - isqrt3 * A[4];
3288 [ + - ]: 792 : h_obj[2][1][0] = tisqrt3 * A[4];
3289 [ + - ]: 792 : h_obj[3][1][0] = A[8];
3290 : :
3291 [ + - ]: 792 : h_obj[1][0][1] = -A[1] - isqrt3 * A[5] - A[9];
3292 [ + - ]: 792 : h_obj[4][0][1] = A[1] - isqrt3 * A[5];
3293 [ + - ]: 792 : h_obj[5][1][0] = tisqrt3 * A[5];
3294 [ + - ]: 792 : h_obj[6][1][0] = A[9];
3295 : :
3296 [ + - ]: 792 : h_obj[2][0][1] = -A[2] - isqrt3 * A[6] - A[10];
3297 [ + - ]: 792 : h_obj[5][0][1] = A[2] - isqrt3 * A[6];
3298 [ + - ]: 792 : h_obj[7][0][1] = tisqrt3 * A[6];
3299 [ + - ]: 792 : h_obj[8][1][0] = A[10];
3300 : :
3301 [ + - ]: 792 : h_obj[3][0][1] = -A[3] - isqrt3 * A[7] - A[11];
3302 [ + - ]: 792 : h_obj[6][0][1] = A[3] - isqrt3 * A[7];
3303 [ + - ]: 792 : h_obj[8][0][1] = tisqrt3 * A[7];
3304 [ + - ]: 792 : h_obj[9][0][1] = A[11];
3305 : :
3306 : : /* Second off-diagonal block */
3307 : 792 : loc2 = matr[5] * loc1;
3308 : 792 : J_C[1] -= loc2;
3309 : 792 : J_C[3] += loc2;
3310 : :
3311 : 792 : loc2 = matr[4] * loc1;
3312 : 792 : J_C[2] += loc2;
3313 : 792 : J_C[6] -= loc2;
3314 : :
3315 : 792 : loc2 = matr[3] * loc1;
3316 : 792 : J_C[5] -= loc2;
3317 : 792 : J_C[7] += loc2;
3318 : :
3319 : 792 : A[0] = -J_C[0] - isqrt3 * J_C[1] - J_C[2];
3320 : 792 : A[1] = J_C[0] - isqrt3 * J_C[1];
3321 : 792 : A[2] = tisqrt3 * J_C[1];
3322 : 792 : A[3] = J_C[2];
3323 : :
3324 : 792 : A[4] = -J_C[3] - isqrt3 * J_C[4] - J_C[5];
3325 : 792 : A[5] = J_C[3] - isqrt3 * J_C[4];
3326 : 792 : A[6] = tisqrt3 * J_C[4];
3327 : 792 : A[7] = J_C[5];
3328 : :
3329 : 792 : A[8] = -J_C[6] - isqrt3 * J_C[7] - J_C[8];
3330 : 792 : A[9] = J_C[6] - isqrt3 * J_C[7];
3331 : 792 : A[10] = tisqrt3 * J_C[7];
3332 : 792 : A[11] = J_C[8];
3333 : :
3334 [ + - ]: 792 : h_obj[0][0][2] = -A[0] - isqrt3 * A[4] - A[8];
3335 [ + - ]: 792 : h_obj[1][2][0] = A[0] - isqrt3 * A[4];
3336 [ + - ]: 792 : h_obj[2][2][0] = tisqrt3 * A[4];
3337 [ + - ]: 792 : h_obj[3][2][0] = A[8];
3338 : :
3339 [ + - ]: 792 : h_obj[1][0][2] = -A[1] - isqrt3 * A[5] - A[9];
3340 [ + - ]: 792 : h_obj[4][0][2] = A[1] - isqrt3 * A[5];
3341 [ + - ]: 792 : h_obj[5][2][0] = tisqrt3 * A[5];
3342 [ + - ]: 792 : h_obj[6][2][0] = A[9];
3343 : :
3344 [ + - ]: 792 : h_obj[2][0][2] = -A[2] - isqrt3 * A[6] - A[10];
3345 [ + - ]: 792 : h_obj[5][0][2] = A[2] - isqrt3 * A[6];
3346 [ + - ]: 792 : h_obj[7][0][2] = tisqrt3 * A[6];
3347 [ + - ]: 792 : h_obj[8][2][0] = A[10];
3348 : :
3349 [ + - ]: 792 : h_obj[3][0][2] = -A[3] - isqrt3 * A[7] - A[11];
3350 [ + - ]: 792 : h_obj[6][0][2] = A[3] - isqrt3 * A[7];
3351 [ + - ]: 792 : h_obj[8][0][2] = tisqrt3 * A[7];
3352 [ + - ]: 792 : h_obj[9][0][2] = A[11];
3353 : :
3354 : : /* Second block of rows */
3355 : 792 : loc3 = matr[3] * f + dg[3] * cross;
3356 : 792 : loc4 = dg[3] * g + matr[3] * cross;
3357 : :
3358 : 792 : J_A[0] = loc0 + loc3 * matr[3] + loc4 * dg[3];
3359 : 792 : J_A[1] = loc3 * matr[4] + loc4 * dg[4];
3360 : 792 : J_A[2] = loc3 * matr[5] + loc4 * dg[5];
3361 : 792 : J_B[0] = loc3 * matr[6] + loc4 * dg[6];
3362 : 792 : J_B[1] = loc3 * matr[7] + loc4 * dg[7];
3363 : 792 : J_B[2] = loc3 * matr[8] + loc4 * dg[8];
3364 : :
3365 : 792 : loc3 = matr[4] * f + dg[4] * cross;
3366 : 792 : loc4 = dg[4] * g + matr[4] * cross;
3367 : :
3368 : 792 : J_A[3] = loc0 + loc3 * matr[4] + loc4 * dg[4];
3369 : 792 : J_A[4] = loc3 * matr[5] + loc4 * dg[5];
3370 : 792 : J_B[3] = loc3 * matr[6] + loc4 * dg[6];
3371 : 792 : J_B[4] = loc3 * matr[7] + loc4 * dg[7];
3372 : 792 : J_B[5] = loc3 * matr[8] + loc4 * dg[8];
3373 : :
3374 : 792 : loc3 = matr[5] * f + dg[5] * cross;
3375 : 792 : loc4 = dg[5] * g + matr[5] * cross;
3376 : :
3377 : 792 : J_A[5] = loc0 + loc3 * matr[5] + loc4 * dg[5];
3378 : 792 : J_B[6] = loc3 * matr[6] + loc4 * dg[6];
3379 : 792 : J_B[7] = loc3 * matr[7] + loc4 * dg[7];
3380 : 792 : J_B[8] = loc3 * matr[8] + loc4 * dg[8];
3381 : :
3382 : : /* Second diagonal block */
3383 : 792 : A[0] = -J_A[0] - isqrt3 * J_A[1] - J_A[2];
3384 : 792 : A[1] = J_A[0] - isqrt3 * J_A[1];
3385 : : // A[2] = tisqrt3*J_A[1];
3386 : : // A[3] = J_A[2];
3387 : :
3388 : 792 : A[4] = -J_A[1] - isqrt3 * J_A[3] - J_A[4];
3389 : 792 : A[5] = J_A[1] - isqrt3 * J_A[3];
3390 : 792 : A[6] = tisqrt3 * J_A[3];
3391 : : // A[7] = J_A[4];
3392 : :
3393 : 792 : A[8] = -J_A[2] - isqrt3 * J_A[4] - J_A[5];
3394 : 792 : A[9] = J_A[2] - isqrt3 * J_A[4];
3395 : 792 : A[10] = tisqrt3 * J_A[4];
3396 : 792 : A[11] = J_A[5];
3397 : :
3398 [ + - ]: 792 : h_obj[0][1][1] = -A[0] - isqrt3 * A[4] - A[8];
3399 [ + - ]: 792 : h_obj[1][1][1] = A[0] - isqrt3 * A[4];
3400 [ + - ]: 792 : h_obj[2][1][1] = tisqrt3 * A[4];
3401 [ + - ]: 792 : h_obj[3][1][1] = A[8];
3402 : :
3403 [ + - ]: 792 : h_obj[4][1][1] = A[1] - isqrt3 * A[5];
3404 [ + - ]: 792 : h_obj[5][1][1] = tisqrt3 * A[5];
3405 [ + - ]: 792 : h_obj[6][1][1] = A[9];
3406 : :
3407 [ + - ]: 792 : h_obj[7][1][1] = tisqrt3 * A[6];
3408 [ + - ]: 792 : h_obj[8][1][1] = A[10];
3409 : :
3410 [ + - ]: 792 : h_obj[9][1][1] = A[11];
3411 : :
3412 : : /* Third off-diagonal block */
3413 : 792 : loc2 = matr[2] * loc1;
3414 : 792 : J_B[1] += loc2;
3415 : 792 : J_B[3] -= loc2;
3416 : :
3417 : 792 : loc2 = matr[1] * loc1;
3418 : 792 : J_B[2] -= loc2;
3419 : 792 : J_B[6] += loc2;
3420 : :
3421 : 792 : loc2 = matr[0] * loc1;
3422 : 792 : J_B[5] += loc2;
3423 : 792 : J_B[7] -= loc2;
3424 : :
3425 : 792 : A[0] = -J_B[0] - isqrt3 * J_B[1] - J_B[2];
3426 : 792 : A[1] = J_B[0] - isqrt3 * J_B[1];
3427 : 792 : A[2] = tisqrt3 * J_B[1];
3428 : 792 : A[3] = J_B[2];
3429 : :
3430 : 792 : A[4] = -J_B[3] - isqrt3 * J_B[4] - J_B[5];
3431 : 792 : A[5] = J_B[3] - isqrt3 * J_B[4];
3432 : 792 : A[6] = tisqrt3 * J_B[4];
3433 : 792 : A[7] = J_B[5];
3434 : :
3435 : 792 : A[8] = -J_B[6] - isqrt3 * J_B[7] - J_B[8];
3436 : 792 : A[9] = J_B[6] - isqrt3 * J_B[7];
3437 : 792 : A[10] = tisqrt3 * J_B[7];
3438 : 792 : A[11] = J_B[8];
3439 : :
3440 [ + - ]: 792 : h_obj[0][1][2] = -A[0] - isqrt3 * A[4] - A[8];
3441 [ + - ]: 792 : h_obj[1][2][1] = A[0] - isqrt3 * A[4];
3442 [ + - ]: 792 : h_obj[2][2][1] = tisqrt3 * A[4];
3443 [ + - ]: 792 : h_obj[3][2][1] = A[8];
3444 : :
3445 [ + - ]: 792 : h_obj[1][1][2] = -A[1] - isqrt3 * A[5] - A[9];
3446 [ + - ]: 792 : h_obj[4][1][2] = A[1] - isqrt3 * A[5];
3447 [ + - ]: 792 : h_obj[5][2][1] = tisqrt3 * A[5];
3448 [ + - ]: 792 : h_obj[6][2][1] = A[9];
3449 : :
3450 [ + - ]: 792 : h_obj[2][1][2] = -A[2] - isqrt3 * A[6] - A[10];
3451 [ + - ]: 792 : h_obj[5][1][2] = A[2] - isqrt3 * A[6];
3452 [ + - ]: 792 : h_obj[7][1][2] = tisqrt3 * A[6];
3453 [ + - ]: 792 : h_obj[8][2][1] = A[10];
3454 : :
3455 [ + - ]: 792 : h_obj[3][1][2] = -A[3] - isqrt3 * A[7] - A[11];
3456 [ + - ]: 792 : h_obj[6][1][2] = A[3] - isqrt3 * A[7];
3457 [ + - ]: 792 : h_obj[8][1][2] = tisqrt3 * A[7];
3458 [ + - ]: 792 : h_obj[9][1][2] = A[11];
3459 : :
3460 : : /* Third block of rows */
3461 : 792 : loc3 = matr[6] * f + dg[6] * cross;
3462 : 792 : loc4 = dg[6] * g + matr[6] * cross;
3463 : :
3464 : 792 : J_A[0] = loc0 + loc3 * matr[6] + loc4 * dg[6];
3465 : 792 : J_A[1] = loc3 * matr[7] + loc4 * dg[7];
3466 : 792 : J_A[2] = loc3 * matr[8] + loc4 * dg[8];
3467 : :
3468 : 792 : loc3 = matr[7] * f + dg[7] * cross;
3469 : 792 : loc4 = dg[7] * g + matr[7] * cross;
3470 : :
3471 : 792 : J_A[3] = loc0 + loc3 * matr[7] + loc4 * dg[7];
3472 : 792 : J_A[4] = loc3 * matr[8] + loc4 * dg[8];
3473 : :
3474 : 792 : loc3 = matr[8] * f + dg[8] * cross;
3475 : 792 : loc4 = dg[8] * g + matr[8] * cross;
3476 : :
3477 : 792 : J_A[5] = loc0 + loc3 * matr[8] + loc4 * dg[8];
3478 : :
3479 : : /* Third diagonal block */
3480 : 792 : A[0] = -J_A[0] - isqrt3 * J_A[1] - J_A[2];
3481 : 792 : A[1] = J_A[0] - isqrt3 * J_A[1];
3482 : : // A[2] = tisqrt3*J_A[1];
3483 : : // A[3] = J_A[2];
3484 : :
3485 : 792 : A[4] = -J_A[1] - isqrt3 * J_A[3] - J_A[4];
3486 : 792 : A[5] = J_A[1] - isqrt3 * J_A[3];
3487 : 792 : A[6] = tisqrt3 * J_A[3];
3488 : : // A[7] = J_A[4];
3489 : :
3490 : 792 : A[8] = -J_A[2] - isqrt3 * J_A[4] - J_A[5];
3491 : 792 : A[9] = J_A[2] - isqrt3 * J_A[4];
3492 : 792 : A[10] = tisqrt3 * J_A[4];
3493 : 792 : A[11] = J_A[5];
3494 : :
3495 [ + - ]: 792 : h_obj[0][2][2] = -A[0] - isqrt3 * A[4] - A[8];
3496 [ + - ]: 792 : h_obj[1][2][2] = A[0] - isqrt3 * A[4];
3497 [ + - ]: 792 : h_obj[2][2][2] = tisqrt3 * A[4];
3498 [ + - ]: 792 : h_obj[3][2][2] = A[8];
3499 : :
3500 [ + - ]: 792 : h_obj[4][2][2] = A[1] - isqrt3 * A[5];
3501 [ + - ]: 792 : h_obj[5][2][2] = tisqrt3 * A[5];
3502 [ + - ]: 792 : h_obj[6][2][2] = A[9];
3503 : :
3504 [ + - ]: 792 : h_obj[7][2][2] = tisqrt3 * A[6];
3505 [ + - ]: 792 : h_obj[8][2][2] = A[10];
3506 : :
3507 [ + - ]: 792 : h_obj[9][2][2] = A[11];
3508 : :
3509 : : // completes diagonal blocks.
3510 [ + - ]: 792 : h_obj[0].fill_lower_triangle();
3511 [ + - ]: 792 : h_obj[4].fill_lower_triangle();
3512 [ + - ]: 792 : h_obj[7].fill_lower_triangle();
3513 [ + - ]: 792 : h_obj[9].fill_lower_triangle();
3514 : :
3515 : 792 : return true;
3516 : : }
3517 : :
3518 : : } // namespace MBMesquite
3519 : :
3520 : : #endif
|
