Menu
▾
▴
[Teem-commit] teem/src/gage gage.h,1.62,1.63 miscGage.c,1.11,1.12 vecGage.c,1.17,1.18
|
From: Gordon K. <kin...@us...> - 2004-04-01 07:19:31
|
Update of /cvsroot/teem/teem/src/gage In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv7913 Modified Files: gage.h miscGage.c vecGage.c Log Message: passing vecGage baton to xavier Index: gage.h =================================================================== RCS file: /cvsroot/teem/teem/src/gage/gage.h,v retrieving revision 1.62 retrieving revision 1.63 diff -C2 -d -r1.62 -r1.63 *** gage.h 28 Feb 2004 17:01:53 -0000 1.62 --- gage.h 1 Apr 2004 07:07:35 -0000 1.63 *************** *** 274,278 **** ** ** The strings gives one of the gageVec airEnum identifiers, and GT[x] ! ** says how many scalars are associated with this scalar. */ enum { --- 274,278 ---- ** ** The strings gives one of the gageVec airEnum identifiers, and GT[x] ! ** says how many scalars are associated with this value. */ enum { *************** *** 287,300 **** gageVecDivergence, /* 4: "d", divergence (based on Jacobian): GT[1] */ gageVecCurl, /* 5: "c", curl (based on Jacobian): GT[3] */ ! gageVecGradient0, /* 6: "g0", gradient of 1st component of vector: GT[3] */ ! gageVecGradient1, /* 7: "g1", gradient of 2nd component of vector: GT[3] */ ! gageVecGradient2, /* 8: "g2", gradient of 3rd component of vector: GT[3] */ ! gageVecMultiGrad, /* 9: "mg", sum of outer products of gradients: GT[9] */ ! gageVecMGFrob, /* 10: "mgfrob", frob norm of multi-gradient: GT[1] */ ! gageVecMGEval, /* 11: "mgeval", eigenvalues of multi-gradient: GT[3] */ ! gageVecMGEvec, /* 12: "mgevec", eigenvectors of multi-gradient: GT[9] */ gageVecLast }; ! #define GAGE_VEC_ITEM_MAX 12 struct gageKind_t; /* dumb forward declaraction, ignore */ --- 287,308 ---- 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] */ ! gageVecCurlGradient,/*10: "cg", curl gradient: GT[9] */ ! gageVecHelGradient,/* 11: "hg", helicity gradient: GT[3] */ ! gageVecDirHelDeriv,/* 12: "dhd", directional derivative of helicity: GT[1] */ ! gageVecProjHelGradient,/* 13: "phg", projected helicity gradient: GT[3] */ ! gageVecGradient0, /* 14: "g0", gradient of 1st component of vector: GT[3] */ ! gageVecGradient1, /* 15: "g1", gradient of 2nd component of vector: GT[3] */ ! gageVecGradient2, /* 16: "g2", gradient of 3rd component of vector: GT[3] */ ! gageVecMultiGrad, /* 17: "mg", sum of outer products of gradients: GT[9] */ ! gageVecMGFrob, /* 18: "mgfrob", frob norm of multi-gradient: GT[1] */ ! gageVecMGEval, /* 19: "mgeval", eigenvalues of multi-gradient: GT[3] */ ! gageVecMGEvec, /* 20: "mgevec", eigenvectors of multi-gradient: GT[9] */ gageVecLast }; ! #define GAGE_VEC_ITEM_MAX 20 struct gageKind_t; /* dumb forward declaraction, ignore */ Index: vecGage.c =================================================================== RCS file: /cvsroot/teem/teem/src/gage/vecGage.c,v retrieving revision 1.17 retrieving revision 1.18 diff -C2 -d -r1.17 -r1.18 *** vecGage.c 30 Mar 2004 04:26:01 -0000 1.17 --- vecGage.c 1 Apr 2004 07:07:35 -0000 1.18 *************** *** 30,33 **** --- 30,43 ---- {gageVecDivergence, 1, 1, {gageVecJacobian, -1, -1, -1, -1}, -1, -1}, {gageVecCurl, 3, 1, {gageVecJacobian, -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}, + {gageVecHessian, 27, 2, {-1, -1, -1, -1, -1}, -1, -1}, + {gageVecCurlGradient, 9, 2, {gageVecHessian, -1, -1, -1, -1}, -1, -1}, + {gageVecHelGradient, 3, 2, {gageVecVector, gageVecJacobian, gageVecCurl, + gageVecCurlGradient, -1}, -1, -1}, + {gageVecDirHelDeriv, 1, 2, {gageVecNormalized, gageVecHelGradient, -1, -1, -1}, -1, -1}, + {gageVecProjHelGradient, 3,2, {gageVecNormalized, gageVecHelGradient, gageVecDirHelDeriv, + -1, -1}, -1, -1}, /* HEY: these should change to sub-items!!! */ {gageVecGradient0, 3, 1, {gageVecJacobian, -1, -1, -1, -1}, -1, -1}, *************** *** 43,47 **** _gageVecFilter (gageContext *ctx, gagePerVolume *pvl) { char me[]="_gageVecFilter"; ! gage_t *fw00, *fw11, *fw22, *vec, *jac; int fd; --- 53,57 ---- _gageVecFilter (gageContext *ctx, gagePerVolume *pvl) { char me[]="_gageVecFilter"; ! gage_t *fw00, *fw11, *fw22, *vec, *jac, *hes; int fd; *************** *** 49,52 **** --- 59,63 ---- 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); *************** *** 62,67 **** gageScl3PFilter2(pvl->iv3 + J*8, pvl->iv2 + J*4, pvl->iv1 + J*2, \ fw00, fw11, fw22, \ ! vec + J, jac + J*3, NULL, \ ! pvl->needD[0], pvl->needD[1], AIR_FALSE) DOIT_2(0); DOIT_2(1); DOIT_2(2); break; --- 73,79 ---- gageScl3PFilter2(pvl->iv3 + J*8, pvl->iv2 + J*4, pvl->iv1 + J*2, \ fw00, fw11, fw22, \ ! vec + J, jac + J*3, hes + J*9, \ ! pvl->needD[0], pvl->needD[1], pvl->needD[2]) ! /* 2nd order derivative computation with this scheme is no good idea */ DOIT_2(0); DOIT_2(1); DOIT_2(2); break; *************** *** 70,75 **** gageScl3PFilter4(pvl->iv3 + J*64, pvl->iv2 + J*16, pvl->iv1 + J*4, \ fw00, fw11, fw22, \ ! vec + J, jac + J*3, NULL, \ ! pvl->needD[0], pvl->needD[1], AIR_FALSE) DOIT_4(0); DOIT_4(1); DOIT_4(2); break; --- 82,87 ---- gageScl3PFilter4(pvl->iv3 + J*64, pvl->iv2 + J*16, pvl->iv1 + J*4, \ fw00, fw11, fw22, \ ! vec + J, jac + J*3, hes + J*9, \ ! pvl->needD[0], pvl->needD[1], pvl->needD[2]) DOIT_4(0); DOIT_4(1); DOIT_4(2); break; *************** *** 80,85 **** pvl->iv2 + J*fd*fd, pvl->iv1 + J*fd, \ fw00, fw11, fw22, \ ! vec + J, jac + J*3, NULL, \ ! pvl->needD[0], pvl->needD[1], AIR_FALSE) DOIT_N(0); DOIT_N(1); DOIT_N(2); break; --- 92,97 ---- pvl->iv2 + J*fd*fd, pvl->iv1 + J*fd, \ fw00, fw11, fw22, \ ! vec + J, jac + J*3, hes + J*9, \ ! pvl->needD[0], pvl->needD[1], pvl->needD[2]) DOIT_N(0); DOIT_N(1); DOIT_N(2); break; *************** *** 92,101 **** _gageVecAnswer (gageContext *ctx, gagePerVolume *pvl) { char me[]="_gageVecAnswer"; ! double tmpMat[9], mgevec[9], mgeval[3]; ! gage_t *vecAns, *normAns, *jacAns; - vecAns = pvl->directAnswer[gageVecVector]; - jacAns = pvl->directAnswer[gageVecJacobian]; - normAns = pvl->directAnswer[gageVecNormalized]; if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecVector)) { /* done if doV */ --- 104,122 ---- _gageVecAnswer (gageContext *ctx, gagePerVolume *pvl) { char me[]="_gageVecAnswer"; ! double cmag, tmpMat[9], mgevec[9], mgeval[3]; ! double symm[9], asym[9], tran[9], eval[3]; ! gage_t *vecAns, *normAns, *jacAns, *curlAns, *hesAns, *curlGradAns, ! *helGradAns, *dirHelDirAns; ! int asw; ! ! vecAns = pvl->directAnswer[gageVecVector]; ! normAns = pvl->directAnswer[gageVecNormalized]; ! jacAns = pvl->directAnswer[gageVecJacobian]; ! curlAns = pvl->directAnswer[gageVecCurl]; ! hesAns = pvl->directAnswer[gageVecHessian]; ! curlGradAns = pvl->directAnswer[gageVecCurlGradient]; ! helGradAns = pvl->directAnswer[gageVecHelGradient]; ! dirHelDirAns = pvl->directAnswer[gageVecDirHelDeriv]; if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecVector)) { /* done if doV */ *************** *** 109,113 **** } if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecNormalized)) { ! if (pvl->directAnswer[gageVecLength]) { ELL_3V_SCALE(normAns, 1.0/pvl->directAnswer[gageVecLength][0], vecAns); } else { --- 130,134 ---- } if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecNormalized)) { ! if (pvl->directAnswer[gageVecLength][0]) { ELL_3V_SCALE(normAns, 1.0/pvl->directAnswer[gageVecLength][0], vecAns); } else { *************** *** 116,120 **** } if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecJacobian)) { ! /* done if doV1 */ /* 0:dv_x/dx 1:dv_x/dy 2:dv_x/dz --- 137,141 ---- } if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecJacobian)) { ! /* done if doD1 */ /* 0:dv_x/dx 1:dv_x/dy 2:dv_x/dz *************** *** 136,144 **** } if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecCurl)) { ! ELL_3V_SET(pvl->directAnswer[gageVecCurl], jacAns[7] - jacAns[5], jacAns[2] - jacAns[6], jacAns[3] - jacAns[1]); } if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecGradient0)) { ELL_3V_SET(pvl->directAnswer[gageVecGradient0], --- 157,274 ---- } if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecCurl)) { ! ELL_3V_SET(curlAns, jacAns[7] - jacAns[5], jacAns[2] - jacAns[6], jacAns[3] - jacAns[1]); } + if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecHelicity)) { + pvl->directAnswer[gageVecHelicity][0] = + ELL_3V_DOT(vecAns, curlAns); + } + if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecNormHelicity)) { + cmag = ELL_3V_LEN(curlAns); + pvl->directAnswer[gageVecNormHelicity][0] = + cmag ? ELL_3V_DOT(normAns, curlAns)/cmag : 0; + } + if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecLambda2)) { + ELL_3M_TRANSPOSE(tran, jacAns); + // symmetric part + ELL_3M_SCALE_ADD2(symm, 0.5, jacAns, 0.5, tran); + // antisymmetric part + ELL_3M_SCALE_ADD2(asym, 0.5, jacAns, -0.5, tran); + // square symmetric part + ELL_3M_MUL(tmpMat, symm, symm); + ELL_3M_COPY(symm, tmpMat); + // square antisymmetric part + ELL_3M_MUL(tmpMat, asym, asym); + // sum of both + ELL_3M_ADD2(symm, symm, tmpMat); + asw = ell_3m_eigenvalues_d(eval, symm, AIR_TRUE); + if (asw == ell_cubic_root_triple) + pvl->directAnswer[gageVecLambda2][0] = eval[0]; + else if (asw == ell_cubic_root_single_double) + pvl->directAnswer[gageVecLambda2][0] = eval[1]; + else if (asw == ell_cubic_root_three) + { + if (eval[0]<=eval[1]) + if (eval[1]<=eval[2]) + pvl->directAnswer[gageVecLambda2][0] = eval[1]; + else + if (eval[0]<=eval[2]) + pvl->directAnswer[gageVecLambda2][0] = eval[2]; + else + pvl->directAnswer[gageVecLambda2][0] = eval[0]; + else + if (eval[2]<=eval[1]) + pvl->directAnswer[gageVecLambda2][0] = eval[1]; + else + if (eval[0]<=eval[2]) + pvl->directAnswer[gageVecLambda2][0] = eval[0]; + else + pvl->directAnswer[gageVecLambda2][0] = eval[2]; + } + } + /* 2nd order vector derivative continued */ + if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecHessian)) { + /* done if doD2 */ + /* the ordering is induced by the scalar hessian computation : + 0:d2v_x/dxdx 1:d2v_x/dxdy 2:d2v_x/dxdz + 3:d2v_x/dydx 4:d2v_x/dydy 5:d2v_x/dydz + 6:d2v_x/dzdx 7:d2v_x/dzdy 8:d2v_x/dzdz + 9:d2v_y/dxdx [...] + [...] + 24:dv2_z/dzdx 25:d2v_z/dzdy 26:d2v_z/dzdz + */ + if (ctx->verbose) { + fprintf(stderr, "%s: hes = \n", me); + ell_3m_PRINT(stderr, hesAns); /* ?? */ + } + } + if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecCurlGradient)) { + pvl->directAnswer[gageVecCurlGradient][0] = hesAns[21]-hesAns[15]; + pvl->directAnswer[gageVecCurlGradient][1] = hesAns[22]-hesAns[16]; + pvl->directAnswer[gageVecCurlGradient][2] = hesAns[23]-hesAns[17]; + pvl->directAnswer[gageVecCurlGradient][3] = hesAns[ 6]-hesAns[18]; + pvl->directAnswer[gageVecCurlGradient][4] = hesAns[ 7]-hesAns[19]; + pvl->directAnswer[gageVecCurlGradient][5] = hesAns[ 8]-hesAns[20]; + pvl->directAnswer[gageVecCurlGradient][6] = hesAns[ 9]-hesAns[ 1]; + pvl->directAnswer[gageVecCurlGradient][7] = hesAns[10]-hesAns[ 2]; + pvl->directAnswer[gageVecCurlGradient][8] = hesAns[11]-hesAns[ 3]; + } + if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecHelGradient)) { + pvl->directAnswer[gageVecHelGradient][0] = + jacAns[0]*curlAns[0]+ + jacAns[3]*curlAns[1]+ + jacAns[6]*curlAns[2]+ + curlGradAns[0]*vecAns[0]+ + curlGradAns[3]*vecAns[1]+ + curlGradAns[6]*vecAns[2]; + pvl->directAnswer[gageVecHelGradient][1] = + jacAns[1]*curlAns[0]+ + jacAns[4]*curlAns[1]+ + jacAns[7]*curlAns[2]+ + curlGradAns[1]*vecAns[0]+ + curlGradAns[4]*vecAns[1]+ + curlGradAns[7]*vecAns[2]; + pvl->directAnswer[gageVecHelGradient][0] = + jacAns[2]*curlAns[0]+ + jacAns[5]*curlAns[1]+ + jacAns[8]*curlAns[2]+ + curlGradAns[2]*vecAns[0]+ + curlGradAns[5]*vecAns[1]+ + curlGradAns[8]*vecAns[2]; + } + if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecDirHelDeriv)) { + pvl->directAnswer[gageVecDirHelDeriv][0] = + ELL_3V_DOT(normAns, helGradAns); + } + if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecProjHelGradient)) { + pvl->directAnswer[gageVecDirHelDeriv][0] = + helGradAns[0]-dirHelDirAns[0]*normAns[0]; + pvl->directAnswer[gageVecDirHelDeriv][1] = + helGradAns[1]-dirHelDirAns[0]*normAns[1]; + pvl->directAnswer[gageVecDirHelDeriv][2] = + helGradAns[2]-dirHelDirAns[0]*normAns[2]; + } if (GAGE_QUERY_ITEM_TEST(pvl->query, gageVecGradient0)) { ELL_3V_SET(pvl->directAnswer[gageVecGradient0], *************** *** 197,200 **** --- 327,338 ---- "divergence", "curl", + "helicity", + "normalized helicity", + "lambda2", + "vector hessian", + "curl gradient", + "helicity gradient", + "directional helicity derivative", + "projected helicity gradient", "gradient0", "gradient1", *************** *** 215,218 **** --- 353,364 ---- "divergence", "curl", + "helicity: dot(vector,curl)", + "normalized helicity", + "lambda2 value for vortex characterization", + "3x3x3 second-order vector derivative", + "3x3 derivative of curl", + "gradient of helicity", + "directional derivative of helicity along flow", + "projection of the helicity gradient onto plane orthogonal to flow", "gradient of 1st component of vector", "gradient of 2nd component of vector", *************** *** 233,236 **** --- 379,390 ---- gageVecDivergence, gageVecCurl, + gageVecHelicity, + gageVecNormHelicity, + gageVecLambda2, + gageVecHessian, + gageVecCurlGradient, + gageVecHelGradient, + gageVecDirHelDeriv, + gageVecProjHelGradient, gageVecGradient0, gageVecGradient1, *************** *** 248,251 **** --- 402,413 ---- #define GV_D gageVecDivergence #define GV_C gageVecCurl + #define GV_H gageVecHelicity + #define GV_NH gageVecNormHelicity + #define GV_LB gageVecLambda2 + #define GV_VH gageVecHessian + #define GV_CG gageVecCurlGradient + #define GV_HG gageVecHelGradient + #define GV_DH gageVecDirHelDeriv + #define GV_PH gageVecProjHelGradient #define GV_G0 gageVecGradient0 #define GV_G1 gageVecGradient1 *************** *** 264,267 **** --- 426,437 ---- "divergence", "div", "d", "curl", "c", + "h", "hel", "hell", + "nh", "nhel", "normhel", "normhell", + "lbda2", "lambda2", + "vh", "vhes", "vhessian", "vector hessian", + "cg", "curlgrad", "curlg", "curljac", "curl gradient", + "hg", "helg", "helgrad", "helicity gradient", + "dirhelderiv", "dhd", "ddh", "directional helicity derivative", + "phg", "projhel", "projhelgrad", "projected helicity gradient", "g0", "grad0", "gradient0", "g1", "grad1", "gradient1", *************** *** 282,285 **** --- 452,463 ---- GV_D, GV_D, GV_D, GV_C, GV_C, + GV_H, GV_H, GV_H, + GV_NH, GV_NH, GV_NH, GV_NH, + GV_LB, GV_LB, + GV_VH, GV_VH, GV_VH, GV_VH, + GV_CG, GV_CG, GV_CG, GV_CG, GV_CG, + GV_HG, GV_HG, GV_HG, GV_HG, + GV_DH, GV_DH, GV_DH, GV_DH, + GV_PH, GV_PH, GV_PH, GV_PH, GV_G0, GV_G0, GV_G0, GV_G1, GV_G1, GV_G1, Index: miscGage.c =================================================================== RCS file: /cvsroot/teem/teem/src/gage/miscGage.c,v retrieving revision 1.11 retrieving revision 1.12 diff -C2 -d -r1.11 -r1.12 *** miscGage.c 13 Mar 2004 20:03:09 -0000 1.11 --- miscGage.c 1 Apr 2004 07:07:35 -0000 1.12 *************** *** 36,40 **** */ gage_t ! gageZeroNormal[3] = {1,0,0}; char --- 36,40 ---- */ gage_t ! gageZeroNormal[3] = {0,0,0}; char |
