From 4c6dae39ca5a6649657ba80e65da9b6b2f56b804 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?C=C3=A9lestin=20Marot?= <celestin.marot@uclouvain.be>
Date: Tue, 5 May 2020 14:26:16 +0200
Subject: [PATCH] only mesh points that are actually part of the mesh

---
 contrib/hxt/tetMesh/src/hxt_tetRefine.c | 42 ++++++++++++++++++++++---
 1 file changed, 37 insertions(+), 5 deletions(-)

diff --git a/contrib/hxt/tetMesh/src/hxt_tetRefine.c b/contrib/hxt/tetMesh/src/hxt_tetRefine.c
index 6f2c1faa16..f10175b0ee 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);
     }
-- 
GitLab