From 9f1c3552e72c45b83fbf9314c9b448d0cfe7bbbf Mon Sep 17 00:00:00 2001 From: Qianqian Fang Date: Sun, 17 Nov 2024 23:50:54 -0500 Subject: [PATCH] [feat] accelerate facenb using std::vector, fix initelem bug --- src/mmc_highorder.cpp | 43 ++++++++++++++----------------------------- src/mmc_mesh.c | 41 +++++++++++++++++------------------------ 2 files changed, 31 insertions(+), 53 deletions(-) diff --git a/src/mmc_highorder.cpp b/src/mmc_highorder.cpp index b3c85ab..1d1e852 100644 --- a/src/mmc_highorder.cpp +++ b/src/mmc_highorder.cpp @@ -100,47 +100,33 @@ void mesh_10nodetet(tetmesh* mesh, mcconfig* cfg) { for (it = edgelist.begin(); it != edgelist.end(); it++) { for (int i = 0; i < 3; i++) { ((float*)(&mesh->node[oldnn + pos]))[i] = - (((float*)(&mesh->node[(*it).first]))[i] + ((float*)(&mesh->node[(*it).second]))[i]) * 0.5f; + (((float*)(&mesh->node[it->first]))[i] + ((float*)(&mesh->node[it->second]))[i]) * 0.5f; } pos++; } } -std::string face_to_string(int a, int b, int c) { - int ar[3]; - ar[0] = MIN(MIN(a, b), c); - ar[2] = MAX(MAX(a, b), c); - ar[1] = - ar[0] - ar[2] + a + b + c; - - std::ostringstream strkey; - strkey << ar[0] << "," << ar[1] << "," << ar[2]; - - return std::string(strkey.str()); -} #ifdef __cplusplus extern "C" #endif void mesh_getfacenb(tetmesh* mesh, mcconfig* cfg) { - std::vector faces; - std::map< std::string, std::pair > facenb; - - faces.reserve(mesh->ne << 2); + std::map< std::vector, std::pair > facenb; + std::map< std::vector, std::pair > ::iterator it; for (int i = 0; i < mesh->ne; i++) { int* ee = mesh->elem + i * mesh->elemlen; for (int j = 0; j < 4; j++) { - faces.push_back(face_to_string(ee[facelist[j][0]], ee[facelist[j][1]], ee[facelist[j][2]])); - } - } + std::vector facevec = { ee[facelist[j][0]], ee[facelist[j][1]], ee[facelist[j][2]] }; + std::sort (facevec.begin(), facevec.end()); - for (unsigned int i = 0; i < faces.size(); i++) { - if (facenb.count(faces[i])) { - facenb[faces[i]].second = (i >> 2) + 1; - } else { - facenb[faces[i]] = std::make_pair((i >> 2) + 1, 0); + if (facenb.count(facevec)) { + facenb[facevec].second = (i << 2) + j + 1; + } else { + facenb[facevec] = std::make_pair((i << 2) + j, 0); + } } } @@ -150,11 +136,10 @@ void mesh_getfacenb(tetmesh* mesh, mcconfig* cfg) { mesh->facenb = (int*)calloc(sizeof(int) * mesh->elemlen, mesh->ne); - for (unsigned int i = 0; i < faces.size(); i++) { - std::pair item = facenb[faces[i]]; - - if (item.second != 0) { - mesh->facenb[i] = ((i >> 2) + 1 == item.first ? item.second : item.first); + for (it = facenb.begin(); it != facenb.end(); it++) { + if (it->second.second != 0) { + mesh->facenb[it->second.first] = ((it->second.second - 1) >> 2) + 1; + mesh->facenb[it->second.second - 1] = (it->second.first >> 2) + 1; } } } diff --git a/src/mmc_mesh.c b/src/mmc_mesh.c index b2e1812..295ef71 100644 --- a/src/mmc_mesh.c +++ b/src/mmc_mesh.c @@ -731,15 +731,11 @@ void mesh_loadelemvol(tetmesh* mesh, mcconfig* cfg) { int tmp, len, i, j, *ee; char fvelem[MAX_FULL_PATH]; - if (cfg->elem && cfg->elemnum && cfg->node && cfg->nodenum) { - mesh_getvolume(mesh, cfg); - return; - } - mesh_filenames("velem_%s.dat", fvelem, cfg); - if ((fp = fopen(fvelem, "rt")) == NULL) { - MESH_ERROR("can not open element volume file"); + if ((cfg->elem && cfg->elemnum && cfg->node && cfg->nodenum) || (fp = fopen(fvelem, "rt")) == NULL) { + mesh_getvolume(mesh, cfg); + return; } len = fscanf(fp, "%d %d", &tmp, &(mesh->ne)); @@ -784,15 +780,11 @@ void mesh_loadfaceneighbor(tetmesh* mesh, mcconfig* cfg) { int* pe; char ffacenb[MAX_FULL_PATH]; - if (cfg->elem && cfg->elemnum) { - mesh_getfacenb(mesh, cfg); - return; - } - mesh_filenames("facenb_%s.dat", ffacenb, cfg); - if ((fp = fopen(ffacenb, "rt")) == NULL) { - MESH_ERROR("can not open face-neighbor list file"); + if ((cfg->elem && cfg->elemnum) || (fp = fopen(ffacenb, "rt")) == NULL) { + mesh_getfacenb(mesh, cfg); + return; } len = fscanf(fp, "%d %d", &(mesh->elemlen), &(mesh->ne)); @@ -927,22 +919,23 @@ int mesh_initelem(tetmesh* mesh, mcconfig* cfg) { 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}; + double pmin[3] = {VERY_BIG, VERY_BIG, VERY_BIG}, pmax[3] = {-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); + pmin[0] = MIN(nodes[elems[j] - 1].x, pmin[0]); + pmin[1] = MIN(nodes[elems[j] - 1].y, pmin[1]); + pmin[2] = MIN(nodes[elems[j] - 1].z, pmin[2]); - 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); + pmax[0] = MAX(nodes[elems[j] - 1].x, pmax[0]); + pmax[1] = MAX(nodes[elems[j] - 1].y, pmax[1]); + pmax[2] = MAX(nodes[elems[j] - 1].z, pmax[2]); } - 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 (cfg->srcpos.x <= pmax[0] && cfg->srcpos.x >= pmin[0] && + cfg->srcpos.y <= pmax[1] && cfg->srcpos.y >= pmin[1] && + cfg->srcpos.z <= pmax[2] && cfg->srcpos.z >= pmin[2]) { + if (mesh_barycentric(i + 1, &(cfg->bary0.x), (FLOAT3*) & (cfg->srcpos), mesh) == 0) { cfg->e0 = i + 1; return 0;