From add7ec0d43e734f8816d664e622693c4ba91622f Mon Sep 17 00:00:00 2001 From: Qianqian Fang Date: Sun, 17 Nov 2024 11:03:23 -0500 Subject: [PATCH] [feat] automatically find initial element id, partially done --- src/mmc_mesh.c | 131 ++++++++++++++++++++++++++++++++++++------------- src/mmc_mesh.h | 2 + 2 files changed, 99 insertions(+), 34 deletions(-) diff --git a/src/mmc_mesh.c b/src/mmc_mesh.c index 2797a4e..b2e1812 100644 --- a/src/mmc_mesh.c +++ b/src/mmc_mesh.c @@ -913,6 +913,96 @@ void mesh_loadseedfile(tetmesh* mesh, mcconfig* cfg) { #endif +/** + * @brief Scan all tetrahedral elements to find the one enclosing the source + * + * @param[in] mesh: the mesh object + * @param[in] cfg: the simulation configuration structure + * + * @return 0 if an enclosing element is found, 1 if not found + */ + +int mesh_initelem(tetmesh* mesh, mcconfig* cfg) { + FLOAT3* nodes = mesh->node; + int i, j; + + for (i = 0; i < mesh->ne; i++) { + float3 pmin = {VERY_BIG, VERY_BIG, VERY_BIG}, pmax = {-VERY_BIG, -VERY_BIG, -VERY_BIG}; + int* elems = (int*)(mesh->elem + i * mesh->elemlen); // convert int4* to int* + + for (j = 0; j < mesh->elemlen; j++) { + pmin.x = MIN(nodes[elems[0] - 1].x, pmin.x); + pmin.y = MIN(nodes[elems[0] - 1].y, pmin.y); + pmin.z = MIN(nodes[elems[0] - 1].z, pmin.z); + + pmax.x = MAX(nodes[elems[0] - 1].x, pmax.x); + pmax.y = MAX(nodes[elems[0] - 1].y, pmax.y); + pmax.z = MAX(nodes[elems[0] - 1].z, pmax.z); + } + + if (cfg->srcpos.x <= pmax.x && cfg->srcpos.x >= pmin.x && + cfg->srcpos.y <= pmax.y && cfg->srcpos.y >= pmin.y && + cfg->srcpos.z <= pmax.z && cfg->srcpos.z >= pmin.z) { + if (mesh_barycentric(i + 1, &(cfg->bary0.x), (FLOAT3*) & (cfg->srcpos), mesh) == 0) { + cfg->e0 = i + 1; + return 0; + } + } + } + + return 1; +} + + +/** + * @brief Compute the barycentric coordinate of the source in the initial element + * + * @param[in] e0: the index (start from 1) of the initial element + * @param[in] bary: the index (start from 1) of the initial element + * @param[in] srcpos: the source position + * @param[in] mesh: the mesh of the domain + * + * @return 0 if all barycentric coordinates are computed and all positive, or 1 any of those is negative + */ + +int mesh_barycentric(int e0, float* bary, FLOAT3* srcpos, tetmesh* mesh) { + int eid = e0 - 1, i; + FLOAT3 vecS = {0.f}, vecAB, vecAC, vecN; + FLOAT3* nodes = mesh->node; + int ea, eb, ec; + float s = 0.f; + int* elems = (int*)(mesh->elem + eid * mesh->elemlen); // convert int4* to int* + + if (eid >= mesh->ne) { + MESH_ERROR("initial element index exceeds total element count"); + } + + for (i = 0; i < 4; i++) { + ea = elems[out[i][0]] - 1; + eb = elems[out[i][1]] - 1; + ec = elems[out[i][2]] - 1; + vec_diff3(&nodes[ea], &nodes[eb], &vecAB); + vec_diff3(&nodes[ea], &nodes[ec], &vecAC); + vec_diff3(&nodes[ea], srcpos, &vecS); + vec_cross3(&vecAB, &vecAC, &vecN); + bary[facemap[i]] = -vec_dot3(&vecS, &vecN); + } + + for (i = 0; i < 4; i++) { + if (bary[i] < 0.f) { + return 1; + } + + s += bary[i]; + } + + for (i = 0; i < 4; i++) { + bary[i] /= s; + } + + return 0; +} + /** * @brief Initialize a data structure storing all pre-computed ray-tracing related data * @@ -935,6 +1025,7 @@ void tracer_init(raytracer* tracer, tetmesh* pmesh, char methodid) { tracer_build(tracer); } + /** * @brief Preparing for the ray-tracing calculations * @@ -954,44 +1045,16 @@ void tracer_prep(raytracer* tracer, mcconfig* cfg) { } else { MESH_ERROR("tracer is not associated with a mesh"); } - } else if ( (cfg->srctype == stPencil || cfg->srctype == stIsotropic || cfg->srctype == stCone || cfg->srctype == stArcSin) && cfg->e0 > 0) { - int eid = cfg->e0 - 1; - FLOAT3 vecS = {0.f}, vecAB, vecAC, vecN; - FLOAT3* nodes = tracer->mesh->node; - int ea, eb, ec; - float s = 0.f, *bary = &(cfg->bary0.x); - int* elems = (int*)(tracer->mesh->elem + eid * tracer->mesh->elemlen); // convert int4* to int* - - if (eid >= tracer->mesh->ne) { - MESH_ERROR("initial element index exceeds total element count"); - } - - for (i = 0; i < 4; i++) { - ea = elems[out[i][0]] - 1; - eb = elems[out[i][1]] - 1; - ec = elems[out[i][2]] - 1; - vec_diff3(&nodes[ea], &nodes[eb], &vecAB); - vec_diff3(&nodes[ea], &nodes[ec], &vecAC); - vec_diff3(&nodes[ea], (FLOAT3*) & (cfg->srcpos), &vecS); - vec_cross3(&vecAB, &vecAC, &vecN); - bary[facemap[i]] = -vec_dot3(&vecS, &vecN); - } - - if (cfg->debuglevel & dlWeight) - fprintf(cfg->flog, "initial bary-centric volumes [%e %e %e %e]\n", - bary[0] / 6., bary[1] / 6., bary[2] / 6., bary[3] / 6.); - - for (i = 0; i < 4; i++) { - if (bary[i] < 0.f) { + } else if ( (cfg->srctype == stPencil || cfg->srctype == stIsotropic || cfg->srctype == stCone || cfg->srctype == stArcSin) ) { + if (cfg->e0 <= 0 || mesh_barycentric(cfg->e0, &cfg->bary0.x, (FLOAT3*) & (cfg->srcpos), tracer->mesh)) { + if (mesh_initelem(tracer->mesh, cfg)) { MESH_ERROR("initial element does not enclose the source!"); } - - s += bary[i]; } - for (i = 0; i < 4; i++) { - bary[i] /= s; - } + if (cfg->debuglevel & dlWeight) + fprintf(cfg->flog, "initial bary-centric volumes [%e %e %e %e]\n", + cfg->bary0.x / 6., cfg->bary0.y / 6., cfg->bary0.z / 6., cfg->bary0.w / 6.); } // TODO: this is a partial fix to the normalization bug described in https://github.com/fangq/mmc/issues/82 diff --git a/src/mmc_mesh.h b/src/mmc_mesh.h index 2da5ca1..6ed26a1 100644 --- a/src/mmc_mesh.h +++ b/src/mmc_mesh.h @@ -158,6 +158,8 @@ void mesh_createdualmesh(tetmesh* mesh, mcconfig* cfg); void mesh_loadroi(tetmesh* mesh, mcconfig* cfg); double mesh_getreff_approx(double n_in, double n_out); double mesh_getreff(double n_in, double n_out); +int mesh_barycentric(int e0, float* bary, FLOAT3* srcpos, tetmesh* mesh); +int mesh_initelem(tetmesh* mesh, mcconfig* cfg); void tracer_init(raytracer* tracer, tetmesh* mesh, char methodid); void tracer_build(raytracer* tracer);