diff --git a/contrib/hxt/tetMesh/exe/Delaunay_CLI.c b/contrib/hxt/tetMesh/exe/Delaunay_CLI.c
index 61a316299eed2eac30e86b65194922fb458f31fb..b61f9f1f38f06f23724d3517158d9203775f2ee1 100644
--- a/contrib/hxt/tetMesh/exe/Delaunay_CLI.c
+++ b/contrib/hxt/tetMesh/exe/Delaunay_CLI.c
@@ -157,6 +157,7 @@ int main(int argc, char** argv) {
.bbox = NULL,
.nodalSizes = NULL,
.numVerticesInMesh = 0,
+ .insertionFirst = 0,
.partitionability = 0,
.verbosity = verbosity,
.reproducible = reproducible,
diff --git a/contrib/hxt/tetMesh/include/hxt_tetDelaunay.h b/contrib/hxt/tetMesh/include/hxt_tetDelaunay.h
index 70093ed7a6a9d0d7e6d7b396ef88d9c57536d5c2..da9357b078d49dcb3e536fbe52b3be4f9445cd63 100644
--- a/contrib/hxt/tetMesh/include/hxt_tetDelaunay.h
+++ b/contrib/hxt/tetMesh/include/hxt_tetDelaunay.h
@@ -44,7 +44,8 @@ typedef struct {
* \warning a segmentation fault will occur if a vertex
* doesn't have a corresponding mesh size */
- uint32_t numVerticesInMesh; /**< The number of vertices in the mesh (all vertices below this point are not (re-)inserted */
+ uint32_t numVerticesInMesh; /**< the number of vertices in the mesh. This value gets incremented by hxtDelaunay */
+ uint32_t insertionFirst; /**< all vertices with an index below insertionFirst are not (re-)inserted and not reordered */
double partitionability; /**< a number between 0 and 1 telling if this mesh is good for making partitions.
Generally, put 0 for an empty mesh, 1-(1/2)^n for a mesh refined n time */
@@ -82,9 +83,10 @@ typedef struct {
* - hxtNodeInfo[i].hilbertDist will change
* - mesh->tetrahedra.* will change
* - mesh->vertices.coord[4*i+3] will change
+ * - options->numVerticesInMesh will change
*
* \param mesh: a valid Delaunay mesh
- * \param options: options to give to the Delaunay algorithm \ref HXTDelaunayOptions (here options->numVerticesInMesh is useless)
+ * \param options: options to give to the Delaunay algorithm \ref HXTDelaunayOptions
* \param[in, out] nodeInfo: the indices of the vertices to insert in the tetrahedral mesh.
* \param nToInsert: the number of element in nodeInfo, hence the number of vertices to insert.
*/
@@ -94,14 +96,14 @@ HXTStatus hxtDelaunaySteadyVertices(HXTMesh* mesh, HXTDelaunayOptions* options,
/**
* \brief Delaunay of a set of vertices
* \details This perform the insertion of the vertices
- * from numVerticesInMesh to mesh->vertices.num\n
+ * from insertionFirst to mesh->vertices.num\n
*
* \warning
- * - the order of mesh->vertices will change
- * - hxtNodeInfo[i].hilbertDist will change
+ * - the order of mesh->vertices will change, except for those below options->insertionFirst
* - mesh->tetrahedra.* will change
* - mesh->vertices.coord[4*i+3] will change
* - vertices that could not be inserted are deleted from mesh->vertices !
+ * - options->numVerticesInMesh will change
*
* \param mesh: a valid Delaunay mesh
* \param options: options to give to the Delaunay algorithm \ref HXTDelaunayOptions
@@ -109,7 +111,7 @@ HXTStatus hxtDelaunaySteadyVertices(HXTMesh* mesh, HXTDelaunayOptions* options,
HXTStatus hxtDelaunay(HXTMesh* mesh, HXTDelaunayOptions* options);
-/**
+/** ! THIS PART IS NOT IMPLEMENTED YET ! (I keep this for the future)
* \brief just add one vertex in the mesh
* \details This perform the single insertion of the vertex vta (vertex to add) in the mesh given by mesh
* To avoid repetitive allocations and effort to directly eliminate deleted tetrahedra from the mesh,
@@ -120,7 +122,7 @@ HXTStatus hxtDelaunay(HXTMesh* mesh, HXTDelaunayOptions* options);
* => options->reproducible and options->delaunayThreads are useless
*
* \param mesh: a valid Delaunay mesh
- * \param options: options to give to the Delaunay algorithm \ref HXTDelaunayOptions (here options->numVerticesInMesh is useless)
+ * \param options: options to give to the Delaunay algorithm \ref HXTDelaunayOptions (here options->insertionFirst is useless)
* \param vtaNodeInfo[in, out]: the index of the vertex to insert in the tetrahedral mesh
* (see \ref hxtDelaunaySteadyVertices to understand how it gets filled)
* \param deleted: a valid pointer to an array of deleted tetrahedra (that must be marked as deleted as well => see hxt_tetFlag.c)
diff --git a/contrib/hxt/tetMesh/src/hxt_tetDelaunay.c b/contrib/hxt/tetMesh/src/hxt_tetDelaunay.c
index 8816d9e94023370cf0bfb9d94f0eb28d7a73b433..14a97ff5abbf6a96a5a5b4a5e116076573acca02 100644
--- a/contrib/hxt/tetMesh/src/hxt_tetDelaunay.c
+++ b/contrib/hxt/tetMesh/src/hxt_tetDelaunay.c
@@ -1429,6 +1429,7 @@ static HXTStatus parallelDelaunay3D(HXTMesh* mesh,
Initializations and allocations
******************************************************/
if(mesh->tetrahedra.num<5){
+ HXT_ASSERT_MSG(options->numVerticesInMesh==0 && mesh->tetrahedra.num==0, "no ghosts or unvalid mesh");
HXT_INFO_COND(options->verbosity>0,
"Initialization of tet. mesh");
HXT_CHECK( hxtTetrahedraInit(mesh, nodeInfo, nToInsert, options->verbosity) );
@@ -1949,40 +1950,14 @@ static HXTStatus DelaunayOptionsInit(HXTMesh* mesh,
HXTDelaunayOptions* options,
HXTBbox* bbox){
HXT_ASSERT(mesh!=NULL);
-HXT_ASSERT_MSG(options->numVerticesInMesh==0 || mesh->tetrahedra.num!=0,
- "Told there was 0 vertices in the mesh, which is incompatible with the number of tetrahedra (>0)");
+HXT_ASSERT_MSG((options->numVerticesInMesh==0 && mesh->tetrahedra.num==0) ||
+ (options->numVerticesInMesh>=4 && mesh->tetrahedra.num>0),
+ "Number of vertices in the mesh and number of tetrahedra cannot match");
const uint32_t nVert = mesh->vertices.num;
- if(options->numVerticesInMesh==0 && mesh->tetrahedra.num!=0) {
-
- HXTVertex* vertices = (HXTVertex*) mesh->vertices.coord;
-
- // count the number of vertices in the mesh
- #pragma omp parallel for
- for (uint32_t i=0; i<nVert; i++) {
- vertices[i].padding.index = 0;
- }
-
- #pragma omp parallel for
- for (uint64_t i=0; i<mesh->tetrahedra.num; i++) {
- vertices[mesh->tetrahedra.node[4*i+0]].padding.index = 1;
- vertices[mesh->tetrahedra.node[4*i+1]].padding.index = 1;
- vertices[mesh->tetrahedra.node[4*i+2]].padding.index = 1;
- if(mesh->tetrahedra.node[4*i+3]!=HXT_GHOST_VERTEX)
- vertices[mesh->tetrahedra.node[4*i+3]].padding.index = 1;
- }
-
- uint32_t numVerticesInMesh = 0;
- #pragma omp parallel for reduction(+:numVerticesInMesh)
- for (uint32_t i=0; i<nVert; i++) {
- numVerticesInMesh += vertices[i].padding.index;
- }
-
- options->numVerticesInMesh = numVerticesInMesh;
- }
-
HXT_ASSERT(options->numVerticesInMesh <= nVert);
+HXT_ASSERT(options->insertionFirst <= nVert);
if(options->bbox==NULL){
options->bbox = bbox;
@@ -2011,30 +1986,21 @@ HXT_ASSERT(options->numVerticesInMesh <= nVert);
* parallel Delaunay
* see header for a complete description
****************************************/
-HXTStatus hxtDelaunay(HXTMesh* mesh, HXTDelaunayOptions* userOptions){
- HXTDelaunayOptions options;
- if(userOptions!=NULL)
- options = *userOptions;
- else
- options = (HXTDelaunayOptions) {.bbox=NULL,
- .nodalSizes=NULL,
- .numVerticesInMesh=0,
- .verbosity=1,
- .delaunayThreads=0,
- .reproducible=0};
-
+HXTStatus hxtDelaunay(HXTMesh* mesh, HXTDelaunayOptions* options){
+HXT_ASSERT(mesh!=NULL);
+HXT_ASSERT(options!=NULL);
HXTBbox bbox;
- HXT_CHECK( DelaunayOptionsInit(mesh, &options, &bbox) );
+ HXT_CHECK( DelaunayOptionsInit(mesh, options, &bbox) );
- const uint32_t nToInsert = mesh->vertices.num - options.numVerticesInMesh;
+ const uint32_t nToInsert = mesh->vertices.num - options->insertionFirst;
- if(options.reproducible && nToInsert<16554) // not worth launching threads and having to reorder tets after...
- options.delaunayThreads = 1;
+ // if(options.reproducible && nToInsert<16554) // not worth launching threads and having to reorder tets after...
+ // options.delaunayThreads = 1;
hxtNodeInfo* nodeInfo;
HXT_CHECK( hxtAlignedMalloc(&nodeInfo, nToInsert*sizeof(hxtNodeInfo)) );
- HXT_CHECK( parallelDelaunay3D(mesh, &options, nodeInfo, nToInsert, 0) );
+ HXT_CHECK( parallelDelaunay3D(mesh, options, nodeInfo, nToInsert, 0) );
HXT_CHECK( hxtAlignedFree(&nodeInfo) );
@@ -2046,29 +2012,19 @@ HXTStatus hxtDelaunay(HXTMesh* mesh, HXTDelaunayOptions* userOptions){
* parallel Delaunay without moving the vertices
* see header for a complete description
***********************************************/
-HXTStatus hxtDelaunaySteadyVertices(HXTMesh* mesh, HXTDelaunayOptions* userOptions, hxtNodeInfo* nodeInfo, uint64_t nToInsert){
+HXTStatus hxtDelaunaySteadyVertices(HXTMesh* mesh, HXTDelaunayOptions* options, hxtNodeInfo* nodeInfo, uint64_t nToInsert){
+HXT_ASSERT(mesh!=NULL);
+HXT_ASSERT(options!=NULL);
HXT_ASSERT(nodeInfo!=NULL);
-
- HXTDelaunayOptions options;
- if(userOptions!=NULL)
- options = *userOptions;
- else
- options = (HXTDelaunayOptions) {.bbox=NULL,
- .nodalSizes=NULL,
- .numVerticesInMesh=0,
- .verbosity=1,
- .delaunayThreads=0,
- .reproducible=0};
-
HXTBbox bbox;
- HXT_CHECK( DelaunayOptionsInit(mesh, &options, &bbox) );
+ HXT_CHECK( DelaunayOptionsInit(mesh, options, &bbox) );
- if(options.reproducible && nToInsert<16554) // not worth launching threads and having to reorder tets after...
- options.delaunayThreads = 1;
+ // if(options.reproducible && nToInsert<16554) // not worth launching threads and having to reorder tets after...
+ // options.delaunayThreads = 1;
-HXT_ASSERT(options.numVerticesInMesh+nToInsert <= mesh->vertices.num);
+ HXT_ASSERT(options->insertionFirst+nToInsert <= mesh->vertices.num);
- HXT_CHECK( parallelDelaunay3D(mesh, &options, nodeInfo, nToInsert, 1) );
+ HXT_CHECK( parallelDelaunay3D(mesh, options, nodeInfo, nToInsert, 1) );
return HXT_STATUS_OK;
}
diff --git a/contrib/hxt/tetMesh/src/hxt_tetMesh.c b/contrib/hxt/tetMesh/src/hxt_tetMesh.c
index 22544f7f95d04dbd1050db04495caf8bc0b5ce5e..950691d63c208c859dcf3a0324e9b80fa9e0317b 100644
--- a/contrib/hxt/tetMesh/src/hxt_tetMesh.c
+++ b/contrib/hxt/tetMesh/src/hxt_tetMesh.c
@@ -108,11 +108,12 @@ HXTStatus hxtTetMesh(HXTMesh* mesh,
HXT_INFO("Recovering %" HXTu64 " missing edge(s)",
nbMissingLines);
+ uint32_t oldNumVertices = mesh->vertices.num;
HXT_CHECK(options->recoveryFun(mesh, options->recoveryData));
- if(delOptions.numVerticesInMesh < mesh->vertices.num) {
+ if(oldNumVertices < mesh->vertices.num) {
HXT_INFO("Steiner(s) point(s) were inserted");
- delOptions.numVerticesInMesh = mesh->vertices.num;
+ delOptions.numVerticesInMesh += mesh->vertices.num - oldNumVertices;
}
t[3] = omp_get_wtime();
diff --git a/contrib/hxt/tetMesh/src/hxt_tetRefine.c b/contrib/hxt/tetMesh/src/hxt_tetRefine.c
index f10175b0eed714b5fc5b2547c8cc5e4071237861..925c5d4da08a0822beea06af71e24575263f7fc2 100644
--- a/contrib/hxt/tetMesh/src/hxt_tetRefine.c
+++ b/contrib/hxt/tetMesh/src/hxt_tetRefine.c
@@ -62,7 +62,6 @@ HXTStatus hxtEmptyMesh(HXTMesh* mesh, HXTDelaunayOptions* delOptions)
}
HXT_CHECK( hxtDelaunaySteadyVertices(mesh, delOptions, nodeInfo, numToInsert) );
- delOptions->numVerticesInMesh = mesh->vertices.num;
#ifdef DEBUG
#pragma omp parallel for simd aligned(nodeInfo:SIMD_ALIGN)
@@ -335,6 +334,7 @@ HXTStatus hxtRefineTetrahedra(HXTMesh* mesh, HXTDelaunayOptions* delOptions,
}
+ delOptions->insertionFirst = mesh->vertices.num; // prepare for the next insertion
uint32_t add = 0;
HXTStatus status = HXT_STATUS_OK;
@@ -484,13 +484,11 @@ HXTStatus hxtRefineTetrahedra(HXTMesh* mesh, HXTDelaunayOptions* delOptions,
mesh->vertices.num += add;
delOptions->partitionability = 1.0 - pow(0.5, iter);
+ uint32_t oldNumVerticesInMesh = delOptions->numVerticesInMesh;
HXT_CHECK(hxtDelaunay(mesh, delOptions));
- uint32_t numAdd = mesh->vertices.num - delOptions->numVerticesInMesh;
- delOptions->numVerticesInMesh = mesh->vertices.num;
-
- if (numAdd == 0) break;
+ if (delOptions->numVerticesInMesh == oldNumVerticesInMesh) break;
}
HXT_CHECK( hxtFree(&startIndex) );