Menu
▾
▴
[Teem-commit] teem/src/gage gage.h,1.64,1.65 vecGage.c,1.19,1.20
|
From: xav <xa...@us...> - 2004-04-11 02:45:16
|
Update of /cvsroot/teem/teem/src/gage In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv7404 Modified Files: gage.h vecGage.c Log Message: added new functions in vecGage.c, especially useful in the context of vortical flows. Index: gage.h =================================================================== RCS file: /cvsroot/teem/teem/src/gage/gage.h,v retrieving revision 1.64 retrieving revision 1.65 diff -C2 -d -r1.64 -r1.65 *** gage.h 1 Apr 2004 22:37:26 -0000 1.64 --- gage.h 11 Apr 2004 02:31:41 -0000 1.65 *************** *** 287,294 **** gageVecDivergence, /* 4: "d", divergence (based on Jacobian): GT[1] */ gageVecCurl, /* 5: "c", curl (based on Jacobian): GT[3] */ ! gageVecHelicity, /* 6: "h", helicity: vec . curl: GT[1] */ ! gageVecNormHelicity,/* 7: "nh", normalized helicity: GT[1] */ ! gageVecLambda2, /* 8: "lmabda2", lambda2 criterion: GT[1] */ ! gageVecHessian, /* 9: "vh", second-order derivative: GT[27] 0:d2v_x/dxdx 1:d2v_x/dxdy 2:d2v_x/dxdz 3:d2v_x/dydx 4:d2v_x/dydy 5:d2v_x/dydz --- 287,296 ---- gageVecDivergence, /* 4: "d", divergence (based on Jacobian): GT[1] */ gageVecCurl, /* 5: "c", curl (based on Jacobian): GT[3] */ ! gageVecCurlNorm, /* 6: "cm", curl magnitude: GT[1] */ ! gageVecHelicity, /* 7: "h", helicity: vec . curl: GT[1] */ ! gageVecNormHelicity,/* 8: "nh", normalized helicity: GT[1] */ ! gageVecLambda2, /* 9: "lmbda2", lambda2 criterion: GT[1] */ ! gageVecImaginaryPart,/* 10: "imag", imag. part of jacobian's eigenv: GT[1] */ ! gageVecHessian, /* 11: "vh", second-order derivative: GT[27] 0:d2v_x/dxdx 1:d2v_x/dxdy 2:d2v_x/dxdz 3:d2v_x/dydx 4:d2v_x/dydy 5:d2v_x/dydz *************** *** 297,316 **** [...] 24:dv2_z/dzdx 25:d2v_z/dzdy 26:d2v_z/dzdz */ ! gageVecCurlGradient,/*10: "cg", curl gradient: GT[9] */ ! gageVecCurlNormGrad,/*11: "cng", curl norm gradient: GT[3] */ ! gageVecNCurlNormGrad, /* 12: "ncng", normalized curl norm gradient: GT[3] */ ! gageVecHelGradient,/* 13: "hg", helicity gradient: GT[3] */ ! gageVecDirHelDeriv,/* 14: "dhd", directional derivative of helicity: GT[1] */ ! gageVecProjHelGradient, /* 15: "phg", projected helicity gradient: GT[3] */ ! gageVecGradient0, /* 16: "g0", gradient of 1st component of vector: GT[3] */ ! gageVecGradient1, /* 17: "g1", gradient of 2nd component of vector: GT[3] */ ! gageVecGradient2, /* 18: "g2", gradient of 3rd component of vector: GT[3] */ ! gageVecMultiGrad, /* 19: "mg", sum of outer products of gradients: GT[9] */ ! gageVecMGFrob, /* 20: "mgfrob", frob norm of multi-gradient: GT[1] */ ! gageVecMGEval, /* 21: "mgeval", eigenvalues of multi-gradient: GT[3] */ ! gageVecMGEvec, /* 22: "mgevec", eigenvectors of multi-gradient: GT[9] */ gageVecLast }; ! #define GAGE_VEC_ITEM_MAX 22 struct gageKind_t; /* dumb forward declaraction, ignore */ --- 299,319 ---- [...] 24:dv2_z/dzdx 25:d2v_z/dzdy 26:d2v_z/dzdz */ ! gageVecDivGradient,/* 12: "dg", divergence gradient: GT[3] */ ! gageVecCurlGradient,/*13: "cg", curl gradient: GT[9] */ ! gageVecCurlNormGrad,/*14: "cng", curl norm gradient: GT[3] */ ! gageVecNCurlNormGrad, /* 15: "ncng", normalized curl norm gradient: GT[3] */ ! gageVecHelGradient,/* 16: "hg", helicity gradient: GT[3] */ ! gageVecDirHelDeriv,/* 17: "dhd", directional derivative of helicity: GT[1] */ ! gageVecProjHelGradient, /* 18: "phg", projected helicity gradient: GT[3] */ ! gageVecGradient0, /* 19: "g0", gradient of 1st component of vector: GT[3] */ ! gageVecGradient1, /* 20: "g1", gradient of 2nd component of vector: GT[3] */ ! gageVecGradient2, /* 21: "g2", gradient of 3rd component of vector: GT[3] */ ! gageVecMultiGrad, /* 22: "mg", sum of outer products of gradients: GT[9] */ ! gageVecMGFrob, /* 23: "mgfrob", frob norm of multi-gradient: GT[1] */ ! gageVecMGEval, /* 24: "mgeval", eigenvalues of multi-gradient: GT[3] */ ! gageVecMGEvec, /* 25: "mgevec", eigenvectors of multi-gradient: GT[9] */ gageVecLast }; ! #define GAGE_VEC_ITEM_MAX 25 struct gageKind_t; /* dumb forward declaraction, ignore */ Index: vecGage.c =================================================================== RCS file: /cvsroot/teem/teem/src/gage/vecGage.c,v retrieving revision 1.19 retrieving revision 1.20 diff -C2 -d -r1.19 -r1.20 *** vecGage.c 1 Apr 2004 22:37:27 -0000 1.19 --- vecGage.c 11 Apr 2004 02:31:41 -0000 1.20 *************** *** 21,24 **** --- 21,61 ---- #include "privateGage.h" + /* + * highly inefficient computation of the imaginary part of complex + * conjugate eigenvalues of a 3x3 non-symmetric matrix + */ + double + gage_imaginary_part_eigenvalues( gage_t *M ) + { + double A, B, C, scale, frob, m[9], _eval[3]; + double beta, gamma; + int roots; + + frob = ELL_3M_FROB(M); + scale = frob > 10 ? 10.0/frob : 1.0; + ELL_3M_SCALE(m, scale, M); + /* + ** from gordon with mathematica; these are the coefficients of the + ** cubic polynomial in x: det(x*I - M). The full cubic is + ** x^3 + A*x^2 + B*x + C. + */ + A = -m[0] - m[4] - m[8]; + B = m[0]*m[4] - m[3]*m[1] + + m[0]*m[8] - m[6]*m[2] + + m[4]*m[8] - m[7]*m[5]; + C = (m[6]*m[4] - m[3]*m[7])*m[2] + + (m[0]*m[7] - m[6]*m[1])*m[5] + + (m[3]*m[1] - m[0]*m[4])*m[8]; + roots = ell_cubic(_eval, A, B, C, AIR_TRUE); + if ( roots != ell_cubic_root_single ) + return 0.; + + /* 2 complex conjuguate eigenvalues */ + beta = A + _eval[0]; + gamma = -C/_eval[0]; + return sqrt( 4.*gamma - beta*beta ); + } + + gageItemEntry _gageVecTable[GAGE_VEC_ITEM_MAX+1] = { *************** *** 30,37 **** --- 67,77 ---- {gageVecDivergence, 1, 1, {gageVecJacobian, -1, -1, -1, -1}, -1, -1}, {gageVecCurl, 3, 1, {gageVecJacobian, -1, -1, -1, -1}, -1, -1}, + {gageVecCurlNorm, 1, 1, {gageVecCurl, -1, -1, -1, -1}, -1, -1}, {gageVecHelicity, 1, 1, {gageVecVector, gageVecCurl, -1, -1, -1}, -1, -1}, {gageVecNormHelicity, 1, 1, {gageVecNormalized, gageVecCurl, -1, -1, -1}, -1, -1}, {gageVecLambda2, 1, 1, {gageVecJacobian, -1, -1, -1, -1}, -1, -1}, + {gageVecImaginaryPart, 1, 1, {gageVecJacobian, -1, -1, -1, -1}, -1, -1}, {gageVecHessian, 27, 2, {-1, -1, -1, -1, -1}, -1, -1}, + {gageVecDivGradient, 3, 1, {gageVecHessian, -1, -1, -1, -1}, -1, -1}, {gageVecCurlGradient, 9, 2, {gageVecHessian, -1, -1, -1, -1}, -1, -1}, {gageVecCurlNormGrad, 3, 2, {gageVecHessian, gageVecCurl, -1, -1, -1}, -1, -1}, *************** *** 59,65 **** fd = GAGE_FD(ctx); ! vec = pvl->directAnswer[gageVecVector]; ! jac = pvl->directAnswer[gageVecJacobian]; ! hes = pvl->directAnswer[gageVecHessian]; if (!ctx->parm.k3pack) { fprintf(stderr, "!%s: sorry, 6pack filtering not implemented\n", me); --- 99,105 ---- fd = GAGE_FD(ctx); ! vec = pvl->directAnswer[gageVecVector]; ! jac = pvl->directAnswer[gageVecJacobian]; ! hes = pvl->directAnswer[gageVecHessian]; if (!ctx->parm.k3pack) { fprintf(stderr, "!%s: sorry, 6pack filtering not implemented\n", me); *************** *** 165,168 **** --- 205,212 ---- jacAns[3] - jacAns[1]); } + if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecCurlNorm)) { + pvl->directAnswer[gageVecCurlNorm][0] = + ELL_3V_LEN( curlAns ); + } if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecHelicity)) { pvl->directAnswer[gageVecHelicity][0] = *************** *** 191,194 **** --- 235,242 ---- pvl->directAnswer[gageVecLambda2][0] = eval[1]; } + if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecImaginaryPart)) { + pvl->directAnswer[gageVecImaginaryPart][0] = + gage_imaginary_part_eigenvalues( jacAns ); + } /* 2nd order vector derivative continued */ if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecHessian)) { *************** *** 207,210 **** --- 255,263 ---- } } + if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecDivGradient)) { + pvl->directAnswer[gageVecDivGradient][0] = hesAns[0] + hesAns[12] + hesAns[24]; + pvl->directAnswer[gageVecDivGradient][1] = hesAns[1] + hesAns[13] + hesAns[25]; + pvl->directAnswer[gageVecDivGradient][2] = hesAns[2] + hesAns[14] + hesAns[26]; + } if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecCurlGradient)) { pvl->directAnswer[gageVecCurlGradient][0] = hesAns[21]-hesAns[15]; *************** *** 335,342 **** --- 388,397 ---- "divergence", "curl", + "curl norm", "helicity", "normalized helicity", "lambda2", "vector hessian", + "div gradient", "curl gradient", "curl norm gradient", *************** *** 363,370 **** --- 418,427 ---- "divergence", "curl", + "curl magnitude", "helicity: dot(vector,curl)", "normalized helicity", "lambda2 value for vortex characterization", "3x3x3 second-order vector derivative", + "gradient of divergence", "3x3 derivative of curl", "gradient of curl norm", *************** *** 391,398 **** --- 448,458 ---- gageVecDivergence, gageVecCurl, + gageVecCurlNorm, gageVecHelicity, gageVecNormHelicity, gageVecLambda2, + gageVecImaginaryPart, gageVecHessian, + gageVecDivGradient, gageVecCurlGradient, gageVecCurlNormGrad, *************** *** 416,423 **** --- 476,486 ---- #define GV_D gageVecDivergence #define GV_C gageVecCurl + #define GV_CM gageVecCurlNorm #define GV_H gageVecHelicity #define GV_NH gageVecNormHelicity #define GV_LB gageVecLambda2 + #define GV_IM gageVecImaginaryPart #define GV_VH gageVecHessian + #define GV_DG gageVecDivGradient #define GV_CG gageVecCurlGradient #define GV_CNG gageVecCurlNormGrad *************** *** 442,449 **** --- 505,515 ---- "divergence", "div", "d", "curl", "c", + "curlnorm", "curl magnitude", "cm", "h", "hel", "hell", "nh", "nhel", "normhel", "normhell", "lbda2", "lambda2", + "imag", "imagpart", "vh", "vhes", "vhessian", "vector hessian", + "dg", "divgrad", "cg", "curlgrad", "curlg", "curljac", "curl gradient", "cng", "curl norm gradient", *************** *** 470,477 **** --- 536,546 ---- GV_D, GV_D, GV_D, GV_C, GV_C, + GV_CM, GV_CM, GV_CM, GV_H, GV_H, GV_H, GV_NH, GV_NH, GV_NH, GV_NH, GV_LB, GV_LB, + GV_IM, GV_IM, GV_VH, GV_VH, GV_VH, GV_VH, + GV_DG, GV_DG, GV_CG, GV_CG, GV_CG, GV_CG, GV_CG, GV_CNG, GV_CNG, |
