diff --git a/contrib/hxt/tetMesh/src/hxt_tetRefine.c b/contrib/hxt/tetMesh/src/hxt_tetRefine.c index 6f2c1faa16f000f6f840d0b706ff70b657270659..f10175b0eed714b5fc5b2547c8cc5e4071237861 100644 --- a/contrib/hxt/tetMesh/src/hxt_tetRefine.c +++ b/contrib/hxt/tetMesh/src/hxt_tetRefine.c @@ -11,6 +11,33 @@ #include "hxt_tetFlag.h" #include "hxt_omp.h" +// mark all the points which are in mesh->(points | lines | triangles) +static void markMeshPoints(HXTMesh* mesh) +{ + #pragma omp parallel for simd + for(uint32_t i=0; i<mesh->vertices.num; i++) { + mesh->vertices.coord[4*i+3] = 0.0; + } + + for(uint64_t i=0; i<mesh->triangles.num; i++) { + for(int j=0; j<3; j++) { + mesh->vertices.coord[4* mesh->triangles.node[3*i+j] + 3] = 1.0; + } + } + + for(uint64_t i=0; i<mesh->lines.num; i++) { + for(int j=0; j<2; j++) { + mesh->vertices.coord[4* mesh->lines.node[2*i+j] + 3] = 1.0; + } + } + + for(uint64_t i=0; i<mesh->points.num; i++) { + for(int j=0; j<2; j++) { + mesh->vertices.coord[4* mesh->lines.node[2*i+j] + 3] = 1.0; + } + } +} + HXTStatus hxtEmptyMesh(HXTMesh* mesh, HXTDelaunayOptions* delOptions) { @@ -23,18 +50,23 @@ HXTStatus hxtEmptyMesh(HXTMesh* mesh, HXTDelaunayOptions* delOptions) hxtNodeInfo* nodeInfo; HXT_CHECK( hxtAlignedMalloc(&nodeInfo, sizeof(hxtNodeInfo)*mesh->vertices.num) ); - #pragma omp parallel for simd aligned(nodeInfo:SIMD_ALIGN) + markMeshPoints(mesh); + + uint32_t numToInsert = 0; for (uint32_t i=0; i<mesh->vertices.num; i++) { - nodeInfo[i].node = i; - nodeInfo[i].status = HXT_STATUS_TRYAGAIN; + if(mesh->vertices.coord[4 * i + 3]==1.0) { + nodeInfo[numToInsert].node = i; + nodeInfo[numToInsert].status = HXT_STATUS_TRYAGAIN; + numToInsert++; + } } - HXT_CHECK( hxtDelaunaySteadyVertices(mesh, delOptions, nodeInfo, mesh->vertices.num) ); + HXT_CHECK( hxtDelaunaySteadyVertices(mesh, delOptions, nodeInfo, numToInsert) ); delOptions->numVerticesInMesh = mesh->vertices.num; #ifdef DEBUG #pragma omp parallel for simd aligned(nodeInfo:SIMD_ALIGN) - for (uint32_t i=0; i<mesh->vertices.num; i++) { + for (uint32_t i=0; i<numToInsert; i++) { if(nodeInfo[i].status!=HXT_STATUS_TRUE){ HXT_WARNING("point %u of the empty mesh was not inserted\n", nodeInfo[i].node); }