1 |
|
|
/* |
2 |
|
|
Teem: Tools to process and visualize scientific data and images . |
3 |
|
|
Copyright (C) 2013, 2012, 2011, 2010, 2009 University of Chicago |
4 |
|
|
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann |
5 |
|
|
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah |
6 |
|
|
|
7 |
|
|
This library is free software; you can redistribute it and/or |
8 |
|
|
modify it under the terms of the GNU Lesser General Public License |
9 |
|
|
(LGPL) as published by the Free Software Foundation; either |
10 |
|
|
version 2.1 of the License, or (at your option) any later version. |
11 |
|
|
The terms of redistributing and/or modifying this software also |
12 |
|
|
include exceptions to the LGPL that facilitate static linking. |
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 |
|
|
along with this library; if not, write to Free Software Foundation, Inc., |
21 |
|
|
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
22 |
|
|
*/ |
23 |
|
|
|
24 |
|
|
#include "gage.h" |
25 |
|
|
#include "privateGage.h" |
26 |
|
|
|
27 |
|
|
/* |
28 |
|
|
* highly inefficient computation of the imaginary part of complex |
29 |
|
|
* conjugate eigenvalues of a 3x3 non-symmetric matrix |
30 |
|
|
*/ |
31 |
|
|
double |
32 |
|
|
gage_imaginary_part_eigenvalues(double *M ) { |
33 |
|
|
double A, B, C, scale, frob, m[9], _eval[3]; |
34 |
|
|
double beta, _gamma; |
35 |
|
|
int roots; |
36 |
|
|
|
37 |
|
|
frob = ELL_3M_FROB(M); |
38 |
|
|
scale = frob > 10 ? 10.0/frob : 1.0; |
39 |
|
|
ELL_3M_SCALE(m, scale, M); |
40 |
|
|
/* |
41 |
|
|
** from gordon with mathematica; these are the coefficients of the |
42 |
|
|
** cubic polynomial in x: det(x*I - M). The full cubic is |
43 |
|
|
** x^3 + A*x^2 + B*x + C. |
44 |
|
|
*/ |
45 |
|
|
A = -m[0] - m[4] - m[8]; |
46 |
|
|
B = m[0]*m[4] - m[3]*m[1] |
47 |
|
|
+ m[0]*m[8] - m[6]*m[2] |
48 |
|
|
+ m[4]*m[8] - m[7]*m[5]; |
49 |
|
|
C = (m[6]*m[4] - m[3]*m[7])*m[2] |
50 |
|
|
+ (m[0]*m[7] - m[6]*m[1])*m[5] |
51 |
|
|
+ (m[3]*m[1] - m[0]*m[4])*m[8]; |
52 |
|
|
roots = ell_cubic(_eval, A, B, C, AIR_TRUE); |
53 |
|
|
if ( roots != ell_cubic_root_single ) |
54 |
|
|
return 0.; |
55 |
|
|
|
56 |
|
|
/* 2 complex conjuguate eigenvalues */ |
57 |
|
|
beta = A + _eval[0]; |
58 |
|
|
_gamma = -C/_eval[0]; |
59 |
|
|
return sqrt( 4.*_gamma - beta*beta ); |
60 |
|
|
} |
61 |
|
|
|
62 |
|
|
|
63 |
|
|
gageItemEntry |
64 |
|
|
_gageVecTable[GAGE_VEC_ITEM_MAX+1] = { |
65 |
|
|
/* enum value len, deriv, prereqs, parent item, parent index, needData */ |
66 |
|
|
{gageVecUnknown, 0, 0, {0}, 0, 0, AIR_FALSE}, |
67 |
|
|
{gageVecVector, 3, 0, {0}, 0, 0, AIR_FALSE}, |
68 |
|
|
{gageVecVector0, 1, 0, {gageVecVector}, gageVecVector, 0, AIR_FALSE}, |
69 |
|
|
{gageVecVector1, 1, 0, {gageVecVector}, gageVecVector, 1, AIR_FALSE}, |
70 |
|
|
{gageVecVector2, 1, 0, {gageVecVector}, gageVecVector, 2, AIR_FALSE}, |
71 |
|
|
{gageVecLength, 1, 0, {gageVecVector}, 0, 0, AIR_FALSE}, |
72 |
|
|
{gageVecNormalized, 3, 0, {gageVecVector, gageVecLength}, 0, 0, AIR_FALSE}, |
73 |
|
|
{gageVecJacobian, 9, 1, {0}, 0, 0, AIR_FALSE}, |
74 |
|
|
{gageVecStrain, 9, 1, {gageVecJacobian}, 0, 0, AIR_FALSE}, |
75 |
|
|
{gageVecDivergence, 1, 1, {gageVecJacobian}, 0, 0, AIR_FALSE}, |
76 |
|
|
{gageVecCurl, 3, 1, {gageVecJacobian}, 0, 0, AIR_FALSE}, |
77 |
|
|
{gageVecCurlNorm, 1, 1, {gageVecCurl}, 0, 0, AIR_FALSE}, |
78 |
|
|
{gageVecHelicity, 1, 1, {gageVecVector, gageVecCurl}, 0, 0, AIR_FALSE}, |
79 |
|
|
{gageVecNormHelicity, 1, 1, {gageVecNormalized, gageVecCurl}, 0, 0, AIR_FALSE}, |
80 |
|
|
{gageVecSOmega, 9, 1, {gageVecJacobian,gageVecStrain}, 0, 0, AIR_FALSE}, |
81 |
|
|
{gageVecLambda2, 1, 1, {gageVecSOmega}, 0, 0, AIR_FALSE}, |
82 |
|
|
{gageVecImaginaryPart, 1, 1, {gageVecJacobian}, 0, 0, AIR_FALSE}, |
83 |
|
|
{gageVecHessian, 27, 2, {0}, 0, 0, AIR_FALSE}, |
84 |
|
|
{gageVecDivGradient, 3, 1, {gageVecHessian}, 0, 0, AIR_FALSE}, |
85 |
|
|
{gageVecCurlGradient, 9, 2, {gageVecHessian}, 0, 0, AIR_FALSE}, |
86 |
|
|
{gageVecCurlNormGrad, 3, 2, {gageVecHessian, gageVecCurl}, 0, 0, AIR_FALSE}, |
87 |
|
|
{gageVecNCurlNormGrad, 3, 2, {gageVecCurlNormGrad}, 0, 0, AIR_FALSE}, |
88 |
|
|
{gageVecHelGradient, 3, 2, {gageVecVector, gageVecJacobian, gageVecCurl, gageVecCurlGradient}, 0, 0, AIR_FALSE}, |
89 |
|
|
{gageVecDirHelDeriv, 1, 2, {gageVecNormalized, gageVecHelGradient}, 0, 0, AIR_FALSE}, |
90 |
|
|
{gageVecProjHelGradient, 3, 2, {gageVecNormalized, gageVecHelGradient, gageVecDirHelDeriv}, 0, 0, AIR_FALSE}, |
91 |
|
|
/* HEY: these should change to sub-items!!! */ |
92 |
|
|
{gageVecGradient0, 3, 1, {gageVecJacobian}, 0, 0, AIR_FALSE}, |
93 |
|
|
{gageVecGradient1, 3, 1, {gageVecJacobian}, 0, 0, AIR_FALSE}, |
94 |
|
|
{gageVecGradient2, 3, 1, {gageVecJacobian}, 0, 0, AIR_FALSE}, |
95 |
|
|
{gageVecMultiGrad, 9, 1, {gageVecGradient0, gageVecGradient1, gageVecGradient2}, 0, 0, AIR_FALSE}, |
96 |
|
|
{gageVecMGFrob, 1, 1, {gageVecMultiGrad}, 0, 0, AIR_FALSE}, |
97 |
|
|
{gageVecMGEval, 3, 1, {gageVecMultiGrad}, 0, 0, AIR_FALSE}, |
98 |
|
|
{gageVecMGEvec, 9, 1, {gageVecMultiGrad, gageVecMGEval}, 0, 0, AIR_FALSE} |
99 |
|
|
}; |
100 |
|
|
|
101 |
|
|
void |
102 |
|
|
_gageVecFilter(gageContext *ctx, gagePerVolume *pvl) { |
103 |
|
59400 |
char me[]="_gageVecFilter"; |
104 |
|
|
double *fw00, *fw11, *fw22, *vec, *jac, *hes; |
105 |
|
|
int fd; |
106 |
|
29700 |
gageScl3PFilter_t *filter[5] = {NULL, gageScl3PFilter2, gageScl3PFilter4, |
107 |
|
|
gageScl3PFilter6, gageScl3PFilter8}; |
108 |
|
|
unsigned int valIdx; |
109 |
|
|
|
110 |
|
29700 |
fd = 2*ctx->radius; |
111 |
|
29700 |
vec = pvl->directAnswer[gageVecVector]; |
112 |
|
29700 |
jac = pvl->directAnswer[gageVecJacobian]; |
113 |
|
29700 |
hes = pvl->directAnswer[gageVecHessian]; |
114 |
✗✓ |
29700 |
if (!ctx->parm.k3pack) { |
115 |
|
|
fprintf(stderr, "!%s: sorry, 6pack filtering not implemented\n", me); |
116 |
|
|
return; |
117 |
|
|
} |
118 |
|
29700 |
fw00 = ctx->fw + fd*3*gageKernel00; |
119 |
|
29700 |
fw11 = ctx->fw + fd*3*gageKernel11; |
120 |
|
29700 |
fw22 = ctx->fw + fd*3*gageKernel22; |
121 |
|
|
/* perform the filtering */ |
122 |
✓✓ |
29700 |
if (fd <= 8) { |
123 |
✓✓ |
163800 |
for (valIdx=0; valIdx<3; valIdx++) { |
124 |
|
140400 |
filter[ctx->radius](ctx->shape, |
125 |
|
70200 |
pvl->iv3 + valIdx*fd*fd*fd, |
126 |
|
70200 |
pvl->iv2 + valIdx*fd*fd, |
127 |
|
70200 |
pvl->iv1 + valIdx*fd, |
128 |
|
|
fw00, fw11, fw22, |
129 |
|
70200 |
vec + valIdx, jac + valIdx*3, hes + valIdx*9, |
130 |
|
70200 |
pvl->needD); |
131 |
|
|
} |
132 |
|
|
} else { |
133 |
✓✓ |
44100 |
for (valIdx=0; valIdx<3; valIdx++) { |
134 |
|
37800 |
gageScl3PFilterN(ctx->shape, fd, |
135 |
|
18900 |
pvl->iv3 + valIdx*fd*fd*fd, |
136 |
|
18900 |
pvl->iv2 + valIdx*fd*fd, |
137 |
|
18900 |
pvl->iv1 + valIdx*fd, |
138 |
|
|
fw00, fw11, fw22, |
139 |
|
18900 |
vec + valIdx, jac + valIdx*3, hes + valIdx*9, |
140 |
|
18900 |
pvl->needD); |
141 |
|
|
} |
142 |
|
|
} |
143 |
|
|
|
144 |
|
29700 |
return; |
145 |
|
29700 |
} |
146 |
|
|
|
147 |
|
|
void |
148 |
|
|
_gageVecAnswer(gageContext *ctx, gagePerVolume *pvl) { |
149 |
|
59400 |
char me[]="_gageVecAnswer"; |
150 |
|
29700 |
double cmag, tmpMat[9], mgevec[9], mgeval[3]; |
151 |
|
29700 |
double asym[9], tran[9], eval[3], tmpVec[3], norm; |
152 |
|
|
double *vecAns, *normAns, *jacAns, *strainAns, *somegaAns, |
153 |
|
|
*curlAns, *hesAns, *curlGradAns, |
154 |
|
|
*helGradAns, *dirHelDirAns, *curlnormgradAns; |
155 |
|
|
/* int asw; */ |
156 |
|
|
|
157 |
|
29700 |
vecAns = pvl->directAnswer[gageVecVector]; |
158 |
|
29700 |
normAns = pvl->directAnswer[gageVecNormalized]; |
159 |
|
29700 |
jacAns = pvl->directAnswer[gageVecJacobian]; |
160 |
|
29700 |
strainAns = pvl->directAnswer[gageVecStrain]; |
161 |
|
29700 |
somegaAns = pvl->directAnswer[gageVecSOmega]; |
162 |
|
29700 |
curlAns = pvl->directAnswer[gageVecCurl]; |
163 |
|
29700 |
hesAns = pvl->directAnswer[gageVecHessian]; |
164 |
|
29700 |
curlGradAns = pvl->directAnswer[gageVecCurlGradient]; |
165 |
|
29700 |
curlnormgradAns = pvl->directAnswer[gageVecCurlNormGrad]; |
166 |
|
29700 |
helGradAns = pvl->directAnswer[gageVecHelGradient]; |
167 |
|
29700 |
dirHelDirAns = pvl->directAnswer[gageVecDirHelDeriv]; |
168 |
|
|
|
169 |
✓✗ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecVector)) { |
170 |
|
|
/* done if doV */ |
171 |
✗✓ |
29700 |
if (ctx->verbose) { |
172 |
|
|
fprintf(stderr, "vec = "); |
173 |
|
|
ell_3v_print_d(stderr, vecAns); |
174 |
|
|
} |
175 |
|
|
} |
176 |
|
|
/* done if doV |
177 |
|
|
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecVector{0,1,2})) { |
178 |
|
|
} |
179 |
|
|
*/ |
180 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecLength)) { |
181 |
|
|
pvl->directAnswer[gageVecLength][0] = ELL_3V_LEN(vecAns); |
182 |
|
|
} |
183 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecNormalized)) { |
184 |
|
|
if (pvl->directAnswer[gageVecLength][0]) { |
185 |
|
|
ELL_3V_SCALE(normAns, 1.0/pvl->directAnswer[gageVecLength][0], vecAns); |
186 |
|
|
} else { |
187 |
|
|
ELL_3V_COPY(normAns, gageZeroNormal); |
188 |
|
|
} |
189 |
|
|
} |
190 |
✓✗ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecJacobian)) { |
191 |
|
|
/* done if doD1 */ |
192 |
|
|
/* |
193 |
|
|
0:dv_x/dx 1:dv_x/dy 2:dv_x/dz |
194 |
|
|
3:dv_y/dx 4:dv_y/dy 5:dv_y/dz |
195 |
|
|
6:dv_z/dx 7:dv_z/dy 8:dv_z/dz |
196 |
|
|
*/ |
197 |
✗✓ |
29700 |
if (ctx->verbose) { |
198 |
|
|
fprintf(stderr, "%s: jac = \n", me); |
199 |
|
|
ell_3m_print_d(stderr, jacAns); |
200 |
|
|
} |
201 |
|
|
} |
202 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecDivergence)) { |
203 |
|
|
pvl->directAnswer[gageVecDivergence][0] = jacAns[0] + jacAns[4] + jacAns[8]; |
204 |
|
|
if (ctx->verbose) { |
205 |
|
|
fprintf(stderr, "%s: div = %g + %g + %g = %g\n", me, |
206 |
|
|
jacAns[0], jacAns[4], jacAns[8], |
207 |
|
|
pvl->directAnswer[gageVecDivergence][0]); |
208 |
|
|
} |
209 |
|
|
} |
210 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecCurl)) { |
211 |
|
|
ELL_3V_SET(curlAns, |
212 |
|
|
jacAns[7] - jacAns[5], |
213 |
|
|
jacAns[2] - jacAns[6], |
214 |
|
|
jacAns[3] - jacAns[1]); |
215 |
|
|
} |
216 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecCurlNorm)) { |
217 |
|
|
pvl->directAnswer[gageVecCurlNorm][0] = ELL_3V_LEN(curlAns); |
218 |
|
|
} |
219 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecHelicity)) { |
220 |
|
|
pvl->directAnswer[gageVecHelicity][0] = |
221 |
|
|
ELL_3V_DOT(vecAns, curlAns); |
222 |
|
|
} |
223 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecNormHelicity)) { |
224 |
|
|
cmag = ELL_3V_LEN(curlAns); |
225 |
|
|
pvl->directAnswer[gageVecNormHelicity][0] = |
226 |
|
|
cmag ? ELL_3V_DOT(normAns, curlAns)/cmag : 0; |
227 |
|
|
} |
228 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecStrain)) { |
229 |
|
|
ELL_3M_TRANSPOSE(tran, jacAns); |
230 |
|
|
/* symmetric part */ |
231 |
|
|
ELL_3M_SCALE_ADD2(strainAns, 0.5, jacAns, 0.5, tran); |
232 |
|
|
/* nested these ifs to make dependency explicit to the compiler */ |
233 |
|
|
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecSOmega)) { |
234 |
|
|
/* antisymmetric part */ |
235 |
|
|
ELL_3M_SCALE_ADD2(asym, 0.5, jacAns, -0.5, tran); |
236 |
|
|
/* square symmetric part */ |
237 |
|
|
ELL_3M_MUL(tmpMat, strainAns, strainAns); |
238 |
|
|
ELL_3M_COPY(somegaAns, tmpMat); |
239 |
|
|
/* square antisymmetric part */ |
240 |
|
|
ELL_3M_MUL(tmpMat, asym, asym); |
241 |
|
|
/* sum of both */ |
242 |
|
|
ELL_3M_ADD2(somegaAns, somegaAns, tmpMat); |
243 |
|
|
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecLambda2)) { |
244 |
|
|
/* get eigenvalues in sorted order */ |
245 |
|
|
/* asw = */ ell_3m_eigenvalues_d(eval, somegaAns, AIR_TRUE); |
246 |
|
|
pvl->directAnswer[gageVecLambda2][0] = eval[1]; |
247 |
|
|
} |
248 |
|
|
} |
249 |
|
|
} |
250 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecImaginaryPart)) { |
251 |
|
|
pvl->directAnswer[gageVecImaginaryPart][0] = |
252 |
|
|
gage_imaginary_part_eigenvalues(jacAns); |
253 |
|
|
} |
254 |
|
|
/* 2nd order vector derivative continued */ |
255 |
✓✗ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecHessian)) { |
256 |
|
|
/* done if doD2 */ |
257 |
|
|
/* the ordering is induced by the scalar hessian computation : |
258 |
|
|
0:d2v_x/dxdx 1:d2v_x/dxdy 2:d2v_x/dxdz |
259 |
|
|
3:d2v_x/dydx 4:d2v_x/dydy 5:d2v_x/dydz |
260 |
|
|
6:d2v_x/dzdx 7:d2v_x/dzdy 8:d2v_x/dzdz |
261 |
|
|
9:d2v_y/dxdx [. . .] |
262 |
|
|
[. . .] |
263 |
|
|
24:dv2_z/dzdx 25:d2v_z/dzdy 26:d2v_z/dzdz |
264 |
|
|
*/ |
265 |
✗✓ |
29700 |
if (ctx->verbose) { |
266 |
|
|
fprintf(stderr, "%s: hes = \n", me); |
267 |
|
|
ell_3m_print_d(stderr, hesAns); /* ?? */ |
268 |
|
|
} |
269 |
|
|
} |
270 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecDivGradient)) { |
271 |
|
|
pvl->directAnswer[gageVecDivGradient][0] = hesAns[0] + hesAns[12] + hesAns[24]; |
272 |
|
|
pvl->directAnswer[gageVecDivGradient][1] = hesAns[1] + hesAns[13] + hesAns[25]; |
273 |
|
|
pvl->directAnswer[gageVecDivGradient][2] = hesAns[2] + hesAns[14] + hesAns[26]; |
274 |
|
|
} |
275 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecCurlGradient)) { |
276 |
|
|
pvl->directAnswer[gageVecCurlGradient][0] = hesAns[21]-hesAns[15]; |
277 |
|
|
pvl->directAnswer[gageVecCurlGradient][1] = hesAns[22]-hesAns[16]; |
278 |
|
|
pvl->directAnswer[gageVecCurlGradient][2] = hesAns[23]-hesAns[17]; |
279 |
|
|
pvl->directAnswer[gageVecCurlGradient][3] = hesAns[ 6]-hesAns[18]; |
280 |
|
|
pvl->directAnswer[gageVecCurlGradient][4] = hesAns[ 7]-hesAns[19]; |
281 |
|
|
pvl->directAnswer[gageVecCurlGradient][5] = hesAns[ 8]-hesAns[20]; |
282 |
|
|
pvl->directAnswer[gageVecCurlGradient][6] = hesAns[ 9]-hesAns[ 1]; |
283 |
|
|
pvl->directAnswer[gageVecCurlGradient][7] = hesAns[10]-hesAns[ 2]; |
284 |
|
|
pvl->directAnswer[gageVecCurlGradient][8] = hesAns[11]-hesAns[ 3]; |
285 |
|
|
} |
286 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecCurlNormGrad)) { |
287 |
|
|
norm = 1./ELL_3V_LEN(curlAns); |
288 |
|
|
|
289 |
|
|
tmpVec[0] = hesAns[21] - hesAns[15]; |
290 |
|
|
tmpVec[1] = hesAns[ 6] - hesAns[18]; |
291 |
|
|
tmpVec[2] = hesAns[ 9] - hesAns[ 3]; |
292 |
|
|
pvl->directAnswer[gageVecCurlNormGrad][0] = |
293 |
|
|
norm*ELL_3V_DOT(tmpVec, curlAns); |
294 |
|
|
|
295 |
|
|
tmpVec[0] = hesAns[22] - hesAns[16]; |
296 |
|
|
tmpVec[1] = hesAns[ 7] - hesAns[19]; |
297 |
|
|
tmpVec[2] = hesAns[10] - hesAns[ 4]; |
298 |
|
|
pvl->directAnswer[gageVecCurlNormGrad][1]= |
299 |
|
|
norm*ELL_3V_DOT(tmpVec, curlAns); |
300 |
|
|
|
301 |
|
|
tmpVec[0] = hesAns[23] - hesAns[17]; |
302 |
|
|
tmpVec[1] = hesAns[ 8] - hesAns[20]; |
303 |
|
|
tmpVec[2] = hesAns[11] - hesAns[ 5]; |
304 |
|
|
pvl->directAnswer[gageVecCurlNormGrad][2]= |
305 |
|
|
norm*ELL_3V_DOT(tmpVec, curlAns); |
306 |
|
|
} |
307 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecNCurlNormGrad)) { |
308 |
|
|
norm = 1./ELL_3V_LEN(curlnormgradAns); |
309 |
|
|
ELL_3V_SCALE(pvl->directAnswer[gageVecNCurlNormGrad], |
310 |
|
|
norm, pvl->directAnswer[gageVecCurlNormGrad]); |
311 |
|
|
} |
312 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecHelGradient)) { |
313 |
|
|
pvl->directAnswer[gageVecHelGradient][0] = |
314 |
|
|
jacAns[0]*curlAns[0]+ |
315 |
|
|
jacAns[3]*curlAns[1]+ |
316 |
|
|
jacAns[6]*curlAns[2]+ |
317 |
|
|
curlGradAns[0]*vecAns[0]+ |
318 |
|
|
curlGradAns[3]*vecAns[1]+ |
319 |
|
|
curlGradAns[6]*vecAns[2]; |
320 |
|
|
pvl->directAnswer[gageVecHelGradient][1] = |
321 |
|
|
jacAns[1]*curlAns[0]+ |
322 |
|
|
jacAns[4]*curlAns[1]+ |
323 |
|
|
jacAns[7]*curlAns[2]+ |
324 |
|
|
curlGradAns[1]*vecAns[0]+ |
325 |
|
|
curlGradAns[4]*vecAns[1]+ |
326 |
|
|
curlGradAns[7]*vecAns[2]; |
327 |
|
|
pvl->directAnswer[gageVecHelGradient][0] = |
328 |
|
|
jacAns[2]*curlAns[0]+ |
329 |
|
|
jacAns[5]*curlAns[1]+ |
330 |
|
|
jacAns[8]*curlAns[2]+ |
331 |
|
|
curlGradAns[2]*vecAns[0]+ |
332 |
|
|
curlGradAns[5]*vecAns[1]+ |
333 |
|
|
curlGradAns[8]*vecAns[2]; |
334 |
|
|
} |
335 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecDirHelDeriv)) { |
336 |
|
|
pvl->directAnswer[gageVecDirHelDeriv][0] = |
337 |
|
|
ELL_3V_DOT(normAns, helGradAns); |
338 |
|
|
} |
339 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecProjHelGradient)) { |
340 |
|
|
pvl->directAnswer[gageVecDirHelDeriv][0] = |
341 |
|
|
helGradAns[0]-dirHelDirAns[0]*normAns[0]; |
342 |
|
|
pvl->directAnswer[gageVecDirHelDeriv][1] = |
343 |
|
|
helGradAns[1]-dirHelDirAns[0]*normAns[1]; |
344 |
|
|
pvl->directAnswer[gageVecDirHelDeriv][2] = |
345 |
|
|
helGradAns[2]-dirHelDirAns[0]*normAns[2]; |
346 |
|
|
} |
347 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecGradient0)) { |
348 |
|
|
ELL_3V_SET(pvl->directAnswer[gageVecGradient0], |
349 |
|
|
jacAns[0], |
350 |
|
|
jacAns[1], |
351 |
|
|
jacAns[2]); |
352 |
|
|
} |
353 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecGradient1)) { |
354 |
|
|
ELL_3V_SET(pvl->directAnswer[gageVecGradient1], |
355 |
|
|
jacAns[3], |
356 |
|
|
jacAns[4], |
357 |
|
|
jacAns[5]); |
358 |
|
|
} |
359 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecGradient2)) { |
360 |
|
|
ELL_3V_SET(pvl->directAnswer[gageVecGradient2], |
361 |
|
|
jacAns[6], |
362 |
|
|
jacAns[7], |
363 |
|
|
jacAns[8]); |
364 |
|
|
} |
365 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecMultiGrad)) { |
366 |
|
|
ELL_3M_IDENTITY_SET(pvl->directAnswer[gageVecMultiGrad]); |
367 |
|
|
ELL_3MV_OUTER_INCR(pvl->directAnswer[gageVecMultiGrad], |
368 |
|
|
pvl->directAnswer[gageVecGradient0], |
369 |
|
|
pvl->directAnswer[gageVecGradient0]); |
370 |
|
|
ELL_3MV_OUTER_INCR(pvl->directAnswer[gageVecMultiGrad], |
371 |
|
|
pvl->directAnswer[gageVecGradient1], |
372 |
|
|
pvl->directAnswer[gageVecGradient1]); |
373 |
|
|
ELL_3MV_OUTER_INCR(pvl->directAnswer[gageVecMultiGrad], |
374 |
|
|
pvl->directAnswer[gageVecGradient2], |
375 |
|
|
pvl->directAnswer[gageVecGradient2]); |
376 |
|
|
} |
377 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecMGFrob)) { |
378 |
|
|
pvl->directAnswer[gageVecMGFrob][0] |
379 |
|
|
= ELL_3M_FROB(pvl->directAnswer[gageVecMultiGrad]); |
380 |
|
|
} |
381 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecMGEval)) { |
382 |
|
|
ELL_3M_COPY(tmpMat, pvl->directAnswer[gageVecMultiGrad]); |
383 |
|
|
/* HEY: look at the return value for root multiplicity? */ |
384 |
|
|
ell_3m_eigensolve_d(mgeval, mgevec, tmpMat, AIR_TRUE); |
385 |
|
|
ELL_3V_COPY(pvl->directAnswer[gageVecMGEval], mgeval); |
386 |
|
|
} |
387 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecMGEvec)) { |
388 |
|
|
ELL_3M_COPY(pvl->directAnswer[gageVecMGEvec], mgevec); |
389 |
|
|
} |
390 |
|
|
|
391 |
|
|
return; |
392 |
|
29700 |
} |
393 |
|
|
|
394 |
|
|
const char * |
395 |
|
|
_gageVecStr[] = { |
396 |
|
|
"(unknown gageVec)", |
397 |
|
|
"vector", |
398 |
|
|
"vector0", |
399 |
|
|
"vector1", |
400 |
|
|
"vector2", |
401 |
|
|
"length", |
402 |
|
|
"normalized", |
403 |
|
|
"Jacobian", |
404 |
|
|
"strain", |
405 |
|
|
"divergence", |
406 |
|
|
"curl", |
407 |
|
|
"curl norm", |
408 |
|
|
"helicity", |
409 |
|
|
"normalized helicity", |
410 |
|
|
"SOmega", |
411 |
|
|
"lambda2", |
412 |
|
|
"imag", |
413 |
|
|
"vector hessian", |
414 |
|
|
"div gradient", |
415 |
|
|
"curl gradient", |
416 |
|
|
"curl norm gradient", |
417 |
|
|
"normalized curl norm gradient", |
418 |
|
|
"helicity gradient", |
419 |
|
|
"directional helicity derivative", |
420 |
|
|
"projected helicity gradient", |
421 |
|
|
"gradient0", |
422 |
|
|
"gradient1", |
423 |
|
|
"gradient2", |
424 |
|
|
"multigrad", |
425 |
|
|
"frob(multigrad)", |
426 |
|
|
"multigrad eigenvalues", |
427 |
|
|
"multigrad eigenvectors", |
428 |
|
|
}; |
429 |
|
|
|
430 |
|
|
const char * |
431 |
|
|
_gageVecDesc[] = { |
432 |
|
|
"unknown gageVec query", |
433 |
|
|
"component-wise-interpolated vector", |
434 |
|
|
"vector[0]", |
435 |
|
|
"vector[1]", |
436 |
|
|
"vector[2]", |
437 |
|
|
"length of vector", |
438 |
|
|
"normalized vector", |
439 |
|
|
"3x3 Jacobian", |
440 |
|
|
"rate-of-strain tensor", |
441 |
|
|
"divergence", |
442 |
|
|
"curl", |
443 |
|
|
"curl magnitude", |
444 |
|
|
"helicity: dot(vector,curl)", |
445 |
|
|
"normalized helicity", |
446 |
|
|
"S squared plus Omega squared for vortex characterization", |
447 |
|
|
"lambda2 value of SOmega", |
448 |
|
|
"imaginary part of complex-conjugate eigenvalues of Jacobian", |
449 |
|
|
"3x3x3 second-order vector derivative", |
450 |
|
|
"gradient of divergence", |
451 |
|
|
"3x3 derivative of curl", |
452 |
|
|
"gradient of curl norm", |
453 |
|
|
"normalized gradient of curl norm", |
454 |
|
|
"gradient of helicity", |
455 |
|
|
"directional derivative of helicity along flow", |
456 |
|
|
"projection of the helicity gradient onto plane orthogonal to flow", |
457 |
|
|
"gradient of 1st component of vector", |
458 |
|
|
"gradient of 2nd component of vector", |
459 |
|
|
"gradient of 3rd component of vector", |
460 |
|
|
"multi-gradient: sum of outer products of gradients", |
461 |
|
|
"frob norm of multi-gradient", |
462 |
|
|
"eigenvalues of multi-gradient", |
463 |
|
|
"eigenvectors of multi-gradient" |
464 |
|
|
}; |
465 |
|
|
|
466 |
|
|
const int |
467 |
|
|
_gageVecVal[] = { |
468 |
|
|
gageVecUnknown, |
469 |
|
|
gageVecVector, |
470 |
|
|
gageVecVector0, |
471 |
|
|
gageVecVector1, |
472 |
|
|
gageVecVector2, |
473 |
|
|
gageVecLength, |
474 |
|
|
gageVecNormalized, |
475 |
|
|
gageVecJacobian, |
476 |
|
|
gageVecStrain, |
477 |
|
|
gageVecDivergence, |
478 |
|
|
gageVecCurl, |
479 |
|
|
gageVecCurlNorm, |
480 |
|
|
gageVecHelicity, |
481 |
|
|
gageVecNormHelicity, |
482 |
|
|
gageVecSOmega, |
483 |
|
|
gageVecLambda2, |
484 |
|
|
gageVecImaginaryPart, |
485 |
|
|
gageVecHessian, |
486 |
|
|
gageVecDivGradient, |
487 |
|
|
gageVecCurlGradient, |
488 |
|
|
gageVecCurlNormGrad, |
489 |
|
|
gageVecNCurlNormGrad, |
490 |
|
|
gageVecHelGradient, |
491 |
|
|
gageVecDirHelDeriv, |
492 |
|
|
gageVecProjHelGradient, |
493 |
|
|
gageVecGradient0, |
494 |
|
|
gageVecGradient1, |
495 |
|
|
gageVecGradient2, |
496 |
|
|
gageVecMultiGrad, |
497 |
|
|
gageVecMGFrob, |
498 |
|
|
gageVecMGEval, |
499 |
|
|
gageVecMGEvec, |
500 |
|
|
}; |
501 |
|
|
|
502 |
|
|
#define GV_V gageVecVector |
503 |
|
|
#define GV_V0 gageVecVector0 |
504 |
|
|
#define GV_V1 gageVecVector1 |
505 |
|
|
#define GV_V2 gageVecVector2 |
506 |
|
|
#define GV_L gageVecLength |
507 |
|
|
#define GV_N gageVecNormalized |
508 |
|
|
#define GV_J gageVecJacobian |
509 |
|
|
#define GV_S gageVecStrain |
510 |
|
|
#define GV_D gageVecDivergence |
511 |
|
|
#define GV_C gageVecCurl |
512 |
|
|
#define GV_CM gageVecCurlNorm |
513 |
|
|
#define GV_H gageVecHelicity |
514 |
|
|
#define GV_NH gageVecNormHelicity |
515 |
|
|
#define GV_SO gageVecSOmega |
516 |
|
|
#define GV_LB gageVecLambda2 |
517 |
|
|
#define GV_IM gageVecImaginaryPart |
518 |
|
|
#define GV_VH gageVecHessian |
519 |
|
|
#define GV_DG gageVecDivGradient |
520 |
|
|
#define GV_CG gageVecCurlGradient |
521 |
|
|
#define GV_CNG gageVecCurlNormGrad |
522 |
|
|
#define GV_NCG gageVecNCurlNormGrad |
523 |
|
|
#define GV_HG gageVecHelGradient |
524 |
|
|
#define GV_DH gageVecDirHelDeriv |
525 |
|
|
#define GV_PH gageVecProjHelGradient |
526 |
|
|
#define GV_G0 gageVecGradient0 |
527 |
|
|
#define GV_G1 gageVecGradient1 |
528 |
|
|
#define GV_G2 gageVecGradient2 |
529 |
|
|
#define GV_MG gageVecMultiGrad |
530 |
|
|
#define GV_MF gageVecMGFrob |
531 |
|
|
#define GV_ML gageVecMGEval |
532 |
|
|
#define GV_MC gageVecMGEvec |
533 |
|
|
|
534 |
|
|
const char * |
535 |
|
|
_gageVecStrEqv[] = { |
536 |
|
|
"v", "vector", "vec", |
537 |
|
|
"v0", "vector0", "vec0", |
538 |
|
|
"v1", "vector1", "vec1", |
539 |
|
|
"v2", "vector2", "vec2", |
540 |
|
|
"l", "length", "len", |
541 |
|
|
"n", "normalized", "normalized vector", |
542 |
|
|
"jacobian", "jac", "j", |
543 |
|
|
"strain", "S", |
544 |
|
|
"divergence", "div", "d", |
545 |
|
|
"curl", "c", |
546 |
|
|
"curlnorm", "curl norm", "curl magnitude", "cm", |
547 |
|
|
"h", "hel", "hell", "helicity", |
548 |
|
|
"nh", "nhel", "normhel", "normhell", "normalized helicity", |
549 |
|
|
"SOmega", |
550 |
|
|
"lbda2", "lambda2", |
551 |
|
|
"imag", "imagpart", |
552 |
|
|
"vh", "vhes", "vhessian", "vector hessian", |
553 |
|
|
"dg", "divgrad", "div gradient", |
554 |
|
|
"cg", "curlgrad", "curlg", "curljac", "curl gradient", |
555 |
|
|
"cng", "curl norm gradient", |
556 |
|
|
"ncng", "norm curl norm gradient", "normalized curl norm gradient", |
557 |
|
|
"hg", "helg", "helgrad", "helicity gradient", |
558 |
|
|
"dirhelderiv", "dhd", "ddh", "directional helicity derivative", |
559 |
|
|
"phg", "projhel", "projhelgrad", "projected helicity gradient", |
560 |
|
|
"g0", "grad0", "gradient0", |
561 |
|
|
"g1", "grad1", "gradient1", |
562 |
|
|
"g2", "grad2", "gradient2", |
563 |
|
|
"mg", "multigrad", |
564 |
|
|
"mgfrob", "frob(multigrad)", |
565 |
|
|
"mgeval", "mg eval", "multigrad eigenvalues", |
566 |
|
|
"mgevec", "mg evec", "multigrad eigenvectors", |
567 |
|
|
"" |
568 |
|
|
}; |
569 |
|
|
|
570 |
|
|
const int |
571 |
|
|
_gageVecValEqv[] = { |
572 |
|
|
GV_V, GV_V, GV_V, |
573 |
|
|
GV_V0, GV_V0, GV_V0, |
574 |
|
|
GV_V1, GV_V1, GV_V1, |
575 |
|
|
GV_V2, GV_V2, GV_V2, |
576 |
|
|
GV_L, GV_L, GV_L, |
577 |
|
|
GV_N, GV_N, GV_N, |
578 |
|
|
GV_J, GV_J, GV_J, |
579 |
|
|
GV_S, GV_S, |
580 |
|
|
GV_D, GV_D, GV_D, |
581 |
|
|
GV_C, GV_C, |
582 |
|
|
GV_CM, GV_CM, GV_CM, GV_CM, |
583 |
|
|
GV_H, GV_H, GV_H, GV_H, |
584 |
|
|
GV_NH, GV_NH, GV_NH, GV_NH, GV_NH, |
585 |
|
|
GV_SO, |
586 |
|
|
GV_LB, GV_LB, |
587 |
|
|
GV_IM, GV_IM, |
588 |
|
|
GV_VH, GV_VH, GV_VH, GV_VH, |
589 |
|
|
GV_DG, GV_DG, GV_DG, |
590 |
|
|
GV_CG, GV_CG, GV_CG, GV_CG, GV_CG, |
591 |
|
|
GV_CNG, GV_CNG, |
592 |
|
|
GV_NCG, GV_NCG, GV_NCG, |
593 |
|
|
GV_HG, GV_HG, GV_HG, GV_HG, |
594 |
|
|
GV_DH, GV_DH, GV_DH, GV_DH, |
595 |
|
|
GV_PH, GV_PH, GV_PH, GV_PH, |
596 |
|
|
GV_G0, GV_G0, GV_G0, |
597 |
|
|
GV_G1, GV_G1, GV_G1, |
598 |
|
|
GV_G2, GV_G2, GV_G2, |
599 |
|
|
GV_MG, GV_MG, |
600 |
|
|
GV_MF, GV_MF, |
601 |
|
|
GV_ML, GV_ML, GV_ML, |
602 |
|
|
GV_MC, GV_MC, GV_MC |
603 |
|
|
}; |
604 |
|
|
|
605 |
|
|
const airEnum |
606 |
|
|
_gageVec = { |
607 |
|
|
"gageVec", |
608 |
|
|
GAGE_VEC_ITEM_MAX, |
609 |
|
|
_gageVecStr, _gageVecVal, |
610 |
|
|
_gageVecDesc, |
611 |
|
|
_gageVecStrEqv, _gageVecValEqv, |
612 |
|
|
AIR_FALSE |
613 |
|
|
}; |
614 |
|
|
const airEnum *const |
615 |
|
|
gageVec = &_gageVec; |
616 |
|
|
|
617 |
|
|
gageKind |
618 |
|
|
_gageKindVec = { |
619 |
|
|
AIR_FALSE, /* statically allocated */ |
620 |
|
|
"vector", |
621 |
|
|
&_gageVec, |
622 |
|
|
1, /* baseDim */ |
623 |
|
|
3, /* valLen */ |
624 |
|
|
GAGE_VEC_ITEM_MAX, |
625 |
|
|
_gageVecTable, |
626 |
|
|
_gageVecIv3Print, |
627 |
|
|
_gageVecFilter, |
628 |
|
|
_gageVecAnswer, |
629 |
|
|
NULL, NULL, NULL, NULL, |
630 |
|
|
NULL |
631 |
|
|
}; |
632 |
|
|
gageKind *const |
633 |
|
|
gageKindVec = &_gageKindVec; |
634 |
|
|
|